Почему встроенная функция lm настолько медленна в R?

Я всегда считал, что функция lm была чрезвычайно быстрой в R, но, как показывает этот пример, закрытое решение, рассчитанное с использованием функции solve, выполняется быстрее.

data<-data.frame(y=rnorm(1000),x1=rnorm(1000),x2=rnorm(1000))
X = cbind(1,data$x1,data$x2)

library(microbenchmark)
microbenchmark(
solve(t(X) %*% X, t(X) %*% data$y),
lm(y ~ .,data=data))

Может кто-нибудь объяснить мне, является ли этот пример игрушкой плохим примером или это так, что lm на самом деле медленный?

EDIT: Как предложил Дирк Эддельбуэттель, поскольку lm нужно решить формулу, сравнение несправедливо, поэтому лучше использовать lm.fit, который не нуждается в разрешении формулы

microbenchmark(
solve(t(X) %*% X, t(X) %*% data$y),
lm.fit(X,data$y))


Unit: microseconds
                           expr     min      lq     mean   median       uq     max neval cld
solve(t(X) %*% X, t(X) %*% data$y)  99.083 108.754 125.1398 118.0305 131.2545 236.060   100  a 
                      lm.fit(X, y) 125.136 136.978 151.4656 143.4915 156.7155 262.114   100   b

Ответ 1

Вы не замечаете, что

solve() возвращает только ваши параметры
lm() возвращает вам (очень богатый) объект со многими компонентами для последующего анализа, вывода, графиков,...
Основная стоимость вашего вызова lm() не, но разрешение формулы y ~ ., из которой должна быть построена матрица модели.

Чтобы проиллюстрировать Rcpp, мы написали несколько вариантов функции fastLm(), выполняя больше того, что lm() делает (т.е. немного больше, чем lm.fit() из базы R) и измеряет его. См. этот тест script, который ясно показывает, что доминирующая стоимость для меньших наборов данных заключается в разборе формулы и построении матрицы модели.

Короче говоря, вы делаете правильную вещь, используя бенчмаркинг, но вы делаете это не так, чтобы попытаться сравнить то, что в основном несравнимо: подмножество с гораздо большей задачей.

Ответ 2

Что-то не так с вашим бенчмаркингом

Неудивительно, что никто этого не заметил!

Вы использовали t(X) %*% X внутри solve(). Вы должны использовать crossprod(X), поскольку X'X - симметричная матрица. crossprod() вычисляет только половину матрицы при копировании остальных. %*% заставляет вычислить все. Таким образом, crossprod() будет в два раза быстрее. Это объясняет, почему в вашем бенчмаркинге у вас примерно одинаковое время между solve() и lm.fit().

На моем старом Intel Nahalem (2008 Intel Core 2 Duo) у меня есть:

X <- matrix(runif(1000*1000),1000)
system.time(t(X)%*%X)
#    user  system elapsed 
#   2.320   0.000   2.079 
system.time(crossprod(X))
#    user  system elapsed 
#    1.22    0.00    0.94

Если ваш компьютер работает быстрее, попробуйте вместо него использовать X <- matrix(runif(2000*2000),2000).

В дальнейшем я объясню детали вычисления, связанные со всеми методами подгонки.

QR-факторизация v.s. Факторизация Cholesky

lm()/lm.fit() основан на QR, а solve() - на основе Холески. Расчетные затраты QR-декомпозиции 2 * n * p^2, тогда как метод Холецкого n * p^2 + p^3 (n * p^2 для вычисления матричного кросс-продукта, p^3 для разложения Холецкого). Таким образом, вы можете видеть, что когда n намного больше, чем p, метод Холески примерно в 2 раза быстрее QR-метода. Таким образом, здесь нет никакой необходимости оценивать здесь. (в случае, если вы не знаете, n - это число данных, p - количество параметров.)

LINPACK QR v.s. LAPACK QR

Как правило, lm.fit() использует (модифицированный) LINPACK QR алгоритм факторизации, а не LAPACK QR-алгоритм факторизации. Может быть, вы не очень хорошо знакомы с BLAS/LINPACK/LAPACK; это код FORTRAN, обеспечивающий ядерные математические вычисления ядра. LINPACK вызывает BLAS уровня 1, а LAPACK вызывает уровень-3 BLAS с использованием блок-алгоритмов. В среднем LAPACK QR в 1,6 раза быстрее, чем LINPACK QR. Критическая причина, по которой lm.fit() не использует версию LAPACK, - это необходимость частичного поворота столбца. LAPACK версия выполняет полный поворот столбцов, что затрудняет использование summary.lm() R матричный коэффициент QR-декомпозиции для получения F-статистики и теста ANOVA.

поворот v.s. no pivoting

fastLm() из RcppEigen пакет использует LAPACK невращающуюся QR-факторизацию. Опять же, вы можете быть неясно о алгоритме факторизации QR и о поворотах. Вам нужно только знать, что QR-факторизация с поворотом имеет только 50% доли уровня 3 BLAS, тогда как QR-факторизация без поворота имеет 100% долю уровня 3 BLAS. В связи с этим отказ от поворота ускорит процесс факторизации QR. Разумеется, конечный результат отличается, и когда модельная матрица имеет ранг недостаточной, никакой поворот не дает опасного результата.

Существует хороший вопрос, связанный с другим результатом, который вы получаете от fastLM: Почему fastLm() возвращает результаты при выполнении регрессии с одним наблюдением?, @BenBolker, @DirkEddelbuettel, и у меня была очень короткая дискуссия в комментариях к ответу Бена.

Вывод: вам нужна скорость или численная стабильность?

В терминах численной устойчивости существует:

LINPACK pivoted QR > LAPACK pivoted QR > pivoted Cholesky > LAPACK non-pivoted QR

В терминах скорости есть:

LINPACK pivoted QR < LAPACK pivoted QR < pivoted Cholesky < LAPACK non-pivoted QR

Как сказал Дирк,

FWIW пакет RcppEigen имеет более полный набор разложений в примере fastLm(). Но вот как Бен так красноречиво заявил: "Это часть цены, которую вы платите за скорость". Мы даем вам достаточно веревки, чтобы повесить себя. Если вы хотите защитить себя от себя, придерживайтесь lm() или lm.fit(), или создайте гибридную "быструю, но безопасную" версию.

Быстрая и стабильная версия

Отметьте мой ответ здесь.