为啥 lm 内存不足而矩阵乘法适用于系数？
Posted
技术标签:
【中文标题】为啥 lm 内存不足而矩阵乘法适用于系数？【英文标题】：Why does lm run out of memory while matrix multiplication works fine for coefficients?为什么 lm 内存不足而矩阵乘法适用于系数？ 【发布时间】：2012-05-06 19:09:18 【问题描述】：我正在尝试使用 R 进行固定效应线性回归。我的数据看起来像
dte yr id v1 v2
. . . . .
. . . . .
. . . . .
然后我决定简单地通过将yr 作为一个因素并使用lm 来做到这一点：
lm(v1 ~ factor(yr) + v2 - 1, data = df)
但是，这似乎内存不足。我的因子中有 20 个级别，df 是 1400 万行，大约需要 2GB 来存储，我在一台有 22GB 专用于这个过程的机器上运行它。
然后我决定尝试老式的方法：为我的每一年创建虚拟变量 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 是什么导致它耗尽内存？谢谢。
【问题讨论】：
为什么不尝试lm.fit 而不是lm 来缩小问题范围？ lm.fit 只是通过 QR 分解或多或少地进行“原始”线性模型拟合——没有关于模型矩阵创建等的无关内容。如果您还遇到 lm.fit 的内存问题，那么 @Jake 的答案似乎是问题所在（QR 与正规方程）。
【参考方案1】：
您可能需要考虑使用biglm 包。它通过使用较小的数据块来拟合 lm 模型。
【讨论】：
【参考方案2】：到目前为止，没有一个答案指向正确的方向。
@idr 接受的答案混淆了lm 和summary.lm。 lm 根本不计算诊断统计数据；相反，summary.lm 会。所以他说的是summary.lm。
@Jake 的答案是关于 QR 分解和 LU/Choleksy 分解的数值稳定性的事实。 @ 987654323@ 的答案通过指出这两个操作背后的浮点运算量来扩展这一点（尽管正如他所说，他没有计入计算矩阵叉积的成本）。但是，不要将 FLOP 计数与内存成本混淆。实际上这两种方法在 LINPACK / LAPACK 中的内存使用情况相同。具体来说，他认为 QR 方法需要更多 RAM 来存储 Q 因子的论点是虚假的。 lm(): What is qraux returned by QR decomposition in LINPACK / LAPACK 中解释的压缩存储阐明了如何计算和存储 QR 分解。 QR vs.的速度问题Chol 在我的回答中有详细说明：Why the built-in lm function is so slow in R?，我在faster lm 上的回答提供了一个小程序lm.chol 使用 Choleksy 方法，比 QR 方法快 3 倍。
@Greg 对biglm 的回答/建议很好，但没有回答问题。由于提到了biglm，我要指出QR decomposition differs in lm and biglm。 biglm 计算户主反射，因此得到的 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 将
-
生成模型框架；
生成模型矩阵；
致电
lm.fit 进行二维码分解；
返回 QR 分解的结果以及 lmObject 中的模型框架。
您说您的输入数据框有 5 列，需要 2 GB 来存储。对于 20 个因子级别，生成的模型矩阵有大约 25 列，占用 10 GB 存储空间。现在让我们看看当我们调用lm 时内存使用量是如何增长的。
lmenvrionment] 然后将其复制到模型框架中，花费 2 GB；
[lm environment]然后生成模型矩阵，消耗 10 GB；
[lm.fit environment] 复制模型矩阵，然后通过 QR 分解覆盖，花费 10 GB；
[lm环境]返回lm.fit的结果，占用10GB；
[全局环境]lm.fit 的结果被lm 进一步返回，又消耗了 10 GB；
[全球环境]模型框架由lm返回，消耗2GB。
因此，总共需要 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 退出时，它会被复制。所以总的来说，我们将拥有X 的 3 个“副本”：
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相比，它们是如此之小，因此内存使用量根本不会增加。
也许您担心在调用crossprod 时会生成X 的本地副本？一点也不！与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 的内存成本是solve(crossprod(X), crossprod(X,y)) 的 3 ~ 6 倍。永远不会达到下限 3，而当模型矩阵与模型框架相同时，会达到上限 6。当没有因子变量或“因子相似”术语时，例如 bs() 和 poly() 等，就会出现这种情况。
【讨论】：
【参考方案3】：详细说明杰克的观点。假设您的回归试图解决：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 调用会耗尽内存。
【讨论】：
【参考方案4】：除了 idris 所说的之外，还值得指出的是 lm() 不会像您在问题中说明的那样使用正规方程来求解参数，而是使用 QR 分解，这种分解效率较低但往往会产生更精确的数值解。
【讨论】：
【参考方案5】：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
【讨论】：
是的，这是真的。但这些其他活动都不会导致 lm 内存不足以上是关于为啥 lm 内存不足而矩阵乘法适用于系数？的主要内容，如果未能解决你的问题，请参考以下文章