---
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);
}
```