Eine Kleinste Quadrate (unterbestimmtes System) schnell lösen (R)

Question

Ich habe ein Programm in R, das eine große Anzahl von Lösungen der kleinsten Quadrate berechnet (> 10.000: typischerweise 100.000+) und nach dem Profiling sind die aktuellen Engpässe für das Programm. Ich habe eine Matrix A mit Spaltenvektoren, die Spanning-Vektoren und eine Lösung b entsprechen. Ich versuche, nach der kleinsten Quadrate Lösung x von Ax=b zu lösen. Die Matrizen sind typischerweise 4xj groß - viele von ihnen sind nicht quadratisch (j < 4) und so sind allgemeine Lösungen für unterbestimmte Systeme das, wonach ich suche.

Die Hauptfrage: Was ist der schnellste Weg, ein unterbestimmtes System in R zu lösen? Ich habe viele Lösungen, die die Normal Equation verwenden, aber ich suche nach einer Routine in R, die schneller als jede der folgenden Methoden ist.

Zum Beispiel: das System für x durch Ax = b gegeben Lösen die folgenden Bedingungen gegeben:

• Das System ist nicht notwendig bestimmt [Regel unter bestimmten] (NcoI (A) < = Länge (b) hält immer). Somit funktioniert solve(A,b) nicht, weil Solve eine quadratische Matrix benötigt.
• Sie können davon ausgehen, dass t(A) %*% A (entspricht crossprod(A)) ist nicht singulär - geprüft wird, früher im Programm
• Sie jedes Paket in R frei verfügbar verwenden können
• Die Lösung ist nicht schön sein muss - es muss nur selten schnell
• Eine obere Schranke für Größe A ist ziemlich 10x10 und Null-Elemente auftreten werden - A ist in der Regel ziemlich dicht

Zwei Zufallsmatrizen zum Testen ...

A = matrix(runif(12), nrow = 4) 
b = matrix(runif(4), nrow = 4) 

Alle folgenden Funktionen wurden profiliert. Sie sind hier wiedergegeben:

f1 = function(A,b) 
{ 
    solve(t(A) %*% A, t(A) %*% b) 
} 
f2 = function(A,b) 
{ 
    solve(crossprod(A), crossprod(A, b)) 
} 
f3 = function(A,b) 
{ 
    ginv(crossprod(A)) %*% crossprod(A,b) # From the `MASS` package 
} 
f4 = function(A,b) 
{ 
    matrix.inverse(crossprod(A)) %*% crossprod(A,b) # From the `matrixcalc` package 
} 
f5 = function(A,b) 
{ 
    qr.solve(crossprod(A), crossprod(A,b)) 
} 
f6 = function(A,b) 
{ 
    svd.inverse(crossprod(A)) %*% crossprod(A,b) 
} 
f7 = function(A,b) 
{ 
    qr.solve(A,b) 
} 
f8 = function(A,b) 
{ 
    Solve(A,b) # From the `limSolve` package 
} 

Nach dem Test f2 ist der aktuelle Gewinner. Ich habe auch lineare Modellmethoden getestet - sie waren lächerlich langsam angesichts all der anderen Informationen, die sie produzieren. Der Code wurde mit der folgenden Profil:

library(ggplot2) 
library(microbenchmark) 

all.equal(
    f1(A,b), 
    f2(A,b), 
    f3(A,b), 
    f4(A,b), 
    f5(A,b), 
    f6(A,b), 
    f7(A,b), 
    f8(A,b), 
     ) 

compare = microbenchmark(
    f1(A,b), 
    f2(A,b), 
    f3(A,b), 
    f4(A,b), 
    f5(A,b), 
    f6(A,b), 
    f7(A,b), 
    f8(A,b), 
    times = 1000) 

autoplot(compare) 

Answer 1

Wie wäre es Rcpp?

library(Rcpp) 
cppFunction(depends='RcppArmadillo', code=' 
      arma::mat fRcpp (arma::mat A, arma::mat b) { 
      arma::mat betahat ; 
      betahat = (A.t() * A).i() * A.t() * b ; 
      return(betahat) ; 
      }         
      ') 

all.equal(f1(A, b), f2(A, b), fRcpp(A, b)) 
#[1] TRUE 
microbenchmark(f1(A, b), f2(A, b), fRcpp(A, b)) 
#Unit: microseconds 
#  expr min  lq  mean median  uq  max neval 
# f1(A, b) 55.110 57.136 67.42110 59.5680 63.0120 160.873 100 
# f2(A, b) 34.444 37.685 43.86145 39.7120 41.9405 117.920 100 
# fRcpp(A, b) 3.242 4.457 7.67109 8.1045 8.9150 39.307 100