r语言 - 如何快速求解最小二乘法(欠定系统)

Question

我在R中有一个程序，它正在计算大量的最小二乘解（>10，000：通常为100，000+），在分析之后，这些是程序的当前瓶颈。我有一个矩阵A，其中包含对应于跨越向量的列向量和一个解决方案b.我正在尝试求解Ax=b的最小二乘解x。矩阵的大小通常为 4xj - 其中许多不是正方形（j <4），因此我正在寻找对不确定系统的通用解决方案。

主要问题：在 R 中求解未确定系统的最快方法是什么？我有许多利用正态方程的解决方案，但我正在寻找比以下任何方法都快的 R 例程。

例如：求解给定以下约束条件的Ax = b给出的x系统：

• 系统不是必要的确定[通常未确定]（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)

Answer 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