Дуальные числа рулят!
Наверное, все знают о комплексных числах с их "мнимой единицей" 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.
Всё как обычно.
Умножение мы можем вывести:
\cdot(b_0+b_1s)=a_0b_0+a_1b_0s+a_0b_1s+\cancel{a_1b_1s^2}=)
s)
Оно похоже на умножение комплексных чисел, но одним слагаемым меньше. И можно заметить важный факт: дуальная ("мнимая") часть выражения вообще не оказывает никакого влияния на действительную. Если в комплексных числах они могли переходить одна в другую при умножении, то здесь никакого возврата нет!
И мы можем вспомнить ситуацию, когда происходит похожее: работа с "бесконечно малыми", например с dx. Добавление бесконечно малой величины не меняет результат как таковой (только другой бесконечно малой добавкой), а при возведении dx в квадрат мы его обычно выкидываем, поскольку интересуемся только "первым порядком".
Каждому дуальному числу можно приписать сопряжённое:
Но если умножение комплексного числа на сопряжённое, равно как умножение кватерниона на сопряжённый давали нам сумму квадратов компонентов (квадрат длины вектора), то здесь ситуация другая:
(a_0-a_1s)=a_0^2-\cancel{a_1^2s^2}=a_0^2)
И снова, дуальная часть по сути "игнорируется".
Чтобы поделить одно дуальное число на другое, мы помножаем числитель и знаменатель на число, сопряжённое к знаменателю:
(b_0-b_1s)}{b_0^2}=\frac{a_0b_0+(a_1b_0-a_0b_1)s}{b_0^2}=)
Как и раньше, действительная часть ровно такая же, как если бы дуальной части и вовсе не было бы. Дуальная часть же подозрительно напоминает выражение для производной от частного: производная числителя умножить на знаменатель, минус числитель умноженный на производную знаменателя, и поделить на знаменатель в квадрате. Разумеется, это неспроста.
С дуальными числами можно проводить и более сложные операции, например, взять экспоненту:
}=1+s+\cancel{\frac{s^2}{2}}+\cancel{\frac{s^3}{6}}+\cancel{\frac{s^4}{24}}+\dots=1+s)
Зная, чему равна экспонента от "дуальной единицы", можно найти экспоненту от любого дуального числа:
}=\exp{(a_0)}\exp{(a_1s)}=\exp{(a_0)}(1+a_1s)=)
}+a_1\exp{(a_0)}s)
Примерно такой же трюк можно провернуть с синусом и косинусом:
Откуда:
}=\sin{a_0}\cos{a_1s}+\cos{a_0}\sin{a_1s}=\sin{a_0}+(a_1\cos{a_0})s\,,)
}=\cos{a_0}\cos{a_1s}-\sin{a_0}\sin{a_1s}=\cos{a_0}-(a_1\sin{a_0})s\,.)
Снова мы отчётливо видим, как у нас "вылезают" формулы для взятия производных.
По сути, если бы мы вместо "дуальной единицы" 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, выражающий координаты точки на матрице, если все неизвестные равны нулю:
где
b_{nz}}{(b_{nx}+t_{0x})^2}\,,)
b_{ny}}{b_{nx}+t_{0x}})\,,)
^2}\,,)
b_{nz}}{b_{nx}+t_{0x}})\,,)
b_{ny}}{(b_{nx}+t_{0x})^2}\,,)
^2}\,,)
Дальше уже дело техники. Находим симметричную положительно определённую матрицу 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!
По сути, в процессе выполнения у нас и получатся ровно те формулы, которые мы выводили "ручками", но только компьютер точно не ошибётся при их "выводе" :)
Разумеется, мы могли бы оптимизировать вычисления - придумать, как задействовать повторяющиеся части (например, сразу найти величину, обратную к знаменателю, и домножать на неё). Небольшую оптимизацию мы уже сделали, заметив, что два элемента нулевых, а ещё один повторяется на двух строках.
Но конкретно в моей задаче я куда больше переживаю за объём памяти, который займёт код, и лёгкость его отладки, да и в современных компьютерах, как мне кажется, такой подход через дуальные числа, даже требуя чуть больше вычислений, может оказаться быстрее и экономичнее, поскольку не так страшны вычисления, как доступ к оперативной памяти.
И даже если мы решим, что считать таким способом производные "в рантайме" - слишком расточительно, всё равно приятно иметь "отладочную версию" программы, где мы сделаем это попросту для подтверждения, что не ошиблись, когда выводили производные "ручками". Конечно, при работе с плавающей точкой можно попробовать дать малое приращение и затем разность поделить на него, но эта операция остаётся численно неустойчивой - нужно раз за разом подбирать правильное приращение, при котором не будет катастрофического падения точности, но которое относительно "локальное", чтобы хорошо аппроксимировать касательную. А здесь можно о таком не задумываться - ответ будет посчитан В ТОЧНОСТИ по формулам для нахождения производных. И всяко это проще, чем припахивать символьные вычисления.
Кроме дуальных чисел, существуют дуальные комплексные числа и дуальные кватернионы, которые позволяют единообразно описывать произвольное движение (вращение и параллельный перенос) на плоскости и в пространстве, соответственно. Есть надежда, что и до них мы в конце концов доберёмся :)
Ваш покорный слуга немножечко рассказал и о кватернионах (и ещё расскажет, когда немного освободится) с тремя "мнимыми единицами", которые очень здорово описывают вращение в пространстве и, несмотря на забвение почти в один век, сейчас они вовсю используются в космической технике, в компьютерной графике (в т.ч в играх) и робототехнике.
Но этим разные интересные "системы чисел" не ограничиваются, и сегодня немного расскажу о дуальных числах. Они внешне похожи на комплексные, только "мнимую единицу" в них обычно называют 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!
По сути, в процессе выполнения у нас и получатся ровно те формулы, которые мы выводили "ручками", но только компьютер точно не ошибётся при их "выводе" :)
Разумеется, мы могли бы оптимизировать вычисления - придумать, как задействовать повторяющиеся части (например, сразу найти величину, обратную к знаменателю, и домножать на неё). Небольшую оптимизацию мы уже сделали, заметив, что два элемента нулевых, а ещё один повторяется на двух строках.
Но конкретно в моей задаче я куда больше переживаю за объём памяти, который займёт код, и лёгкость его отладки, да и в современных компьютерах, как мне кажется, такой подход через дуальные числа, даже требуя чуть больше вычислений, может оказаться быстрее и экономичнее, поскольку не так страшны вычисления, как доступ к оперативной памяти.
И даже если мы решим, что считать таким способом производные "в рантайме" - слишком расточительно, всё равно приятно иметь "отладочную версию" программы, где мы сделаем это попросту для подтверждения, что не ошиблись, когда выводили производные "ручками". Конечно, при работе с плавающей точкой можно попробовать дать малое приращение и затем разность поделить на него, но эта операция остаётся численно неустойчивой - нужно раз за разом подбирать правильное приращение, при котором не будет катастрофического падения точности, но которое относительно "локальное", чтобы хорошо аппроксимировать касательную. А здесь можно о таком не задумываться - ответ будет посчитан В ТОЧНОСТИ по формулам для нахождения производных. И всяко это проще, чем припахивать символьные вычисления.
Кроме дуальных чисел, существуют дуальные комплексные числа и дуальные кватернионы, которые позволяют единообразно описывать произвольное движение (вращение и параллельный перенос) на плоскости и в пространстве, соответственно. Есть надежда, что и до них мы в конце концов доберёмся :)