gpt4 book ai didi

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

转载 作者:行者123 更新时间:2023-12-04 01:05:34 27 4
gpt4 key购买 nike

我在 R 中有一个程序,它正在计算大量的最小二乘解(> 10,000:通常为 100,000+),并且在分析之后,这些是程序当前的瓶颈。我有一个矩阵 A具有对应于生成向量的列向量和解决方案 b .我正在尝试求解最小二乘解 xAx=b .矩阵的大小通常为 4xj - 其中许多不是正方形(j < 4),因此我正在寻找针对欠定系统的通用解决方案。

主要问题:在 R 中解决欠定系统的最快方法是什么?我有许多利用 Normal Equation 的解决方案但我正在 R 中寻找比以下任何方法都快的例程。

例如:
解决系统为xAx = b 提供鉴于以下限制:

  • 系统没有必要确定[通常欠确定](ncol (A) <= length(b) 总是成立)。因此solve(A,b)不起作用,因为 solve 需要一个方阵。
  • 您可以假设 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)

    最佳答案

    怎么样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

    关于r - 如何快速求解最小二乘法(欠定系统)?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/27674866/

    27 4 0
    Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
    广告合作:1813099741@qq.com 6ren.com