Cinco maneras de hacer mínimos cuadrados (con antorcha)

Cinco maneras de hacer minimos cuadrados con antorcha

Nota: Esta publicación es una versión resumida de un capítulo de la tercera parte del próximo libro, Deep Learning and Scientific Computing with R torch. La tercera parte está dedicada a la computación científica más allá del aprendizaje profundo. A lo largo del libro, me concentro en los conceptos subyacentes, esforzándome por explicarlos de la manera más “verbal” que puedo. Esto no significa saltarse las ecuaciones; significa cuidarse de explicar por qué son como son.

¿Cómo se calcula la regresión lineal de mínimos cuadrados? En R, usando lm(); en torchhay linalg_lstsq().

Donde R, a veces, oculta la complejidad del usuario, marcos de computación de alto rendimiento como torch tienden a pedir un poco más de esfuerzo por adelantado, ya sea leyendo cuidadosamente la documentación, o jugando con algunos, o ambos. Por ejemplo, aquí está la pieza central de documentación para linalg_lstsq()profundizando en la driver parámetro a la función:

`driver` chooses the LAPACK/MAGMA function that will be used.
For CPU inputs the valid values are 'gels', 'gelsy', 'gelsd, 'gelss'.
For CUDA input, the only valid driver is 'gels', which assumes that A is full-rank.
To choose the best driver on CPU consider:
  -   If A is well-conditioned (its condition number is not too large), or you do not mind some precision loss:
     -   For a general matrix: 'gelsy' (QR with pivoting) (default)
     -   If A is full-rank: 'gels' (QR)
  -   If A is not well-conditioned:
     -   'gelsd' (tridiagonal reduction and SVD)
     -   But if you run into memory issues: 'gelss' (full SVD).

Si necesitará saber esto dependerá del problema que esté resolviendo. Pero si lo hace, sin duda ayudará a tener una idea de lo que se alude allí, aunque solo sea de una manera de alto nivel.

En nuestro problema de ejemplo a continuación, vamos a tener suerte. Todos los controladores devolverán el mismo resultado, pero solo una vez que hayamos aplicado una especie de «truco». El libro analiza por qué eso funciona; No haré eso aquí, para mantener la publicación razonablemente breve. Lo que haremos en su lugar es profundizar en los diversos métodos utilizados por linalg_lstsq()así como algunos otros de uso común.

El plan

La forma en que organizaremos esta exploración es resolviendo un problema de mínimos cuadrados desde cero, haciendo uso de varias factorizaciones de matrices. Concretamente, abordaremos la tarea:

  1. Por medio de los llamados ecuaciones normalesla forma más directa, en el sentido de que resulta inmediatamente de un enunciado matemático del problema.

  2. De nuevo, partiendo de las ecuaciones normales, pero haciendo uso de Factorización de Cholesky en resolverlos.

  3. Una vez más, tomando las ecuaciones normales como punto de partida, pero procediendo por medio de LU descomposición.

  4. A continuación, empleando otro tipo de factorización: código QR – que, junto con el último, da cuenta de la gran mayoría de las descomposiciones aplicadas “en el mundo real”. Con la descomposición QR, el algoritmo de solución no parte de las ecuaciones normales.

  5. Y, por último, haciendo uso de Valor singular de descomposición (SVD). Aquí tampoco se necesitan las ecuaciones normales.

Regresión para la predicción del tiempo

El conjunto de datos que usaremos está disponible en el Repositorio de aprendizaje automático de UCI.

