Qué es SVD, la que sostiene medio machine learning

El problema que arregla 🌸

Los autovalores del capítulo 10 tienen una limitación grande: solo existen para matrices cuadradas. Y tu tabla de datos casi nunca lo es. La mía tiene 1.793 filas y 5 columnas.

Por eso en el capítulo anterior no le sacamos los autovectores a la tabla: se los sacamos a su matriz de covarianza, que sí es cuadrada. Fue un rodeo.

La SVD quita el rodeo. Funciona con cualquier matriz, del tamaño que sea.

Abrimos el archivo con más columnas

Este bloque necesita una tabla más ancha. Con las tres de siempre no hay nada que descomponer, porque están casi sin relación entre ellas. Así que meto dos más: el descuento y el monto final facturado.

import pandas as pd
import numpy as np

COLS = ['unidades', 'monto', 'descuento', 'satisfaccion',
        'monto_final_facturado']

df = pd.read_csv('ventas-miss-yera.csv')
for c in COLS:
    df[c] = pd.to_numeric(df[c], errors='coerce')
M = df[COLS].dropna().to_numpy()

Z = (M - M.mean(axis=0)) / M.std(axis=0)
print('M.shape:', M.shape)

M.shape: (1793, 5)

Esa Z es la tabla escalada: a cada columna le resté su promedio y la dividí por su desviación. Después de eso todas tienen media cero y desviación uno, o sea que ninguna manda por tener números más grandes. En el capítulo 7 viste por qué eso hace falta.

Las tres piezas

U, S, Vt = np.linalg.svd(Z, full_matrices=False)
print('U ', U.shape)
print('S ', S.shape)
print('Vt', Vt.shape)
print('valores singulares:', np.round(S, 4))

U  (1793, 5)
S  (5,)
Vt (5, 5)
valores singulares: [55.9718 43.1322 42.2936 41.5137 21.4392]

Girar, estirar, girar. Cualquier matriz del mundo se puede escribir así, y eso es lo que la hace tan útil.

Fíjate en que S viene como un vector de 5 números y no como una matriz. NumPy te da solo la diagonal porque el resto son ceros y sería un desperdicio guardarlos.

Que sí reconstruye

print(np.allclose(U @ np.diag(S) @ Vt, Z))

True

Nada se perdió. Es la misma tabla escrita de otra manera.

Y ahora la conexión con PCA

C = np.cov(Z, rowvar=False)
vals = np.sort(np.linalg.eigh(C)[0])[::-1]
print('S al cuadrado entre n-1:', np.round(S ** 2 / (Z.shape[0] - 1), 4))
print('autovalores de C       :', np.round(vals, 4))

S al cuadrado entre n-1: [1.7482 1.0382 0.9982 0.9617 0.2565]
autovalores de C       : [1.7482 1.0382 0.9982 0.9617 0.2565]

Los mismos números. No parecidos: los mismos.

SVD y PCA son la misma cosa calculada de dos formas. Y las librerías de verdad usan SVD, porque hace las cuentas sobre tus datos directos sin tener que construir la covarianza primero, y esa construcción es justamente donde se pierde precisión.

vals, vecs = np.linalg.eigh(C)
orden = np.argsort(vals)[::-1]
print('primera fila de Vt:', np.round(Vt[0], 4))
print('primer autovector :', np.round(vecs[:, orden][:, 0], 4))

primera fila de Vt: [ 0.     -0.6996  0.0969 -0.0631 -0.7051]
primer autovector : [-0.      0.6996 -0.0969  0.0631  0.7051]

Los mismos números con el signo cambiado. Y eso no es un error: una dirección y su contraria son la misma dirección. Si un autovector apunta al noreste, el mismo con signo cambiado apunta al suroeste, y la recta que definen es idéntica.

Te lo digo porque va a pasar: corres el mismo PCA dos veces, o en dos máquinas, y los signos salen cambiados. No te has vuelto loca 🙂

