Szybkie rozwiązywanie problemu z najmniejszymi kwadratami (niedookreślonym systemem) (R)

Question

Mam program w R, który oblicza dużą liczbę rozwiązań najmniejszych kwadratów (> 10 000: zwykle 100 000+), a po profilowaniu to są obecne wąskie gardła dla program. Mam matrycę A z wektorami kolumnowymi, które odpowiadają wektorom spinającym i rozwiązaniom b. Próbuję rozwiązać dla rozwiązania najmniejszych kwadratów x z Ax=b. Macierze mają zwykle rozmiar 4xj - wiele z nich nie jest kwadratowych (j < 4), więc szukam ogólnych rozwiązań dla niedookreślonych systemów.

Główne pytanie: Jaki jest najszybszy sposób rozwiązania niedocenianego systemu w R? Mam wiele rozwiązań, które używają Normal Equation, ale szukam rutyny w R, która jest szybsza niż którąkolwiek z poniższych metod.

Na przykład: rozwiązuje system x podanej przez Ax = b otrzymuje następujące ograniczenia:

• systemu nie jest to konieczne określana [zazwyczaj poniżej ustalonej] (Ncol (A) < = długość (B) zawsze trzyma). Zatem solve(A,b) nie działa, ponieważ rozwiązanie wymaga kwadratowej macierzy.
• Można założyć, że t(A) %*% A (odpowiednik crossprod(A)) jest nieosobliwe - to sprawdzone wcześniej w programie
• Można użyć dowolnego pakietu swobodnie dostępne w R
• Rozwiązanie nie musi być ładny - po prostu potrzebuje być szybki
• górna związana od wielkości A jest rozsądnie 10x10 i zerowe elementy występują rzadko - A jest zwykle dość gęsty

Dwie macierze losowe do testowania ...

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

Wszystkie poniższe funkcje zostały profilowane. Są tu przytoczony:

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 
} 

Po testach, f2 jest obecny zwycięzca. Przetestowałem również metody modeli liniowych - były absurdalnie powolne, biorąc pod uwagę wszystkie inne informacje, które produkują. Kod został wyprofilowany przy użyciu następującego:

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

Co powiecie na 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