Perché lm esaurisce la memoria mentre la moltiplicazione della matrice funziona correttamente per i coefficienti?

Question

Io sto cercando di fare effetti fissi regressione lineare con R. I miei dati sembra

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

Ho quindi deciso di fare semplicemente questo facendo yr un fattore e utilizzare lm:

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

Tuttavia, questo sembra a corto di memoria. Ho 20 livelli nel mio fattore e df è di 14 milioni di righe che richiede circa 2 GB per l'archiviazione, sto eseguendo questo su una macchina con 22 GB dedicati a questo processo.

Ho quindi deciso di provare le cose alla vecchia maniera: creare variabili dummy per ciascuno dei miei anni t1 a t20 facendo:

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

e semplicemente calcolare:

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

Questo viene eseguito senza un problema e produce la risposta quasi subito.

Sono particolarmente curioso di cosa si tratta lm che lo rende a corto di memoria quando posso calcolare bene i coefficienti? Grazie.

Answer 1

Nessuna delle risposte finora ha indicato la giusta direzione.

La risposta accettata da @idr crea confusione tra lm e summary.lm. lm non calcola statistiche diagnostiche su tutti; invece, lo fa summary.lm. Quindi sta parlando di summary.lm.

@Jake La risposta è un dato di fatto sulla stabilità numerica della fattorizzazione QR e della fattorizzazione LU/Choleksy. La risposta di Aravindakshan la espande, sottolineando la quantità di operazioni in virgola mobile dietro entrambe le operazioni (sebbene, come ha detto, non ha considerato i costi per il calcolo del prodotto a matrice incrociata). Ma non confondere i conteggi FLOP con i costi di memoria. In realtà entrambi i metodi hanno lo stesso utilizzo della memoria in LINPACK/LAPACK. In particolare, il suo argomento secondo cui il metodo QR costa più RAM per memorizzare il fattore Q è falso. La memoria compattata come spiegato in lm(): What is qraux returned by QR decomposition in LINPACK/LAPACK chiarisce in che modo la fattorizzazione QR viene calcolata e archiviata. Problema di velocità di QR v.s. Chol è dettagliato nella mia risposta: Why the built-in lm function is so slow in R?, e la mia risposta su faster lm fornisce una piccola routine lm.chol utilizzando il metodo Choleksy, che è 3 volte più veloce del metodo QR.

@Greg La risposta/suggerimento per biglm è buona, ma non risponde alla domanda. Dal momento che si parla di biglm, vorrei sottolineare che QR decomposition differs in lm and biglm. biglm calcola la riflessione della famiglia in modo che il risultante fattore R abbia diagonali positive. Vedi Cholesky factor via QR factorization per i dettagli. Il motivo per cui biglm esegue questa operazione, è che il risultante R sarà uguale al fattore Cholesky, vedere QR decomposition and Choleski decomposition in R per informazioni. Inoltre, a parte biglm, è possibile utilizzare mgcv. Leggi la mia risposta: biglm predict unable to allocate a vector of size xx.x MB per ulteriori informazioni.

---

Dopo una sintesi, è il momento di inviare la mia risposta.

Per montare un modello lineare, lm sarà

1. genera un modello di telaio;
2. genera una matrice modello;
3. chiamata lm.fit per fattorizzazione QR;
4. restituisce il risultato della fattorizzazione QR nonché il frame del modello in lmObject.

Hai detto che il frame dati di input con 5 colonne costa 2 GB per l'archiviazione. Con 20 livelli di fattore la matrice del modello risultante ha circa 25 colonne che occupano 10 GB di spazio. Ora vediamo come cresce l'utilizzo della memoria quando chiamiamo lm.

• [ambiente globale] inizialmente avete 2 GB di storage per il frame di dati;
• [lm envrionment] quindi viene copiato in un frame del modello, che costa 2 GB;
• [lm ambiente] quindi viene generata una matrice modello, che costa 10 GB;
• [lm.fit ambiente] una copia della matrice del modello viene quindi sovrascritta dalla fattorizzazione QR, che costa 10 GB;
• [lm ambiente] viene restituito il risultato di lm.fit, che costa 10 GB;
• [ambiente globale] il risultato di lm.fit è ulteriormente restituito da lm, costano altri 10 GB;
• [ambiente globale] il frame del modello viene restituito da lm, che costa 2 GB.

Quindi, è richiesto un totale di 46 GB di RAM, molto più grande della RAM disponibile di 22 GB.

