Унитарная матрица преобразования в 3D

[ilghiz]

Добрый день,

пусть у меня имеется обычное трехмерное пространство и в нем задан набор точек $X = [x_1, \dots, x_k] \in {\mathbb R}^{3 \times k}, ~ 3 \le \le 4$, которые посредством матрицы трансформации и вектора линейного переноса переводятся в $Y = [y_1, \dots, y_k] \in {\mathbb R}^{3 \times k}$.

Запишем такую трансформацию в виде $A X + v u^T = Y$, $A \in {\mathbb R}^{3 \times 3}$, $v \in {\mathbb R}^3$, $u$ - вектор, состоящий из $k$ единиц.

В общем случае, чтобы найти $A$ и $v$ нам необходимо иметь хотя бы четыре точки ($k=4$), тогда, решая линейную систему

$$[X^T, u] [A, v]^T = Y^T$$

с тремя правыми частями мы получаем искомые $[A, v]^T$.

В моем же случае, есть два неприятных нюанса:

1. $X$ и $Y$ заданы с некоторой погрешностью.
2. мне необходимо получить $A$ унитарной, более того, без инверсий.

Понимая, что унитарная $A$ задана только тремя степенями свободы, и имея дополнительно три степени свободы от вектора перемещения $v$, хочется решить эту задачу для $k=3$, но, что-то подсказывает мне, что тут будут какие-то неоднозначности.

Также, если таки $k=4$ и $X$ и $Y$ заданы с некоторой погрешностью, хочется как-то это регуляризовать, чтобы затянуть $A$ в подкласс унитарных.

Скажите, пожалуйста,

1. есть ли такое решение для $k=3$, и, если есть, наставьте, пожалуйста, на путь истинный,
2. если нет, то как для $k>3$ правильнее регуляризовать с учетом унитарности $A$?

Спасибо!

[StaticZero]

Если использовать нормальную нотацию, то вы хотите узнать матрицу $A$ и вектор $\mathbf b$ такие, что $A \mathbf x_i + \mathbf b = \mathbf y_i$ для некоторого заданного набора $\mathbf x_i$, $\mathbf y_i$.

Преобразование общего вида, конечно, требует четырёх пар. Касательно остальных вопросов: надо фиксировать какой-то из векторов, скажем, $\mathbf x_1$ и перейти к набору векторов $\mathbf x'_i = \mathbf x_i - \mathbf x_1$, $i = 2, 3, 4$, и преобразвание станет

$$ A \mathbf x'_i + \mathbf b = A \mathbf x_i + \mathbf b - A \mathbf x_1 = \mathbf y_i - (\mathbf y_1 - \mathbf b) = \mathbf y'_i + \mathbf b, $$

так что для штрихованных векторов получаем $A \mathbf x'_i = \mathbf y'_i$ для $i = 2, 3, 4$. Три векторных уравнения определяют девять компонент матрицы $A$; мы успешно выгнали из уравнений трансляцию $\mathbf b$.

Теперь ваши вопросы.
1)

ilghiz в сообщении #1502600 писал(а):

1. есть ли такое решение [дающее $A \in \mathrm{SO}(3)$ -- добавлено мной] для $k=3$,

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

Если расширить класс искомых матриц до $\mathrm O(3)$, то есть разрешить $\det = -1$, то и тут не всё гладко: такое преобразование всё ещё будет сохранять длины, то есть векторы-результаты $\mathbf y'_i$ должны этому требованию соответствовать.

Возьмите ваши ладони в руки и сделайте из указательного и большого пальцев "пистолетики", оттопырив их (а остальные согнув). Поверните обе ладони кверху. Пусть большой палец на правой руке это вектор, скажем, $\mathbf x'_2$, указательный это $\mathbf x'_3$; на левой руке, соответственно, $\mathbf y'_2$ будет большой и $\mathbf y'_3$ указательный. Эти вы можете совместить обыкновенным поворотом вокруг одного их пальцев, но у вас возникнут проблемы, если вы будете пытаться большой палец отображать в указательный.

2)

