Почему мой расчет SVD отличается от расчета SVD этой матрицы в numpy?

Я пытаюсь вручную вычислить SVD матрицы A, определенной ниже, но у меня возникают некоторые проблемы. Вычисление вручную и с помощью метода svd в numpy дает два разных результата.

Вычислено вручную ниже:

import numpy as np
A = np.array([[3,2,2], [2,3,-2]])
V = np.linalg.eig(A.T @ A)[1]
U = np.linalg.eig(A @ A.T)[1]
S = np.c_[np.diag(np.sqrt(np.linalg.eig(A @ A.T)[0])), [0,0]]
print(A)
print(U @ S @ V.T)

И вычисляется с помощью метода svd numpy:

X,Y,Z = np.linalg.svd(A)
Y = np.c_[np.diag(Y), [0,0]]
print(A)
print(X @ Y @ Z)

Когда эти два кода запущены. Ручной расчет не равен методу svd. Почему существует расхождение между этими двумя расчетами?


python numpy svd
person Oliver G    schedule 24.09.2018    source источник
comment
Хм, я разобрал задачу поэлементно и обнаружил, что если сравнивать только X с U = np.linalg.eig(A @ A.T)[1], получается не та же матрица (знаки несколько другие). Даже Z и V отличаются в том смысле, что 2-я строка Z похожа на 3-ю строку V и наоборот с некоторыми противоположными знаками. Кажется, может быть проблема с numpy linalg. Проверьте второй ответ здесь   -  person Sheldore    schedule 24.09.2018
comment
Я вижу, что ряды и/или знаки либо перепутаны, но я не понимаю, почему это так. U @ S @ VT тоже не равно A, чего я не понимаю.   -  person Oliver G    schedule 24.09.2018
comment
Ответ заключается в том, что U и V изначально вычисляются неправильно. Для этого просто возьмите матрицу 2,2 и вычислите A с помощью какого-нибудь онлайн-инструмента, а затем сравните, даст ли вам eig linalg тот же ответ. Если нет, то с ним что-то не так   -  person Sheldore    schedule 24.09.2018

Ответы (1)


arrow_upward
1
arrow_downward

Посмотрите на собственные значения, возвращаемые np.linalg.eig(A.T @ A):

In [57]: evals, evecs = np.linalg.eig(A.T @ A)

In [58]: evals
Out[58]: array([2.50000000e+01, 3.61082692e-15, 9.00000000e+00])

Таким образом (игнорируя обычную неточность с плавающей запятой), он вычислил [25, 0, 9]. Собственные векторы, связанные с этими собственными значениями, находятся в столбцах evecs в том же порядке. Но ваша конструкция S не соответствует этому порядку; вот твой S:

In [60]: S
Out[60]: 
array([[5., 0., 0.],
       [0., 3., 0.]])

Когда вы вычисляете U @ S @ V.T, значения в S @ V.T выровнены неправильно.

В качестве быстрого исправления вы можете перезапустить свой код с явно установленным S следующим образом:

S = np.array([[5, 0, 0],
              [0, 0, 3]])

С этим изменением ваш код выводит

[[ 3  2  2]
 [ 2  3 -2]]
[[-3. -2. -2.]
 [-2. -3.  2.]]

Так лучше, но почему знаки неправильные? Теперь проблема в том, что вы независимо вычислили U и V. Собственные векторы не уникальны; они являются базой собственного пространства, и такая база не уникальна. Если собственное значение простое и если вектор нормализован до длины, равной единице (что и делает numpy.linalg.eig), все равно остается выбор знака. То есть, если v является собственным вектором, то и -v тоже. Выбор, сделанный eig при вычислении U и V, не обязательно приведет к восстановлению знака A при вычислении U @ S @ V.T.

Оказывается, вы можете получить ожидаемый результат, просто поменяв местами все знаки либо в U, либо в V. Вот модифицированная версия вашего скрипта, которая генерирует ожидаемый результат:

import numpy as np

A = np.array([[3,  2,  2],
              [2,  3, -2]])

U = np.linalg.eig(A @ A.T)[1]
V = -np.linalg.eig(A.T @ A)[1]
#S = np.c_[np.diag(np.sqrt(np.linalg.eig(A @ A.T)[0])), [0,0]]
S = np.array([[5, 0, 0],
              [0, 0, 3]])

print(A)
print(U @ S @ V.T)

Выход:

[[ 3  2  2]
 [ 2  3 -2]]
[[ 3.  2.  2.]
 [ 2.  3. -2.]]
person Warren Weckesser    schedule 24.09.2018
comment
Но разве для СВД матрица S не должна иметь сингулярные значения по диагонали в порядке убывания? - person Oliver G; 25.09.2018
comment
Нет, но матрицы в SVD должны быть построены так, чтобы сингулярные значения и соответствующие сингулярные векторы были выровнены надлежащим образом. Тот факт, что numpy.linalg.svd возвращает единичные значения в порядке убывания, является деталью реализации, а не математическим требованием. - person Warren Weckesser; 25.09.2018