Warum ist die integrierte LM-Funktion in R so langsam?

Question

Ich dachte immer, dass die lm Funktion war extrem schnell in R, aber wie dieses Beispiel würde vorschlagen, die geschlossene Lösung berechnet mit der solve Funktion ist viel schneller.

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)) 

Kann mir jemand erklären, wenn dieses Spielzeug Beispiel ein schlechtes Beispiel ist, oder es ist der Fall, dass lm eigentlich langsam?

EDIT: Wie von Dirk Eddelbuettel vorgeschlagen, wie lm Bedürfnisse der Formel zu lösen, ist der Vergleich unfair, so besser lm.fit zu verwenden, die nicht die Formel

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 

Answer 1

zu lösen braucht

Sie sind mit Blick auf die

• solve() liefert nur Ihre Parameter
• lm() kehren Sie ein (sehr reich) Objekt mit viele Komponenten für die nachfolgende Analyse, Inferenz, Plots, ...
• die Hauptkosten Ihres lm() Anruf nicht der Vorsprung aber die Auflösung der Formel y ~ ., aus dem die Modellmatrix aufgebaut werden muss.

Zur Veranschaulichung RCPP wir einige Varianten einer Funktion fastLm() mehr zu tun, was geschrieben hat lm() tut (dh etwas mehr als lm.fit() von der Basis R) und gemessen. Siehe z.B. this benchmark script was deutlich zeigt, dass die dominierende Kosten für kleinere Datensätze in der Analyse der Formel und Erstellen der Modellmatrix ist.

Kurz gesagt, Sie tun die Richtige Sache mit Benchmarking, aber Sie tun es nicht so richtig beim Versuch zu vergleichen, was meist unvergleichbar ist: eine Teilmenge mit einer viel größeren Aufgabe.

Answer 2

etwas falsch mit Ihrem Benchmarking

Es ist wirklich erstaunlich, dass niemand dies beobachtet!

Sie haben t(X) %*% X innerhalb solve() verwendet. Sie sollten crossprod(X) verwenden, da X'X eine symmetrische Matrix ist. crossprod() berechnet nur die Hälfte der Matrix, während der Rest kopiert wird. %*% erzwingt die Berechnung aller. So wird crossprod() zweimal schneller sein. Dies erklärt, warum Sie in Ihrem Benchmarking zwischen solve() und lm.fit() ungefähr das gleiche Timing haben.

Auf meiner alten Intel Nahalem (2008 Intel Core 2 Duo), ich habe:

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 

Wenn Ihr Computer schneller, versuchen X <- matrix(runif(2000*2000),2000) stattdessen verwenden.

Im Folgenden werde ich die Rechendetails bei allen Anpassungsmethoden erläutern.

---

QR-Faktorisierung vs. Cholesky Faktorisierung

lm()/lm.fit() ist QR basiert, während solve() ist Cholesky basiert. Die rechnerischen Kosten der QR-Zerlegung sind 2 * n * p^2, während die Cholesky-Methode n * p^2 + p^3 (n * p^2 für die Berechnung des Matrix-Kreuzprodukts, p^3 für die Cholesky-Zerlegung) ist. So können Sie sehen, dass, wenn n ist viel viel größer als p, Cholesky Methode ist etwa 2-mal schneller als QR-Methode. Es besteht also keine Notwendigkeit, hier zu benchmarken. (, falls Sie nicht wissen, n ist die Anzahl der Daten, p die Anzahl der Parameter ist.)

---

LINPACK QR v.s. LAPACK QR

Typischerweise lm.fit() Nutzungen (modifiziert) LINPACK QR-Faktorisierung Algorithmus, anstatt LAPACK QR-Faktorisierung Algorithmus. Vielleicht sind Sie nicht sehr vertraut mit BLAS/LINPACK/LAPACK; Dies sind FORTRAN-Codes, die kerntechnische Matrixberechnungen liefern. LINPACK ruft level-1 BLAS auf, während LAPACK level-3 BLAS unter Verwendung von Blockalgorithmen aufruft. Im Durchschnitt ist LAPACK QR 1,6 Mal schneller als LINPACK QR. Der kritische Grund, dass lm.fit()LAPACK Version nicht verwendet, ist die Notwendigkeit für partielle Spaltenschwenkung.LAPACK Version macht volle Spalte schwenken, so dass summary.lm() schwieriger R Matrixfaktor der QR-Zerlegung zu F-Statistik und ANOVA Test zu produzieren.

---

schwenken vs. kein Schwenken

fastLm() von RcppEigen Paket verwendet LAPACK nicht drehbar QR-Faktorisierung. Auch hier können Sie sich über den QR-Faktorisierungsalgorithmus und die Schwenkprobleme nicht im Klaren sein. Sie müssen nur wissen, dass QR-Faktorisierung mit Pivotierung nur 50% Anteil von Level-3 BLAS, während QR-Faktorisierung ohne Pivotierung 100% Anteil von Level-3 BLAS. In diesem Zusammenhang wird das Aufgeben der Schwenkbewegung den QR-Faktorisierungsprozess beschleunigen. Sicherlich ist das Endergebnis anders, und wenn die Modellmatrix einen Rangmangel aufweist, führt kein Pivotieren zu einem gefährlichen Ergebnis.

Es gibt eine gute Frage zu dem anderen Ergebnis, das Sie von fastLM erhalten: Why does fastLm() return results when I run a regression with one observation?. @BenBolker, @DirkEddelbuettel und ich hatten eine sehr kurze Diskussion in Kommentaren von Bens Antwort.

---

Fazit: Haben Sie Geschwindigkeit wollen, oder numerische Stabilität?

Im Hinblick auf die numerische Stabilität gibt es:

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

In Bezug auf die Geschwindigkeit, gibt es:

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

Wie Dirk sagte

FWIW die RcppEigen Paket hat eine umfangreichere Reihe von Zerlegungen in seinem fastLm() Beispiel. Aber hier ist es, wie Ben so eloquent sagte: "Dies ist ein Teil des Preises, den Sie für die Geschwindigkeit zahlen." Wir geben dir genug Seil, um dich aufzuhängen. Wenn Sie sich vor sich selbst schützen möchten, bleiben Sie bei lm() oder lm.fit() oder erstellen Sie eine hybride "Fast-but-Safe" -Version.

---

Schnelle und stabile Version

prüfen my answer Here.