ilghiz в сообщении #1502600 писал(а):

более того, без инверсий.

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

Этот пример иллюстрирует тот факт, что изучение всех возможных преобразований таких, которые конкретный вектор $\mathbf x_i$ переводят в конкретный вектор $\mathbf y_i$, может быть слишком ограничительным. Более содержательным может оказаться разговор о том, какие преобразования переводят $\mathbf x_i$ в какой-нибудь любой другой вектор $\mathbf y_j$.

3)

ilghiz в сообщении #1502600 писал(а):

если нет, то как для $k>3$ правильнее регуляризовать с учетом унитарности $A$?

МНК здесь может помочь, но сперва надо уточнить, чего именно вы хотите: преобразований "жёстких", то есть $\mathbf x_i \mapsto \mathbf y_i$, или "свободных" (преобразование точек как множества), $\mathbf x_i \mapsto \mathbf y_j$?

[arseniiv]

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

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

[ilghiz]

Спасибо большое, StaticZero, за ответ!

Не, мне только надо $x_i$ в $y_i$!

Позвольте уточнить и возразить, в первом случае, для $k=3$, я конечно же предполагал, что
$\forall i \ne j=1,...,3: ||x_i-x_j||_2 \simeq ||y_i - y_j||_2$, а, следовательно, такое преобразование с учетом унитарности $A$ найти можно.

Пусть $x_1,x_2,x_3$ - координаты вершин треугольника, и $y_1,y_2,y_3$ соответствующие координаты вершин такого же треугольника после преобразования. Очевидно, что трансляцию можно вытащить из $y_1-x_1$, а, далее нам потребуется две степени свободы для совмещения $x_2$ и $y_2$ и еще одна степень свободы для совмещения $x_3$ и $y_3$.

Вот будет ли получаемые $A$ и $b$ - единственными - мне не понятно, и, даже, если и будут, то будет ли численно устойчиво нахождение $A$ и $b$ в зависимости от малых возмущений по $x_1,x_2,x_3$ и $y_1,y_2,y_3$ - очень не очевидно.

Мой второй вопрос был как раз в том, на сколько малые возмущения по $x$ и $y$ могут изменить сами значения в $A$ и $b$. Мне понятно нужна унитарность по $A$, но я могу ее дополнительно требовать или нет. Если ее не требовать, то все сводится к решению линейной системы с матрицей $4 \times 4$ и тремя правыми частями, но если потребовать унитарность, то хоть и неизвестных становится не много и градиенты выписываются, решать такую нелинейность очень бы не хотелось бы.

В моем случае ошибка сидит в $x$ и $y$, следовательно наименьшие квадраты надо писать, если по-честному, только по их отклонениям.

----

Сейчас подумал и понял, что мне нужна вот такая задача:

у меня есть тройки точек $x_1,x_2,x_3 \in {\mathbb R}^3$ и $y_1,y_2,y_3 \in {\mathbb R}^3$, про которые известно, что $x_i \to y_i, \forall i=1,2,3$, и известно, что расстояния внутри каждой тройки более-менее сохраняются (численная погрешность конечно имеет место быть). Мне надобно найти для каждой такой пары троек матрицу трансформаций и вектор переноса причем так, чтобы значения в этой матрице и векторе были численно устойчивы.

То есть пусть у меня есть $A$ и $b$, и я подействовал на сотню ($N$) каких-то точек $X=[x_1,\dots,x_N]$ и получил соотвествующий набор $Y=[y_1,\dots\y_N]$. Я хочу, чтобы любая комбинация трех точек из $X$ и соответсвующая ей комбинация точек из $Y$ давала бы мне одни и те же $A$ и $b$.

