SVD-разложение && направление собственных векторов

[allchemist]

Доброй ночи.
Пытаюсь реализовать алгоритм SVD-разложения матрицы X на U, S и V:
U - матрица собственных векторов матрицы $XX^T$, V - матрица собственных векторов $X^TX$.
Проверяю правильность: $UU^T=I, VV^T=I$ - все отлично. $USV^T!=X$
Стал искать ошибку в процедуре нахождения собственных векторов. Сначала использовал метод Якоби, и он не дал правильного SVD-разложения, хотя работал правильно. Написал QR-разложение, убедился в правильности, нашел с его помощью собственные вектора, но и он не дал верного SVD-разложения.
Потом заметил закопанную собаку в том, что собственные вектора инвариантны к изменению направления (знака), т.к. равенство Ax=lx при этом не нарушается.
Два метода поиска собств. векторов дали один результат с точностью до знака некоторых векторов.
Забил ту же матрицу в scilab - и он нашел правильное разложение. Матрица S совпадала с предыдущими, U и V отличались некоторыми знаками.
Понятно, почему SVD-разложения, найденные двумя способами, отличались друг от друга и от правильного разложение.

Если кто может, объясните, пожалуйста,
* почему разные методы диагонализации дают разные направления собственных векторов?
* есть ли какое-то "правильное" направление?
* как в сложившихся условиях найти правильное SVD-разложение?

Спасибо!

[ewert]

allchemist в сообщении #242545 писал(а):

есть ли какое-то "правильное" направление?

Конечно, нет. Однако непонятно, как конкретно Вы искали собственные векторы $\vec u_k$ для матрицы $X\,X^T$ и $\vec v_k$ для матрицы $X^TX$. Эти наборы ведь должны быть не абы какими, а согласованными друг с другом: $A\vec v_k=s_k\vec u_k$ (и, так сказать, наоборот), где $s_k$ -- соответствующее сингулярное число. Т.е. находить надо только матрицу $V^T$ (из столбцов $\vec v_k$), а матрицу $U$ вычислять уже по ней (или наоборот). Если же вычислять эти матрицы независимо друг от друга -- естественно, на согласованность надеяться не приходится.

[Nik_Svan]

allchemist в сообщении #242545 писал(а):

...Пытаюсь реализовать алгоритм SVD-разложения матрицы X ...

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

[ewert]

Nik_Svan в сообщении #242581 писал(а):

Сингулярный анализ - дело тонкое.

А что в нём такого тонкого (если, конечно, хотя бы одна из матриц -- $X\,X^T$ или $X^TX$ -- невырожденна и не плохо обусловлена, но это общая проблема для любых задач)?

[allchemist]

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

Цитата:

$A\vec v_k=s_k \vec u_k$

До этого я не догадался, спасибо.

Пишу не лиспе, кстати

[Nik_Svan]

ewert в сообщении #242585 писал(а):

Nik_Svan в сообщении #242581 писал(а):

Сингулярный анализ - дело тонкое.

А что в нём такого тонкого (если, конечно, хотя бы одна из матриц -- $X\,X^T$ или $X^TX$ -- невырожденна и не плохо обусловлена, но это общая проблема для любых задач)?

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

[allchemist]

$\mathbb{V}^T=\mathbb{S}^{-1}\mathbb{U}^T\mathbb{A}$!
Как же я до этой очевидной вещи сам не додумался...

Кстати, а что имеется в виду под грамотной программной реализацией?
Устойчивость? Работа с плохо обусловленными матрицами?

Если считать, что svd сводится к диагонализации $\mathbb{XX}^T$ с последующим вычислением $\mathbb{V}^T=\mathbb{S}^{-1}\mathbb{U}^T\mathbb{A}$, то грамотность реализации svd зависит от способа нахождения спектра действительной квадратной симметричной матрицы, а эту задачу мы решать умеем (например, методом Хаусхолдера)

Или Вы имеете ввиду нахождение svd через преобразование Карунена-Лоева (pca)?

[ewert]

allchemist в сообщении #242798 писал(а):

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

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

Впрочем, подобных проблем не возникнет, если:

1) матрица хорошо обусловлена; или:

2) матрица пусть и вырождена, но её ненулевые сингулярные числа не безумно сильно отличаются друг от друга.

А что имел в виду Nik_Svan под "грамотной реализацией" -- я не знаю. Проблема-то имеет объективный характер. Для плохо обусловленных матриц задача численно неустойчива, и тут уж никуда не денешься, как ни крути. Единственно -- можно привлечь какую-то дополнительную информацию.

[allchemist]

Цитата:

Единственно -- можно привлечь какую-то дополнительную информацию.

Я хотел использовать svd как еще один способ реализации метода главных компонент. К счастью, есть другой замечательный способ реализации pca - с помощью нейронной сети, обучаемый по обобщенному правилу Хебба (GHA). Он может найти собственные числа и вектора матрицы корелляции $\mathbb{XX}^T$ без ее построения, сразу по $\mathbb{X}$.

Цитата:

Это-то мы умеем, но там проблема -- что делать, если матрица вырождена или плохо обусловлена

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

[MrD]

Вброшу свои две копейки.
Пару лет назад успешно пользовался этим:
http://www.uic.nnov.ru/~zny/nl/doc/html/
- это если нужно "ехать" а не "шашечки". Я в детали не вдавался, мне важен был результат.

Впрочем, там и теория неплохо изложена, в частности, есть вот такой документ:
http://www.uic.nnov.ru/~zny/nl/BookSomePages.pdf

Исходники здесь:
http://code.google.com/p/nl-software/downloads/list

p.s. в LAPACK'е это все тоже есть, но ковыряться в фортрановских текстах - это не для ленивых.