Почему встроенная функция 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 на самом деле медленно?

РЕДАКТИРОВАТЬ: Как предложил Дирк Эддельбюттель, как 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
Источник

2 ответа

Решение

Вы пропускаете это

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

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

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

Источник

Что-то не так с вашим тестом

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

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

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

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-факторизация против холесского факторизации

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 против 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 тестовое задание.

поворот против не поворот

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() или создайте гибридную "быструю, но безопасную" версию.

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

Проверьте мой ответ здесь.

Источник
Другие вопросы по тегам