What is a typical use case of iterative solvers for linear systems in R?

Question

I wonder if iterative solvers are a faster way to solve linear systems (non-sparse, symmetric, positive definite).

I tried conjugate gradient methods from the R packages Rlinsolve and cPCG, but both seem to be not very accurate and slower compared to the non-iterative base::solve().

library(Rlinsolve)
library(cPCG)
library(microbenchmark)

n <- 2000
A <- tcrossprod(matrix(rnorm(n^2),nrow=n) + diag(rep(10,n)))
x <- rnorm(n)
b <- A%*%x

mean(abs(x - solve(A,b)))
## [1] 1.158205e-08
mean(abs(x - lsolve.cg(A,b)$x))
## [1] 0.03836865
mean(abs(x - cgsolve(A,b)))
## [1] 0.02642611
mean(abs(x - pcgsolve(A,b)))
## [1] 0.02638013

microbenchmark(solve(A, b), lsolve.cg(A, b),
               cgsolve(A, b), pcgsolve(A, b), times=5)
## Unit: milliseconds
##            expr       min        lq      mean    median         uq       max neval  cld
##     solve(A, b)  183.3039  188.6678  189.7362  188.8665   189.8514  197.9914     5 a   
##  lsolve.cg(A, b) 7178.7477 7784.7646 7934.8406 8114.5838 8218.7356 8377.3714     5    d
##   cgsolve(A, b) 1907.0940 2020.8368 2226.0513 2121.2917  2483.1947 2597.8393     5  b  
##  pcgsolve(A, b) 4059.5856 4109.0319 4203.4093 4242.7750  4275.9537 4329.7005     5   c 

(R version 3.6.1 with OpenBLAS and 4 cores.)

Am I missing something? What is a typical use-case for such iterative methods?

What is a good R example for non-sparse linear systems showing the advantages of iterative solvers?

Answer 1

As a creator of Rlinsolve package, I kind of disagree with what this question implicitly argues. If you have dense matrix A, all the benefits of storing sparse matrix with tailored computation disappear at once. I've seen some decent usages of sparse solvers when covariance structure under Gaussian model is banded but such literature is extremely thin.

Please keep in mind that no single tool is designed to solve every problem. If you have symmetric, positive-definite matrix, cholesky or EVD-based solution is a great tool.

FYI, I've seen Rlinsolve to be used in a statistical computing paper that compares the performance of EM-based iterative solver, which is their creation, against methods delivered from my package. I believe that serves a good role somehow.