11

我正在尝试使用 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 是什么导致内存不足？谢谢。

4

5 回答 5

10

除了 idris 所说的之外，还值得指出的是 lm() 不会像您在问题中说明的那样使用正规方程来求解参数，而是使用 QR 分解，这种分解效率较低，但往往会产生更准确的数值解决方案。

于 2012-04-26T08:12:39.867 回答
10

到目前为止，没有一个答案指向正确的方向。

@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 vs Chol的速度问题在我的回答中有详细说明：为什么内置的 lm 函数在 R 中这么慢？，而我对fasterlm的回答提供了一个使用Choleksy方法的小程序lm.chol，比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 的数据框存储空间；
  • [lm环境]然后将其复制到模型框架中，花费2 GB；
  • [lm环境]然后生成模型矩阵，消耗10GB；
  • [lm.fit环境]复制模型矩阵，然后通过 QR 分解覆盖，花费 10 GB；
  • [lm环境]返回结果lm.fit，消耗10GB；
  • [全局环境]的结果由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执行 read 和 write to不同X，crossprod只执行 reads 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()等）时就是这种情况。
于 2016-09-04T17:43:49.643 回答
9

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

于 2012-04-26T13:33:08.433 回答
7

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
于 2012-04-26T04:48:23.020 回答
5

详细说明杰克的观点。假设您的回归试图解决：（y = AxA 是设计矩阵）。具有 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 会很有趣。以及为什么lm调用内存不足。

于 2012-08-12T19:58:21.487 回答