Decompositions, factorisations, inverses and equation solvers (dense matrices) • cpp11armadillo

This vignette is adapted from the official Armadillo documentation.

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

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

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.

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

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.

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

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.

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

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.

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

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.

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

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.

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

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.

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

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.

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

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.

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

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.

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

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

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

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

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

[[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 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 Sanderson and Curtin (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).

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

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

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

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

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

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

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

Sanderson, Conrad, and Ryan Curtin. 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. https://doi.org/10.1109/ICSPCS.2017.8270510.