# 10 - Decompositions, factorisations, inverses and equation solvers (dense matrices)

Adapted from the official Armadillo documentation.

# 1 Cholesky decomposition of symmetric matrix

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:

- 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()`.

## 1.1 Examples

    [[cpp4r::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;
    }

# 2 Eigen decomposition of a symmetric matrix

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:

- `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.

## 2.1 Examples

    [[cpp4r::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;
    }

# 3 Eigen decomposition of a general matrix

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:

- `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.

## 3.1 Examples

    [[cpp4r::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_complexes(eigval);
      out[2] = as_complexes_matrix(eigvec);

      return out;
    }

# 4 Eigen decomposition for a pair of general matrices

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:

- `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.

## 4.1 Examples

    [[cpp4r::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_complexes(eigval);
      out[2] = as_complexes_matrix(eigvec);

      return out;
    }

# 5 Upper Hessenberg decomposition

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:

- `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.

## 5.1 Examples

    [[cpp4r::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;
    }

# 6 Inverse of a general matrix

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:

- `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.

## 6.1 Examples

    [[cpp4r::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;
    }

# 7 Inverse of a symmetric positive definite matrix

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:

- `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.

## 7.1 Examples

    [[cpp4r::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;
    }

# 8 Lowe-upper decomposition with partial pivoting

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:

- `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.

## 8.1 Examples

    [[cpp4r::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;
    }

# 9 Find the orthonormal basis of the null space of a matrix

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:

- \(\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.

## 9.1 Examples

    [[cpp4r::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;
    }

# 10 Find the orthonormal basis of the range space of a matrix

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:

- \(\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.

## 10.1 Examples

    [[cpp4r::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;
    }

# 11 Moore-Penrose pseudo-inverse

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:

- \(\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.

## 11.1 Examples

    [[cpp4r::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;
    }

# 12 QR decomposition

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:

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

## 12.1 Examples

    [[cpp4r::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;
    }

# 13 Economic QR decomposition

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

## 13.1 Examples

    [[cpp4r::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;
    }

# 14 Generalized Schur decomposition

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

## 14.1 Examples

    [[cpp4r::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;
    }

# 15 Schur decomposition of a square matrix

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.

## 15.1 Examples

    [[cpp4r::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;
    }

# 16 Solve a system of linear equations

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 SandersonCurtin 2017.

\(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:

    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:

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

## 16.1 Examples

    [[cpp4r::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);
    }

# 17 Singular value decomposition

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

## 17.1 Examples

    [[cpp4r::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;
    }

# 18 Economic singular value decomposition

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

## 18.1 Examples

    [[cpp4r::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;
    }

# 19 Sylvester equation solver

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

## 19.1 Examples

    [[cpp4r::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);
    }

# 20 References

Sanderson C, Curtin R (2017). “An open source C++ implementation of multi-threaded Gaussian mixture models, k-means and expectation maximisation.” In 2017 11th International Conference on Signal Processing and Communication Systems (ICSPCS), 1–8. doi:10.1109/ICSPCS.2017.8270510. https://ieeexplore.ieee.org/abstract/document/8270510?casa_token=OezPt_eRZ1gAAAAA:zVbHP2fNhLNvrTZ9pjSfTrPSH-DtNUQgfnONPinr3NAGBbdmLkUz7M1NCC1fM4o1_-Ad355j5g.