Rows: 7,588
Columns: 25
$ station           <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11,…
$ Date              <date> 2013-06-30, 2013-06-30,…
$ Present_Tmax      <dbl> 28.7, 31.9, 31.6, 32.0, 31.4, 31.9,…
$ Present_Tmin      <dbl> 21.4, 21.6, 23.3, 23.4, 21.9, 23.5,…
$ LDAPS_RHmin       <dbl> 58.25569, 52.26340, 48.69048,…
$ LDAPS_RHmax       <dbl> 91.11636, 90.60472, 83.97359,…
$ LDAPS_Tmax_lapse  <dbl> 28.07410, 29.85069, 30.09129,…
$ LDAPS_Tmin_lapse  <dbl> 23.00694, 24.03501, 24.56563,…
$ LDAPS_WS          <dbl> 6.818887, 5.691890, 6.138224,…
$ LDAPS_LH          <dbl> 69.45181, 51.93745, 20.57305,…
$ LDAPS_CC1         <dbl> 0.2339475, 0.2255082, 0.2093437,…
$ LDAPS_CC2         <dbl> 0.2038957, 0.2517714, 0.2574694,…
$ LDAPS_CC3         <dbl> 0.1616969, 0.1594441, 0.2040915,…
$ LDAPS_CC4         <dbl> 0.1309282, 0.1277273, 0.1421253,…
$ LDAPS_PPT1        <dbl> 0.0000000, 0.0000000, 0.0000000,…
$ LDAPS_PPT2        <dbl> 0.000000, 0.000000, 0.000000,…
$ LDAPS_PPT3        <dbl> 0.0000000, 0.0000000, 0.0000000,…
$ LDAPS_PPT4        <dbl> 0.0000000, 0.0000000, 0.0000000,…
$ lat               <dbl> 37.6046, 37.6046, 37.5776, 37.6450,…
$ lon               <dbl> 126.991, 127.032, 127.058, 127.022,…
$ DEM               <dbl> 212.3350, 44.7624, 33.3068, 45.7160,…
$ Slope             <dbl> 2.7850, 0.5141, 0.2661, 2.5348,…
$ `Solar radiation` <dbl> 5992.896, 5869.312, 5863.556,…
$ Next_Tmax         <dbl> 29.1, 30.5, 31.1, 31.7, 31.2, 31.5,…
$ Next_Tmin         <dbl> 21.2, 22.5, 23.9, 24.3, 22.5, 24.0,…

Por la forma en que enmarcamos la tarea, casi todo en el conjunto de datos sirve como predictor. Como objetivo, usaremos Next_Tmax, la temperatura máxima alcanzada en el día siguiente. Esto significa que debemos eliminar Next_Tmin del conjunto de predictores, ya que sería una pista demasiado poderosa. Haremos lo mismo para stationla identificación de la estación meteorológica y Date. Esto nos deja con veintiún predictores, incluidas las mediciones de la temperatura real (Present_Tmax, Present_Tmin), modelos de pronósticos de varias variables (LDAPS_*) e información auxiliar (lat, lony `Solar radiation`entre otros).

Tenga en cuenta cómo, arriba, he agregado una línea a estandarizar los predictores. Este es el «truco» al que me refería anteriormente. Para ver qué sucede sin estandarización, consulte el libro. (La conclusión es: tendrías que llamar linalg_lstsq() con argumentos no predeterminados).

Para torchdividimos los datos en dos tensores: una matriz Aque contiene todos los predictores y un vector b que sostiene el objetivo.

weather <- torch_tensor(weather_df %>% as.matrix())
A <- weather[ , 1:-2]
b <- weather[ , -1]

dim(A)
[1] 7588   21

Ahora, primero determinemos el resultado esperado.

Establecer expectativas con lm()

Si hay una implementación de mínimos cuadrados en la que «creemos», seguramente debe ser lm().

fit <- lm(Next_Tmax ~ . , data = weather_df)
fit %>% summary()
Call:
lm(formula = Next_Tmax ~ ., data = weather_df)

Residuals:
     Min       1Q   Median       3Q      Max