In realtà se può essere "in linea" in lm, potremmo risparmiare 20 GB di costi. Ma non c'è modo di allineare una funzione R in un'altra funzione R.

Forse possiamo fare un piccolo esempio per vedere cosa succede intorno 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 

così effettivamente, X viene copiato quando passò in lm.fit. Diamo un'occhiata a ciò che ha 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 

nota che la compatta a matrice QR qrfit$qr$qr è grande come matrice modello X. Viene creato all'interno di lm.fit, ma all'uscita di lm.fit, viene copiato. Quindi, in totale, avremo 3 "copie" di X:

• quello originale in ambiente globale;
• quello copiato in lm.fit, sovrascritto dalla fattorizzazione QR;
• quello restituito da lm.fit.

Nel suo caso, X è 10 GB, quindi i costi di memoria associati lm.fit sola è già 30 GB. Per non parlare degli altri costi associati a lm.

---

D'altra parte, diamo un'occhiata a

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

X prende 10 GB, ma crossprod(X) è solo una matrice 25 * 25 e crossprod(X,y) è solo un vettore lunghezza-25. Sono così piccoli rispetto a X, quindi l'utilizzo della memoria non aumenta affatto.

Forse si è preoccupati che venga eseguita una copia locale di X quando viene chiamato crossprod? Affatto! A differenza di lm.fit che esegue sia la lettura che la scrittura su X, crossprod legge solo X, quindi non viene eseguita alcuna copia. Siamo in grado di verificare questo con la nostra matrice giocattolo X da:

tracemem(X) 
crossprod(X) 

Si vedrà alcun messaggio copia!

---

Se si desidera un breve riassunto per tutti sopra, qui è:

• costi di memoria per lm.fit(X, y) (o anche .lm.fit(X, y)) è tre volte più grande di quella per solve(crossprod(X), crossprod(X,y));
• A seconda di quanto più grande è la matrice del modello rispetto al frame del modello, i costi di memoria per lm sono 3 ~ 6 volte più grandi di quello per solve(crossprod(X), crossprod(X,y)). Il limite inferiore 3 non viene mai raggiunto, mentre il limite superiore 6 viene raggiunto quando la matrice del modello è uguale a quella del modello. Questo è il caso in cui non ci sono variabili fattore o termini "fattore-alike", come bs() e poly(), ecc

Answer 2

lm fa molto di più che trovare i coefficienti per le funzionalità di input. Ad esempio, fornisce statistiche diagnostiche che ti dicono di più sui coefficienti delle tue variabili indipendenti incluso l'errore standard e il valore t di ciascuna delle tue variabili indipendenti.

Penso che la comprensione di queste statistiche diagnostiche sia importante quando si eseguono regressioni per capire quanto sia valida la regressione.

Questi calcoli aggiuntivi fanno sì che lm sia più lento del semplice risolvere le equazioni di matrice per la regressione.

Ad esempio, utilizzando il mtcars set di dati:

>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

In aggiunta a quanto Idris ha detto, è anche la pena sottolineare che lm() non risolve per i parametri utilizzando le equazioni normali come te illustrato nella tua domanda, ma utilizza piuttosto la decomposizione QR, che è meno efficiente ma tende a produrre più soluzioni numericamente accurate.

Answer 4

Si potrebbe prendere in considerazione l'utilizzo del pacchetto biglm. Si adatta a modelli lm usando blocchi di dati più piccoli.

Answer 5

Per elaborare il punto di Jake. Diciamo che la tua regressione sta cercando di risolvere: y = Ax (A è la matrice di progettazione). Con m osservazioni e n variabili indipendenti A è una matrice mxn. Quindi il costo del QR è ~ m*n^2. Nel tuo caso sembra m = 14x10^6 e n = 20. Quindi m*n^2 = 14*10^6*400 è un costo significativo.

Tuttavia con le normali equazioni si sta tentando di invertire A'A ('indica la trasposizione), che è quadrata e della dimensione nxn. La soluzione viene solitamente eseguita utilizzando LU che costa n^3 = 8000. Questo è molto più piccolo del costo computazionale del QR. Ovviamente questo non include il costo della moltiplicazione della matrice.

Inoltre, se la routine QR tenta di memorizzare la matrice Q di dimensioni mxm=14^2*10^12 (!), La memoria sarà insufficiente. È possibile scrivere QR per non avere questo problema. Sarebbe interessante sapere quale versione di QR sia effettivamente utilizzata. E perché esattamente la chiamata lm esaurisce la memoria.