6

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?

Nairolf
  • 2,418
  • 20
  • 34
  • 1
    There is some useful discussion of this question at [When to use iterative methods for solving systems of linear equation](https://math.stackexchange.com/questions/1299343/when-to-use-iterative-methods-for-solving-systems-of-linear-equation). In brief, very large systems (much larger than n=2000). Also, if the system is sparse and in a pattern, some iterative methods enable one to avoid storing many of the zeroes in the matrix. – Simon Oct 03 '19 at 04:42
  • 1
    For the accuracy side, you have to set the tolerance to a lower value (argument `reltol` of `lsolve.cg` and `tol` for `cgsolve`). – nicola Oct 08 '19 at 14:34
  • @nicola good point. But this would slowdown computations even more. It would be nice to see an example where the iterative methods are faster than `base::solve()`. – Nairolf Oct 08 '19 at 16:38
  • When I tried the least squares conjugate gradient solver from the Eigen package, https://eigen.tuxfamily.org/dox/classEigen_1_1LeastSquaresConjugateGradient.html, via Rcpp and with a sparse (non-square) covariate matrix, it was orders of magnitude faster than base::solve(). So for me that was a good use case... – Tom Wenseleers Sep 01 '22 at 13:18

1 Answers1

4

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.

Kisung You
  • 66
  • 1
  • Thanks for your reply. To goal of this question is to understand the practical use of iterative solvers. For which applications would you recommend iterative solvers? Are those applications limited to overdetermined system? – Nairolf Oct 09 '19 at 15:15
  • I understand from your answer that for a non-sparse symmetric, positive-definite matrices `A` you do not recommend iterative solvers. Do you recommend iterative solvers for a sparse symmetric, positive-definite matrix `A`? Can you provide a R example, where you compare it against solves from the R package **spam** or **Matrix**? – Nairolf Oct 09 '19 at 15:19
  • Of course I agree that "no single tool is designed to solve every problem". That is exactly the reason I want to understand for what applications iterative solver are good for. – Nairolf Oct 09 '19 at 15:27
  • Can you give the DOI of the mentioned article? – Nairolf Oct 09 '19 at 15:32
  • I guess I have to be a bit careful to the author of `spam` package, I become a bit careful. A limited number of examples in statistical applications would, well in fact mostly, include sparse design matrix (which is mostly not symmetric nor positive definite but many seem to believe its Gram matrix be sparse) and algorithms for sampling from high-dimensional Gaussian random variables with sparse precision such as banded type. – Kisung You Oct 09 '19 at 17:46
  • I'm running on my desktop for solving 1-dimensional Poisson equation with varying number grid points, and it gains computational supremacy once `n` reaches about 10000. For much larger ones, `solve` can't handle the problem as `Error in asMethod(object) : Cholmod error 'problem too large' at file ../Core/cholmod_dense.c, line 105` message to discourage our efforts, while `Rlinsolve` does its job. – Kisung You Oct 09 '19 at 17:49
  • Also, **Rlinsolve** is not supertailored with much degrees of freedom. I added many checks such as detecting `Inf` or `NaN` values, checking symmetry upon argument input and so on. Therefore, direct comparison is, at this moment, not viable to see the sole effectiveness of algorithms implemented. For that part, I feel regretful not to be able to provide low-level functionalities, but this is **R** and some of emails I've received seemed to say users were satisfied with decent performance level. – Kisung You Oct 09 '19 at 17:54
  • Comparing with **spam** or **Matrix** package would be an intriguing example. Once I update the package to remove extra add-on's, let's try that part as time permits. It may take a while though as I'm in my 4th year of doctoral program - you would understand what it's like.. yeah.. – Kisung You Oct 09 '19 at 17:56
  • Finally, the paper I mentioned is a working paper which I received personally to go over. Considering they do not prefer to use arXiv, it will only be available to the public once it's published somewhere. – Kisung You Oct 09 '19 at 17:58
  • Finally again, what about making a small consent on acquiring an **inverse of sparse covariance/precision matrix** as a model problem. – Kisung You Oct 09 '19 at 17:59
  • Yes, please edit your answer and add an illustrative example with R code, which demonstrates a good use case of iterative methods. If the example features sparse matrices this is perfectly fine. – Nairolf Oct 09 '19 at 19:03