Descomposición en valores singulares

Para ilustrar la descomposición en valores singulares de una matriz vamos a hacer uso del Ejemplo VI.1.8 de [I. Ojeda, J. Gago, Métodos matemáticos para la Estadística].

Otro buen ejemplo de cálculo de esta descomposición paso a paso se puede encontrar aquí.

Vamos a utilizar singular_value_decomposition que apareció en la versión 1.8 de sympy. Si nuestra versión es anterior, podemos ejecutar pip install "sympy>=1.8" (luego tenemos que reiniciar el núcleo).

from sympy import Matrix

A=Matrix([(2,0,1),(3,-1,1),(-2,4,1),(1,1,1)])
A

\(\displaystyle \left[\begin{matrix}2 & 0 & 1\\3 & -1 & 1\\-2 & 4 & 1\\1 & 1 & 1\end{matrix}\right]\)

U,D,V=A.singular_value_decomposition()
U

\(\displaystyle \left[\begin{matrix}\frac{1}{2} & - \frac{\sqrt{14}}{14}\\\frac{1}{2} & - \frac{\sqrt{14}}{7}\\\frac{1}{2} & \frac{3 \sqrt{14}}{14}\\\frac{1}{2} & 0\end{matrix}\right]\)

D

\(\displaystyle \left[\begin{matrix}2 \sqrt{3} & 0\\0 & 2 \sqrt{7}\end{matrix}\right]\)

V

\(\displaystyle \left[\begin{matrix}\frac{\sqrt{3}}{3} & - \frac{\sqrt{2}}{2}\\\frac{\sqrt{3}}{3} & \frac{\sqrt{2}}{2}\\\frac{\sqrt{3}}{3} & 0\end{matrix}\right]\)

U*D*V.T==A

True

Proceso paso a paso a partir de los valores propios de \(A^tA\)

Vamos a hacer el proceso paso a paso. Como sabemos, los valores singulares son las raíces cuadradas de los valores propios de \(A^tA\).

from sympy import sqrt, eye, GramSchmidt

A.T*A

\(\displaystyle \left[\begin{matrix}18 & -10 & 4\\-10 & 18 & 4\\4 & 4 & 4\end{matrix}\right]\)

ev=(A.T*A).eigenvals()
ev

{28: 1, 12: 1, 0: 1}

Tomamos las raíces cuadradas de los valores no nulos (queremos una descomposición corta).

d=[sqrt(a) for a in ev.keys() if a!=0]
d

[2*sqrt(7), 2*sqrt(3)]

La matriz \(V\) está formada por los vectores propios correspondientes a 28 y 12.

V28=(A.T*A - 28*eye(3)).nullspace()
V28

[Matrix([
 [-1],
 [ 1],
 [ 0]])]

V12=(A.T*A - 12*eye(3)).nullspace()
V12

[Matrix([
 [1],
 [1],
 [1]])]

Juntamos las bases, ortogonalizamos por Gram-Schmidt y las ponemos como columnas en una matriz.

V=Matrix.hstack(*GramSchmidt(V28+V12,True))
V

\(\displaystyle \left[\begin{matrix}- \frac{\sqrt{2}}{2} & \frac{\sqrt{3}}{3}\\\frac{\sqrt{2}}{2} & \frac{\sqrt{3}}{3}\\0 & \frac{\sqrt{3}}{3}\end{matrix}\right]\)

D=Matrix.diag(d)
D

\(\displaystyle \left[\begin{matrix}2 \sqrt{7} & 0\\0 & 2 \sqrt{3}\end{matrix}\right]\)

Como \(A=UDV^t\), tenemos \(AV=UD\), luego \(U\) se calcula como sigue.

U=A*V*D.inv()
U

\(\displaystyle \left[\begin{matrix}- \frac{\sqrt{14}}{14} & \frac{1}{2}\\- \frac{\sqrt{14}}{7} & \frac{1}{2}\\\frac{3 \sqrt{14}}{14} & \frac{1}{2}\\0 & \frac{1}{2}\end{matrix}\right]\)

Comprobemos que es efectivamente una descomposición de \(A\).

U*D*V.T==A

True