Por que não consigo obter um SVD válido de X por decomposição de autovalor de XX 'e X'X?

9

Estou tentando fazer SVD manualmente:

m<-matrix(c(1,0,1,2,1,1,1,0,0),byrow=TRUE,nrow=3)

U=eigen(m%*%t(m))$vector
V=eigen(t(m)%*%m)$vector
D=sqrt(diag(eigen(m%*%t(m))$values))

U1=svd(m)$u
V1=svd(m)$v
D1=diag(svd(m)$d)

U1%*%D1%*%t(V1)
U%*%D%*%t(V)

Mas a última linha não retorna m. Por quê? Parece ter algo a ver com os sinais desses autovetores ... Ou entendi mal o procedimento?

failstatistician
fonte
É-me dito repetidamente que o sinal não importa nos SVDs ... assim
falha no estatístico
@Amoeba Obrigado por esclarecer isso. Eu estava focando mais na pergunta em inglês do que no código. Estatístico com falha: veja o que D=diag(c(-1,1,1)*sqrt(eigen(m%*%t(m))$values))faz e lembre-se de que a raiz quadrada (assim como qualquer vetor próprio normalizado) é definida apenas para assinar. Para obter mais informações, altere mpara m <- matrix(-2,1,1)e inclua ,1,1)no final de cada uma das chamadas para diag. Este é um exemplo que cria o mesmo problema - mas é tão simples que a natureza do problema se tornará completamente óbvia. 1×1
whuber
11
Entendi. Obrigado! Você tem uma regra geral para determinar o vetor c (-1, 1, 1)? Ou como os sinais das duas decomposições devem ser ligados?
#Re
11
Observe que o truque do @ whuber com c(-1,1,1)funciona, mas Ddefinido assim não está fornecendo valores singulares. Todos os valores singulares devem ser positivos por definição. A questão de como vincular os sinais de Ue Vé boa, e eu não tenho uma resposta. Por que você simplesmente não faz um SVD? :-)
ameba

Respostas:

13

Análise do Problema

O SVD de uma matriz nunca é único. Deixe a matriz ter dimensões deixe que seu SVD sejan × kAn×k

A=UDV

para um matriz com colunas ortonormais, uma diagonal matriz com entradas não-negativos, e uma matriz com colunas ortonormais.U p × p D k × p Vn×pUp×pDk×pV

Agora escolha, arbitrariamente , qualquer diagonal matriz tendo s na diagonal, para que é a identidade . Entãop×p± 1 S 2 = I p × p I pS±1S2=Ip×pIp

A=UDV=UIDIV=U(S2)D(S2)V=(US)(SDS)(VS)

também é um SVD de porque demonstra que tem colunas ortonormais e um cálculo semelhante demonstra que o tem colunas ortonormais. Além disso, como e são diagonais, eles comutam, de onde mostra que ainda possui entradas não negativas.( U S ) ' ( U S ) = S ' L ' L S = S ' I P S = S ' S = S 2 = I P L S V S S D S D S = D S 2 = D DA

(US)(US)=SUUS=SIpS=SS=S2=Ip
USVSSD
SDS=DS2=D
D

O método implementado no código para encontrar um SVD encontra um que diagonaliza e, da mesma forma, um que diagonaliza Ele procede ao cálculo de em termos dos valores próprios encontrados em . O problema é que isto não assegurar uma correspondência consistente das colunas de com as colunas de .A A = ( U D V ) ( U D V ) = U D V V D U = U D 2 U V A ' A = V D 2 V . D D 2 U VU

AA=(UDV)(UDV)=UDVVDU=UD2U
V
AA=VD2V.
DD2UV

Uma solução

Em vez disso, depois de encontrar um e um , use-os para calcularVUV

UAV=U(UDV)V=(UU)D(VV)=D

direta e eficientemente. Os valores diagonais deste não são necessariamente positivos. D (Isso ocorre porque não há nada sobre o processo de diagonalizar ou que garantirá que, uma vez que esses dois processos foram realizados separadamente.) Torne-os positivos escolhendo as entradas ao longo da diagonal de igualar os sinais das entradas de , para que tenha todos os valores positivos. Para compensar isso, multiplique por :AAAASDSDUS

A=UDV=(US)(SD)V.

Isso é um SVD.

Exemplo

Seja com . Um SVD én=p=k=1A=(2)

(2)=(1)(2)(1)

com , e .U=(1)D=(2)V=(1)

Se você diagonalizar , naturalmente escolheria e . Da mesma forma, se você diagonalizar , escolheria . Infelizmente, Em vez disso, calcule Como isso é negativo, defina . Isso ajusta para e para . Você obteve que é um dos dois SVDs possíveis (mas não o mesmo que o original!).AA=(4)U=(1)D=(4)=(2)AA=(4)V=(1)

