¿Por qué no puedo obtener una SVD válida de X a través de la descomposición de valores propios de XX 'y X'X?

Estoy tratando de hacer SVD a mano:

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)

Pero la última línea no regresa m. ¿Por qué? Parece que tiene algo que ver con los signos de estos vectores propios ... ¿O entendí mal el procedimiento?

Análisis del problema

La SVD de una matriz nunca es única. Deje que la matriz tenga dimensiones y deje que su SVD sean × $A$$n\times k$

$$A = U D V^\prime$$

para una matriz con columnas ortonormales, una diagonal matriz con entradas no negativas y una matriz con columnas ortonormales.U p × p D k × p $n\times p$$U$$p\times p$$D$$k\times p$$V$

Ahora elija, arbitrariamente , cualquier diagonal matriz tenga s en la diagonal, de modo que sea ​​la identidad . EntoncesS ± 1 S 2 = I p × p I $p\times p$$S$$\pm 1$$S^2 = I$$p\times p$$I_p$

$$A = U D V^\prime = U I D I V^\prime = U (S^2) D (S^2) V^\prime = (US) (SDS) (VS)^\prime$$

también es una SVD de porque demuestra que tiene columnas ortonormales y un cálculo similar demuestra que tiene columnas ortonormales. Además, dado que y son diagonales, conmutan, de donde muestra que todavía tiene entradas no negativas.$A$

$$(US)^\prime(US) = S^\prime U^\prime U S = S^\prime I_p S = S^\prime S = S^2 = I_p$$ $US$$VS$$S$$D$ $$S D S = DS^2 = D$$ $D$

El método implementado en el código para encontrar un SVD encuentra una que diagonaliza y, de manera similar, una que diagonaliza Se procede a calcular en términos de los valores propios encontrados en . El problema es esto no asegura un juego consistente de las columnas de con las columnas de .$U$

$$AA^\prime = (UDV^\prime)(UDV^\prime)^\prime = UDV^\prime V D^\prime U^\prime = UD^2 U^\prime$$ $V$ $$A^\prime A = VD^2V^\prime.$$ $D$$D^2$$U$$V$

Una solución

En cambio, después de encontrar una y una , úselas para calcular$U$$V$

$$U^\prime A V = U^\prime (U D V^\prime) V = (U^\prime U) D (V^\prime V) = D$$

directa y eficientemente Los valores diagonales de esta no son necesariamente positivos. $D$ (Esto se debe a que no hay nada en el proceso de diagonalización de o que garantice eso, ya que esos dos procesos se llevaron a cabo por separado). Hazlos positivos eligiendo las entradas a lo largo de la diagonal de para igualar los signos de las entradas de , de modo que tenga todos los valores positivos. Compensar esto multiplicando a la derecha por :$A^\prime A$$AA^\prime$$S$$D$$SD$$U$$S$

$$A = U D V^\prime = (US) (SD) V^\prime.$$

Esa es una SVD.

Ejemplo

Sea con . Una SVD es$n=p=k=1$$A=(-2)$

$$(-2) = (1)(2)(-1)$$

con , y .$U=(1)$$D=(2)$$V=(-1)$

Si diagonaliza , naturalmente elegiría y . Del mismo modo, si diagonaliza , elegiría . Desafortunadamente, En cambio, calcule Como esto es negativo, establezca . Esto ajusta a y a . Ha obtenido que es uno de los dos SVD posibles (¡pero no el mismo que el original!).$A^\prime A = (4)$$U=(1)$$D=(\sqrt{4})=(2)$$AA^\prime=(4)$$V=(1)$

$$UDV^\prime = (1)(2)(1) = (2) \ne A.$$ $$D=U^\prime A V = (1)^\prime (-2) (1) = (-2).$$ $S=(-1)$$U$$US = (1)(-1)=(-1)$$D$$SD = (-1)(-2)=(2)$ $$A = (-1)(2)(1),$$

Código

Aquí hay un código modificado. Su salida confirma

1. El método se recrea mcorrectamente.
2. $U$ y realmente todavía son ortonormales.$V$
3. Pero el resultado no es el mismo SVD devuelto por svd. (Ambos son 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

Como describí en un comentario a la respuesta de @ whuber, este método para calcular la SVD no funciona para todas las matrices . El problema no se limita a los signos.

El problema es que puede haber valores propios repetidos, y en este caso la descomposición propia de y no es única y no todas las opciones de y se pueden utilizar para recuperar el factor diagonal de la SVD. Por ejemplo, si toma una matriz ortogonal no diagonal (por ejemplo, ), entonces . Entre todas las opciones posibles para la matriz de vector propio de , devolverá , por lo tanto, en este caso, no es diagonal.$A'A$$AA'$$U$$V$$A=\begin{bmatrix}3/5&4/5\\-4/5&3/5\end{bmatrix}$$AA'=A'A=I$$I$eigen$U=V=I$$U'AV=A$

Intuitivamente, esta es otra manifestación del mismo problema que @whuber describe, que tiene que haber una "coincidencia" entre las columnas de y , y calcular dos descomposiciones propias por separado no lo asegura.$U$$V$

Si todos los valores singulares de son distintos, entonces la descomposición propia es única (hasta escala / signos) y el método funciona. Observación: todavía no es una buena idea usarlo en el código de producción en una computadora con aritmética de coma flotante, porque cuando se forman los productos y el resultado calculado puede verse perturbado por una cantidad del orden de , donde es la precisión de la máquina. Si las magnitudes de los valores singulares difieren mucho (de más de , aproximadamente), esto es perjudicial para la precisión numérica de los más pequeños.A ′ A A A ′ ‖ A ‖ 2 u u ≈ 2 × 10 - 16 10 - $A$$A'A$$AA'$$\|A\|^2u$$u \approx 2\times 10^{-16}$$10^{-8}$

Calcular el SVD a partir de las dos descomposiciones propias es un gran ejemplo de aprendizaje, pero en las aplicaciones de la vida real siempre se usa la svdfunción de R para calcular la descomposición de valores singulares.