Efficient Triangular Backsubstitution in R

Question

I am interested in solving the linear system of equations Ax=b where A is a lower-triangular matrix (n × n) and b is a (n × 1) vector where n ≈ 600k.

I coded up backsubstitution in R and it works fast for matrices of size upto 1000, but is really slow for larger n (≈ 600k). I know that the naive backsubstitution is O(n^2).

My R function is below; does anyone know of a more efficient (vectorized, parallelized etc.) way of doing it, which scales to large n?

Backsubstitution

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

How about the R function backsolve? It calls the Level 3 BLAS routine dtrsm which is probably what you want to be doing. In general, you won't beat BLAS/LAPACK linear algebra routines: they're insanely optimized.