Para qué se usa: quedarse con lo grande

k = 2
Zk = U[:, :k] @ np.diag(S[:k]) @ Vt[:k]
print('rango de Z :', np.linalg.matrix_rank(Z))
print('rango de Zk:', np.linalg.matrix_rank(Zk))
print('error medio:', round(float(np.abs(Z - Zk).mean()), 4))
print('numeros de la tabla entera:', Z.size)
print('numeros que hacen falta   :', U[:, :k].size + k + Vt[:k].size)

rango de Z : 5
rango de Zk: 2
error medio: 0.4905
numeros de la tabla entera: 8965
numeros que hacen falta   : 3598

Me quedé con los dos valores singulares más grandes y tiré los tres chicos. El resultado es una tabla de rango 2 que se guarda con 3.598 números en vez de 8.965.

Eso se llama aproximación de rango bajo, y hay un teorema que dice que de todas las matrices de rango 2 que existen, esta es la más parecida a la original. No es una heurística: es lo mejor que se puede hacer.

Con mis datos el error es 0.49 y no compensa, por lo mismo de siempre: acá no había redundancia que exprimir. Pero cámbiale los datos y la historia cambia.

Dónde la vas a encontrar

• Sistemas de recomendación. La tabla de usuarios por productos está casi vacía y es enorme. Quedarse con los k valores singulares grandes es lo que hay detrás de "a quien vio esto también le gustó aquello".
• Compresión de imágenes. Una foto es una matriz de píxeles. Con el 10% de los valores singulares se ve casi igual.
• Buscar temas en textos. Palabras por documentos, misma receta.
• La pseudoinversa. Cuando no hay inversa porque la matriz es singular, como en el capítulo 6, la SVD da lo más parecido que existe.

En todos los casos la idea es la misma: casi toda la información está en unas pocas direcciones, y las demás son ruido 🌟

Lo que te llevas

• La SVD parte cualquier matriz en giro, estirado y giro.
• Funciona con matrices que no son cuadradas, y los autovalores no.
• S al cuadrado entre n menos uno son los autovalores de la covarianza.
• Las librerías calculan PCA con SVD, y por eso.
• Los signos pueden salir cambiados y da igual.
• Quedarse con los k más grandes es la mejor aproximación de rango k que existe.

Comprueba que se entendió

Ejercicios

U, S, Vt = np.linalg.svd(Z, full_matrices=False)
print('S                :', np.round(S, 4))
print('ya venia ordenado:', bool(np.all(np.diff(S) <= 0)))
print('todos positivos  :', bool(np.all(S > 0)))

S                : [55.9718 43.1322 42.2936 41.5137 21.4392]
ya venia ordenado: True
todos positivos  : True

U, S, Vt = np.linalg.svd(Z, full_matrices=False)
for k in range(1, 6):
    Zk = U[:, :k] @ np.diag(S[:k]) @ Vt[:k]
    print(k, '-> error', round(float(np.abs(Z - Zk).mean()), 4),
          '| guarda', U[:, :k].size + k + Vt[:k].size, 'numeros')

1 -> error 0.5975 | guarda 1799 numeros
2 -> error 0.4905 | guarda 3598 numeros
3 -> error 0.3578 | guarda 5397 numeros
4 -> error 0.0966 | guarda 7196 numeros
5 -> error 0.0 | guarda 8995 numeros

Zd = np.column_stack([Z, Z[:, 0] * 3])
_, Sd, _ = np.linalg.svd(Zd, full_matrices=False)
print('valores singulares:', np.round(Sd, 6))
print('cuantos no son cero:', int((Sd > 1e-10).sum()))
print('matrix_rank dice   :', np.linalg.matrix_rank(Zd))

valores singulares: [133.913495  55.971759  42.440309  42.136853  21.439308   0.      ]
cuantos no son cero: 5
matrix_rank dice   : 5