nk027 / sanic

Fast routines for solving large systems of linear equations in R. Makes Eigen Cholesky-, LU-, QR-, and iterative (Conjugate Gradient, BiCGSTAB) solvers for both dense and sparse problems available.

sanic: Solving Ax = b Nimbly in C++

Routines for solving large systems of linear equations in R. Direct and iterative solvers from the Eigen C++ library are made available. Solvers include Cholesky, LU, QR, and Krylov subspace methods (Conjugate Gradient, BiCGSTAB). Both dense and sparse problems are supported.

Installation

sanic is available on CRAN. The development version can be installed from GitHub.

install.packages("sanic")
devtools::install_github("nk027/sanic")

Usage

To solve a linear system of equations, use the solve2() function for automatic dispatch to a specific solver or access the LU, QR, Cholesky or Conjugate Gradient solvers directly.

To solve an eigenproblem, use the eigen2() or svd2() functions, or the arnoldi() function.

Solvers

Solver	Function	Notes	Sparse	Reference
LU decomposition	solve_lu()	Partial pivoting, full pivoting	Yes	1, 2, 3
Householder QR decomposition	solve_qr()	Column pivoting, full pivoting, no pivoting	Yes	1, 2, 3, 4
Cholesky decomposition	solve_chol()	LDLT for semidefinite problems, LLT for positive definite problems	Yes	1, 2 3, 4
Conjugate Gradient (CG)	solve_cg()	Biconjugate gradient stabilised (BiCGTAB) for square problems, least squares (LSCG) for rectangular problems, classic CG for symmetric positive definite problems, preconditioners	Always	1, 2, 3

Eigenproblems

Solver	Function	Notes	Sparse	Reference
Spectral decomposition	eigen2()	Square and symmetric problems	No	1, 2
Singular value decomposition	svd2()	Bidiagonal Divide and Conquer SVD for large and Jacobi SVD for small problems	No	1, 2
Arnoldi iteration	arnoldi()	Square problems using an iteratively constructed Hessenberg matrix	Always	1
Lanczos algorithm	lanczos()	Symmetric problems using an iteratively constructed tridiagonal matrix	Always	1

sanic's Issues

Add support for more (sparse) matrix formats

Could start with more from Matrix, but also possibly around spam?

Return failed CG guess

Currently, CG estimates that fail to converge lead to an error. It might make sense to change this to a warning, and return the latest guess -- especially if the number of iterations is capped. Happens here:

sanic/src/solve_iterative.cpp

Line 39 in d64bf0b

Rcpp::stop("Solving failed.");

Patch non-symmetric matrices

QR and SVD should work anyways, but it seems buggy atm.
Consider checking symmetry.

Consider supporting eigenvalue problems

Also iteratively, Arnoldi and Lanzcos should be implementable.

Allow more interfacing

Atm a matrix is fed in and a result is fed out. Would be interesting to cache decompositions and interface more directly with Eigen.