-1.94439 -0.27097  0.01407  0.28931  2.04015

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)        2.605e-15  5.390e-03   0.000 1.000000    
Present_Tmax       1.456e-01  9.049e-03  16.089  < 2e-16 ***
Present_Tmin       4.029e-03  9.587e-03   0.420 0.674312    
LDAPS_RHmin        1.166e-01  1.364e-02   8.547  < 2e-16 ***
LDAPS_RHmax       -8.872e-03  8.045e-03  -1.103 0.270154    
LDAPS_Tmax_lapse   5.908e-01  1.480e-02  39.905  < 2e-16 ***
LDAPS_Tmin_lapse   8.376e-02  1.463e-02   5.726 1.07e-08 ***
LDAPS_WS          -1.018e-01  6.046e-03 -16.836  < 2e-16 ***
LDAPS_LH           8.010e-02  6.651e-03  12.043  < 2e-16 ***
LDAPS_CC1         -9.478e-02  1.009e-02  -9.397  < 2e-16 ***
LDAPS_CC2         -5.988e-02  1.230e-02  -4.868 1.15e-06 ***
LDAPS_CC3         -6.079e-02  1.237e-02  -4.913 9.15e-07 ***
LDAPS_CC4         -9.948e-02  9.329e-03 -10.663  < 2e-16 ***
LDAPS_PPT1        -3.970e-03  6.412e-03  -0.619 0.535766    
LDAPS_PPT2         7.534e-02  6.513e-03  11.568  < 2e-16 ***
LDAPS_PPT3        -1.131e-02  6.058e-03  -1.866 0.062056 .  
LDAPS_PPT4        -1.361e-03  6.073e-03  -0.224 0.822706    
lat               -2.181e-02  5.875e-03  -3.713 0.000207 ***
lon               -4.688e-02  5.825e-03  -8.048 9.74e-16 ***
DEM               -9.480e-02  9.153e-03 -10.357  < 2e-16 ***
Slope              9.402e-02  9.100e-03  10.331  < 2e-16 ***
`Solar radiation`  1.145e-02  5.986e-03   1.913 0.055746 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4695 on 7566 degrees of freedom
Multiple R-squared:  0.7802,    Adjusted R-squared:  0.7796
F-statistic:  1279 on 21 and 7566 DF,  p-value: < 2.2e-16

Con una varianza explicada del 78%, el pronóstico está funcionando bastante bien. Esta es la línea de base con la que queremos verificar todos los demás métodos. Con ese fin, almacenaremos las predicciones respectivas y los errores de predicción (este último se operacionaliza como error cuadrático medio, RMSE). Por ahora, solo tenemos entradas para lm():

rmse <- function(y_true, y_pred) 

all_preds <- data.frame(
  b = weather_df$Next_Tmax,
  lm = fit$fitted.values
)
all_errs <- data.frame(lm = rmse(all_preds$b, all_preds$lm))
all_errs
       lm
1 40.8369

Usando torchla forma rápida: linalg_lstsq()

Ahora, por un momento, supongamos que no se trata de explorar diferentes enfoques, sino de obtener un resultado rápido. En torchtenemos linalg_lstsq(), una función dedicada específicamente a resolver problemas de mínimos cuadrados. (Esta es la función cuya documentación estaba citando arriba). Al igual que hicimos con lm()probablemente seguiríamos adelante y lo llamaríamos, haciendo uso de la configuración predeterminada:

x_lstsq <- linalg_lstsq(A, b)$solution

all_preds$lstsq <- as.matrix(A$matmul(x_lstsq))
all_errs$lstsq <- rmse(all_preds$b, all_preds$lstsq)

tail(all_preds)
              b         lm      lstsq
7583 -1.1380931 -1.3544620 -1.3544616
7584 -0.8488721 -0.9040997 -0.9040993
7585 -0.7203294 -0.9675286 -0.9675281
7586 -0.6239224 -0.9044044 -0.9044040
7587 -0.5275154 -0.8738639 -0.8738635
7588 -0.7846007 -0.8725795 -0.8725792

Las predicciones se parecen a las de lm() muy de cerca; de hecho, tan de cerca que podemos suponer que esas pequeñas diferencias se deben solo a errores numéricos que surgen desde el fondo de las respectivas pilas de llamadas. RMSE, por lo tanto, también debería ser igual:

       lm    lstsq
