nabbla (nabbla1) wrote,
nabbla
nabbla1

Categories:

Дуальные числа рулят!

Наверное, все знают о комплексных числах с их "мнимой единицей" i (а у электронщиков: j, чтобы никто не догадался), "корнем из минус единицы", и о том, как здорово они "прижились" во всех областях физики.

Ваш покорный слуга немножечко рассказал и о кватернионах (и ещё расскажет, когда немного освободится) с тремя "мнимыми единицами", которые очень здорово описывают вращение в пространстве и, несмотря на забвение почти в один век, сейчас они вовсю используются в космической технике, в компьютерной графике (в т.ч в играх) и робототехнике.

Но этим разные интересные "системы чисел" не ограничиваются, и сегодня немного расскажу о дуальных числах. Они внешне похожи на комплексные, только "мнимую единицу" в них обычно называют s, а то и вовсе ε, и её свойство немножко другое: s2=0.

Расскажу не просто так, а применительно к задаче, которую они позволяют решить очень элегантно...


Итак, у нас есть дуальные числа A=a0+a1s и B=b0+b1s.

Мы можем их складывать и вычитать:
A+B = (a0+b0) + (a1+b1)s,
A-B = (a0-b0) + (a1-b1)s.

Всё как обычно.

Умножение мы можем вывести:



Оно похоже на умножение комплексных чисел, но одним слагаемым меньше. И можно заметить важный факт: дуальная ("мнимая") часть выражения вообще не оказывает никакого влияния на действительную. Если в комплексных числах они могли переходить одна в другую при умножении, то здесь никакого возврата нет!

И мы можем вспомнить ситуацию, когда происходит похожее: работа с "бесконечно малыми", например с dx. Добавление бесконечно малой величины не меняет результат как таковой (только другой бесконечно малой добавкой), а при возведении dx в квадрат мы его обычно выкидываем, поскольку интересуемся только "первым порядком".

Каждому дуальному числу можно приписать сопряжённое:


Но если умножение комплексного числа на сопряжённое, равно как умножение кватерниона на сопряжённый давали нам сумму квадратов компонентов (квадрат длины вектора), то здесь ситуация другая:


И снова, дуальная часть по сути "игнорируется".

Чтобы поделить одно дуальное число на другое, мы помножаем числитель и знаменатель на число, сопряжённое к знаменателю:



Как и раньше, действительная часть ровно такая же, как если бы дуальной части и вовсе не было бы. Дуальная часть же подозрительно напоминает выражение для производной от частного: производная числителя умножить на знаменатель, минус числитель умноженный на производную знаменателя, и поделить на знаменатель в квадрате. Разумеется, это неспроста.

С дуальными числами можно проводить и более сложные операции, например, взять экспоненту:


Зная, чему равна экспонента от "дуальной единицы", можно найти экспоненту от любого дуального числа:



Примерно такой же трюк можно провернуть с синусом и косинусом:



Откуда:



Снова мы отчётливо видим, как у нас "вылезают" формулы для взятия производных.
По сути, если бы мы вместо "дуальной единицы" s поставили бы dx, и стали бы точно так же преобразовывать выражения, то получили бы такие же результаты.

Дуальное число - это как бы совокупность значения выражения в некоторой точки (действительная часть) и производной от выражения в этой же точке (дуальная часть).

Возьмём как пример синус (см. выражение выше). Мы знаем значение аргумента (a0) и его производную (a1). Когда мы берём от него синус, мы получаем значение функции sin(a0) и её производной, которая (как и требуют правила взятия производной от сложной функции) равна производной аргумента, умноженной на производную самой функции, т.е a1cos(a0).

Собственно, это одно из главных применений дуальных чисел - автоматическое дифференцирование выражений!

Обратимся к конкретной задаче, которая у меня возникла.

Есть набор точек Bn (bnx, bny, nnz) - расположение точек некоторого объекта в пространстве. Но мы наблюдаем их повёрнутыми на неизвестный малый угол (φx, φy, φz), и ещё смещёнными в сторону на "большой", известный нам вектор (tx, ty, tz) и на малый, неизвестный вектор (Δtx, Δty, Δtz).

Иными словами, мы наблюдаем точки Cn:


Но и их мы не можем наблюдать непосредственно, как набор из трёх координат. Вместо этого, мы видим их изображение на фотоприёмной матрице, которое проецируется следующим образом:



В действительности мы наблюдаем точки (fny, fnz) и должны решить оптимизационную задачу:

найти такой угол φ и вектор Δt, чтобы dn как можно лучше совпадали с fn.

Под наилучшим совпадением можно считать обеспечение минимума суммы квадратов разностей. Кроме того, учитывая "малость" искомых величин, мы можем линеаризовать выражения для dny и dnz.

Приходим к задаче линейной безусловной оптимизации - вполне милая вещь, решаемая точно. Разумеется, мы начали с задачи нелинейной условной оптимизации, когда нужно найти кватернион поворота и вектора параллельного переноса (7 чисел, на 4 из которых наложено условие единичной нормы), при которых мы увидим ровно такую картинку, но уже примерно зная, чего ожидать с прошлого кадра, мы заранее повернули точки объекта на "большой" кватернион с прошлого шага (собственно поэтому наши точки называются Bn, исходные называются An) и добавили к ним "большой" параллельный перенос, оставив только малые поправки, которые и должны сейчас посчитать.

