Risolvere rapidamente un minimo (sistema sottodeterminato) (R)

Question

Ho un programma in R che sta calcolando una grande quantità di soluzioni dei minimi quadrati (> 10.000: tipicamente 100.000+) e, dopo la profilazione, questi sono gli attuali colli di bottiglia per il programma. Ho una matrice A con i vettori di colonna che corrispondono ai vettori di espansione e una soluzione b. Sto tentando di risolvere la soluzione dei minimi quadrati x di Ax=b. Le matrici sono in genere di dimensioni 4xj, molte delle quali non sono quadrate (j < 4) e quindi le soluzioni generali per i sistemi sottodeterminati sono ciò che sto cercando.

La domanda principale: qual è il modo più veloce per risolvere un sistema sottodeterminato in R? Ho molte soluzioni che utilizzano il Normal Equation ma sto cercando una routine in R che sia più veloce di uno qualsiasi dei metodi seguenti.

Ad esempio: risolvere il sistema di x in Ax = b dato i seguenti vincoli:

• Il sistema non è necessario determinare [solitamente sotto-determinato] (ncol (A) < = lunghezza (b) tiene sempre). Pertanto, solve(A,b) non funziona perché la risoluzione richiede una matrice quadrata.
• si può supporre che t(A) %*% A (equivalente a crossprod(A)) è non singolare - si è verificata in precedenza nel programma
• È possibile utilizzare qualsiasi pacchetto liberamente disponibili in R
• La soluzione non deve essere abbastanza - ha solo bisogno per essere veloce
• un limite superiore alla dimensione del A è ragionevolmente 10x10 e zero elementi si verificano di rado - A è di solito abbastanza densa

due matrici casuali per i test ...

01.235.164,106174 millions

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

Tutte le funzioni seguenti sono state profilate. Essi sono riprodotti qui:

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 
} 

Dopo la prova, f2 è il vincitore corrente. Ho anche provato i metodi del modello lineare - erano ridicolmente lenti a causa di tutte le altre informazioni che producevano. Il codice è stato profilato utilizzando la seguente:

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

Come su 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