sanic (version 0.0.2)

solve_chol: Solve Systems of Equations Using Direct Methods

Description

Direct solvers using Cholesky, LU, or QR decompositions for systems of equations \(Ax = b\). Dense or sparse methods are used depending on the format of the input matrix (see sparsify).

Usage

solve_chol(a, b, pivot = 1L, ordering = 0L)
solve_lu(a, b, pivot = 1L, ordering = 1L)
solve_qr(a, b, pivot = 1L, ordering = 1L)

Value

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

Arguments

a

Square numeric matrix with the coefficients of the linear system. Both dense and sparse matrices are supported (see sparsify).

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.

pivot

Integer scalar indicating the pivoting scheme to be used. Defaults to partial pivoting. See the Details for further information.

ordering

Integer scalar indicating the ordering scheme to be used. See the Details for further information.

Details

Pivoting schemes for dense matrices can be set to 0 for no pivoting, 1 (default) for partial pivoting, or 2 for full pivoting. Not all schemes are available for every decomposition:

* solve_chol() The default is pivot = 1 for the robust LDLT decomposition of \(A\), such that \(A = P'LDL^*P\). For the LDLT \(A\) needs to be positive or negative semidefinite. Set pivot = 0 for the plain LLT decomposition of \(A\), such that \(A = LL^* = U^*U\). For the LLT \(A\) needs to be positive definite and preferably numerically stable. * solve_lu() The default is pivot = 1 for the partial pivoting LU decomposition of \(A\), such that \(A = PLU\). For this scheme \(A\) needs to be invertible and preferably numerically stable. Set pivot = 2 for the complete pivoting LU decomposition of \(A\), such that \(A = P^{-1}LUQ^{-1}\). This scheme is applicable to square matrices, rank-revealing, and stable. solve_qr() The default is pivot = 1 for the column pivoting Householder QR decomposition of \(A\), such that \(AP = QR\). This scheme is generally applicable, rank-revealing, and stable. Set pivot = 2 for the full pivoting Householder QR decomposition of \(A\), such that \(PAP' = QR\). This scheme is generally applicable, rank-revealing, and optimally stable. Set pivot = 0 for an unpivoted Householder QR decomposition of \(A\), such that \(A = QR\). This scheme is generally applicable, but not as stable as pivoted variants.

Ordering schemes for sparse matrices can be set to 0 for approximate minimum degree (AMD) ordering, 1 for column approximate minimum degree (COLAMD) ordering, or 2 for natural ordering. Not all orderings are available for every decomposition:

* solve_chol() The default is ordering = 0 for AMD ordering. Set ordering = 2 for natural ordering. * solve_lu() The default is ordering = 1 for COLAMD ordering. Set ordering = 0 for AMD or ordering = 2 for natural ordering. * solve_qr() The default is ordering = 1 for COLAMD ordering. Set ordering = 0 for AMD or ordering = 2 for natural ordering.

Examples

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

# Solve via QR for general matrices
A <- matrix(rnorm(12), nrow = 4, ncol = 3)
b <- A %*% x
norm(solve_qr(A, b) - x)

# Solve via LU for square matrices
A <- matrix(rnorm(9), nrow = 3, ncol = 3)
b <- A %*% x
norm(solve_lu(A, b) - x)

# Solve via Cholesky for symmetric matrices
A <- crossprod(matrix(rnorm(9), nrow = 3, ncol = 3))
b <- A %*% x
norm(solve_chol(A, b) - x)

# Sparse methods are available for the 'dgCMatrix' class from Matrix
A <- crossprod(matrix(rnorm(9), nrow = 3, ncol = 3))
b <- A %*% x
norm(solve_qr(sparsify(A), b))
norm(solve_lu(sparsify(A), b))
norm(solve_chol(sparsify(A), b))