1 40.8369 40.8369

Es; y este es un resultado satisfactorio. Sin embargo, en realidad solo surgió debido a ese “truco”: la normalización. (Nuevamente, tengo que pedirle que consulte el libro para obtener más detalles).

Ahora, exploremos lo que podemos hacer sin usar linalg_lstsq().

Mínimos cuadrados (I): Las ecuaciones normales

Empezamos indicando el objetivo. Dada una matriz, (mathbf)que contiene características en sus columnas y observaciones en sus filas, y un vector de resultados observados, (mathbf)queremos encontrar coeficientes de regresión, uno para cada característica, que nos permitan aproximar (mathbf) tan bien como sea posible. Llame al vector de coeficientes de regresión (mathbf). Para obtenerlo, necesitamos resolver un sistema de ecuaciones simultáneas, que en notación matricial se presenta como

[
mathbf = mathbf
]

Si (mathbf) Si fuera una matriz cuadrada e invertible, la solución podría calcularse directamente como (mathbf = mathbf^mathbf). Sin embargo, esto casi nunca será posible; (con suerte) siempre tendremos más observaciones que predictores. Se necesita otro enfoque. Comienza directamente desde el enunciado del problema.

Cuando usamos las columnas de (mathbf) por (mathbf) para aproximar (mathbf)esa aproximación es necesariamente en el espacio columna de (mathbf). (mathbf), por otro lado, normalmente no lo será. Queremos que esos dos estén lo más cerca posible. En otras palabras, queremos minimizar la distancia entre ellos. Al elegir la norma 2 para la distancia, esto produce el objetivo

[
minimize ||mathbf-mathbf||^2
]

Esta distancia es la longitud (al cuadrado) del vector de errores de predicción. Ese vector es necesariamente ortogonal a (mathbf) sí mismo. Es decir, cuando lo multiplicamos por (mathbf)obtenemos el vector cero:

[
mathbf^T(mathbf – mathbf) = mathbf
]

Un reordenamiento de esta ecuación produce el llamado ecuaciones normales:

[
mathbf^T mathbf mathbf = mathbf^T mathbf
]

Estos pueden ser resueltos por (mathbf)calculando el inverso de (mathbf^Tmathbf):

[
mathbf = (mathbf^T mathbf)^ mathbf^T mathbf
]

(mathbf^Tmathbf) es una matriz cuadrada. Todavía podría no ser invertible, en cuyo caso se calcularía el llamado pseudoinverso. En nuestro caso, esto no será necesario; Ya sabemos (mathbf) tiene rango completo, y también (mathbf^Tmathbf).

Por tanto, a partir de las ecuaciones normales hemos derivado una receta para calcular (mathbf). Pongámoslo en uso y comparemos con lo que obtuvimos de lm() y linalg_lstsq().

AtA <- A$t()$matmul(A)
Atb <- A$t()$matmul(b)
inv <- linalg_inv(AtA)
x <- inv$matmul(Atb)

all_preds$neq <- as.matrix(A$matmul(x))
all_errs$neq <- rmse(all_preds$b, all_preds$neq)

all_errs
       lm   lstsq     neq
1 40.8369 40.8369 40.8369

Habiendo confirmado que la forma directa funciona, podemos permitirnos algo de sofisticación. Aparecerán cuatro factorizaciones de matrices diferentes: Cholesky, LU, QR y Descomposición de valores singulares. El objetivo, en todos los casos, es evitar el costoso cálculo de la (pseudo-) inversa. Eso es lo que todos los métodos tienen en común. Sin embargo, no difieren “solo” en la forma en que se factoriza la matriz, sino también, en cuales matriz es. Esto tiene que ver con las restricciones que imponen los diversos métodos. En términos generales, el orden en que se enumeran arriba refleja una pendiente descendente de condiciones previas, o dicho de otra manera, una pendiente ascendente de generalidad. Debido a las limitaciones involucradas, los dos primeros (Cholesky, así como la descomposición LU) se realizarán en (mathbf^Tmathbf)mientras que los dos últimos (QR y SVD) operan en (mathbf) directamente. Con ellos, nunca hay necesidad de calcular (mathbf^Tmathbf).

