爲什麼lm內存不足，而矩陣乘法對係數工作正常？

Question

我試圖做固定效應與R.線性迴歸我的數據看起來像

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

我再決定簡單地通過使yr一個因素做到這一點，使用lm：

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

然而，這似乎耗盡內存。我在我的因素中有20個級別，而df是1400萬行，這需要大約2GB的存儲空間，我在22 GB專用於此過程的計算機上運行此操作。

然後我決定嘗試的東西的老式方法：通過做t1到t20創建虛擬變量，我的每個歲：

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

，簡單地計算：

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

此運行而不一個問題，並幾乎馬上產生答案。

我特別好奇，當我可以計算係數就好了時，它會使內存耗盡嗎？謝謝。

Answer 1

迄今沒有任何答案指出了正確的方向。

@idr接受的答案正在混淆lm和summary.lm。 lm根本不計算診斷統計量;相反，summary.lm呢。所以他在談論summary.lm。

@Jake的答案是關於QR因子分解和LU/Choleksy因子分解的數值穩定性的一個事實。 Aravindakshan的答案通過指出兩個操作背後浮點操作的數量來擴展這個問題（儘管如他所說，他沒有計算計算矩陣交叉乘積的成本）。但是，不要混淆FLOP計數和內存成本。實際上這兩種方法在LINPACK/LAPACK中具有相同的內存使用情況。具體來說，他認爲QR方法花費更多內存來存儲Q因子是一個虛假的因素。如lm(): What is qraux returned by QR decomposition in LINPACK/LAPACK中所解釋的壓縮存儲器闡明瞭如何計算和存儲QR分解。 QR v.s.的速度問題Chol在我的回答中有詳細說明：Why the built-in lm function is so slow in R?，而我在faster lm上的回答提供了一個使用Choleksy方法的小例程lm.chol，它比QR方法快3倍。

@Greg對biglm的回答/建議很好，但它沒有回答這個問題。由於提到了biglm，我會指出QR decomposition differs in lm and biglm。 biglm計算householder反射，以便生成的R因子具有正對角線。詳情請參閱Cholesky factor via QR factorization。 biglm這樣做的原因是由此產生的R將與Cholesky因子相同，請參閱QR decomposition and Choleski decomposition in R以獲取信息。此外，除了biglm，您可以使用mgcv。閱讀我的答案：biglm predict unable to allocate a vector of size xx.x MB瞭解更多信息。

---

經過總結，是時候發佈我的回答。

爲了適合的線性模型，lm將

1. 生成模型框架;
2. 生成模型矩陣;
3. 致電lm.fit QR分解;
4. 返回QR分解的結果以及lmObject中的模型框架。

你說你的5列輸入數據幀需要2 GB的存儲空間。在20個因子水平下，最終的模型矩陣大約有25列，存儲容量爲10 GB。現在讓我們看看當我們調用lm時內存使用情況如何增長。

• [全球環境]最初，您有2 GB的數據存儲空間;
• [lm envrionment]然後它被複制到模型框架，花費2 GB;
• [lm環境]然後生成模型矩陣，成本爲10GB;
• [lm.fit環境]複製模型矩陣，然後用QR分解重寫，成本10 GB;
• [lm環境]返回結果lm.fit，成本10 GB;
• [全球環境]lm.fit的結果進一步由lm返回，花費另外10 GB;
• [全球環境]模型框架由lm返回，成本爲2 GB。

因此，總共需要46 GB RAM，遠遠大於可用的22 GB RAM。

其實如果lm.fit可以「內聯」爲lm，我們可以節省20 GB的成本。但是無法在另一個R函數中內聯R函數。

也許我們可以舉一個小例子，看看周圍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 

所以事實上，當傳遞到lm.fitX被複制。讓我們來看看什麼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 

注意，緊湊QR矩陣qrfit$qr$qr是模型矩陣X一樣大。它在lm.fit內部創建，但在lm.fit的退出時被複制。所以總共我們將有3個「副本」X：

• 原來在全球環境中的一個;
• 一個複製到lm.fit，被QR分解覆蓋;
• 由lm.fit返回的那一個。

在你的情況下，X是10 GB，所以與lm.fit單獨相關的內存成本已經是30 GB。更Let論與lm相關的其他費用。

---

在另一方面，讓我們來看看

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

X需要10 GB，但crossprod(X)只是25 * 25矩陣，crossprod(X,y)只是一個長度爲25載體。與X相比，它們非常小，因此內存使用量根本不增加。

也許你擔心當crossprod被調用時會產生本地副本X？一點也不！與執行X的讀取和寫入的lm.fit不同，crossprod僅讀取X，因此不進行復制。我們可以通過我們的玩具矩陣X驗證這一點：

tracemem(X) 
crossprod(X) 

你會看到沒有複製的消息！

---

如果你想上面所有的簡短總結，那就是：爲lm.fit(X, y)（甚至.lm.fit(X, y)）

• 內存成本是對於solve(crossprod(X), crossprod(X,y))三倍大;
• 根據模型矩陣比模型框架大多少，lm的存儲器成本是solve(crossprod(X), crossprod(X,y))的3〜6倍。下限3永遠不會達到，而當模型矩陣與模型框架相同時達到上限6。出現這種情況時，有沒有因子變量或「因子都」而言，像bs()和poly()等

Answer 2

lm不僅僅是爲您的輸入功能找到係數。例如，它提供的診斷統計數據可以告訴您更多關於自變量的係數，包括每個自變量的標準誤差和t值。

我認爲理解這些診斷統計信息在運行迴歸分析以瞭解迴歸的有效性時很重要。

這些額外的計算結果會導致lm比簡單地求解迴歸的矩陣方程要慢。

例如，使用mtcars數據集：

>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

除了什麼伊德里斯說，這也是值得指出的是LM（）就像你在說明沒有解決使用正規方程的參數你問題，而是使用QR分解，效率較低但傾向於產生更多數值精確的解決方案。

Answer 4

您可能要考慮使用biglm包。它通過使用更小的數據塊適合lm模型。

Answer 5

詳細說明傑克的觀點。假設您的迴歸試圖解決：y = Ax（A是設計矩陣）。有m個觀測值和n個獨立變量A是一個m×n矩陣。然後QR的費用是〜m*n^2。在你的情況下，它看起來像m = 14x10^6和n = 20。所以m*n^2 = 14*10^6*400是一個重要的成本。

但是，使用正規方程式，您正嘗試將A'A（'表示轉置），它是正方形，大小爲nxn。解決方案通常使用成本爲n^3 = 8000的LU完成。這比QR的計算成本小得多。 當然這不包括矩陣乘法的代價。

此外，如果QR例程試圖存儲大小爲mxm=14^2*10^12（！）的Q矩陣，那麼您的內存將不足。可以寫QR來解決這個問題。知道什麼版本的QR實際上正在使用會很有趣。爲什麼lm呼叫的內存不足。