Type:	Package
Title:	Solving Ax = b Nimbly in C++
Version:	0.0.2
Date:	2023-08-22
Author:	Nikolas Kuschnig [aut, cre], Lukas Vashold [ctb], Yixuan Qiu [ctb]
Maintainer:	Nikolas Kuschnig <nikolas.kuschnig@wu.ac.at>
Description:	Routines for solving large systems of linear equations and eigenproblems 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). Dense and sparse problems are supported.
URL:	https://github.com/nk027/sanic
BugReports:	https://github.com/nk027/sanic/issues
Depends:	R (≥ 3.3.0)
Imports:	Rcpp (≥ 1.0.5), Matrix, methods
License:	GPL-3 | file LICENSE
Encoding:	UTF-8
LinkingTo:	Rcpp, RcppEigen
RoxygenNote:	7.2.3
NeedsCompilation:	yes
Packaged:	2023-08-23 10:23:09 UTC; nikolas
Repository:	CRAN
Date/Publication:	2023-08-23 16:00:07 UTC

Krylov Subspace Spectral Decomposition

Description

Arnoldi iteration and Lanczos method to iteratively approximate the Hessenberg or tridiagonal form of a matrix A and find its eigenvalues.

Usage

arnoldi(
  a,
  b,
  symmetric,
  iter = nrow(a),
  tol = .Machine$double.eps,
  eigen = TRUE,
  orthogonalise = TRUE
)

lanczos(
  a,
  b,
  iter = nrow(a),
  tol = .Machine$double.eps,
  eigen = TRUE,
  orthogonalise = TRUE
)

Arguments

Value

Returns a list with slots "H" for the Hessenberg form of 'a' or slots "diagonal" and "subdiagonal" for its triangular form, slot "Q" with the orthonormal basis, and, if requested, eigenvalues in the slot "values".

Examples

set.seed(42)
# Compute Hessenberg of a square matrix
A <- matrix(rnorm(9), nrow = 3, ncol = 3)
ks <- arnoldi(A, symmetric = FALSE)

# Compute tridiagonal of a symmetric matrix
A <- crossprod(matrix(rnorm(9), nrow = 3, ncol = 3))
ks <- lanczos(A)
ks <- arnoldi(A, symmetric = TRUE) # Short-hand

Spectral Decomposition

Description

Solvers for eigenproblems around the matrix A. Compute eigenvalues \lambda and eigenvectors v of A, such that Av = \lambda v.

Usage

eigen2(a, symmetric, vectors = TRUE)

Arguments

Value

Solves the eigenproblem and returns a list with eigenvalues in the "values" slot and, if requested, eigenvectors in the slot "vectors".

Examples

set.seed(42)
# Compute eigenvalues and eigenvectors for a square matrix
A <- matrix(rnorm(9), nrow = 3, ncol = 3)
ev <- eigen2(A, symmetric = FALSE)

# Compute eigenvalues and eigenvectors for a symmetric matrix
A <- crossprod(matrix(rnorm(9), nrow = 3, ncol = 3))
ev <- eigen2(A, symmetric = TRUE)
# Check reconstruction
norm(A %*% ev$vectors - ev$vectors %*% diag(ev$values))

Solving Ax = b Nimbly in C++

Description

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.

Solve Systems of Equations

Description

Solve systems of equations Ax = b using an automatically chosen direct method (see solve_chol). Methods are chosen for speed at reasonable accuracy. Please choose a suitable method manually if numerical stability is the main consideration.

Usage

solve2(a, b, ...)

Arguments

Value

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

Examples

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

# Solve using a general matrix
A <- matrix(rnorm(9), nrow = 3, ncol = 3)
b <- A %*% x
norm(solve2(A, b) - x)

# Solve using a symmetric matrix
A <- crossprod(matrix(rnorm(9), nrow = 3, ncol = 3))
b <- A %*% x
norm(solve2(A, b) - x)

# Solve using a square matrix
A <- matrix(rnorm(12), nrow = 4, ncol = 3)
b <- A %*% x
norm(solve2(A, b) - x)

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
)

Arguments

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.

Value

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

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)

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)

Arguments

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.

Value

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

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))

Transform a Matrix to Be Sparse.

Description

Concise function to transform dense to sparse matrices of class dgCMatrix (see sparseMatrix).

Usage

sparsify(x)

Arguments

Value

Returns 'x' as dgCMatrix.

Examples

sparsify(matrix(rnorm(9L), 3L))

Singular Value Decomposition

Description

Solvers for generalized eigenproblems around the matrix A. Compute singular values \Sigma, left singular vectors U and right singular vectors V of A, such that A = U \Sigma V^*. Two different types are available: (1) bidiagonal divide and conquer strategy (BDC) SVD, and (2) two-sided Jacobi SVD for small matrices (<16) and high accuracy.

Usage

svd2(a, type = c("BDC", "Jacobi"), vectors = TRUE, thin = TRUE)

Arguments

Value

Solves the generalised eigenproblem and returns a list with singular values in the "d" component and, if requested, singular vectors in the components "u" and "v".

Examples

set.seed(42)
# Compute singular values and vectors using BDC
A <- matrix(rnorm(9), nrow = 3, ncol = 3)
sv <- svd2(A)

# Compute singular values using Jacobi
A <- matrix(rnorm(9), nrow = 3, ncol = 3)
sv <- svd2(A, type = "J", vectors = FALSE)

# Compute singular values and full vectors using BDC
A <- matrix(rnorm(12), nrow = 4, ncol = 3)
sv <- svd2(A, type = "B", thin = FALSE)
A <- matrix(rnorm(12), nrow = 3, ncol = 4)
sv <- svd2(A, type = "B", thin = FALSE)