2015-01-05 11 views
15

Рассмотрим разложение особых значений M = USV *. Тогда разложение на собственные значения M * M дает M * M = V (S * S) V * = VS * U * USV *. Я хотел бы проверить это равенство с NumPy, показав, что собственные векторы, возвращаемые eigh функции такие же, как те, возвращаемый svd функции:Собственные векторы, вычисляемые с помощью eigh и svd numpy, не соответствуют

import numpy as np 
np.random.seed(42) 
# create mean centered data 
A=np.random.randn(50,20) 
M= A-np.array(A.mean(0),ndmin=2) 

# svd 
U1,S1,V1=np.linalg.svd(M) 
S1=np.square(S1) 
V1=V1.T 

# eig 
S2,V2=np.linalg.eigh(np.dot(M.T,M)) 
indx=np.argsort(S2)[::-1] 
S2=S2[indx] 
V2=V2[:,indx] 

# both Vs are in orthonormal form 
assert np.all(np.isclose(np.linalg.norm(V1,axis=1), np.ones(V1.shape[0]))) 
assert np.all(np.isclose(np.linalg.norm(V1,axis=0), np.ones(V1.shape[1]))) 
assert np.all(np.isclose(np.linalg.norm(V2,axis=1), np.ones(V2.shape[0]))) 
assert np.all(np.isclose(np.linalg.norm(V2,axis=0), np.ones(V2.shape[1]))) 

assert np.all(np.isclose(S1,S2)) 
assert np.all(np.isclose(V1,V2)) 

Последнее утверждение не удается. Зачем?

+0

Вы можете добавить положительное число ко всем диагональным элементам, т. е. сделать M2 = M + a * I, где a достаточно велико, чтобы сделать M2 положительным полудефектом. Тогда SVD и eigh должны лучше соглашаться. –

ответ

14

Просто играйте с небольшими номерами, чтобы отладить вашу проблему.

Начните с A=np.random.randn(3,2) вместо вашей гораздо большей матрицы с размером (50,20)

В моем случайном случае, я считаю, что

v1 = array([[-0.33872745, 0.94088454], 
    [-0.94088454, -0.33872745]]) 

и v2:

v2 = array([[ 0.33872745, -0.94088454], 
    [ 0.94088454, 0.33872745]]) 

они отличаются только знак и, очевидно, даже если нормализовать единичный модуль, вектор может отличаться для знака.

Теперь, если вы попытаетесь Хитрость

assert np.all(np.isclose(V1,-1*V2)) 

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

Правильный способ проверить равенство между векторами:

assert allclose(abs((V1*V2).sum(0)),1.) 

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

(V1*V2).sum(0) 

, что на самом деле есть либо +1 или -1 в зависимости от вектора:

array([ 1., -1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 
    1., -1., 1., 1., 1., -1., -1.]) 

EDIT: Это произойдет в большинстве случаев, особенно если вы начинаете с случайной матрицы. Обратите внимание, однако, что этот тест будет, вероятно, потерпеть неудачу, если один или несколько собственных значений имеет собственное подпространство размерности больше, чем 1, как отметил @Sven Marnach в своем комментарии ниже:

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

+0

@matus ОК, я потерян :) Но я верю в ваше суждение, и поэтому я удалю свои комментарии, чтобы не путать будущих читателей. Ура! – BartoszKP

+0

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

+0

@SvenMarnach, это очень актуальная точка. Я отредактирую сообщение, чтобы указать на это оговорку – gg349

 Смежные вопросы

  • Нет связанных вопросов^_^