Как быстро решить наименьших квадратов (недоопределенная система)?

У меня есть программа на R, которая вычисляет большое количество решений наименьших квадратов (> 10000: обычно 100, 000+), и после профилирования это текущие узкие места для программы. У меня есть матрица A с векторами столбцов, которые соответствуют охватывающим векторам и решению b. Я пытаюсь найти решение для наименьших квадратов x из Ax=b. Матрицы типично имеют размер 4xj - многие из них не квадратные (j <4), и поэтому я ищу общие решения для недостаточно определенных систем.

Главный вопрос: Какой самый быстрый способ решить недостаточно определенную систему в R? У меня есть много решений, в которых используется нормальное уравнение, но я ищу подпрограмму в R, которая работает быстрее, чем любой из методов ниже.

Например: Решите систему для x заданного Ax = b учетом следующих ограничений:

  • Система необязательно определяется [обычно ncol (A) <= length(b) ] (ncol (A) <= length(b) всегда выполняется). Таким образом, solve(A,b) не работает, потому что для решения требуется квадратная матрица.
  • Вы можете предположить, что t(A) %*% A (эквивалентно crossprod(A)) не является единичным - это проверено ранее в программе
  • Вы можете использовать любой пакет, свободно доступный в R
  • Решение не должно быть красивым - оно должно быть быстрым
  • Верхняя граница размера A разумно равна 10x10, и нулевые элементы встречаются редко - A обычно довольно плотный

Две случайные матрицы для тестирования...

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

Все функции ниже были профилированы. Они воспроизведены здесь:

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
}

После тестирования f2 - текущий победитель. Я также проверил методы линейной модели - они были смехотворно медленными, учитывая всю другую информацию, которую они производят. Код был профилирован с использованием следующего:

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)

Ответы

Ответ 1

Как насчет 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