X <- matrix(c( 10, -1, -2, 2,
4, 8, -1, 3,
-2, 5, -3, 4,
-7, -8, 6, -2,
-2 , 5, 4, -4
),
byrow = T, nrow = 5, ncol = 4,
dimnames = list(c("P1", "P2", "P3", "P4","P5"),
c("V1", "V2", "V3","V4")))
X
#> V1 V2 V3 V4
#> P1 10 -1 -2 2
#> P2 4 8 -1 3
#> P3 -2 5 -3 4
#> P4 -7 -8 6 -2
#> P5 -2 5 4 -4Vuelta a la facultad. SVD
2023
R
Álgebra
Unos compis del trabajo están haciendo cosas muy chulas utilizando SVD (Descomposición en valores singulares) y me ha recordado a los tiempos de la universidad.
La SVD es una factorización de una matriz rectangular tal que así.
\[ \mathbf{X} = \mathbf{U} \mathbf{\Sigma} \mathbf{V^T} \]
Veamos un ejemplo, dónde tenemos una matriz de 5x4
Hacemos la svd
dvs <- svd(X)
dvs
#> $d
#> [1] 17.033392 10.982941 6.873957 2.997880
#>
#> $u
#> [,1] [,2] [,3] [,4]
#> [1,] -0.41759401 0.66860794 0.28085132 -0.28960722
#> [2,] -0.50584653 -0.30661800 0.13614448 -0.62790360
#> [3,] -0.23820970 -0.40000698 -0.61431429 -0.21204097
#> [4,] 0.71395758 0.08354597 -0.00214809 -0.68888308
#> [5,] 0.05705485 -0.54033629 0.72470722 0.04840934
#>
#> $v
#> [,1] [,2] [,3] [,4]
#> [1,] -0.6360870 0.61508773 0.4578647 -0.08614145
#> [2,] -0.6015603 -0.77316598 0.2003873 -0.01358129
#> [3,] 0.3855737 -0.13572378 0.5864215 -0.69930092
#> [4,] -0.2912926 0.07389525 -0.6374282 -0.70948815Y vemos que efectivamente
dvs$u %*% diag(dvs$d) %*% t(dvs$v)
#> [,1] [,2] [,3] [,4]
#> [1,] 10 -1 -2 2
#> [2,] 4 8 -1 3
#> [3,] -2 5 -3 4
#> [4,] -7 -8 6 -2
#> [5,] -2 5 4 -4Ahora, la proyección de las filas en el espacio vectorial definido por la descomposión en valores singulares sería
\[ Proj = \mathbf{X} \mathbf{V} \]
De forma que si tenemos unos nuevos datos y queremos proyectarlos sobre la estructura definida por la DVS ya calculada se haría así.
nuevaX <- matrix(c( 20, -2, -3, 2.1,
5, 2, -1, 3
),
byrow = T, nrow = 2, ncol = 4,
dimnames = list(c("P6", "P7"),
c("V1", "V2", "V3","V4")))
nuevaX
#> V1 V2 V3 V4
#> P6 20 -2 -3 2.1
#> P7 5 2 -1 3.0
(proyeccion = nuevaX %*% dvs$v)
#> [,1] [,2] [,3] [,4]
#> P6 -13.287055 14.410438 5.6586547 -1.087689
#> P7 -5.643007 1.886516 0.1913919 -1.887033Y para proyectar los valores de X originales, podemos hacer lo mismo o como sabemos que
\[ \mathbf{X} = \mathbf{U} \mathbf{\Sigma} \mathbf{V^T} \]
entonces
\[ Proj = \mathbf{X} \mathbf{V} = \mathbf{U} \mathbf{\Sigma} \mathbf{V^T} \mathbf{V} = \mathbf{U} \mathbf{\Sigma} \mathbf{I} = \mathbf{U} \mathbf{\Sigma} \]
Es decir, la proyección de las filas son los vectores singulares izquierdos por la matriz diagonal de los valores singulares.
¿Y para qué sirve esto? Pues una cosa interesante de la SVD es que nos sirve para reducir la dimensión de los datos, quedándonos con los k primeros valores singulares y los respectivos k primeros vectores singulares.
Supongamos que queremos quedarnos con k = 3
dvs_red <- svd(X, nu = 3, nv = 3)
dvs_red
#> $d
#> [1] 17.033392 10.982941 6.873957 2.997880
#>
#> $u
#> [,1] [,2] [,3]
#> [1,] -0.41759401 0.66860794 0.28085132
#> [2,] -0.50584653 -0.30661800 0.13614448
#> [3,] -0.23820970 -0.40000698 -0.61431429
#> [4,] 0.71395758 0.08354597 -0.00214809
#> [5,] 0.05705485 -0.54033629 0.72470722
#>
#> $v
#> [,1] [,2] [,3]
#> [1,] -0.6360870 0.61508773 0.4578647
#> [2,] -0.6015603 -0.77316598 0.2003873
#> [3,] 0.3855737 -0.13572378 0.5864215
#> [4,] -0.2912926 0.07389525 -0.6374282Podemos reconstruir la matriz original tal que así
# matriz original
X
#> V1 V2 V3 V4
#> P1 10 -1 -2 2
#> P2 4 8 -1 3
#> P3 -2 5 -3 4
#> P4 -7 -8 6 -2
#> P5 -2 5 4 -4
# reconstrucción
dvs_red$u %*% diag(dvs_red$d[1:3]) %*% t(dvs_red$v)
#> [,1] [,2] [,3] [,4]
#> [1,] 9.925211 -1.011791 -2.607138 1.384017
#> [2,] 3.837849 7.974435 -2.316350 1.664474
#> [3,] -2.054758 4.991367 -3.444527 3.548997
#> [4,] -7.177898 -8.028048 4.555812 -3.465227
#> [5,] -1.987499 5.001971 4.101486 -3.897035Y así almacenando U, d y V podemos tener una aproximación de X ocupando menos espacio.
Una utilidad de esto es que si pensamos en las filas como observaciones y las columnnas como variables, la SVD nos sirve para hacer una reducción de la dimensionalidad, de hecho los valores singulares al cuadrado son los autovalores de \(\mathbf{X^T} \mathbf{X}\)
Z = t(X) %*% X
eigen(Z)
#> eigen() decomposition
#> $values
#> [1] 290.136446 120.624987 47.251282 8.987284
#>
#> $vectors
#> [,1] [,2] [,3] [,4]
#> [1,] 0.6360870 0.61508773 -0.4578647 0.08614145
#> [2,] 0.6015603 -0.77316598 -0.2003873 0.01358129
#> [3,] -0.3855737 -0.13572378 -0.5864215 0.69930092
#> [4,] 0.2912926 0.07389525 0.6374282 0.70948815
dvs_red$d^2
#> [1] 290.136446 120.624987 47.251282 8.987284Si nos quedamos con los k primeros singulares podríamos querer tener por ejemplo la proyección de las filas en ese subespacio, y usar esas coordenadas para ver distancias entre filas. Una propiedad deseable de ese subespacio es que las distancias entre los individuos sean parecidas o al menos ordenen de la misma forma que en el espacio original.
Por ejemplo las distancias euclídeas entre las filas de X
dist(X)
#> P1 P2 P3 P4
#> P2 10.908712
#> P3 13.601471 7.071068
#> P4 20.445048 21.236761 17.635192
#> P5 15.874508 10.908712 10.630146 14.212670Vemos que P5 está más cerca de P3 que de P2.
Si usamos sólo la matriz \(\mathbf{U_{n \times k}}\) para ver distancias tenemos que no se mantiene esa relación entre las distancias
dist(dvs_red$u)
#> 1 2 3 4
#> 2 0.9898455
#> 3 1.4055028 0.8022090
#> 4 1.3049119 1.2881284 1.2309319
#> 5 1.3725326 0.8472829 1.3783512 1.1614943En cambio si utilizamos la proyección tenemos que las distancias son del mismo orden que en los datos originales y que se mantiene la relación entre las diferentes filas.
proyeccion <- (dvs_red$u %*% diag(dvs_red$d[1:3]))
dist(proyeccion)
#> 1 2 3 4
#> 2 10.861467
#> 3 13.599483 6.960296
#> 4 20.409979 21.235974 17.577158
#> 5 15.842132 10.718639 10.601432 14.039748Por eso, si queremos utilizar SVD para hacer una reducción de dimensionalidad y posteriormente utilizar técnicas que utilicen distancias puede tener más sentido utilizar la proyección \(\mathbf{U_{n\times k}} \mathbf{\Sigma_{k \times k}}\) que sólo \(\mathbf{U_{n\times k}}\)