[[10. 0. 0.] [ 0. 5. 0.] [ 0. 0. 1.]] Criterio satisfecho en 7 iteraciones para ϵ = 1e-10
Métodos potenciales de prospección, FCAG, 2024.
Macros $\LaTeX$ $\newcommand{\vec}[1]{\mathbf{#1}}$ $\newcommand{\SVD}{U\Sigma V^T}$ $\newcommand{\SVDdecon}{\sigma_1\vec{u}_1\vec{v}_1^T+\ldots+\sigma_r\vec{u}_r\vec{v}_r^T}$ $\newcommand{\R}{\mathbb R}$ $\newcommand{\norm}[2]{||\mathbf{#1}_{#2}||}$
Con este algoritmo podemos indagar en una forma para calcular la SVD. Algoritmos mucho más eficientes y robustos existen. El método tiene ramificaciones que se emplean diariamente en la labor científica actual (randomized linear algebra, matrix completion, etc.)
Lo interesante del método es que podemos aproximar los vectores singulares y valores singulares por medio de las potencias de la matriz $A$ aplicadas a un vector de números aleatorios. En particular, las potencias de $A^TA$ aplicadas a un vector aleatorio deben producir un vector en la misma dirección y sentido que $\vec{v}_1$. Luego de normalizar, se obtiene el valor singular $\sigma_1$ y el vector singular $\vec{u}_1=A\vec{v}_1$.
Si consideramos la primer potencia vista desde la descomposición SVD tenemos $A^TA=V\Sigma^2V^T=\sigma_1^2\vec{v}_1\vec{v}_1^T + \ldots$. Las sucesivas iteraciones potencian aún más al primer valor singular respecto del resto, lo que se traduce en obtener un vector prácticamente en la dirección de $\vec{v}_1$ al aumentar las potencias y evaluar en algún $\vec{x}$: $A^TA\vec{x} = [\sigma_1^2(\vec{v}_1^T\vec{x})]\vec{v}_1 + \ldots$. Más formalmente, si $\vec{x}=\Sigma_j c_j \vec{v}_j$, entonces $(A^TA)^p\vec{x}=\Sigma_j c_j \sigma_j^{2p} \vec{v}_j$.
Alcanza con evaluar las potencias en un vector que contenga al menos una componente en la dirección de $\vec{v}_1$ para ir bien encaminados; de allí que se considere un vector aleatorio para $\vec{x}$, ya que es poco probable que $\vec{x}$ no contenga ninguna componente en la dirección del vector singular.
Dada una matriz $A \in \R^{n\times m}$, para obtener la descomposición SVD $A=\SVD=\SVDdecon$:
El método funciona mientras los valores singulares consecutivos tengan valores bien distintos.
Criterio de corte. Una forma de detener las iteraciones es considerar el ángulo formado por dos vectores consecutivos en la power iteration. Si estos vectores tienen aproximadamente la misma direccion y apuntan en el mismo sentido (el producto escalar normalizado es aproximadamente 1), paramos. Esto es:
Implementamos el algoritmo directamente como primer aproximación sobre la matriz
\begin{equation} A= \begin{bmatrix} 10 &0 &0 \\ 0 &5 &0 \\ 0& 0& 1 \end{bmatrix}. \end{equation}La descomposición SVD es sencillamente $A = I A I $. Claramente los valores singulares son $\{10,5,1\}$.
A: [[10. 0. 0.] [ 0. 5. 0.] [ 0. 0. 1.]] Iteración 1 u1: [[ 1.00000000e+00] [-6.92772059e-13] [-7.50518520e-44]],σ1: 10.0, v1^T: [[ 1.00000000e+00 -1.38554412e-12 -7.50518520e-43]] Iteración 2 u2: [[-2.77108824e-12] [-1.00000000e+00] [-2.91742256e-52]],σ2: 5.0, v2^T: [[-1.38554412e-12 -1.00000000e+00 -1.45871128e-51]] Iteración 3 u3: [[ 7.50518520e-42] [-7.29875578e-51] [ 1.00000000e+00]],σ3: 1.0, v3^T: [[ 7.50518520e-43 -1.45975116e-51 1.00000000e+00]]
Analizamos detener el número de iteraciones si los vectores son prácticamente colineales (producto escalar cerca de 1):
[[10. 0. 0.] [ 0. 5. 0.] [ 0. 0. 1.]] Criterio satisfecho en 7 iteraciones para ϵ = 1e-10
Probamos la subrutina incorporando el criterio de corte. Construimos un ejemplo donde $A \in \R^{3,3}$ es diagonal:
\begin{equation} A= \begin{bmatrix} 10 &0 &0 \\ 0 &5 &0 \\ 0& 0& 1 \end{bmatrix} \end{equation}Claramente los valores singulares son $\{10,5,1\}$.
SNR (algoritmo)= 0.0 SNR (numpy) = 0.0 A: [[10. 0. 0.] [ 0. 5. 0.] [ 0. 0. 1.]] S: [10. 5. 1.] U: [[ 1.00000000e+00 -2.39004491e-09 -9.79913124e-30] [ 9.56017963e-09 1.00000000e+00 -2.31416515e-11] [-5.53095862e-18 -5.78541286e-10 -1.00000000e+00]] VT: [[ 1.00000000e+00 4.78008981e-09 -5.53095862e-19] [-4.78008981e-09 1.00000000e+00 -1.15708257e-10] [-9.79913124e-29 -1.15708257e-10 -1.00000000e+00]] Comparación con numpy.linalg.svd(): S: [10. 5. 1.] U: [[1. 0. 0.] [0. 1. 0.] [0. 0. 1.]] VT: [[1. 0. 0.] [0. 1. 0.] [0. 0. 1.]]
Ahora un ejemplo con una matriz de números aleatorios $\R^{10,10}$.
SNR (algoritmo)= 1.45 SNR (numpy) = 0.0
Finalmente, analizamos un caso donde el método falla. La matriz es ($3\times 2$) de rango 1, debe entonces tener un valor singular 0.
σ2: 0 Resultado truncado hasta u1,v1^T A: [[10. 20.] [ 1. 2.] [ 1. 2.]] S: [22.58317958] U: [[0.99014754] [0.09901475] [0.09901475]] VT: [[0.4472136 0.89442719]] Comparación con numpy.linalg.svd(): S: [2.25831796e+01 5.06714263e-17] U: [[-0.99014754 0.14002801] [-0.09901475 -0.70014004] [-0.09901475 -0.70014004]] VT: [[-0.4472136 -0.89442719] [ 0.89442719 -0.4472136 ]]
Nota. El resultado del algoritmo obtenido es similar al producido por la SVD "económica", donde $U$ no contiene su complemento ortogonal:
A = numpy.linalg.svd(A,full_matrices=False)
Cuando se utiliza:
A = numpy.linalg.svd(A,full_matrices=True)
obtenemos además de los vectores columna en $U$ dados por el método ilustrado aquí, productos vectoriales entre los vectores de $U$ hasta llenar la dimensión. Notamos que la SVD puede "diferir" en los resultados de sus vectores singulares, traslando el signo de un vector de un lado a otro.
A continuación vemos un ejemplo para comprender la diferencia entre Full_matrices = False | True
A: [[10. 0.] [ 1. 1.] [ 1. 1.]] Full_matrices = False ===================== S: [10.10148425 1.40000571] U: [[-0.98975576 0.14277093] [-0.1009543 -0.69986301] [-0.1009543 -0.69986301]] VT: [[-0.99980022 -0.01998801] [ 0.01998801 -0.99980022]] Cross product between vector in U given by full_matrices=False (minus the crossproduct in fact): [[-0. ] [-0.70710678] [ 0.70710678]] Full_matrices = True ===================== [10.10148425 1.40000571] [[-9.89755758e-01 1.42770935e-01 -1.38094227e-17] [-1.00954296e-01 -6.99863008e-01 -7.07106781e-01] [-1.00954296e-01 -6.99863008e-01 7.07106781e-01]] [[-0.99980022 -0.01998801] [ 0.01998801 -0.99980022]] The last column vector in U is: [[-1.38094227e-17] [-7.07106781e-01] [ 7.07106781e-01]]
Eso es todo por hoy.