С четырьмя точками - понятно, все будет работать, но хочется с тремя, ибо при большом $N$ все комбинации четырех точек и все комбинации из трех точек - это две очень большие разницы :(

-- 25.01.2021, 01:38 --

Спасибо большое, arseniiv, за ответ!

arseniiv в сообщении #1502613 писал(а):

Кажется, там какое-то выражение с псевдообратными должно быть (они связаны с МНК, но я забыл как), и тогда есть полная надежда его красиво, быстро и точно считать всякими библиотеками численной линейной алгебры.

да я довел это до конца, формулы не сильно страшные, но всяко три раза сингулярное разложение, и корни методом деления отрезка пополам, выкладки пока не осилил написать, на бумаге больше 3-х страниц получились. Мне как-то это показалось слишком громоздким, а хочется чего-то красивого. Если эти выкладки имеют смысл, завтра постораюсь их здесь записать.

[StaticZero]

ilghiz в сообщении #1502615 писал(а):

С четырьмя точками - понятно, все будет работать, но хочется с тремя, ибо при большом $N$ все комбинации четырех точек и все комбинации из трех точек - это две очень большие разницы :(

Ну, выгоните отсюда трансляцию, если вы хотите единый $\mathbf b$ для всех, будут вам тройки точек. Правда это всё равно $O(N^3)$. Может быть, можно лучше? Например, сделать пристрелку по выпуклой оболочке, а потом уточнить. Надо дальше формализовывать и думать.

[arseniiv]

Да, я навёл себе справки, и ну вот: если нам надо «решить» несовместную систему $A\mathbf v = \mathbf b$ или получить решение совместной, имеющее наименьший скалярный квадрат, то мы берём $\mathbf v^\circ = A^+ \mathbf b$, и аналогично для $A X = B$ нам нужно взять $X^\circ = A^+ B$, где $A^+$, конечно же, псевдообращение $A$.

И вот этот-то $X$ дальше и раскладывать полярно.

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

И в вашем случае $A$ это стопка иксов и $B$ — стопка игреков, а вот $A^+ B$ можно найти сразу одной функцией, ну смотря что вы используете. Например в NumPy* это функция lstsq (под капотом использующая функции ?gelsd из LAPACK, и наконец эти функции используют SVD), а вот полярного разложения в текущей версии NumPy нет, и придётся опять использовать SVD, уже в явную. А если вы используете другой язык и библиотеку, то может там уже реализовано. И очень надеюсь, что вы не пишете эти алгоритмы руками с нуля — вычислительная математика очень коварная штука с кучей подводных камней, и чтобы иметь уверенность в том, что реализация имеет все нужные свойства, надо обкладываться справочниками. (Или может быть специалисты в численных методах вас успокоят, но боюсь, что могут согласиться.)

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

[ilghiz]

Спасибо большое, StaticZero, что помогаете!

StaticZero в сообщении #1502618 писал(а):

Ну, выгоните отсюда трансляцию, если вы хотите единый $\mathbf b$ для всех, будут вам тройки точек.

так не устойчиво получается, я уже попробовал :(

Грубо говоря, если считать по четверкам, то матрицы и вектора преобразования более-менее хорошо получаются, а если единожды трансляцию выкинуть, то набегает ошибка. Можно правда трансляцию разную подставлять, но это все равно по арифметической сложности как если по четверкам считать.

[arseniiv]

Во, и совет про трансляцию: возьмите вместо исходных столбцов столбцы, дополненные единицей (в компграфике часто говорят «проектвное пространство», но это одна из небольших ошибок, и на самом деле там используется вложение аффинного пространства в линейное), и ищите матрицу порядком на один больше. Последний столбец матрицы даст вектор параллельного переноса, если последний элемент там будет тоже единицей (а последняя строка кроме этой единицы должна быть нулевой). Правда я сомневаюсь, что теперь мы можем использовать МНК с обычным скалярным произведением. Но попробуйте на каких-то тестах, вдруг работает. Если не очень работает, значит надо будет кое-что предпринять, пока лень думать как оно выразится.

-- Пн янв 25, 2021 05:45:35 --

Скорее всего не сработает, но чем чёрт не шутит.

[ilghiz]

Спасибо большое, arseniiv!

arseniiv в сообщении #1502620 писал(а):

И вот этот-то $X$ дальше и раскладывать полярно.

Да, верно, только тут-то как раз собака-то и порылась.

Итак, пусть у нас 3 точки, пусть мы трансляцию выгнали, то есть имеем:

$$A x_1 + b = y_1$$

$$A (x_2-x_1) = y_2-y_1$$

$$A (x_3-x_1) = y_3-y_1$$

Запишем для матрицы $[x_2-x_1, x_3-x_1]$ ее сингулярное разложение $U_xD_xV^T_x$, и перекинем направо правые сингулярные векторы и диагональ.
Тогда имеем $A U_x = H$, где $H=[y_2-y_1,y_3-y_1]V_x D_x^{-1}$.

Очевидно, что $A = H U_x^T + v w^T$, где $U_x^T v = 0$, а $w$ - произвольное неизвестное.

Чтобы потребовать унитарности $A$, потребуем

$$||AA^T - I||_2^2 = ... = ||T + x x^T||_2^2, T = H H^T - I,$$

по построению, $T$ - симметричная, то есть $T = VDV^T$, и поиск сводится к поиску неизвестного $z$, который минимизирует $||D - zz^T||_2^2$, с диагональной $D$, а это не имеет разумного решения, если у $D$ на диагонали больше одного не нулевого элемента.

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

arseniiv в сообщении #1502620 писал(а):

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

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

[arseniiv]

ilghiz в сообщении #1502622 писал(а):

так не устойчиво получается, я уже попробовал :(

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

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

-- Пн янв 25, 2021 05:55:57 --

Ага, вот попробуйте вычитать барицентр — вдруг стабильность улучшится. Я не могу быть уверен, потому что численные алгоритмы для меня довольно тёмная вещь, но будет приятно если догадки оправдаются хоть насколько-то.

-- Пн янв 25, 2021 05:56:36 --

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

[ilghiz]

arseniiv в сообщении #1502623 писал(а):

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

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

Хочется как раз все замутить с тройками.

[arseniiv]

ilghiz в сообщении #1502626 писал(а):

Так я в самом первом моем сообщении написал, что так и делал.

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

А барицентр отнимать попробуйте — если и он будет себя вести плохо, будет уже интересная история!

[StaticZero]

(Оффтоп)

arseniiv в сообщении #1502625 писал(а):

дисперсия в корень из числа точек раз больше дисперсии тех векторчиков

Или меньше?

[ilghiz]

Спасибо большое, StaticZero и arseniiv за советы!!!

Похоже с тройками тоже все делается, надо вначале вычислить $A$ по возможности максимально устойчиво. Если наложить ограничения, что исходные и конечные расстояния в треугольнике из трех точек равны, и углы в треугольнике не близки к 0 или 180, то можно гарантировать численно устойчивое построение $A$ в виде унитарной матрицы. А дальше, по совету arseniiv, использовать барицентр для поиска вектора трансформации $b$, тут вроде с устойчивостью тоже все в порядке получается, так как $A$ уже унитарная и любые средние как раз минимизируют евклидову норму.

Появился правда сопутствующий вопрос... Матрица трансформаций задана тремя углами, хотя вычисляется в виде девяти чисел. Мне далее ее надо сравнивать, и хотелось бы уменьшить размерность. Сразу видно, что можно матрицу пересчитать в углы и для каждого угла записать его синус и косинус. Тогда будет только шесть чисел. Углы сами сравнивать не удобно, так как при переходе от $-\pi$ к $+\pi$ реальное расстояние между углами может быть маленькое, а разница в $l_2$-норме далеко не такой маленькой.

Скажите, пожалуйста, а есть ли какое-то представление матрицы поворотов с не большим числом переменных $g$, чтобы любое расстояние векторов $g$ по $l_2$-норме не сильно отличалось от сравнения либо векторов, составленных из элементов самих матриц, или из пар синус-косинус?

[StaticZero]

ilghiz, я ничего не понял про норму, но поворот задаётся осью и углом, а ось задаётся единичным вектором вдоль неё. Всего 4 компонента. три, если вектор оси единичный

PS Чтобы вставить никнейм, надо нажать прямо на него: он от этого скопируется в форму ответа.