Боль оптимизационных задач
Не от хорошей жизни я всякую дичь пробую, сначала это были кватернионы в анизотропном пространстве ( раз, два), а недавно - иерархические математические структуры в QuatCore.
Всё дело в том, что алгоритм сопровождения, основанный на решении нелинейной задачи условной оптимизации, упорно не хочет вписываться в 16-битные целые числа.
Сегодня, чтобы на интуитивном уровне понять, что с этим делать, решил рассмотреть простейшую "игрушечную задачу", снизив число степеней свободы с шести ДО ОДНОЙ, и число точек ДО ДВУХ. Как ни странно, даже тогда ФУНДАМЕНТАЛЬНАЯ ПРОБЛЕМА сохраняется, хотя здесь с ней ещё как-то можно справиться...
Итак, у нас есть неизвестное число x, которое мы хотим измерить. У нас есть два зашумлённых входа h1, h2 которые реагируют на изменение числа x, но чувствительность различается:
Здесь ξ η - случайные величины с нулевым математическим ожиданием (мы считаем, что всё хорошо откалибровали и всю систематику убрали) и дисперсией, для простоты, равной ЕДИНИЦЕ.
К примеру, h1, h2 измеряются в вольтах, а величина x - в метрах. Тогда "чувствительности" k1, k2 имеют размерность В/м.
Задача: как лучше всего на основании h1, h2 оценить x, т.е найти
, чтобы результат был несмещённым, а дисперсия - как можно ниже?
Исходя из размерности, вроде как "чувствительности" должны оказаться в знаменателе, чтобы вольты снова превратить в метры. По крайней мере, если бы вход был лишь один, h1, то нам бы ничего другого не оставалось, как записать
Но стоит начать оптимизацию - и дела немного усложняются...
Выводится всё очень легко. Ответ должен быть взвешенной суммой h1 и h2:
Подставим сюда выражения для h1, h2:
x+\alpha\xi+\beta\eta)
Чтобы выражение было несмещённым (математическое ожидание равнялось в точности x), скобочка должна обращаться в единицу:
У нас осталась одна степень свободы в выборе α, β. Поэтому выпишем выражение для дисперсии:
=\alpha^2D(\xi)+\beta^2D(\eta)=\alpha^2+\beta^2\rightarrow{min})
Но оптимизировать значения α, β надо с учётом ограничения выше, на несмещённость. Тут можно было бы просто β выразить через α, но это некрасиво, давайте уж воспользуемся множителем Лагранжа для условной оптимизации:
\rightarrow{min})
Чтобы найти экстремум, надо взять частные производные по α, β и λ и приравнять к нулю:
Удобно взять чуть другую, смасштабированную "лямбду" (внести в неё знак минуса и деление на два, чтобы лишнего не писать):
Теперь α и β, выраженные через λ' подставляем в третью строчку:
=1\Leftrightarrow\lambda%27=\frac{1}{k_1^2+k_2^2}.)
И наконец, находим α и β
Оценка x запишется так:
По размерности всё верно. То, что мы для начала УМНОЖАЕМ каждый вход на его чувствительность, вместо того, чтобы делить - имеет ясную интерпретацию, ЧЕМ ЧУВСТВИТЕЛЬНЕЕ ВХОД (при том же шуме, что и другие) - ТЕМ ОН ПОЛЕЗНЕЕ, поэтому именно он будет иметь бОльший вес! Можно рассмотреть вырожденные случаи:
1. h2 вообще не реагирует на измеряемую величину, т.е k2=0. Тогда он выкидывается из выражения, остаётся попросту h1/k1.
2. h2 СУЩЕСТВЕННО более чувствителен, т.е k2>>k1. Тогда, поделив числитель и знаменатель на его квадрат и вычеркнув "бесконечно малые", придём к выражению h2/k2, т.е не видим никакого смысла использовать вход h1.
Подобное выражение весьма универсально, в более сложных задачах оно остаётся практически неизменным, только числа меняются на матрицы. В моём алгоритме выражение следующее:
^{-1}(M_1h_1+M_2h_2+\dots+M_nh_n))
Здесь hn - вектора 2×1, координаты точки на фотоприёмной матрице, Mn - матрицы 6×2, выражающие "чувствительность" каждой координаты каждой точки к 6 параметрам сближения (3 угла и 3 смещения). Как можно заметить, правая скобка становится вектором 6×1, левая скобка становится симметричной матрицей 6×6, при обращении она, разумеется, сохраняет свой размер и симметричность. И наконец, умножая матрицу 6×6 на вектор 6×1, снова получаем вектор 6×1 - это наши параметры сближения (точнее, их коррекция к предыдущей итерации, т.к изначально задача нелинейная, но мы её линеаризовали вблизи предыдущего приближения).
По сути, ТО ЖЕ САМОЕ.
Так вот, с точки зрения целочисленных вычислений, это полный песец!
Точнее так: полный песец наступает, если чувствительности могут колебаться в широком диапазоне. Тогда динамический диапазон, необходимый при решении этой задачи, ещё примерно УТРАИВАЕТСЯ ПО ЧИСЛУ ПОРЯДКОВ.
Наибольшую боль мне доставляет измерение дальности. Дальность может меняться от 0,5 до 300 метров. Измеряем мы, по сути, угловой размер фиксированной базы, т.е где-то на ответной части стоит два уголковых отражателя, между которыми фиксированное расстояние D метров. А мы измеряем угол φ≈D/L, где L - дальность. "Чувствительность" к изменению дальности данного угла равна производной этого выражения по L:
Квадрат в знаменателе - это БОЛЬ. Он означает, что чувствительность меняется в 6002=360 000 раз! Это уже монструозно, да просто чтобы это число представить, требуется 19 бит (тогда можно надеяться, что чувствительность будет меняться от целого значения 1 до значения 360 000, хотя понятно, какова погрешность в этой единице, когда плюс-минус 100% и никак иначе её отрегулировать нельзя!)
Но если вычислять это выражение в лоб, это будет означать: сначала мы должны умножить на этот коэффициент (который может скакать в 360 000 раз), а потом мы ещё обязаны будем поделить на ЕГО КВАДРАТ! А это и вовсе разброс в 130 трлн раз (!!!) Чтобы это число записать, нужно 37 бит, т.е тут уже и 32-битные значения нас не устраивают...
Поэтому о вычислениях "в лоб" никогда и мысли не было... Довольно быстро решено было искомый вектор параллельного переноса (как изделие расположено относительно оптической мишени) представлять "с общей экспонентой", т.е есть 16-битные значения x,y,z (первое беззнаковое, от 0 до 2, остальные знаковые, от -1 до +1, считаем, что ровно один разряд перед запятой) домножаются на 2E, где E принимает значения от -1 до 8. Ещё точнее, домножаться они должны на 2E-1, чтобы E принимало более удобные значения от 0 до 9.
Значение чувствительности по горизонтали на изменение дальности (т.е значения x вектора параллельного переноса) записывается так:
^2})
Взяв вектор (x,y,z) с общей экспонентой E, мы первым делом приводили к той же экспоненте точки (bnx,bny,bnz) - координаты элементов оптической мишени, в метрах. Как только мы это делали, в выражении для mn14 экспонента выносится за скобки, частично сокращается, и остаётся значение 2(-(E-1)). Но это если мы хотим найти чувствительность в "пикселях на метр". На больших дальностях эта чувствительность будет чрезвычайно малой. Как мы знаем по размерности, "входные значения" (невязки в пикселях) должны в итоге поделиться на эти чувствительности, чтобы получить метры. Поэтому экспонента снова выползет вверх, давая 2E-1. Это хорошо: если просто эту экспоненту выкинуть, мы будем попросту искать мантиссу (x,y,z) при той же экспоненте E. Потом, игнорируя экспоненты, мы просто скорректируем мантиссу (x,y,z), и только если она сильно при этом поменяет масштаб - подкорректируем его на одну ступень, то ли в плюс, то ли в минус. Так-то мы знаем, что шажки маленькие, движения все неторопливые, сразу на несколько ступеней прыгнуть невозможно.
Теперь всё хорошо?
Ситуация, безусловно, улучшилась. Теперь у нас в знаменателе будет получаться, грубо говоря, от 1 до 4 (когда x принимает значения от 1 до 2). В числителе у нас есть координата точки мишени, порядка метра, но гораздо большее значение может принимать y. Поле зрения прибора: 24 градусов по горизонтали и вертикали, значит, если мишень оказалась в самом краю поля зрения (12 градусов от оптической оси), получим:
Так что максимальное значение mn14 будет около 0,21, или 6965, если использовать формат Q1.15. Когда мы это число возведём в квадрат "для знаменателя", получится 0,045, или 1480. В знаменателе эти квадраты будут складываться, мы ожидаем 8 примерно одинаковых значений, что даст 11843, при максимальном значении 32767 для знакового числа или 65535 для беззнакового. Маловато, но выглядит адекватно.
Но вот в чём проблема: нас больше интересует ситуация y=0, и тогда в числителе остаётся только координата точки мишени, её величина: порядка метра. Когда мы работаем с дальности 300 метров, координаты мишени приводятся к масштабу вектора параллельного переноса, т.е делятся на 256. Получаем число 1/256 в числителе, и (300/256)2 в знаменателе, а всё вместе: 0,0028, или 93 в формате Q1.15. Если теперь это число возвести в квадрат, получается 8·10-6 или "0,26" в формате Q1.15, которое округлится до нуля. У нас произошла полная потеря точности!
Простейшее решение, как хоть что-то выцарапать - повысить масштаб значения mn14. Мы видели, что даже когда мы складывали 8 квадратов, у нас "по наихудшему случаю" получалось 11843, есть куда расти.
Именно тогда мне очень захотелось ввести анизотропию, т.е поперечные компоненты умножить на 2. В смысле, снизить цену деления вдвое, поэтому те же точки будут представлены увеличенными вдвое координатами. В итоге, без грубого вмешательства ручками у нас mn14 удваивается, а его квадрат - учетверяется. Так что максимальное значение на диагоналях матрицы 6×6 (на которых суммируются квадраты чувствительностей по горизонтали и вертикали по всем обнаруженным точкам) - порядка 47375. Как видно, нужно для диагональных элементов симметричной матрицы использовать беззнаковые числа, но пока всё получается.
Теперь при работе с 300 метров строго по центру, мы получим квадрат, равный 1,06, что округлится до 1. На бровях избежали полной потери точности... Практика показала: из-за ошибок вычислений, которых здесь хватает, округление получилось до нуля... А это всё, матрица выходит с нулевым определителем, уравнение не имеет решения.
Понятно, можно тут ещё костылей понаставить - при маленьких y дополнительно смасштабировать этот коэффициент, а при больших y этого не делать, чтобы не вызвать переполнения. Но тогда вдруг оказывается, что мы всё ближе и ближе к использованию плавающей точки, но вместо того, чтобы поручить вычисления процессору, мы что-то такое придумываем "ручками", рискуя где-нибудь выйти за диапазон, ошибиться на степени двойки и так далее. Наверное, это осуществимо, но лично мне кажется уродством...
Обидно, имея полную свободу действий, в итоге приходить вот к такому! Вот потому и думаю перейти с 16 бит на что-нибудь покрупнее, то ли 32 бита фиксированных, либо 16 бит мантиссы + экспонента. К сожалению, в нынешней реализации QuatCore делать такой код - сущее мучение. Это, в конце концов, TTA - Transport Triggered Architecture, т.е мы говорим, откуда взять 16 бит на шине, и куда их в итоге запихнуть. Стоит размер чисел увеличить - это надо будет удваивать все-все операции, да и адресация вся летит псу под хвост. Просто процессор переделать в 32-битный - тогда он в ПЛИС 5576ХС6Т уже не влезет, причём не по ЛЭ, и даже не по памяти, а по интерконнекторам. Не знаю точно, почему, но у него не получается прошвырнуть через весь кристалл полноценную 32-битную шину, он рушится именно на этапе Place&Route. Т.е сам по себе Analysis&Synthesis завершается успешно, не так уж много ЛЭ, ну допустим 10% от кристалла, но объединить их в общую схему не выходит. Грубо говоря, там всего "два слоя дорожек", и их не хватает.
Впрочем, появилась ещё одна идейка, как это всё решить малой кровью, оставшись в честных 16 битах. Про неё напишу завтра. Снова о пользе геометрии - пока чисто алгебраически размышлял о 6 степенях свободы - в упор не видел. А нарисовал примитивный частный случай - 2 степени свободы - тут же осенило.
Всё дело в том, что алгоритм сопровождения, основанный на решении нелинейной задачи условной оптимизации, упорно не хочет вписываться в 16-битные целые числа.
Сегодня, чтобы на интуитивном уровне понять, что с этим делать, решил рассмотреть простейшую "игрушечную задачу", снизив число степеней свободы с шести ДО ОДНОЙ, и число точек ДО ДВУХ. Как ни странно, даже тогда ФУНДАМЕНТАЛЬНАЯ ПРОБЛЕМА сохраняется, хотя здесь с ней ещё как-то можно справиться...
Итак, у нас есть неизвестное число x, которое мы хотим измерить. У нас есть два зашумлённых входа h1, h2 которые реагируют на изменение числа x, но чувствительность различается:
Здесь ξ η - случайные величины с нулевым математическим ожиданием (мы считаем, что всё хорошо откалибровали и всю систематику убрали) и дисперсией, для простоты, равной ЕДИНИЦЕ.
К примеру, h1, h2 измеряются в вольтах, а величина x - в метрах. Тогда "чувствительности" k1, k2 имеют размерность В/м.
Задача: как лучше всего на основании h1, h2 оценить x, т.е найти
, чтобы результат был несмещённым, а дисперсия - как можно ниже?
Исходя из размерности, вроде как "чувствительности" должны оказаться в знаменателе, чтобы вольты снова превратить в метры. По крайней мере, если бы вход был лишь один, h1, то нам бы ничего другого не оставалось, как записать
Но стоит начать оптимизацию - и дела немного усложняются...
Выводится всё очень легко. Ответ должен быть взвешенной суммой h1 и h2:
Подставим сюда выражения для h1, h2:
Чтобы выражение было несмещённым (математическое ожидание равнялось в точности x), скобочка должна обращаться в единицу:
У нас осталась одна степень свободы в выборе α, β. Поэтому выпишем выражение для дисперсии:
Но оптимизировать значения α, β надо с учётом ограничения выше, на несмещённость. Тут можно было бы просто β выразить через α, но это некрасиво, давайте уж воспользуемся множителем Лагранжа для условной оптимизации:
Чтобы найти экстремум, надо взять частные производные по α, β и λ и приравнять к нулю:
Удобно взять чуть другую, смасштабированную "лямбду" (внести в неё знак минуса и деление на два, чтобы лишнего не писать):
Теперь α и β, выраженные через λ' подставляем в третью строчку:
И наконец, находим α и β
Оценка x запишется так:
По размерности всё верно. То, что мы для начала УМНОЖАЕМ каждый вход на его чувствительность, вместо того, чтобы делить - имеет ясную интерпретацию, ЧЕМ ЧУВСТВИТЕЛЬНЕЕ ВХОД (при том же шуме, что и другие) - ТЕМ ОН ПОЛЕЗНЕЕ, поэтому именно он будет иметь бОльший вес! Можно рассмотреть вырожденные случаи:
1. h2 вообще не реагирует на измеряемую величину, т.е k2=0. Тогда он выкидывается из выражения, остаётся попросту h1/k1.
2. h2 СУЩЕСТВЕННО более чувствителен, т.е k2>>k1. Тогда, поделив числитель и знаменатель на его квадрат и вычеркнув "бесконечно малые", придём к выражению h2/k2, т.е не видим никакого смысла использовать вход h1.
Подобное выражение весьма универсально, в более сложных задачах оно остаётся практически неизменным, только числа меняются на матрицы. В моём алгоритме выражение следующее:
Здесь hn - вектора 2×1, координаты точки на фотоприёмной матрице, Mn - матрицы 6×2, выражающие "чувствительность" каждой координаты каждой точки к 6 параметрам сближения (3 угла и 3 смещения). Как можно заметить, правая скобка становится вектором 6×1, левая скобка становится симметричной матрицей 6×6, при обращении она, разумеется, сохраняет свой размер и симметричность. И наконец, умножая матрицу 6×6 на вектор 6×1, снова получаем вектор 6×1 - это наши параметры сближения (точнее, их коррекция к предыдущей итерации, т.к изначально задача нелинейная, но мы её линеаризовали вблизи предыдущего приближения).
По сути, ТО ЖЕ САМОЕ.
Так вот, с точки зрения целочисленных вычислений, это полный песец!
Точнее так: полный песец наступает, если чувствительности могут колебаться в широком диапазоне. Тогда динамический диапазон, необходимый при решении этой задачи, ещё примерно УТРАИВАЕТСЯ ПО ЧИСЛУ ПОРЯДКОВ.
Наибольшую боль мне доставляет измерение дальности. Дальность может меняться от 0,5 до 300 метров. Измеряем мы, по сути, угловой размер фиксированной базы, т.е где-то на ответной части стоит два уголковых отражателя, между которыми фиксированное расстояние D метров. А мы измеряем угол φ≈D/L, где L - дальность. "Чувствительность" к изменению дальности данного угла равна производной этого выражения по L:
Квадрат в знаменателе - это БОЛЬ. Он означает, что чувствительность меняется в 6002=360 000 раз! Это уже монструозно, да просто чтобы это число представить, требуется 19 бит (тогда можно надеяться, что чувствительность будет меняться от целого значения 1 до значения 360 000, хотя понятно, какова погрешность в этой единице, когда плюс-минус 100% и никак иначе её отрегулировать нельзя!)
Но если вычислять это выражение в лоб, это будет означать: сначала мы должны умножить на этот коэффициент (который может скакать в 360 000 раз), а потом мы ещё обязаны будем поделить на ЕГО КВАДРАТ! А это и вовсе разброс в 130 трлн раз (!!!) Чтобы это число записать, нужно 37 бит, т.е тут уже и 32-битные значения нас не устраивают...
Поэтому о вычислениях "в лоб" никогда и мысли не было... Довольно быстро решено было искомый вектор параллельного переноса (как изделие расположено относительно оптической мишени) представлять "с общей экспонентой", т.е есть 16-битные значения x,y,z (первое беззнаковое, от 0 до 2, остальные знаковые, от -1 до +1, считаем, что ровно один разряд перед запятой) домножаются на 2E, где E принимает значения от -1 до 8. Ещё точнее, домножаться они должны на 2E-1, чтобы E принимало более удобные значения от 0 до 9.
Значение чувствительности по горизонтали на изменение дальности (т.е значения x вектора параллельного переноса) записывается так:
Взяв вектор (x,y,z) с общей экспонентой E, мы первым делом приводили к той же экспоненте точки (bnx,bny,bnz) - координаты элементов оптической мишени, в метрах. Как только мы это делали, в выражении для mn14 экспонента выносится за скобки, частично сокращается, и остаётся значение 2(-(E-1)). Но это если мы хотим найти чувствительность в "пикселях на метр". На больших дальностях эта чувствительность будет чрезвычайно малой. Как мы знаем по размерности, "входные значения" (невязки в пикселях) должны в итоге поделиться на эти чувствительности, чтобы получить метры. Поэтому экспонента снова выползет вверх, давая 2E-1. Это хорошо: если просто эту экспоненту выкинуть, мы будем попросту искать мантиссу (x,y,z) при той же экспоненте E. Потом, игнорируя экспоненты, мы просто скорректируем мантиссу (x,y,z), и только если она сильно при этом поменяет масштаб - подкорректируем его на одну ступень, то ли в плюс, то ли в минус. Так-то мы знаем, что шажки маленькие, движения все неторопливые, сразу на несколько ступеней прыгнуть невозможно.
Теперь всё хорошо?
Ситуация, безусловно, улучшилась. Теперь у нас в знаменателе будет получаться, грубо говоря, от 1 до 4 (когда x принимает значения от 1 до 2). В числителе у нас есть координата точки мишени, порядка метра, но гораздо большее значение может принимать y. Поле зрения прибора: 24 градусов по горизонтали и вертикали, значит, если мишень оказалась в самом краю поля зрения (12 градусов от оптической оси), получим:
Так что максимальное значение mn14 будет около 0,21, или 6965, если использовать формат Q1.15. Когда мы это число возведём в квадрат "для знаменателя", получится 0,045, или 1480. В знаменателе эти квадраты будут складываться, мы ожидаем 8 примерно одинаковых значений, что даст 11843, при максимальном значении 32767 для знакового числа или 65535 для беззнакового. Маловато, но выглядит адекватно.
Но вот в чём проблема: нас больше интересует ситуация y=0, и тогда в числителе остаётся только координата точки мишени, её величина: порядка метра. Когда мы работаем с дальности 300 метров, координаты мишени приводятся к масштабу вектора параллельного переноса, т.е делятся на 256. Получаем число 1/256 в числителе, и (300/256)2 в знаменателе, а всё вместе: 0,0028, или 93 в формате Q1.15. Если теперь это число возвести в квадрат, получается 8·10-6 или "0,26" в формате Q1.15, которое округлится до нуля. У нас произошла полная потеря точности!
Простейшее решение, как хоть что-то выцарапать - повысить масштаб значения mn14. Мы видели, что даже когда мы складывали 8 квадратов, у нас "по наихудшему случаю" получалось 11843, есть куда расти.
Именно тогда мне очень захотелось ввести анизотропию, т.е поперечные компоненты умножить на 2. В смысле, снизить цену деления вдвое, поэтому те же точки будут представлены увеличенными вдвое координатами. В итоге, без грубого вмешательства ручками у нас mn14 удваивается, а его квадрат - учетверяется. Так что максимальное значение на диагоналях матрицы 6×6 (на которых суммируются квадраты чувствительностей по горизонтали и вертикали по всем обнаруженным точкам) - порядка 47375. Как видно, нужно для диагональных элементов симметричной матрицы использовать беззнаковые числа, но пока всё получается.
Теперь при работе с 300 метров строго по центру, мы получим квадрат, равный 1,06, что округлится до 1. На бровях избежали полной потери точности... Практика показала: из-за ошибок вычислений, которых здесь хватает, округление получилось до нуля... А это всё, матрица выходит с нулевым определителем, уравнение не имеет решения.
Понятно, можно тут ещё костылей понаставить - при маленьких y дополнительно смасштабировать этот коэффициент, а при больших y этого не делать, чтобы не вызвать переполнения. Но тогда вдруг оказывается, что мы всё ближе и ближе к использованию плавающей точки, но вместо того, чтобы поручить вычисления процессору, мы что-то такое придумываем "ручками", рискуя где-нибудь выйти за диапазон, ошибиться на степени двойки и так далее. Наверное, это осуществимо, но лично мне кажется уродством...
Обидно, имея полную свободу действий, в итоге приходить вот к такому! Вот потому и думаю перейти с 16 бит на что-нибудь покрупнее, то ли 32 бита фиксированных, либо 16 бит мантиссы + экспонента. К сожалению, в нынешней реализации QuatCore делать такой код - сущее мучение. Это, в конце концов, TTA - Transport Triggered Architecture, т.е мы говорим, откуда взять 16 бит на шине, и куда их в итоге запихнуть. Стоит размер чисел увеличить - это надо будет удваивать все-все операции, да и адресация вся летит псу под хвост. Просто процессор переделать в 32-битный - тогда он в ПЛИС 5576ХС6Т уже не влезет, причём не по ЛЭ, и даже не по памяти, а по интерконнекторам. Не знаю точно, почему, но у него не получается прошвырнуть через весь кристалл полноценную 32-битную шину, он рушится именно на этапе Place&Route. Т.е сам по себе Analysis&Synthesis завершается успешно, не так уж много ЛЭ, ну допустим 10% от кристалла, но объединить их в общую схему не выходит. Грубо говоря, там всего "два слоя дорожек", и их не хватает.
Впрочем, появилась ещё одна идейка, как это всё решить малой кровью, оставшись в честных 16 битах. Про неё напишу завтра. Снова о пользе геометрии - пока чисто алгебраически размышлял о 6 степенях свободы - в упор не видел. А нарисовал примитивный частный случай - 2 степени свободы - тут же осенило.