r - R 中的高效三角反代换

Question

我对求解线性方程组 Ax=b 很感兴趣，其中 A 是下三角矩阵 (n × n)，b 是 (n × 1) 向量，其中 n ≈ 600k。

我在 R 中编写了反代换代码，它对大小高达 1000 的矩阵运行速度很快，但对于较大的 n (≈ 600k) 确实很慢。我知道朴素的反向代入是 O(n^2)。

我的 R 函数如下；有谁知道更有效(矢量化、并行化等)的方法，可以扩展到大 n？

回代

backsub=function(X,y)
{ 
 l=dim(X)
 n=l[1]  
 p=l[2]
 y=as.matrix(y)
    
   for (j in seq(p,1,-1))
    {  
      y[j,1]=y[j,1]/X[j,j]
      if((j-1)>0)
         y[1:(j-1),1]=y[1:(j-1),1]-(y[j,1]*X[1:(j-1),j])
    }  
    return(y)
}

Answer 1

R 函数 backsolve 怎么样？它调用 Level 3 BLAS routine dtrsm这可能是你想要做的。通常，您不会击败 BLAS/LAPACK 线性代数例程:它们经过疯狂优化。