sanic (version 0.0.2)

solve_cg: Solve Systems of Equations Using Iterative Methods

Description

Iterative solvers using the Conjugate Gradient method for sparse systems of equations \(Ax = b\). Three different types are available: (1) stabilized bi-conjugate gradient (BiCGSTAB) for square matrices, (2) conjugate gradient for rectangular least-squares (LSCG), and (3) classic conjugate gradient (CG) for symmetric positive definite matrices.

Usage

solve_cg(
  a,
  b,
  x0,
  type = c("BiCGSTAB", "LSCG", "CG"),
  iter,
  tol,
  precond = 1L,
  verbose = FALSE
)

Value

Solves for \(x\) and returns a numeric matrix with the results.

Arguments

a

Square numeric matrix with the coefficients of the linear system. Dense and sparse matrices are supported, but the format must be sparse (see sparsify). Dense matrices are coerced automatically.

b

Numeric vector or matrix at the right-hand side of the linear system. If missing, 'b' is set to an identity matrix and 'a' is inverted.

x0

Numeric vector or matrix with an initial guess. Must be of the same dimension as 'b'.

type

Character scalar. Whether to use the BiCGSTAB, least squares CG or classic CG method.

iter

Integer scalar with the maximum number of iterations. Defaults to the theoretical maximum, i.e. the number of columns in 'a'.

tol

Numeric scalar with the desired tolerance. Defaults to the machine precision.

precond

Integer scalar indicating the type of preconditioner to be used. Defaults to diagonal preconditioning. See the Details for further information.

verbose

Logical scalar. Whether to print iterations and tolerance.

Details

Preconditioners can be set to 0 for no / identity preconditioning, 1 (default) for Jacobi / diagonal preconditioning, or 2 for incomplete factorisation. Not all schemes are available for every type:

* type = "BiCGSTAB" The default is precond = 1 for diagonal preconditioning. Set precond = 0 for no preconditioning, or precond = 2 for an incomplete LUT preconditioner. * type = "LSCG" The default is precond = 1 for diagonal least squares preconditioning. Set precond = 0 for no preconditioning. * type = "CG" The default is precond = 1 for diagonal preconditioning. Set precond = 0 for no preconditioning, or precond = 2 for an incomplete Cholesky preconditioner.

Examples

set.seed(42)
x <- rnorm(3)

# Solve via BiCGSTAB for square matrices
A <- matrix(rnorm(9), nrow = 3, ncol = 3)
b <- A %*% x
norm(solve_cg(A, b, type = "B") - x)

# Solve via LSCG for rectangular matrices
A <- matrix(rnorm(12), nrow = 4, ncol = 3)
b <- A %*% x
norm(solve_cg(A, b, type = "LS") - x)

# Solve via classic CG for symmetric matrices
A <- crossprod(matrix(rnorm(9), nrow = 3, ncol = 3))
b <- A %*% x
norm(solve_cg(A, b, type = "CG") - x)

# The input matrix A should always be in sparse format
A <- sparsify(crossprod(matrix(rnorm(9), nrow = 3, ncol = 3)))
# The right-hand side should be a dense matrix
b <- as.matrix(A %*% x)

# We can check the speed of convergence and quality directly
solve_cg(A, b, verbose = TRUE)
# And provide guesses as starting value
solve_cg(A, b, x0 = x, verbose = TRUE)