Mínimos cuadrados (II): descomposición de Cholesky

En la descomposición de Cholesky, una matriz se factoriza en dos matrices triangulares del mismo tamaño, siendo una la transpuesta de la otra. Esto comúnmente se escribe ya sea

[
mathbf = mathbf mathbf^T
]
o

[
mathbf = mathbf^Tmathbf
]

Aquí símbolos (mathbf) y (mathbf) denotan matrices triangulares inferiores y triangulares superiores, respectivamente.

Para que la descomposición de Cholesky sea posible, una matriz debe ser simétrica y definida positiva. Estas son condiciones bastante fuertes, que a menudo no se cumplirán en la práctica. En nuestro caso, (mathbf) no es simétrico. Esto implica inmediatamente que tenemos que operar en (mathbf^Tmathbf) en cambio. Y desde (mathbf) ya es positivo definido, sabemos que (mathbf^Tmathbf) es, también.

En torchobtenemos la descomposición de Cholesky de una matriz usando linalg_cholesky(). De forma predeterminada, esta llamada volverá (mathbf)una matriz triangular inferior.

# AtA = L L_t
AtA <- A$t()$matmul(A)
L <- linalg_cholesky(AtA)

Vamos a comprobar que podemos reconstruir (mathbf) desde (mathbf):

LLt <- L$matmul(L$t())
diff <- LLt - AtA
linalg_norm(diff, ord = "fro")
torch_tensor
0.00258896
[ CPUFloatType ]

Aquí, calculé la norma de Frobenius de la diferencia entre la matriz original y su reconstrucción. La norma de Frobenius resume individualmente todas las entradas de la matriz y devuelve la raíz cuadrada. En teoría, nos gustaría ver cero aquí; pero en presencia de errores numéricos, el resultado es suficiente para indicar que la factorización funcionó bien.

Ahora que tenemos (mathbfmathbf^T) en vez de (mathbf^Tmathbf), ¿cómo nos ayuda eso? Es aquí donde ocurre la magia, y encontrarás el mismo tipo de magia en el trabajo en los tres métodos restantes. La idea es que por alguna descomposición surja una forma más eficiente de resolver el sistema de ecuaciones que constituye una tarea dada.

Con (mathbfmathbf^T)el punto es que (mathbf) es triangular, y cuando ese es el caso, el sistema lineal se puede resolver por sustitución simple. Eso es mejor visible con un pequeño ejemplo:

[
begin
1 & 0 & 0\
2 & 3 & 0\
3 & 4 & 1
end
begin
x1\
x2\
x3
end
=
begin
1\
11\
15
end
]

Comenzando en la fila superior, inmediatamente vemos que (x1) es igual (1); y una vez que sabemos ese es sencillo calcular, a partir de la fila dos, que (x2) debe ser (3). La última fila nos dice que (x3) debe ser (0).

En codigo, torch_triangular_solve() se utiliza para calcular de manera eficiente la solución de un sistema lineal de ecuaciones donde la matriz de predictores es triangular inferior o superior. Un requisito adicional es que la matriz sea simétrica, pero esa condición ya la teníamos que cumplir para poder utilizar la factorización de Cholesky.

Por defecto, torch_triangular_solve() espera que la matriz sea triangular superior (no inferior); pero hay un parámetro de función, upper, que nos permite corregir esa expectativa. El valor devuelto es una lista y su primer elemento contiene la solución deseada. Para ilustrar, aquí está torch_triangular_solve()aplicado al ejemplo del juguete que resolvimos manualmente arriba:

some_L <- torch_tensor(
  matrix(c(1, 0, 0, 2, 3, 0, 3, 4, 1), nrow = 3, byrow = TRUE)
)
some_b <- torch_tensor(matrix(c(1, 11, 15), ncol = 1))

x <- torch_triangular_solve(
  some_b,
  some_L,
  upper = FALSE
)[[1]]
x
torch_tensor
 1
 3
 0
[ CPUFloatType ]

Volviendo a nuestro ejemplo actual, las ecuaciones normales ahora se ven así:

[
mathbfmathbf^T mathbf = mathbf^T mathbf
]

Introducimos una nueva variable, (mathbf)reposar durante (mathbf^T mathbf),

[
mathbfmathbf = mathbf^T mathbf
]

y calcular la solución de este sistema:

Atb <- A$t()$matmul(b)

y <- torch_triangular_solve(
  Atb$unsqueeze(2),
  L,
  upper = FALSE
)[[1]]

Ahora que tenemos (y)volvemos a ver cómo se definió:

[
mathbf = mathbf^T mathbf
]

Para determinar (mathbf)por lo que podemos volver a utilizar torch_triangular_solve():

x <- torch_triangular_solve(y, L$t())[[1]]

Y ahí estamos.

Como de costumbre, calculamos el error de predicción:

all_preds$chol <- as.matrix(A$matmul(x))
all_errs$chol <- rmse(all_preds$b, all_preds$chol)

all_errs
       lm   lstsq     neq    chol
1 40.8369 40.8369 40.8369 40.8369

Ahora que ha visto el fundamento detrás de la factorización de Cholesky y, como ya se sugirió, la idea se traslada a todas las demás descomposiciones, es posible que desee ahorrarse algo de trabajo utilizando una función de conveniencia dedicada, torch_cholesky_solve(). Esto hará obsoletas las dos llamadas a torch_triangular_solve().

Las siguientes líneas producen el mismo resultado que el código anterior, pero, por supuesto, hacer ocultar la magia subyacente.

L <- linalg_cholesky(AtA)

x <- torch_cholesky_solve(Atb$unsqueeze(2), L)

all_preds$chol2 <- as.matrix(A$matmul(x))
all_errs$chol2 <- rmse(all_preds$b, all_preds$chol2)
all_errs
       lm   lstsq     neq    chol   chol2
1 40.8369 40.8369 40.8369 40.8369 40.8369

Pasemos al siguiente método, de manera equivalente, a la siguiente factorización.

Mínimos cuadrados (III): factorización LU

La factorización LU recibe su nombre de los dos factores que introduce: una matriz triangular inferior, (mathbf)así como uno triangular superior, (mathbf). En teoría, no hay restricciones en la descomposición LU: siempre que permitamos intercambios de filas, convirtiendo efectivamente (mathbf = mathbfmathbf) en (mathbf = mathbfmathbfmathbf) (dónde (mathbf) es una matriz de permutación), podemos factorizar cualquier matriz.

En la práctica, sin embargo, si queremos hacer uso de torch_triangular_solve() , la matriz de entrada tiene que ser simétrica. Por lo tanto, aquí también tenemos que trabajar con (mathbf^Tmathbf)no (mathbf) directamente. (Y es por eso que estoy mostrando la descomposición LU justo después de Cholesky: son similares en lo que nos hacen hacer, aunque no en absoluto similares en espíritu).

Trabajando con (mathbf^Tmathbf) significa que nuevamente estamos comenzando desde las ecuaciones normales. Nosotros factorizamos (mathbf^Tmathbf), luego resuelve dos sistemas triangulares para llegar a la solución final. Estos son los pasos, incluida la matriz de permutación que no siempre es necesaria (mathbf):

[
begin
mathbf^T mathbf mathbf &= mathbf^T mathbf \
mathbf mathbfmathbf mathbf &= mathbf^T mathbf \
mathbf mathbf &= mathbf^T mathbf^T mathbf \
mathbf &= mathbf mathbf
end
]

