¿Por qué lm se queda sin memoria mientras que la multiplicación de la matriz funciona bien para los coeficientes?

Question

que estoy tratando de hacer los efectos fijos de regresión lineal con R. Mis datos parece

dte yr id v1 v2 
    . . . . . 
    . . . . . 
    . . . . . 

Entonces decidí simplemente hacer esto haciendo yr un factor y utilizar lm:

lm(v1 ~ factor(yr) + v2 - 1, data = df) 

Sin embargo, esto parece quedarse sin memoria. Tengo 20 niveles en mi factor y df son 14 millones de filas que requieren aproximadamente 2 GB para almacenar, estoy ejecutando esto en una máquina con 22 GB dedicados a este proceso.

Entonces me decidí a probar cosas de la manera antigua: crear variables ficticias para cada una de mis años t1 a t20 haciendo:

df$t1 <- 1*(df$yr==1) 
df$t2 <- 1*(df$yr==2) 
df$t3 <- 1*(df$yr==3) 
... 

y simplemente calcular:

solve(crossprod(x), crossprod(x,y)) 

Esto funciona sin un problema y produce la respuesta casi de inmediato.

Tengo una curiosidad particular acerca de qué es lo que hace que la memoria se quede sin memoria cuando puedo calcular los coeficientes bien? Gracias.

Answer 1

Ninguna de las respuestas hasta ahora ha apuntado en la dirección correcta.

La respuesta aceptada por @idr está haciendo confusión entre lm y summary.lm. lm no se calcula ninguna estadística de diagnóstico; en cambio, summary.lm hace. Entonces él está hablando de summary.lm.

@Jake La respuesta es un hecho sobre la estabilidad numérica de la factorización QR y la factorización LU/Choleksy. La respuesta de Aravindakshan amplía esto, señalando la cantidad de operaciones de coma flotante detrás de ambas operaciones (aunque, como dijo, no contabilizó los costos para calcular el producto cruzado de la matriz). Pero no confundas los conteos de FLOP con los costos de memoria. En realidad, ambos métodos tienen el mismo uso de memoria en LINPACK/LAPACK. Específicamente, su argumento de que el método QR cuesta más RAM para almacenar el factor Q es falso. El almacenamiento compactado como se explica en lm(): What is qraux returned by QR decomposition in LINPACK/LAPACK aclara cómo se calcula y almacena la factorización QR. Emisión de velocidad de QR v.s. Chol se detalla en mi respuesta: Why the built-in lm function is so slow in R?, y mi respuesta en faster lm proporciona una pequeña rutina lm.chol utilizando el método Choleksy, que es 3 veces más rápido que el método QR.

@Greg La respuesta/sugerencia del cliente para biglm es buena, pero no responde la pregunta. Como se menciona biglm, señalaría que QR decomposition differs in lm and biglm. biglm calcula la reflexión del propietario para que el factor resultante R tenga diagonales positivas. Vea Cholesky factor via QR factorization para más detalles. La razón por la que biglm hace esto, es que el resultado R será el mismo que el factor Cholesky, consulte QR decomposition and Choleski decomposition in R para obtener información. Además, aparte de biglm, puede usar mgcv. Lea mi respuesta: biglm predict unable to allocate a vector of size xx.x MB para más.

---

Después de un resumen, es el momento de publicar mi respuesta.

Con el fin de ajustar un modelo lineal, lm voluntad

1. genera un marco de modelo;
2. genera una matriz modelo;
3. llamada lm.fit para la factorización QR;
4. devuelve el resultado de la factorización QR, así como el marco del modelo en lmObject.

Dijo que su marco de datos de entrada con 5 columnas cuesta 2 GB para almacenar. Con 20 niveles de factor, la matriz de modelo resultante tiene alrededor de 25 columnas que toman 10 GB de almacenamiento. Ahora veamos cómo crece el uso de la memoria cuando llamamos al lm.

