---
title: "Decompositions, factorisations, inverses and equation solvers (dense matrices)"
output: rmarkdown::html_vignette
bibliography: "references.bib"
vignette: >
%\VignetteIndexEntry{Decompositions, factorisations, inverses and equation solvers (dense matrices)}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
```
This vignette is adapted from the official Armadillo
[documentation](https://arma.sourceforge.net/docs.html).
# 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:
```cpp
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:
- Forms 1 and 2 with `layout = "upper"`: $X = R^T R$.
- Forms 1 and 2 with `layout = "lower"`: $X = R R^T$.
- Form 3 with `layout = "upper"`: $X(P,P) = R^T R$, where $X(P,P)$ is a
non-contiguous view of $X$.
- Form 3 with `layout = "lower"`: $X(P,P) = R * R^T$, where $X(P,P)$ is a
non-contiguous view of $X$.
- Form 4 with `layout = "upper"`: $X = P R^T R P^T$.
- Form 4 with `layout = "lower"`: $X = P R R^T P^T$.
Caveats:
- `R = chol(X)` and `R = chol(X,layout)` reset `R` and throw an error if the
decomposition fails. The other forms reset `R` and `P`, and return a bool set
to false without an error.
- To find the inverse of a positive definite matrix, use `inv_sympd()`.
## Examples
```cpp
[[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:
```cpp
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:
- `eigval = eig_sym(X)` resets `eigval` and throws an error.
- `eig_sym(eigval, X)` resets `eigval` and returns a bool set to false without
an error.
- `eig_sym(eigval, eigvec, X)` resets `eigval` and `eigvec`, and returns a bool
set to false without an error.
## Examples
```cpp
[[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:
```cpp
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:
- `eigval = eig_gen(X)` resets `eigval` and throws an error.
- `eig_gen(eigval, X)` resets `eigval` and returns a bool set to false without
an error.
- `eig_gen(eigval, eigvec, X)` resets `eigval` and `eigvec`, and returns a bool
set to false without an error.
- `eig_gen(eigval, leigvec, reigvec, X)` resets `eigval`, `leigvec`, and
`reigvec`, and returns a bool set to false without an error.
## Examples
```cpp
[[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:
```cpp
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:
- `eigval = eig_pair(A, B)` resets `eigval` and throws an error.
- `eig_pair(eigval, A, B)` resets `eigval` and returns a bool set to false
without an error.
- `eig_pair(eigval, eigvec, A, B)` resets `eigval` and `eigvec`, and returns a
bool set to false without an error.
- `eig_pair(eigval, leigvec, reigvec, A, B)` resets `eigval`, `leigvec`, and
`reigvec`, and returns a bool set to false without an error.
## Examples
```cpp
[[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:
```cpp
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:
- `H = hess(X)` resets `H` and throws an error.
- `hess(H, X)` resets `H` and returns a bool set to false without an error.
- `hess(U, H, X)` resets `U` and `H`, and returns a bool set to false without
an error.
Caveat: in general, upper Hessenberg decomposition is not unique.
## Examples
```cpp
[[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:
```cpp
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:
- `inv_opts::no_ugly`: do not provide inverses for poorly conditioned matrices
(where `rcond < A.n_rows * datum::eps`).
- `inv_opts::allow_approx`: allow approximate inverses for rank deficient or
poorly conditioned matrices.
The reciprocal condition number is optionally calculated and stored in `rcond`:
- `rcond` close to 1 suggests that `A` is well-conditioned.
- `rcond` close to 0 suggests that `A` is badly conditioned.
If `A` is not square sized, an error is thrown.
If `A` appears to be singular:
- `B = inv(A)` resets `B` and throws an error.
- `inv(B, A)` resets `B` and returns a bool set to false without an error.
- `inv(B, rcond, A)` resets `B`, sets `rcond` to zero, and returns a bool set to
false without an error.
Caveats:
- If matrix `A` is known to be symmetric positive definite, `inv_sympd()` is
faster.
- If matrix `A` is known to be diagonal, use `inv(diagmat(A))`.
- If matrix `A` is known to be triangular, use `inv(trimatu(A))` or
`inv(trimatl(A))`.
To solve a system of linear equations, such as `Z = inv(X) * Y`, `solve()` can
be faster and/or more accurate.
## Examples
```cpp
[[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:
```cpp
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:
- `inv_opts::no_ugly`: do not provide inverses for poorly conditioned matrices
(where `rcond < A.n_rows * datum::eps`).
- `inv_opts::allow_approx`: allow approximate inverses for rank deficient or
poorly conditioned matrices.
The reciprocal condition number is optionally calculated and stored in `rcond`:
- `rcond` close to 1 suggests that `A` is well-conditioned.
- `rcond` close to 0 suggests that `A` is badly conditioned.
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:
- `B = inv_sympd(A)` resets `B` and throws an error.
- `inv_sympd(B, A)` resets `B` and returns a bool set to false without an error.
- `inv_sympd(B, rcond, A)` resets `B`, sets `rcond` to zero, and returns a bool
set to false without an error.
Caveats:
- If matrix `A` is known to be symmetric, use `inv_sympd(symmatu(A))` or
`inv_sympd(symmatl(A))`.
- If matrix `A` is known to be diagonal, use `inv(diagmat(A))`.
- If matrix `A` is known to be triangular, use `inv(trimatu(A))` or
`inv(trimatl(A))`.
- To solve a system of linear equations, such as `Z = inv(X) * Y`, `solve()` can
be faster and/or more accurate.
## Examples
```cpp
[[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:
```cpp
// 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:
- `lu(L, U, P, X)` resets `L`, `U`, and `P`, and returns a bool set to false
without an error.
- `lu(L, U, X)` resets `L` and `U`, and returns a bool set to false without an
error.
## Examples
```cpp
[[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:
```cpp
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:
- $\max_{\text{rc}} = \max(n_\text{rows}, n_\text{cols})$
- $\max_{\text{sv}} = \max(\text{maximum singular value of } A)$ (i.e. spectral norm)
- $\epsilon = 1 - \text{least value greater than 1 that is representable}$
The computation is based on singular value decomposition. If the decomposition
fails:
- `B = null(A)` resets `B` and throws an error.
- `null(B, A)` resets `B` and returns a bool set to false without an error.
## Examples
```cpp
[[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:
```cpp
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:
- $\max_{\text{rc}} = \max(n_\text{rows}, n_\text{cols})$
- $\max_{\text{sv}} = \max(\text{maximum singular value of } A)$ (i.e. spectral norm)
- $\epsilon = 1 - \text{least value greater than 1 that is representable}$
The computation is based on singular value decomposition. If the decomposition
fails:
- `B = orth(A)` resets `B` and throws an error.
- `orth(B, A)` resets `B` and returns a bool set to false without an error.
## Examples
```cpp
[[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:
```cpp
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:
- $\max_{\text{rc}} = \max(n_\text{rows}, n_\text{cols})$
- $\max_{\text{sv}} = \max(\text{maximum singular value of } A)$ (i.e. spectral norm)
- $\epsilon = 1 - \text{least value greater than 1 that is representable}$
Any singular values less than `tolerance` are treated as zero.
The `method` argument is optional; `method` is either `"dc"` or `"std"`:
- `"dc"` indicates divide-and-conquer method (default setting).
- `"std"` indicates standard method.
- The divide-and-conquer method provides slightly different results than the
standard method, but is considerably faster for large matrices.
If the decomposition fails:
- `B = pinv(A)` resets `B` and throws an error.
- `pinv(B, A)` resets `B` and returns a bool set to false without an error.
Caveats:
- To find approximate solutions to under-/over-determined or rank deficient
systems of linear equations, `solve()` can be considerably faster and/or more
accurate.
- If the given matrix `A` is square-sized and only occasionally rank deficient,
using `inv()` or `inv_sympd()` with the `inv_opts::allow_approx` option is
faster.
## Examples
```cpp
[[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:
```cpp
// 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:
- Form 1: $Q * R = X$.
- Form 2: $Q * R = Y$, where `Y` is a non-contiguous view of `X` with columns
permuted by `P`, and `P` is a permutation vector with type `uvec`.
- Form 3: $Q * R = X * P$, where `P` is a permutation matrix with type `umat`.
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
```cpp
[[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
```cpp
[[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:
- `"none"`: no ordering (default operation).
- `"lhp"`: left-half-plane: eigenvalues with real part < 0.
- `"rhp"`: right-half-plane: eigenvalues with real part > 0.
- `"iuc"`: inside-unit-circle: eigenvalues with absolute value < 1.
- `"ouc"`: outside-unit-circle: eigenvalues with absolute value > 1.
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
```cpp
[[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:
- `S = schur(X)` resets `S` and throws an error.
- `schur(S, X)` resets `S` and returns a bool set to false (an error is not
thrown).
- `schur(U, S, X)` resets `U` and `S`, and returns a bool set to false (an error
is not thrown).
Caveat: in general, Schur decomposition is not unique.
## Examples
```cpp
[[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:
- `solve_opts::fast`: fast mode: disable determining solution quality via
`rcond`, disable iterative refinement, disable equilibration.
- `solve_opts::refine`: apply iterative refinement to improve solution quality
(matrix `A` must be square).
- `solve_opts::equilibrate`: equilibrate the system before solving (matrix `A`
must be square).
- `solve_opts::likely_sympd`: indicate that matrix `A` is likely
symmetric/hermitian positive definite (sympd).
- `solve_opts::allow_ugly`: keep solutions of systems that are singular to
working precision.
- `solve_opts::no_approx`: do not find approximate solutions for rank deficient
systems.
- `solve_opts::force_sym`: force use of the symmetric/hermitian solver (not
limited to sympd matrices).
- `solve_opts::force_approx`: force use of the approximate solver.
The above settings can be combined using the `+` operator; for example:
```cpp
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:
```cpp
#define ARMA_WARN_LEVEL 1
#include <armadillo>
```
Caveats:
- Using `solve_opts::fast` will speed up finding the solution, but for poorly
conditioned systems the solution may have lower quality.
- Not all sympd matrices are automatically detected; to directly indicate that
matrix `A` is likely sympd, use `solve_opts::likely_sympd`.
- Using `solve_opts::force_approx` is only advised if the system is known to be
rank deficient; the approximate solver is considerably slower.
If no solution is found:
- `X = solve(A,B)` resets `X` and throws an error.
- `solve(X,A,B)` resets `X` and returns a bool set to false (an error is not
thrown).
## Examples
```cpp
[[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"`:
- `"dc"` indicates divide-and-conquer method (default setting).
- `"std"` indicates standard method.
- The divide-and-conquer method provides slightly different results than the
standard method, but is considerably faster for large matrices.
If the decomposition fails:
- `s = svd(X)` resets `s` and throws an error.
- `svd(s, X)` resets `s` and returns a bool set to false (an error is not
thrown).
- `svd(U, s, V, X)` resets `U`, `s`, and `V`, and returns a bool set to false
(an error is not thrown).
## Examples
```cpp
[[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:
- `"both"`: compute both left and right singular vectors (default operation).
- `"left"`: compute only left singular vectors.
- `"right"`: compute only right singular vectors.
The `method` argument is optional; `method` is either `"dc"` or `"std"`:
- `"dc"` indicates divide-and-conquer method (default setting).
- `"std"` indicates standard method.
- The divide-and-conquer method provides slightly different results than the
standard method, but is considerably faster for large matrices.
If the decomposition fails, `U`, `s`, and `V` are reset, and a bool set to false
is returned (an error is not thrown).
## Examples
```cpp
[[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:
- `syl(A,B,C)` resets `X` and throws an error.
- `syl(X,A,B,C)` resets `X` and returns a bool set to false (an error is not
thrown).
## Examples
```cpp
[[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);
}
```