Vemos que cuando (mathbf) es necesario, hay un cálculo adicional: siguiendo la misma estrategia que hicimos con Cholesky, queremos mover (mathbf) de izquierda a derecha. Afortunadamente, lo que puede parecer costoso (calcular el inverso) no lo es: para una matriz de permutación, su transposición invierte la operación.

En cuanto al código, ya estamos familiarizados con la mayor parte de lo que debemos hacer. La única pieza que falta es torch_lu(). torch_lu() devuelve una lista de dos tensores, el primero una representación comprimida de las tres matrices (mathbf), (mathbf)y (mathbf). Podemos descomprimirlo usando torch_lu_unpack() :

lu <- torch_lu(AtA)

c(P, L, U) %<-% torch_lu_unpack(lu[[1]], lu[[2]])

Nos movemos (mathbf) al otro lado:

Todo lo que queda por hacer es resolver dos sistemas triangulares, y hemos terminado:

y <- torch_triangular_solve(
  Atb$unsqueeze(2),
  L,
  upper = FALSE
)[[1]]
x <- torch_triangular_solve(y, U)[[1]]

all_preds$lu <- as.matrix(A$matmul(x))
all_errs$lu <- rmse(all_preds$b, all_preds$lu)
all_errs[1, -5]
       lm   lstsq     neq    chol      lu
1 40.8369 40.8369 40.8369 40.8369 40.8369

Al igual que con la descomposición de Cholesky, podemos ahorrarnos la molestia de llamar torch_triangular_solve() dos veces. torch_lu_solve() toma la descomposición y devuelve directamente la solución final:

lu <- torch_lu(AtA)
x <- torch_lu_solve(Atb$unsqueeze(2), lu[[1]], lu[[2]])

all_preds$lu2 <- as.matrix(A$matmul(x))
all_errs$lu2 <- rmse(all_preds$b, all_preds$lu2)
all_errs[1, -5]
       lm   lstsq     neq    chol      lu      lu
1 40.8369 40.8369 40.8369 40.8369 40.8369 40.8369

Ahora, observamos los dos métodos que no requieren el cálculo de (mathbf^Tmathbf).

Mínimos cuadrados (IV): factorización QR

Cualquier matriz se puede descomponer en una matriz ortogonal, (mathbf)y una matriz triangular superior, (mathbf). La factorización QR es probablemente el enfoque más popular para resolver problemas de mínimos cuadrados; es, de hecho, el método utilizado por R’s lm(). ¿De qué manera, entonces, simplifica la tarea?

En cuanto a (mathbf)ya sabemos su utilidad: en virtud de ser triangular, define un sistema de ecuaciones que puede resolverse paso a paso, mediante mera sustitución. (mathbf) es aún mejor. Una matriz ortogonal es aquella cuyas columnas son ortogonales, es decir, los productos punto mutuos son todos cero, y tienen una norma unitaria; y lo bueno de tal matriz es que su inversa es igual a su transpuesta. En general, la inversa es difícil de calcular; la transposición, sin embargo, es fácil. Ver cómo el cálculo de un inverso – resolver (mathbf=mathbf^mathbf) – es solo la tarea central en mínimos cuadrados, es inmediatamente claro cuán importante es esto.

En comparación con nuestro esquema habitual, esto conduce a una receta ligeramente más corta. No hay una variable «ficticia» (mathbf) ya no. En su lugar, nos movemos directamente (mathbf) al otro lado, calculando la transpuesta (que es el inverso). Todo lo que queda, entonces, es la sustitución hacia atrás. Además, dado que cada matriz tiene una descomposición QR, ahora comenzamos directamente desde (mathbf) en vez de (mathbf^Tmathbf):

[
begin
mathbfmathbf &= mathbf\
mathbfmathbfmathbf &= mathbf\
mathbfmathbf &= mathbf^Tmathbf\
end
]