• [entorno global] inicialmente usted tiene 2 GB de almacenamiento de la trama de datos;
• [lm envrionment] continuación, se copia a un modelo de marco, con un costo de 2 GB;
• [lm entorno] a continuación, se genera una matriz de modelo, con un costo de 10 GB;
• [lm.fit entorno] una copia de la matriz del modelo se sobrescribe luego por factorización QR, con un costo de 10 GB;
• [lm entorno] se devuelve el resultado de lm.fit, con un costo de 10 GB;
• [entorno global] el resultado de lm.fit es devuelto por lm, con un costo de otros 10 GB;
• [entorno global] el modelo de marco es devuelto por lm, con un costo de 2 GB.

Por lo tanto, se requiere un total de 46 GB de RAM, mucho mayor que la RAM de 22 GB disponible.

En realidad, si lm.fit puede ser "incluido" en lm, podríamos ahorrar 20 GB de costos. Pero no hay forma de alinear una función R en otra función R.

Tal vez podemos tomar un pequeño ejemplo para ver lo que sucede a su alrededor lm.fit:

X <- matrix(rnorm(30), 10, 3) # a `10 * 3` model matrix 
y <- rnorm(10) ## response vector 

tracemem(X) 
# [1] "<0xa5e5ed0>" 

qrfit <- lm.fit(X, y) 
# tracemem[0xa5e5ed0 -> 0xa1fba88]: lm.fit 

Así que de hecho, X se copia cuando pasó a lm.fit. Vamos a echar un vistazo a lo que tiene qrfit

str(qrfit) 
#List of 8 
# $ coefficients : Named num [1:3] 0.164 0.716 -0.912 
# ..- attr(*, "names")= chr [1:3] "x1" "x2" "x3" 
# $ residuals : num [1:10] 0.4 -0.251 0.8 -0.966 -0.186 ... 
# $ effects  : Named num [1:10] -1.172 0.169 1.421 -1.307 -0.432 ... 
# ..- attr(*, "names")= chr [1:10] "x1" "x2" "x3" "" ... 
# $ rank   : int 3 
# $ fitted.values: num [1:10] -0.466 -0.449 -0.262 -1.236 0.578 ... 
# $ assign  : NULL 
# $ qr   :List of 5 
# ..$ qr : num [1:10, 1:3] -1.838 -0.23 0.204 -0.199 0.647 ... 
# ..$ qraux: num [1:3] 1.13 1.12 1.4 
# ..$ pivot: int [1:3] 1 2 3 
# ..$ tol : num 1e-07 
# ..$ rank : int 3 
# ..- attr(*, "class")= chr "qr" 
# $ df.residual : int 7 

Tenga en cuenta que la matriz compacta QR qrfit$qr$qr es tan grande como la matriz del modelo X. Se crea dentro de lm.fit, pero al salir de lm.fit, se copia. Entonces, en total, tendremos 3 "copias" de X:

• la original en el entorno global;
• el copiado en lm.fit, sobrescrito por factorización QR;
• la devuelta por lm.fit.

En su caso, X es de 10 GB, por lo que los costes de memoria asociados con lm.fit sola ya es de 30 GB. Menos aún otros costos asociados con lm.

---

Por otra parte, vamos a echar un vistazo a

solve(crossprod(X), crossprod(X,y)) 

X tiene 10 GB, pero crossprod(X) sólo es una matriz 25 * 25, y crossprod(X,y) es sólo un vector talla 25. Son muy pequeños en comparación con X, por lo que el uso de la memoria no aumenta en absoluto.

Tal vez le preocupa que se realice una copia local de X cuando se llama crossprod? ¡De ningún modo! A diferencia de lm.fit que realiza tanto lectura como escritura en X, crossprod solo lee X, por lo que no se realiza ninguna copia. Podemos verificar esto con nuestra matriz de juguete X por:

tracemem(X) 
crossprod(X) 

Usted verá ningún mensaje de copia!

---

Si quieres un breve resumen de todo lo anterior, aquí está:

