cpp11armadillo
An 'Armadillo' Interface

Package index
Search the cpp11armadillo package
Vignettes
Functions
12
Source code
28
Man pages
5
Home
/
CRAN
/
cpp11armadillo
/
Decompositions, factorisations, inverses and equation solvers (dense matrices)

Decompositions, factorisations, inverses and equation solvers (dense matrices)
In cpp11armadillo: An 'Armadillo' Interface 

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

This vignette is adapted from the official Armadillo documentation.

Cholesky decomposition of symmetric matrix {#chol}

Cholesky decomposition of symmetric (or hermitian) matrix X into a triangular matrix R, with an optional permutation vector (or matrix) P. By default, R is upper triangular.

Usage:

chol(R, P, X, layout, output)

// form 1
mat R = chol(X) // chol(X, "upper") or chol(X, "lower") also work

// form 2
chol(R, X) // chol(R, X, "upper") or chol(R, X, "lower") also work

// form 3
chol(R, P, X, "upper", "vector") // chol(R, P, X, "lower", "vector") also work

// form 4
chol(R, P, X, "upper", "matrix") // chol(R, P, X, "lower", "matrix") also work

The optional argument layout is either "upper" or "lower", which specifies whether R is upper or lower triangular.

Forms 1 and 2 require X to be positive definite.

Forms 3 and 4 require X to be positive semi-definite. The pivoted decomposition provides a permutation vector or matrix P with type uvec or umat.

The decomposition has the following form:

Caveats:

Examples

[[cpp11::register]] list chol1_(const doubles_matrix<>& x,
                                const char* layout,
                                const char* output) {
  mat X = as_mat(x);

  mat Y = X.t() * X;

  mat R;
  umat P;

  writable::list out(2);
  bool ok = chol(R, P, Y, layout, output);
  out[0] = writable::logicals({ok});
  out[1] = as_doubles_matrix(R);

  return out;
}

Eigen decomposition of a symmetric matrix {#eig_sym}

Eigen decomposition of dense symmetric (or hermitian) matrix X into eigenvalues eigval and eigenvectors eigvec.

Usage:

vec eigval = eig_sym(X)

eig_sym(eigval, X)
eig_sym(eigval, eigvec, X, "dc") // eig_sym(eigval, eigvec, X, "std") also works

The eigenvalues and corresponding eigenvectors are stored in eigval and eigvec, respectively. The eigenvalues are in ascending order. The eigenvectors are stored as column vectors.

The optional argument method is either "dc" or "std", which specifies the method used for the decomposition. The divide-and-conquer method (dc) provides slightly different results than the standard method (std), but is considerably faster for large matrices.

If X is not square sized, an error is thrown.

If the decomposition fails:

Examples

[[cpp11::register]] list eig_sym1_(const doubles_matrix<>& x,
                                   const char* method) {
  mat X = as_mat(x);

  vec eigval;
  mat eigvec;

  bool ok = eig_sym(eigval, eigvec, X, method);

  writable::list out(3);
  out[0] = writable::logicals({ok});
  out[1] = as_doubles(eigval);
  out[2] = as_doubles_matrix(eigvec);

  return out;
}

Eigen decomposition of a general matrix {#eig_gen}

Eigen decomposition of dense general (non-symmetric/non-hermitian) square matrix X into eigenvalues eigval and eigenvectors eigvec.

Usage:

cx_vec eigval = eig_gen(X, bal)

eig_gen(eigval, X, bal)
eig_gen(eigval, eigvec, X, bal)
eig_gen(eigval, leigvec, reigvec, X, bal)

The eigenvalues and corresponding right eigenvectors are stored in eigval and eigvec, respectively. If both left and right eigenvectors are requested, they are stored in leigvec and reigvec, respectively. The eigenvectors are stored as column vectors.

The optional argument balance is either "balance" or "nobalance", which specifies whether to diagonally scale and permute X to improve conditioning of the eigenvalues. The default operation is "nobalance".

If X is not square sized, an error is thrown.

If the decomposition fails:

Examples

[[cpp11::register]] list eig_gen1_(const doubles_matrix<>& x,
                                   const char* balance) {
  mat X = as_mat(x);

  cx_vec eigval;
  cx_mat eigvec;

  bool ok = eig_gen(eigval, eigvec, X, balance);

  writable::list out(3);
  out[0] = writable::logicals({ok});
  out[1] = as_complex_doubles(eigval);
  out[2] = as_complex_matrix(eigvec);

  return out;
}

Eigen decomposition for a pair of general matrices {#eig_pair}

Eigen decomposition for pair of general dense square matrices A and B of the same size, such that A * eigvec = B * eigvec * diagmat(eigval). The eigenvalues and corresponding right eigenvectors are stored in eigval and eigvec, respectively. If both left and right eigenvectors are requested, they are stored in leigvec and reigvec, respectively. The eigenvectors are stored as column vectors.

Usage:

cx_vec eigval = eig_pair(A, B)

eig_pair(eigval, A, B)
eig_pair(eigval, eigvec, A, B)
eig_pair(eigval, leigvec, reigvec, A, B)

If A or B is not square sized, an error is thrown.

If the decomposition fails:

Examples

[[cpp11::register]] list eig_pair1_(const doubles_matrix<>& a,
                                    const doubles_matrix<>& b) {
  mat A = as_mat(a);
  mat B = as_mat(b);

  cx_vec eigval;
  cx_mat eigvec;

  bool ok = eig_pair(eigval, eigvec, A, B);

  writable::list out(3);
  out[0] = writable::logicals({ok});
  out[1] = as_complex_doubles(eigval);
  out[2] = as_complex_matrix(eigvec);

  return out;
}

Upper Hessenberg decomposition {#hess}

Upper Hessenberg decomposition of square matrix X, such that X = U * H * U.t(). U is a unitary matrix containing the Hessenberg vectors. H is a square matrix known as the upper Hessenberg matrix, with elements below the first subdiagonal set to zero.

Usage:

mat H = hess(X)

hess(U, H, X)
hess(H, X)

If X is not square sized, an error is thrown.

If the decomposition fails:

Caveat: in general, upper Hessenberg decomposition is not unique.

Examples

[[cpp11::register]] list hess1_(const doubles_matrix<>& x) {
  mat X = as_mat(x);

  mat H;
  bool ok = hess(H, X);

  writable::list out(2);
  out[0] = writable::logicals({ok});
  out[1] = as_doubles_matrix(H);

  return out;
}

Inverse of a general matrix {#inv}

Inverse of general square matrix A.

Usage:

mat B = inv(A)
mat B = inv(A, settings)

inv(B, A)
inv(B, A, settings)
inv(B, rcond, A)

The settings argument is optional, it is one of the following:

The reciprocal condition number is optionally calculated and stored in rcond:

If A is not square sized, an error is thrown.

If A appears to be singular:

Caveats:

To solve a system of linear equations, such as Z = inv(X) * Y, solve() can be faster and/or more accurate.

Examples

[[cpp11::register]] list inv1_(const doubles_matrix<>& a) {
  mat A = as_mat(a);

  mat B;
  bool ok = inv(B, A, inv_opts::allow_approx);

  writable::list out(2);
  out[0] = writable::logicals({ok});
  out[1] = as_doubles_matrix(B);

  return out;
}

Inverse of a symmetric positive definite matrix {#inv_sympd}

Inverse of symmetric/hermitian positive definite matrix A.

Usage:

mat B = inv_sympd(A)
mat B = inv_sympd(A, settings)

inv_sympd(B, A)
inv_sympd(B, A, settings)
inv_sympd(B, rcond, A)

The settings argument is optional, it is one of the following:

The reciprocal condition number is optionally calculated and stored in rcond:

If A is not square sized, an error is thrown.

If A is not symmetric positive definite, an error is thrown.

If A appears to be singular:

Caveats:

Examples

[[cpp11::register]] list inv_sympd1_(const doubles_matrix<>& a) {
  mat A = as_mat(a);

  mat B;
  bool ok = inv_sympd(B, A, inv_opts::allow_approx);

  writable::list out(2);
  out[0] = writable::logicals({ok});
  out[1] = as_doubles_matrix(B);

  return out;
}

Lowe-upper decomposition with partial pivoting {#lu}

Lower-upper decomposition (with partial pivoting) of matrix X.

Usage:

// form 1
lu(L, U, P, X)

// form 2
lu(L, U, X)

The first form provides a lower-triangular matrix L, an upper-triangular matrix U, and a permutation matrix P, such that P.t() * L * U = X.

The second form provides permuted L and U, such that L * U = X. Note that in this case L is generally not lower-triangular.

If the decomposition fails:

Examples

[[cpp11::register]] list lu1_(const doubles_matrix<>& x) {
  mat X = as_mat(x);

  mat L, U, P;

  bool ok = lu(L, U, P, X);

  writable::list out(4);
  out[0] = writable::logicals({ok});
  out[1] = as_doubles_matrix(L);
  out[2] = as_doubles_matrix(U);
  out[3] = as_doubles_matrix(P);

  return out;
}

Find the orthonormal basis of the null space of a matrix {#null}

Find the orthonormal basis of the null space of matrix A.

Usage:

mat B = null(A)
B = null(A, tolerance)

null(B, A)
null(B, A, tolerance)

The dimension of the range space is the number of singular values of A not greater than tolerance.

The tolerance argument is optional; by default tolerance = max_rc * max_sv * epsilon, where:

The computation is based on singular value decomposition. If the decomposition fails:

Examples

[[cpp11::register]] list null1_(const doubles_matrix<>& a) {
  mat A = as_mat(a);

  mat B;
  bool ok = null(B, A);

  writable::list out(2);
  out[0] = writable::logicals({ok});
  out[1] = as_doubles_matrix(B);

  return out;
}

Find the orthonormal basis of the range space of a matrix {#range}

Find the orthonormal basis of the range space of matrix A, so that $B^t B \approx I_{r,r}$, where $r = \text{rank}(A)$.

Usage:

mat B = orth(A)
B = orth(A, tolerance)

orth(B, A)
orth(B, A, tolerance)

The dimension of the range space is the number of singular values of A greater than tolerance.

The tolerance argument is optional; by default tolerance = max_rc * max_sv * epsilon, where:

The computation is based on singular value decomposition. If the decomposition fails:

Examples

[[cpp11::register]] list orth1_(const doubles_matrix<>& a) {
  mat A = as_mat(a);

  mat B;
  bool ok = orth(B, A);

  writable::list out(2);
  out[0] = writable::logicals({ok});
  out[1] = as_doubles_matrix(B);

  return out;
}

Moore-Penrose pseudo-inverse {#pinv}

Moore-Penrose pseudo-inverse (generalised inverse) of matrix A based on singular value decomposition.

Usage:

B = pinv(A)
B = pinv(A, tolerance)
B = pinv(A, tolerance, method)

pinv(B, A)
pinv(B, A, tolerance)
pinv(B, A, tolerance, method)

The tolerance argument is optional; by default tolerance = max_rc * max_sv * epsilon, where:

Any singular values less than tolerance are treated as zero.

The method argument is optional; method is either "dc" or "std":

If the decomposition fails:

Caveats:

Examples

[[cpp11::register]] list pinv1_(const doubles_matrix<>& a) {
  mat A = as_mat(a);

  mat B = pinv(A);

  writable::list out(1);
  out[0] = as_doubles_matrix(B);

  return out;
}

QR decomposition {#qr}

QR decomposition of matrix X into an orthogonal matrix Q and a right triangular matrix R, with an optional permutation matrix/vector P.

Usage:

// form 1
qr(Q, R, X)

// form 2
qr(Q, R, P, X, "vector")

// form 3
qr(Q, R, P, X, "matrix")

The decomposition has the following form:

If P is specified, a column pivoting decomposition is used.

The diagonal entries of R are ordered from largest to smallest magnitude.

If the decomposition fails, Q, R and P are reset, and the function returns a bool set to false (an error is not thrown).

Examples

[[cpp11::register]] list qr1_(const doubles_matrix<>& x) {
  mat X = as_mat(x);

  mat Q, R;
  umat P;

  bool ok = qr(Q, R, P, X, "matrix");

  writable::list out(4);
  out[0] = writable::logicals({ok});
  out[1] = as_doubles_matrix(Q);
  out[2] = as_doubles_matrix(R);
  out[3] = as_integers_matrix(P);

  return out;
}

Economic QR decomposition {#qr_econ}

Economical QR decomposition of matrix X into an orthogonal matrix Q and a right triangular matrix R, such that $QR = X$.

If the number of rows of X is greater than the number of columns, only the first n rows of R and the first n columns of Q are calculated, where n is the number of columns of X.

If the decomposition fails, Q and R are reset, and the function returns a bool set to false (an error is not thrown).

Examples

[[cpp11::register]] list qr_econ1_(const doubles_matrix<>& x) {
  mat X = as_mat(x);

  mat Q, R;

  bool ok = qr_econ(Q, R, X);

  writable::list out(3);
  out[0] = writable::logicals({ok});
  out[1] = as_doubles_matrix(Q);
  out[2] = as_doubles_matrix(R);

  return out;
}

Generalized Schur decomposition {#qz}

Generalised Schur decomposition for pair of general square matrices A and B of the same size, such that $A = Q^T AA Z^T$ and $B = Q^T BB Z^T$.

The select argument is optional and specifies the ordering of the top left of the Schur form. It is one of the following:

The left and right Schur vectors are stored in Q and Z, respectively.

In the complex-valued problem, the generalised eigenvalues are found in diagvec(AA) / diagvec(BB). If A or B is not square sized, an error is thrown.

If the decomposition fails, AA, BB, Q and Z are reset, and the function returns a bool set to false (an error is not thrown).

Examples

[[cpp11::register]] list qz1_(const doubles_matrix<>& a, const doubles_matrix<>& b,
                                 const char* select) {
  mat A = as_mat(a);
  mat B = as_mat(b);

  mat AA, BB, Q, Z;

  bool ok = qz(AA, BB, Q, Z, A, B, select);

  writable::list out(5);
  out[0] = writable::logicals({ok});
  out[1] = as_doubles_matrix(AA);
  out[2] = as_doubles_matrix(BB);
  out[3] = as_doubles_matrix(Q);
  out[4] = as_doubles_matrix(Z);

  return out;
}

Schur decomposition of a square matrix {#schur}

Schur decomposition of square matrix X into an orthogonal matrix U and an upper triangular matrix S, such that $X = U S U^T$.

If the decomposition fails:

Caveat: in general, Schur decomposition is not unique.

Examples

[[cpp11::register]] list schur1_(const doubles_matrix<>& x) {
  mat X = as_mat(x);

  mat U, S;

  bool ok = schur(U, S, X);

  writable::list out(3);
  out[0] = writable::logicals({ok});
  out[1] = as_doubles_matrix(U);
  out[2] = as_doubles_matrix(S);

  return out;
}

Solve a system of linear equations {#solve}

Solve a dense system of linear equations, $AX = B$, where $X$ is unknown. This is similar functionality to the \ operator in Matlab/Octave (e.g., X = A \ B). The implementation details are available in @Sanderson2017.

$A$ can be square sized (critically determined system), non-square (under/over-determined system), or rank deficient.

$B$ can be a vector or matrix.

The number of rows in A and B must be the same.

By default, matrix A is analysed to automatically determine whether it is a general matrix, band matrix, diagonal matrix, or symmetric/hermitian positive definite (sympd) matrix. Based on the detected matrix structure, a specialised solver is used for faster execution. If no solution is found, an approximate solver is automatically used as a fallback.

If A is known to be a triangular matrix, the solution can be computed faster by explicitly indicating that A is triangular through trimatu() or trimatl() (see examples below).

The settings argument is optional; it is one of the following, or a combination thereof:

The above settings can be combined using the + operator; for example:

solve_opts::fast + solve_opts::no_approx

If a rank deficient system is detected and the solve_opts::no_approx option is not enabled, a warning is emitted and an approximate solution is attempted. Since Armadillo 10.4, this warning can be disabled by setting ARMA_WARN_LEVEL to 1 before including the armadillo header:

#define ARMA_WARN_LEVEL 1
#include <armadillo>

Caveats:

If no solution is found:

Examples

[[cpp11::register]] doubles_matrix<> solve1_(const doubles_matrix<>& a,
                                             const doubles_matrix<>& b) {
  mat A = as_mat(a);
  mat B = as_mat(b);

  mat X = solve(A, B);

  return as_doubles_matrix(X);
}

Singular value decomposition {#svd}

Singular value decomposition of dense matrix X.

If X is square, it can be reconstructed using $X = U S V^T$, where $S$ is a diagonal matrix containing the singular values.

The singular values are in descending order.

The method argument is optional; method is either "dc" or "std":

If the decomposition fails:

Examples

[[cpp11::register]] list svd1_(const doubles_matrix<>& x) {
  mat X = as_mat(x);

  mat U;
  vec s;
  mat V;

  bool ok = svd(U, s, V, X);

  writable::list out(4);
  out[0] = writable::logicals({ok});
  out[1] = as_doubles_matrix(U);
  out[2] = as_doubles(s);
  out[3] = as_doubles_matrix(V);

  return out;
}

Economic singular value decomposition {#svd_econ}

Economical singular value decomposition of dense matrix X.

The singular values are in descending order.

The mode argument is optional; mode is one of:

The method argument is optional; method is either "dc" or "std":

If the decomposition fails, U, s, and V are reset, and a bool set to false is returned (an error is not thrown).

Examples

[[cpp11::register]] list svd_econ1_(const doubles_matrix<>& x) {
  mat X = as_mat(x);

  mat U;
  vec s;
  mat V;

  svd_econ(U, s, V, X);

  writable::list out(3);
  out[0] = as_doubles_matrix(U);
  out[1] = as_doubles(s);
  out[2] = as_doubles_matrix(V);

  return out;
}

Sylvester equation solver {#syl}

Solve the Sylvester equation, i.e., $AX + XB + C = 0$, where $X$ is unknown.

Matrices A, B, and C must be square sized.

If no solution is found:

Examples

[[cpp11::register]] doubles_matrix<> syl1_(const doubles_matrix<>& a,
                                          const doubles_matrix<>& b,
                                          const doubles_matrix<>& c) {
  mat A = as_mat(a);
  mat B = as_mat(b);
  mat C = as_mat(c);

  mat X = syl(A, B, C);

  return as_doubles_matrix(X);
}


Try the cpp11armadillo package in your browser

Any scripts or data that you put into this service are public.

cpp11armadillo documentation built on June 8, 2025, 9:40 p.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close