Нам необходимо посчитать матрицы частных производных, размером 2х6, и вектор gn, выражающий координаты точки на матрице, если все неизвестные равны нулю:

где















Дальше уже дело техники. Находим симметричную положительно определённую матрицу J (6х6):


(она положительно определённая, если задача имеет однозначное решение. Если мы, к примеру, попытаемся определить наше расстояние до звёзд по картинке звёздного неба, то матрица выйдет неотрицательно определённой, поскольку малое изменение расстояния никак не повлияет на изображение, пока речь идёт о тысячах и миллионах километров)

И также находим вектор K (6x1):



где

-
разность между измеренным положением точки на фотоприёмной матрице и предполагаемым положением.

И наконец, решаем уравнение:


(система из 6 линейных уравнений с 6 неизвестными)

- и золотой ключик у нас в кармане!

Как можно заметить, большинство действий записываются довольно компактно и столь же компактно реализуются в коде: что умножение матриц, что их сложение, транспонирование, решение системы линейных уравнений. Существует замечательный код на Фортране, реализующий всё это, он же портирован на Си. Зная о симметрии и положительной определённости матрицы J, можно не вычислять все 36 её коэффициентов, а только 21 (6*7/2) из них, а затем совершить LDL-разложение, и затем, выполнив две подстановки, получить ответ. Это всё тоже сто лет как реализовано.

Но вот конкретные выражения для матрицы Mn вызывают вселенскую тоску... Как думаешь о том, как будешь вбивать это всё на ассемблере, а потом ещё искать ошибки и исправлять их, и сколько же строчек кода это всё займёт - и становится страшно.

К счастью, можно этого избежать. Мы просто вычисляем dn в дуальных числах, положив, что Bn и t0 имеет нулевую дуальную часть, тогда как φ и Δt имеют нулевую действительную часть, а в дуальной части мы последовательно ставим единичку сначала в φx, на следующей итерации - в φy и так далее. И каждый раз мы попросту считаем значения dny и dnz.

Их действительная часть каждый раз будет неизменна - gny и gnz, тогда как дуальная часть даст нам последовательно все столбцы матрицы Mn!

По сути, в процессе выполнения у нас и получатся ровно те формулы, которые мы выводили "ручками", но только компьютер точно не ошибётся при их "выводе" :)

Разумеется, мы могли бы оптимизировать вычисления - придумать, как задействовать повторяющиеся части (например, сразу найти величину, обратную к знаменателю, и домножать на неё). Небольшую оптимизацию мы уже сделали, заметив, что два элемента нулевых, а ещё один повторяется на двух строках.

Но конкретно в моей задаче я куда больше переживаю за объём памяти, который займёт код, и лёгкость его отладки, да и в современных компьютерах, как мне кажется, такой подход через дуальные числа, даже требуя чуть больше вычислений, может оказаться быстрее и экономичнее, поскольку не так страшны вычисления, как доступ к оперативной памяти.

И даже если мы решим, что считать таким способом производные "в рантайме" - слишком расточительно, всё равно приятно иметь "отладочную версию" программы, где мы сделаем это попросту для подтверждения, что не ошиблись, когда выводили производные "ручками". Конечно, при работе с плавающей точкой можно попробовать дать малое приращение и затем разность поделить на него, но эта операция остаётся численно неустойчивой - нужно раз за разом подбирать правильное приращение, при котором не будет катастрофического падения точности, но которое относительно "локальное", чтобы хорошо аппроксимировать касательную. А здесь можно о таком не задумываться - ответ будет посчитан В ТОЧНОСТИ по формулам для нахождения производных. И всяко это проще, чем припахивать символьные вычисления.


Кроме дуальных чисел, существуют дуальные комплексные числа и дуальные кватернионы, которые позволяют единообразно описывать произвольное движение (вращение и параллельный перенос) на плоскости и в пространстве, соответственно. Есть надежда, что и до них мы в конце концов доберёмся :)
Tags: кватернионы-это просто (том 1), математика, работа, странные девайсы
Subscribe

Recent Posts from This Journal

  • Нахождение двух самых отдалённых точек

    Пока компьютер долго и упорно мучал симуляцию, я пытался написать на ассемблере алгоритм захвата на ближней дистанции. А сейчас на этом коде можно…

  • Слишком общительный счётчик

    Вчера я чуть поторопился отсинтезировать проект,параметры не поменял: RomWidth = 8 вместо 7, RamWidth = 9 вместо 8, и ещё EnableByteAccess=1, чтобы…

  • Балансируем конвейер QuatCore

    В пятницу у нас всё замечательно сработало на симуляции, первые 16 миллисекунд полёт нормальный. А вот прошить весь проект на ПЛИС и попробовать "в…

  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic

    Your IP address will be recorded 

  • 10 comments

Recent Posts from This Journal

  • Нахождение двух самых отдалённых точек

    Пока компьютер долго и упорно мучал симуляцию, я пытался написать на ассемблере алгоритм захвата на ближней дистанции. А сейчас на этом коде можно…

  • Слишком общительный счётчик

    Вчера я чуть поторопился отсинтезировать проект,параметры не поменял: RomWidth = 8 вместо 7, RamWidth = 9 вместо 8, и ещё EnableByteAccess=1, чтобы…

  • Балансируем конвейер QuatCore

    В пятницу у нас всё замечательно сработало на симуляции, первые 16 миллисекунд полёт нормальный. А вот прошить весь проект на ПЛИС и попробовать "в…