• costes de memoria para lm.fit(X, y) (o incluso .lm.fit(X, y)) es tres veces más grande que el de solve(crossprod(X), crossprod(X,y));
• Según cuánto más grande es el modelo de matriz que el modelo de marco, los costos de memoria para lm son 3 ~ 6 veces más grandes que para solve(crossprod(X), crossprod(X,y)). El límite inferior 3 nunca se alcanza, mientras que el límite superior 6 se alcanza cuando la matriz del modelo es igual que el marco del modelo. Este es el caso cuando no hay factores o variables de términos "factor de semejanza", como bs() y poly(), etc.

Answer 2

lm hace mucho más que simplemente encontrar los coeficientes para sus características de entrada. Por ejemplo, proporciona estadísticas de diagnóstico que le informan más sobre los coeficientes en sus variables independientes, incluido el error estándar y el valor t de cada una de sus variables independientes.

Creo que la comprensión de estas estadísticas de diagnóstico es importante cuando se ejecutan regresiones para comprender qué tan válida es la regresión.

Estos cálculos adicionales causan que lm sea más lento que simplemente resolver las ecuaciones matriciales para la regresión.

Por ejemplo, utilizando el mtcars conjunto de datos:

>data(mtcars) 
>lm_cars <- lm(mpg~., data=mtcars) 
>summary(lm_cars) 

Call:               
lm(formula = mpg ~ ., data = mtcars)       

Residuals:              
    Min  1Q Median  3Q  Max      
-3.4506 -1.6044 -0.1196 1.2193 4.6271      

Coefficients:             
      Estimate Std. Error t value Pr(>|t|)    
(Intercept) 12.30337 18.71788 0.657 0.5181    
cyl   -0.11144 1.04502 -0.107 0.9161    
disp   0.01334 0.01786 0.747 0.4635    
hp   -0.02148 0.02177 -0.987 0.3350    
drat   0.78711 1.63537 0.481 0.6353    
wt   -3.71530 1.89441 -1.961 0.0633 .    
qsec   0.82104 0.73084 1.123 0.2739    
vs   0.31776 2.10451 0.151 0.8814    
am   2.52023 2.05665 1.225 0.2340    
gear   0.65541 1.49326 0.439 0.6652    
carb  -0.19942 0.82875 -0.241 0.8122    
---               
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 2.65 on 21 degrees of freedom   
Multiple R-squared: 0.869,  Adjusted R-squared: 0.8066  
F-statistic: 13.93 on 10 and 21 DF, p-value: 3.793e-07  

Answer 3

Además de lo Idris dijo, también vale la pena señalar que lm() no resuelve para los parámetros utilizando las ecuaciones normales como que se ilustra en su pregunta, sino que usa la descomposición QR, que es menos eficiente pero tiende a producir más soluciones numéricamente precisas.

Answer 4

Es posible que desee considerar el uso del paquete biglm. Se adapta a modelos de lm mediante el uso de trozos de datos más pequeños.

Answer 5

Para explicar el punto de Jake. Digamos que su regresión está tratando de resolver: y = Ax (A es la matriz de diseño). Con m observaciones yn variables independientes A es una matriz mxn. Entonces el costo de QR es ~ m*n^2. En tu caso, parece que m = 14x10^6 yn = 20. Entonces m*n^2 = 14*10^6*400 es un costo significativo.

Sin embargo, con las ecuaciones normales está tratando de invertir A'A ('indica transposición), que es cuadrado y de tamaño nxn. La solución generalmente se realiza usando LU que cuesta n^3 = 8000. Esto es mucho más pequeño que el costo computacional de QR. Por supuesto, esto no incluye el costo de la matriz multiplicada.

Además, si la rutina QR intenta almacenar la matriz Q que es del tamaño mxm=14^2*10^12 (!), Entonces su memoria será insuficiente. Es posible escribir QR para no tener este problema. Sería interesante saber qué versión de QR se está usando en realidad. Y por qué exactamente la llamada lm se queda sin memoria.