Почему у lm заканчивается память, а умножение матриц прекрасно работает для коэффициентов?

Я пытаюсь сделать фиксированные эффекты линейной регрессии с R. Мои данные выглядят как

<code>dte   yr   id   v1   v2
  .    .    .    .    .
  .    .    .    .    .
  .    .    .    .    .
</code>

Затем я решил просто сделать это, сделавyr фактор и использованиеlm:

<code>lm(v1 ~ factor(yr) + v2 - 1, data = df)
</code>

Однако, это, кажется, исчерпывает память. У меня есть 20 уровней в моем факторе иdf Это 14 миллионов строк, для хранения которых требуется около 2 ГБ, я запускаю это на машине с 22 ГБ, выделенной для этого процесса.

Затем я решил попробовать вещи старомодным способом: создать фиктивные переменные для каждого из моих летt1 вt20 при выполнении:

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

и просто вычислить:

<code>solve(crossprod(x), crossprod(x,y))
</code>

Это работает без проблем и дает ответ почти сразу.

Мне особенно любопытно, что же такое в lm, что заставляет его исчерпать память, когда я могу просто вычислить коэффициенты? Благодарю.

 Ben Bolker26 апр. 2012 г., 17:54
почему бы тебе не попробоватьlm.fit вместоlm сузить проблему?lm.fit просто делает более-менее «сырой» подгонка линейной модели с помощью QR-разложения - ничего лишнего о создании матрицы модели и т. д. Если у вас также возникают проблемы с памятьюlm.fit, тогда может показаться, что ответ @ Jake является проблемой (QR против нормальных уравнений).

Ответы на вопрос(5)

что сказал idris, также стоит отметить, что lm () не решает для параметров, используя нормальные уравнения, как вы проиллюстрировали в своем вопросе, а скорее использует разложение QR, которое менее эффективно, но имеет тенденцию производить больше численно точные решения.

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       

ваша регрессия пытается решить: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.Of course this doesn't include the cost of the matrix multiply.

Further if the QR routine tries to store the Q matrix which is of size mxm=14^2*10^12 (!), then your memory will be insufficient. Можно написать QR, чтобы не иметь этой проблемы. Было бы интересно узнать, какая версия QR на самом деле используется. И почему именноlm звонок заканчивается из памяти.

Решение Вопроса

Принятый ответ от@idr делает путаницу междуlm а такжеsummary.lm. lm computes no diagnostic statistics at all; вместо,summary.lm делает. Так он говорит оsummary.lm.

@Джейкответ является фактом о числовой стабильности QR-факторизации и LU / Choleksy факторизации.AravindakshanОтвет расширяет это, указывая количество операций с плавающей запятой, стоящих за обеими операциями (хотя, по его словам, он не учитывал затраты на вычисление матричного перекрестного продукта). Но не путайте подсчеты FLOP со стоимостью памяти. На самом деле оба метода имеют одинаковое использование памяти в LINPACK / LAPACK. В частности, его аргумент, что метод QR стоит больше оперативной памяти для храненияQ фактор является поддельным. Уплотненное хранилище, как описано вlm (): что такое qraux, возвращаемое разложением QR в LINPACK / LAPACK разъясняет, как QR-факторизация вычисляется и сохраняется. Скорость выдачи QR v.s. Чол подробно изложен в моем ответе:Почему встроенная функция lm так медленно работает в R?и мой ответ наБыстрееlm обеспечивает небольшую рутинуlm.chol с использованием метода Холекса, который в 3 раза быстрее, чем метод QR.

@Gregответ / предложение дляbiglm это хорошо, но это не отвечает на вопрос. посколькуbiglm упоминается, я хотел бы отметить, чтоQR-разложение отличаетсяlm and biglm. biglm вычисляет отражение домохозяина так, чтобы в результатеR фактор имеет положительные диагонали. УвидетьФактор Холецкого через QR-факторизацию для деталей. Причина того, чтоbiglm это то, что в результатеR будет таким же, как фактор Холецкого, см.QR-разложение и разложение Холецкого в R для информации. Кроме того, кромеbiglm, ты можешь использоватьmgcv, Прочитайте мой ответ:biglm predict unable to allocate a vector of size xx.x MB для большего.

After a summary, it is time to post my answer.

Чтобы соответствовать линейной модели,lm будут

generates a model frame; generates a model matrix; call lm.fit for QR factorization; returns the result of QR factorization as well as the model frame in lmObject.

Вы сказали, что ваш входной фрейм данных с 5 столбцами стоит 2 ГБ для хранения. При 20 факторных уровнях полученная матрица модели имеет около 25 столбцов, занимающих 10 ГБ памяти. Теперь давайте посмотрим, как растет использование памяти, когда мы вызываемlm.

[global environment] initially you have 2 GB storage for the data frame; [lm envrionment] then it is copied to a model frame, costing 2 GB; [lm environment] then a model matrix is generated, costing 10 GB; [lm.fit environment] a copy of model matrix is made then overwritten by QR factorization, costing 10 GB; [lm environment] the result of lm.fit is returned, costing 10 GB; [global environment] the result of lm.fit is further returned by lm, costing another 10 GB; [global environment] the model frame is returned by lm, costing 2 GB.

Таким образом, требуется всего 46 ГБ ОЗУ, что намного больше, чем доступная 22 ГБ ОЗУ.

На самом деле, еслиlm.fit может быть "встроенным" вlmМы могли бы сэкономить 20 ГБ затрат. Но нет способа встроить 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:

the original one in global environment; the one copied into lm.fit, the overwritten by QR factorization; the one returned by lm.fit.

В твоем случае,X составляет 10 ГБ, поэтому затраты памяти связаны сlm.fit в одиночку уже 30 гб. Не говоря уже о других расходах, связанных сlm.

С другой стороны, давайте посмотрим на

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

X занимает 10 гб, ноcrossprod(X) это только25 * 25 матрица иcrossprod(X,y) просто вектор длины 25 Они такие крошечные по сравнению сXпри этом использование памяти вообще не увеличивается.

Может быть, вы беспокоитесь, что локальная копияX будет сделано, когдаcrossprod называется? Не за что! В отличие отlm.fit который выполняет чтение и запись вX, crossprod только читаетX, так что копия не сделана. Мы можем проверить это с нашей игрушечной матрицейX от:

tracemem(X)
crossprod(X)

Вы не увидите копирующего сообщения!

If you want a short summary for all above, here it is:

memory costs for lm.fit(X, y) (or even .lm.fit(X, y)) is three times as large as that for solve(crossprod(X), crossprod(X,y)); Depending on how much larger the model matrix is than the model frame, memory costs for lm is 3 ~ 6 times as large as that for solve(crossprod(X), crossprod(X,y)). The lower bound 3 is never reached, while the upper bound 6 is reached when the model matrix is as same as the model frame. This is the case when there is no factor variables or "factor-alike" terms, like bs() and poly(), etc.

biglm пакет. Он подходит для моделей с использованием меньших порций данных.

Ваш ответ на вопрос