当矩阵乘法对系数工作正常时,为什么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))

这没有问题,几乎立即产生答案.

我特别好奇当我能够很好地计算系数时,lm会使内存耗尽吗？谢谢.

Answer 1

除了idris所说的,还值得指出的是lm()并没有像你在问题中所说明的那样使用正规方程求解参数,而是使用QR分解,效率较低但往往产生更多的数值精确度解决方案.

Answer 2

您可能需要考虑使用biglm包.它通过使用较小的数据块来适应lm模型.

Answer 3

到目前为止,答案都没有指出正确的方向.

@idr接受的答案是在lm和之间混淆summary.lm.lm根本不计算诊断统计数据 ; 相反,summary.lm确实如此.所以他在谈论summary.lm.

@Jake的回答是关于QR分解和LU/Choleksy分解的数值稳定性的事实.Aravindakshan的答案通过指出两个操作背后的浮点运算量来扩展这一点(尽管如他所说,他没有计算计算矩阵交叉产品的成本).但是,不要将FLOP计数与内存成本混淆.实际上,这两种方法在LINPACK/LAPACK中具有相同的内存使用量.具体来说,他认为QR方法花费更多RAM来存储Q因子的说法是假的.压缩存储如lm()中所述:LINPACK/LAPACK中QR分解返回的qraux是如何计算和存储QR分解的.QR和Chol的速度问题详见我的答案:为什么内置的lm功能在R中如此慢？我的回答更快,lm提供了一个lm.chol使用Choleksy方法的小程序,比QR方法快3倍.

@Greg的答案/建议biglm很好,但它没有回答这个问题.既然biglm提到了,我会指出QR分解在lm和biglm.biglm计算住户反映,使得结果R因子具有正对角线.有关详细信息,请参阅QR分解的Cholesky因子.这样biglm做的原因是结果R将与Cholesky因子相同,参见R中的QR分解和Choleski分解以获取信息.此外,除了biglm,你可以使用mgcv.阅读我的回答:biglm预测无法分配大小为xx.x MB的向量.

总结之后,是时候发布我的答案了.

为了适应线性模型,lm将会

1. 生成模型框架;
2. 生成模型矩阵;
3. 要求lm.fitQR分解;
4. 返回QR分解的结果以及模型框架lmObject.

您说您的5列输入数据帧需要2 GB才能存储.使用20个因子级别,得到的模型矩阵有大约25列,占用10 GB存储空间.现在让我们看一下调用时内存使用量的增长情况lm.

• [全局环境]最初你有2 GB的数据存储空间;
• [ lmenvrionment]然后将其复制到模型框架,成本为2 GB;
• [ lm环境]然后生成一个模型矩阵,成本为10 GB;
• [ 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 

确实,X在传入时会被复制lm.fit.让我们来看看什么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.更别说与之相关的其他费用了lm.

另一方面,让我们来看看

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

X需要10 GB,但crossprod(X)只是一个25 * 25矩阵,crossprod(X,y)只是一个长度为25的向量.与它们相比它们非常小X,因此内存使用量根本没有增加.

也许你担心X在crossprod被叫时会制作本地副本？一点也不!与lm.fit执行读取和写入不同X,crossprod只读取X,因此不进行复制.我们可以通过我们的玩具矩阵验证这一点X:

tracemem(X)
crossprod(X)

你将看不到复制信息!

如果您想要上述所有内容的简短摘要,请点击此处:

• lm.fit(X, y)(或甚至.lm.fit(X, y))的记忆成本是其三倍solve(crossprod(X), crossprod(X,y));
• 根据模型矩阵比模型框架大多少,内存成本lm是3~6倍solve(crossprod(X), crossprod(X,y)).从未达到下限3,而当模型矩阵与模型框架相同时,达到上限6.出现这种情况时,有没有因子变量或"因子相似"条款,如bs()和poly()等.

Answer 4

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 5

详细阐述杰克的观点.假设您的回归试图解决:y = Ax(A是设计矩阵).m个观测值和n个独立变量A是mxn矩阵.然后QR的成本是〜m*n^2.在你的情况下,它看起来像m = 14x10 ^ 6和n = 20.因此m*n^2 = 14*10^6*400成本很高.

然而,对于正常的方程,您试图反转A'A('表示转置),这是方形和大小nxn.解决方案通常使用LU来完成n^3 = 8000.这远小于QR的计算成本.当然,这不包括矩阵乘法的成本.

此外,如果QR例程试图存储大小为mxm=14^2*10^12(!)的Q矩阵,那么你的记忆力将不足.可以写QR来解决这个问题.了解实际使用的QR版本会很有趣.为什么这个lm电话会耗尽内存.