Ликбез по кватернионам, Мортари-Маркли берут реванш!
Оглавление "Ликбеза по кватернионам":
[Spoiler (click to open)]
Часть 1 - история вопроса
Часть 2 - основные операции
Часть 3 - запись вращения через кватернионы
Часть 4 - кватернионы и спиноры; порядок перемножения
Часть 5 - практическая реализация поворота
Часть 5 1/2 - введение метрики, "расстояния" между поворотами
часть 5 5/8 - метрика ненормированных кватернионов
часть 5 11/16 - красивая псевдометрика произвольных кватернионов
Часть 5 3/4 - исследуем "пространство поворотов"
Часть 5 7/8 - почти изотропный ёжик
Часть 6 - поворот по кратчайшему пути
Часть 6 1/4 - кратчайший поворот в общем случае
Часть 6 2/4 - поворот, совмещающий два направления
Часть 6 3/4 - кватернион из синуса и косинуса угла
Часть 6 7/8 - "уполовинивание угла" на плоскости
Часть 7 - интегрирование угловых скоростей, углы Эйлера-Крылова
Часть 8 - интегрирование угловых скоростей, матрицы поворота
Часть 8 1/2 - ортонормирование матрицы и уравнения Пуассона
Часть 9 - интегрирование угловых скоростей с помощью кватернионов
Часть 10 - интегрирование угловых скоростей, методы 2-го порядка
Часть 10 1/2 - интегрирование с поддержанием нормы
Часть 11 - интегрирование угловых скоростей, методы высших порядков (в разработке)
Часть 12 - навигационная задача
Часть 13 - Дэвенпорт берёт след!
Часть 14 - линейный метод Мортари-Маркли
Часть 15 - среднее от двух кватернионов
Часть 15 1/2 - проверка и усреднение кватернионов
Часть 16 - разложение кватерниона на повороты
Часть 17 - лидарная задача
Задачка к части 17
Дэвенпорт VS Мортари-Маркли
Мортари-Маркли берут реванш!
Дэвенпорт VS Мортари-Маркли, раунд 3
В пятницу мы попробовали найти поворот из зашумлённых векторов, используя два разных метода - более классический Дэвенпорта и относительно новый (2007 года) Мортари-Маркли. Первый, будучи куда более сложным в реализации (особенно на очень простом "космическом" железе, умеющем складывать и умножать в целых числах, и ничего более), вышел на 11% точнее.
Но это была всего лишь одна выборка случайного шума, по одной выборке делать результаты - совершенно недопустимо! Поэтому решил всё-таки собрать статистику чуточку побольше. Но тут уже нельзя опираться на Wolfram Alpha, который найдёт все собственные значения и собственные вектора - надо научиться делать это самостоятельно, "малой кровью". К счастью, для данной конкретной задачи есть несколько интересных наблюдений, существенно упрощающих нам жизнь!
Степенной метод
Это очень простой итеративный метод поиска собственного вектора, соответствующего максимальному по модулю собственному значению - ровно то, что нам надо!
Особенно в ситуации, когда нам хочется сравнить между собой два метода, т.е кватернион, найденный по методу Мортари-Маркли - очень хорошо сгодится в качестве нулевого приближения.
А дальше, умножаем его на нашу матрицу K (матрицу Дэвенпорта), нормируем на единицу - и получаем более точное значение собственного вектора - вот, собственно, и всё!
Увы, за эту простоту приходится платить очень медленной сходимостью. Примерно так:
(0,70713; 0,00311; 0,00411; -0,70706),
(0,70711; 0,00297; 0,00402; -0,70708),
(0,70713; 0,00292; 0,00392; -0,70707),
(0,70712; 0,00283; 0,00388; -0,70708),
(0,70713; 0,00281; 0,00381; -0,70707)
тогда как прийти оно должно к ответу:
(0,70712; 0,00267; 0,0037; -0,70707)
Да, ползёт безусловно в нужную сторону, но как-то по-черепашьи! Не очень...
Для "хороших" задач собственное значение уже известно!
У нас в прошлый раз получились следующие собственные значения и собственные вектора:
Так вот, если мы сложим произведения длин исходных векторов (обычных, трёхмерных, "эталонных" и "измеренных"):
то получим значение 109666,5895. То есть, оно отличается в бОльшую сторону от нужного нам собственного значения всего на 0,024% (!). Так и будет всегда происходить, когда измеренные вектора не взяты совсем уж "с потолка", т.е вектор ошибки (разности между истинным и измеренным) у нас на порядки меньше самих векторов.
Именно это наблюдение в своё время стало ключевым для разработки алгоритма QUEST (QUaternion ESTimator). Получалось, что вовсе не нужно мучительно составлять характеристическое уравнение для матрицы 4х4 (там выкладки довольно громоздкие), потом его решать через кубическое, прибегая к куче квадратных и кубических корней и, возможно, комплексным числам для промежуточных вычислений (которые должны сократиться в конце). А можно либо просто взять довольно точное "для всех практических целей" приближение, либо, для очистки совести, ещё провести одну итерацию алгоритма Ньютона - и получить совсем шикарный результат. Не сколько для более точного нахождения поворота, а как раз для оценки "невязки" - её эта разница между λ0 и λmax характеризует очень хорошо.
Впрочем, у QUEST векторы полагались единичными, но каждому вектору приписывался вес ωk, чтобы некоторые направления, измеренные более точно (скажем, более яркая звезда, позволяющая с меньшим шумом находить яркостный центр) сильнее влияли на результат. И там "нулевое приближение" для собственного значения было скорее таким:
А мы, наоборот, сейчас все веса задали единичными, а вот вектора единичными не являются в лидарной задаче.
Очевидно, наиболее общее выражение таково:
Ну а зная собственное значение, найти собственный вектор уже не так трудно, это уже просто система линейных уравнений, притом одно уравнение "лишнее", так что вектор будет с точностью до множителя.
В нашем примере мы получим следующую матрицу K-λE:
)
"Невооружённым взглядом" видно, что первая и четвёртая строка практически одинаковые, т.е четвёртую строку можно получить с помощью линейной комбинации первых трёх (с большим весом 1-я, и ещё небольшая добавка второй и третьей, чтобы все коэффициенты совсем занулить). Если так, то можно постановить qz=-1 (четвёртый компонент собственного вектора, он же кватернион), и просто-напросто решить такую систему:
\left(\begin{array}{c}q_0\\q_x\\q_y\end{array}\right)=\left(\begin{array}{c}-100697\\-39\\-221\end{array}\right))
И отсюда получаем ненормированный кватернион
qненорм (1,00004; 0,00377; 0,00523; -1)
и после нормировки:
q (0,70711; 0,00266; 0,00369; -0,70708)
тогда как точный результат:
q (0,70712; 0,00267; 0,00370; -0,70707)
То есть, отличие в 5-м разряде после запятой, и не более чем на единичку. "Расстояние" между этими кватернионами: 7,5 угловых секунд (0,002°), для наших целей хватит.
Набираем статистику
Вот теперь, когда понятно, как можно посчитать кватернион по обоим методам, пора их сравнить чуть более научно.
Будем брать разные выборки случайного шума измерений и смотреть, какие ошибки даёт один и другой метод.
В зависимости от конкретной выборки шума, ошибка измерений скакала в широких пределах: от 0,008° до 0,776° для метода Мортари-Маркли и 0,06°..0,684° для метода Дэвенпорта (100 выборок). Среднее значение ошибки составило 0,282° для Дэвенпорта и 0,302° для Мортари-Маркли, т.е всё же на 6,5% выше.
Попробуем ещё "повысить точность" измерений с 1 мм до 0,1 мм (это уже похоже на реальную точность лазерного трекера). Снова делаем 100 выборок шума измерений. Получаем разброс ошибки от 0,002° до 0,0856° (Дэвенпорт) и от 0,006° до 0,0872° (Мортари-Маркли). Средняя ошибка составила 0,0359° для Дэвенпорта и 0,0370° для Мортари-Маркли, в этот раз на 3% выше.
Что ж, метод Мортари-Маркли частично оправдан: в нашем примере он в среднем даёт ошибку на 6,5% больше при точности измерений в 1 мм, и на 3% больше при точности измерений в 0,1 мм. Хотя статистики до сих пор маловато: мы рассмотрели только один поворот, на 90°, и ровно одну конкретную мишень.
Что интересно, и метод Дэвенпорта мы частично оправдали: оказалось, что и по нему можно находить кватернион "в пол-пинка", поскольку нужное нам собственное значение мы можем узнать по простой формуле через входные данные, с большой точностью. Собственно, так и начинается алгоритм QUEST (QUaternion ESTimator).
Наверное, пора немножко отвлечься от этой математики (есть не один, а целых два алгоритма, которые меня вполне устраивают!) и вернуться к макету. Скоро ехать к заказчику, надо как-то немножко лоску навести на всё это безобразие. Не знаю, правда, за что браться.
[Spoiler (click to open)]
Часть 1 - история вопроса
Часть 2 - основные операции
Часть 3 - запись вращения через кватернионы
Часть 4 - кватернионы и спиноры; порядок перемножения
Часть 5 - практическая реализация поворота
Часть 5 1/2 - введение метрики, "расстояния" между поворотами
часть 5 5/8 - метрика ненормированных кватернионов
часть 5 11/16 - красивая псевдометрика произвольных кватернионов
Часть 5 3/4 - исследуем "пространство поворотов"
Часть 5 7/8 - почти изотропный ёжик
Часть 6 - поворот по кратчайшему пути
Часть 6 1/4 - кратчайший поворот в общем случае
Часть 6 2/4 - поворот, совмещающий два направления
Часть 6 3/4 - кватернион из синуса и косинуса угла
Часть 6 7/8 - "уполовинивание угла" на плоскости
Часть 7 - интегрирование угловых скоростей, углы Эйлера-Крылова
Часть 8 - интегрирование угловых скоростей, матрицы поворота
Часть 8 1/2 - ортонормирование матрицы и уравнения Пуассона
Часть 9 - интегрирование угловых скоростей с помощью кватернионов
Часть 10 - интегрирование угловых скоростей, методы 2-го порядка
Часть 10 1/2 - интегрирование с поддержанием нормы
Часть 11 - интегрирование угловых скоростей, методы высших порядков (в разработке)
Часть 12 - навигационная задача
Часть 13 - Дэвенпорт берёт след!
Часть 14 - линейный метод Мортари-Маркли
Часть 15 - среднее от двух кватернионов
Часть 15 1/2 - проверка и усреднение кватернионов
Часть 16 - разложение кватерниона на повороты
Часть 17 - лидарная задача
Задачка к части 17
Дэвенпорт VS Мортари-Маркли
Мортари-Маркли берут реванш!
Дэвенпорт VS Мортари-Маркли, раунд 3
В пятницу мы попробовали найти поворот из зашумлённых векторов, используя два разных метода - более классический Дэвенпорта и относительно новый (2007 года) Мортари-Маркли. Первый, будучи куда более сложным в реализации (особенно на очень простом "космическом" железе, умеющем складывать и умножать в целых числах, и ничего более), вышел на 11% точнее.
Но это была всего лишь одна выборка случайного шума, по одной выборке делать результаты - совершенно недопустимо! Поэтому решил всё-таки собрать статистику чуточку побольше. Но тут уже нельзя опираться на Wolfram Alpha, который найдёт все собственные значения и собственные вектора - надо научиться делать это самостоятельно, "малой кровью". К счастью, для данной конкретной задачи есть несколько интересных наблюдений, существенно упрощающих нам жизнь!
Степенной метод
Это очень простой итеративный метод поиска собственного вектора, соответствующего максимальному по модулю собственному значению - ровно то, что нам надо!
Особенно в ситуации, когда нам хочется сравнить между собой два метода, т.е кватернион, найденный по методу Мортари-Маркли - очень хорошо сгодится в качестве нулевого приближения.
А дальше, умножаем его на нашу матрицу K (матрицу Дэвенпорта), нормируем на единицу - и получаем более точное значение собственного вектора - вот, собственно, и всё!
Увы, за эту простоту приходится платить очень медленной сходимостью. Примерно так:
(0,70713; 0,00311; 0,00411; -0,70706),
(0,70711; 0,00297; 0,00402; -0,70708),
(0,70713; 0,00292; 0,00392; -0,70707),
(0,70712; 0,00283; 0,00388; -0,70708),
(0,70713; 0,00281; 0,00381; -0,70707)
тогда как прийти оно должно к ответу:
(0,70712; 0,00267; 0,0037; -0,70707)
Да, ползёт безусловно в нужную сторону, но как-то по-черепашьи! Не очень...
Для "хороших" задач собственное значение уже известно!
У нас в прошлый раз получились следующие собственные значения и собственные вектора:
Так вот, если мы сложим произведения длин исходных векторов (обычных, трёхмерных, "эталонных" и "измеренных"):
то получим значение 109666,5895. То есть, оно отличается в бОльшую сторону от нужного нам собственного значения всего на 0,024% (!). Так и будет всегда происходить, когда измеренные вектора не взяты совсем уж "с потолка", т.е вектор ошибки (разности между истинным и измеренным) у нас на порядки меньше самих векторов.
Именно это наблюдение в своё время стало ключевым для разработки алгоритма QUEST (QUaternion ESTimator). Получалось, что вовсе не нужно мучительно составлять характеристическое уравнение для матрицы 4х4 (там выкладки довольно громоздкие), потом его решать через кубическое, прибегая к куче квадратных и кубических корней и, возможно, комплексным числам для промежуточных вычислений (которые должны сократиться в конце). А можно либо просто взять довольно точное "для всех практических целей" приближение, либо, для очистки совести, ещё провести одну итерацию алгоритма Ньютона - и получить совсем шикарный результат. Не сколько для более точного нахождения поворота, а как раз для оценки "невязки" - её эта разница между λ0 и λmax характеризует очень хорошо.
Впрочем, у QUEST векторы полагались единичными, но каждому вектору приписывался вес ωk, чтобы некоторые направления, измеренные более точно (скажем, более яркая звезда, позволяющая с меньшим шумом находить яркостный центр) сильнее влияли на результат. И там "нулевое приближение" для собственного значения было скорее таким:
А мы, наоборот, сейчас все веса задали единичными, а вот вектора единичными не являются в лидарной задаче.
Очевидно, наиболее общее выражение таково:
Ну а зная собственное значение, найти собственный вектор уже не так трудно, это уже просто система линейных уравнений, притом одно уравнение "лишнее", так что вектор будет с точностью до множителя.
В нашем примере мы получим следующую матрицу K-λE:
"Невооружённым взглядом" видно, что первая и четвёртая строка практически одинаковые, т.е четвёртую строку можно получить с помощью линейной комбинации первых трёх (с большим весом 1-я, и ещё небольшая добавка второй и третьей, чтобы все коэффициенты совсем занулить). Если так, то можно постановить qz=-1 (четвёртый компонент собственного вектора, он же кватернион), и просто-напросто решить такую систему:
И отсюда получаем ненормированный кватернион
qненорм (1,00004; 0,00377; 0,00523; -1)
и после нормировки:
q (0,70711; 0,00266; 0,00369; -0,70708)
тогда как точный результат:
q (0,70712; 0,00267; 0,00370; -0,70707)
То есть, отличие в 5-м разряде после запятой, и не более чем на единичку. "Расстояние" между этими кватернионами: 7,5 угловых секунд (0,002°), для наших целей хватит.
Набираем статистику
Вот теперь, когда понятно, как можно посчитать кватернион по обоим методам, пора их сравнить чуть более научно.
Будем брать разные выборки случайного шума измерений и смотреть, какие ошибки даёт один и другой метод.
В зависимости от конкретной выборки шума, ошибка измерений скакала в широких пределах: от 0,008° до 0,776° для метода Мортари-Маркли и 0,06°..0,684° для метода Дэвенпорта (100 выборок). Среднее значение ошибки составило 0,282° для Дэвенпорта и 0,302° для Мортари-Маркли, т.е всё же на 6,5% выше.
Попробуем ещё "повысить точность" измерений с 1 мм до 0,1 мм (это уже похоже на реальную точность лазерного трекера). Снова делаем 100 выборок шума измерений. Получаем разброс ошибки от 0,002° до 0,0856° (Дэвенпорт) и от 0,006° до 0,0872° (Мортари-Маркли). Средняя ошибка составила 0,0359° для Дэвенпорта и 0,0370° для Мортари-Маркли, в этот раз на 3% выше.
Что ж, метод Мортари-Маркли частично оправдан: в нашем примере он в среднем даёт ошибку на 6,5% больше при точности измерений в 1 мм, и на 3% больше при точности измерений в 0,1 мм. Хотя статистики до сих пор маловато: мы рассмотрели только один поворот, на 90°, и ровно одну конкретную мишень.
Что интересно, и метод Дэвенпорта мы частично оправдали: оказалось, что и по нему можно находить кватернион "в пол-пинка", поскольку нужное нам собственное значение мы можем узнать по простой формуле через входные данные, с большой точностью. Собственно, так и начинается алгоритм QUEST (QUaternion ESTimator).
Наверное, пора немножко отвлечься от этой математики (есть не один, а целых два алгоритма, которые меня вполне устраивают!) и вернуться к макету. Скоро ехать к заказчику, надо как-то немножко лоску навести на всё это безобразие. Не знаю, правда, за что браться.