En torch, linalg_qr() nos da las matrices (mathbf) y (mathbf).

c(Q, R) %<-% linalg_qr(A)

En el lado derecho, solíamos tener una «variable de conveniencia» que contenía (mathbf^Tmathbf) ; aquí, nos saltamos ese paso y, en su lugar, hacemos algo «inmediatamente útil»: mover (mathbf) al otro lado.

El único paso restante ahora es resolver el sistema triangular restante.

x <- torch_triangular_solve(Qtb$unsqueeze(2), R)[[1]]

all_preds$qr <- as.matrix(A$matmul(x))
all_errs$qr <- rmse(all_preds$b, all_preds$qr)
all_errs[1, -c(5,7)]
       lm   lstsq     neq    chol      lu      qr
1 40.8369 40.8369 40.8369 40.8369 40.8369 40.8369

A estas alturas, esperará que termine esta sección diciendo «también hay un solucionador dedicado en torch/torch_linalg, a saber…”). Bueno, no literalmente, no; pero efectivamente, sí. si llamas linalg_lstsq() paso driver = "gels"se utilizará la factorización QR.

Mínimos cuadrados (V): descomposición de valores singulares (SVD)

En un verdadero orden culminante, el último método de factorización que discutimos es el más versátil, el más diversamente aplicable y el más significativo desde el punto de vista semántico: Descomposición de valores singulares (SVD). El tercer aspecto, por fascinante que sea, no se relaciona con nuestra tarea actual, por lo que no lo abordaré aquí. Aquí, lo que importa es la aplicabilidad universal: cada matriz se puede componer en componentes al estilo SVD.

Descomposición de valores singulares factoriza una entrada (mathbf) en dos matrices ortogonales, llamadas (mathbf) y (mathbf^T)y una diagonal, denominada (mathbf)tal que (mathbf = mathbf mathbf mathbf^T). Aquí (mathbf) y (mathbf^T) son los izquierda y Vectores de right singulary (mathbf) sostiene el valores singulares.

[
begin
mathbfmathbf &= mathbf\
mathbfmathbfmathbf^Tmathbf &= mathbf\
mathbfmathbf^Tmathbf &= mathbf^Tmathbf\
mathbf^Tmathbf &= mathbf\
end
]

Empezamos por obtener la factorización, usando linalg_svd(). El argumento full_matrices = FALSE dice torch que queremos un (mathbf) de dimensionalidad igual que (mathbf)no ampliado a 7588 x 7588.

c(U, S, Vt) %<-% linalg_svd(A, full_matrices = FALSE)

dim(U)
dim(S)
dim(Vt)
[1] 7588   21
[1] 21
[1] 21 21

Nos movemos (mathbf) al otro lado – una operación barata, gracias a (mathbf) siendo ortogonal.

Con ambos (mathbf^Tmathbf) y (mathbf) siendo vectores de la misma longitud, podemos usar la multiplicación por elementos para hacer lo mismo para (mathbf). Introducimos una variable temporal, ypara mantener el resultado.

Ahora queda con el sistema final para resolver, (mathbf)nuevamente nos beneficiamos de la ortogonalidad – esta vez, de la matriz (mathbf^T).

Para terminar, calculemos las predicciones y el error de predicción:

all_preds$svd <- as.matrix(A$matmul(x))
all_errs$svd <- rmse(all_preds$b, all_preds$svd)

all_errs[1, -c(5, 7)]
       lm   lstsq     neq    chol      lu     qr      svd
1 40.8369 40.8369 40.8369 40.8369 40.8369 40.8369 40.8369

Eso concluye nuestro recorrido por los algoritmos de mínimos cuadrados importantes. La próxima vez, presentaré extractos del capítulo sobre la transformada discreta de Fourier (DFT), reflejando nuevamente el enfoque en comprender de qué se trata. ¡Gracias por leer!

Foto por Pearse O’Halloran sobre Unsplash

Fuente del artículo

Deja un comentario