UDV=(1)(2)(1)=(2)A.
D=UAV=(1)(2)(1)=(2).
S=(1)UUS=(1)(1)=(1)DSD=(1)(2)=(2)
A=(1)(2)(1),

Código

Aqui está o código modificado. Sua saída confirma

  1. O método recria mcorretamente.
  2. U e realmente ainda são ortonormais.V
  3. Mas o resultado não é o mesmo SVD retornado por svd. (Ambos são igualmente válidos.)
m <- matrix(c(1,0,1,2,1,1,1,0,0),byrow=TRUE,nrow=3)

U <- eigen(tcrossprod(m))$vector
V <- eigen(crossprod(m))$vector
D <- diag(zapsmall(diag(t(U) %*% m %*% V)))
s <- diag(sign(diag(D)))  # Find the signs of the eigenvalues
U <- U %*% s              # Adjust the columns of U
D <- s %*% D              # Fix up D.  (D <- abs(D) would be more efficient.)

U1=svd(m)$u
V1=svd(m)$v
D1=diag(svd(m)$d,n,n)

zapsmall(U1 %*% D1 %*% t(V1)) # SVD
zapsmall(U %*% D %*% t(V))    # Hand-rolled SVD
zapsmall(crossprod(U))        # Check that U is orthonormal
zapsmall(tcrossprod(V))       # Check that V' is orthonormal
whuber
fonte
11
+1. Isto é muito claro. Eu apenas acrescentaria que, na prática, é suficiente calcular um Uou outro Ve então obter outra matriz através da multiplicação A. Dessa maneira, é possível executar apenas uma (em vez de duas) composições automáticas, e os sinais sairão corretos.
Ameba
2
@Amoeba É isso mesmo: no espírito da computação manual de um SVD, que evidentemente é um exercício educacional, nenhuma atenção é dada aqui à eficiência.
whuber
2
Obrigado por sua amável ajuda! Acho que entendi esse problema (finalmente).
#Re
3
@Federico Obrigado por esse lembrete. Você está certo - assumi implicitamente que todos os autovalores são distintos, pois esse certamente será o caso em aplicações estatísticas e alguém perde o hábito de considerar as ambiguidades dos espaços auto-degenerados.
whuber
3
Você está correto, este é apenas um caso extremo, e realmente complicado. Em certo sentido, é outra manifestação do mesmo problema que você delinear em sua resposta, que esse método não garante um "matching" entre as colunas de e . Computar o SVD a partir das composições eigend ainda é um ótimo exemplo de aprendizado. UV
Federico Poloni
5

Como descrevi em um comentário à resposta do @ whuber, esse método para calcular o SVD não funciona para todas as matrizes . A questão não se limita a sinais.

O problema é que pode haver autovalores repetidos e, neste caso, a composição automática de e não é única e nem todas as opções de e podem ser usadas para recuperar o fator diagonal do SVD. Por exemplo, se você usar qualquer matriz ortogonal não diagonal (digamos, ), então . Entre todas as opções possíveis para a matriz de vetor próprio de , retornará , portanto, neste caso, não é diagonal.Um Um ' L V A = [ 3 / 5 4 / 5 - 4 / 5 3 / 5 ] Uma Um ' = A ' A = I I U = V = I L ' A V = AAAAAUVA=[3/54/54/53/5]AA=AA=IIeigenU=V=IUAV=A

Intuitivamente, essa é outra manifestação do mesmo problema que o @whuber descreve, que deve haver uma "correspondência" entre as colunas de e , e calcular duas composições independentes separadamente não garante isso.VUV

Se todos os valores singulares de forem distintos, a composição automática é única (até escala / sinais) e o método funciona. Observação: ainda não é uma boa idéia usá-lo no código de produção em um computador com aritmética de ponto flutuante, porque quando você forma os produtos e o resultado calculado pode ser perturbado por uma quantidade da ordem de , onde é a precisão da máquina. Se as magnitudes dos valores singulares diferem bastante (mais de , aproximadamente), isso é prejudicial à precisão numérica dos menores.AAAAAu 2 × 10 - 16 10 - 8A2uu2×1016108

A computação do SVD a partir das duas composições independentes é um ótimo exemplo de aprendizado, mas nas aplicações da vida real sempre use a svdfunção de R para calcular a decomposição de valor singular.

Federico Poloni
fonte
11
Este comentário é um bom conselho. Observe, porém, que esse tópico não está preocupado com a maneira correta de calcular o SVD (e acredito que ninguém argumentaria contra sua recomendação). O OP aceita implicitamente que svdfunciona. De fato, eles o usam como um padrão contra o qual comparar um cálculo manual, cujo objetivo é verificar o entendimento, e não substituir de svdforma alguma.
whuber
@whuber Observação correta; Eu reformulei o último parágrafo.
Federico Poloni