Decompositions, factorisations, inverses and equation solvers (sparse matrices) (original) (raw)

This vignette is adapted from the official Armadillo documentation.

Eigen decomposition of a symmetric matrix

Obtain a limited number of eigenvalues and eigenvectors of sparse symmetric real matrix X.

Usage:

vec eigval = eigs_sym(X, k)
vec eigval = eigs_sym(X, k, form)
vec eigval = eigs_sym(X, k, form, opts)
vec eigval = eigs_sym(X, k, sigma)
vec eigval = eigs_sym(X, k, sigma, opts)

eigs_sym(eigval, X, k)
eigs_sym(eigval, X, k, form)
eigs_sym(eigval, X, k, form, opts)
eigs_sym(eigval, X, k, sigma)
eigs_sym(eigval, X, k, sigma, opts)

eigs_sym(eigval, eigvec, X, k)
eigs_sym(eigval, eigvec, X, k, form)
eigs_sym(eigval, eigvec, X, k, form, opts)
eigs_sym(eigval, eigvec, X, k, sigma)
eigs_sym(eigval, eigvec, X, k, sigma, opts)

k specifies the number of eigenvalues and eigenvectors.

The argument form is optional and is one of the following:

The argument sigma is optional; if sigma is given, eigenvalues closest to sigma are found via shift-invert mode. Note that to use sigma, both ARMA_USE_ARPACK and ARMA_USE_SUPERLU must be enabled in armadillo/config.hpp.

The opts argument is optional; opts is an instance of the eigs_opts structure:

struct eigs_opts
{
  double       tol;     // default: 0
  unsigned int maxiter; // default: 1000
  unsigned int subdim;  // default: max(2*k+1, 20)
};

The eigenvalues and corresponding eigenvectors are stored in eigval and eigvec, respectively.

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

If the decomposition fails:

Caveats:

Examples

[[cpp11::register]] list eig_sym2_(const doubles_matrix<>& x,
                                   const char* method,
                                   const int& k) {
  sp_mat X = as_SpMat(x);

  sp_mat Y = X.t() * X;

  vec eigval;
  mat eigvec;

  eigs_opts opts;
  opts.maxiter = 10000;
  bool ok = eigs_sym(eigval, eigvec, Y, k, method, opts);

  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

Obtain a limited number of eigenvalues and eigenvectors of sparse general (non-symmetric/non-hermitian) square matrix X.

Usage:

cx_vec eigval = eigs_gen(X, k)
cx_vec eigval = eigs_gen(X, k, form)
cx_vec eigval = eigs_gen(X, k, sigma)
cx_vec eigval = eigs_gen(X, k, form, opts)
cx_vec eigval = eigs_gen(X, k, sigma, opts)

eigs_gen(eigval, X, k)
eigs_gen(eigval, X, k, form)
eigs_gen(eigval, X, k, sigma)
eigs_gen(eigval, X, k, form, opts)
eigs_gen(eigval, X, k, sigma, opts)

eigs_gen(eigval, eigvec, X, k)
eigs_gen(eigval, eigvec, X, k, form)
eigs_gen(eigval, eigvec, X, k, sigma)
eigs_gen(eigval, eigvec, X, k, form, opts)
eigs_gen(eigval, eigvec, X, k, sigma, opts)

k specifies the number of eigenvalues and eigenvectors.

The argument form is optional; form is one of the following:

The argument sigma is optional; if sigma is given, eigenvalues closest to sigma are found via shift-invert mode. Note that to use sigma, both ARMA_USE_ARPACK and ARMA_USE_SUPERLU must be enabled in armadillo/config.hpp.

The opts argument is optional; opts is an instance of the eigs_opts structure:

struct eigs_opts
{
  double       tol;     // default: 0
  unsigned int maxiter; // default: 1000
  unsigned int subdim;  // default: max(2*k+1, 20)
};

The eigenvalues and corresponding eigenvectors are stored in eigval and eigvec, respectively.

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

If the decomposition fails:

Caveats:

Examples

[[cpp11::register]] list eig_gen2_(const doubles_matrix<>& x,
                                   const char* method,
                                   const int& k) {
  sp_mat X = as_SpMat(x);

  cx_vec eigval;
  cx_mat eigvec;

  eigs_opts opts;
  opts.maxiter = 10000;

  bool ok = eigs_gen(eigval, eigvec, X, k, method, opts);

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

  return out;
}

Truncated singular value decomposition

Obtain a limited number of singular values and singular vectors (truncated SVD) of a sparse matrix.

The singular values and vectors are calculated via sparse eigen decomposition of:

\[ \begin{bmatrix} 0_{n \times n} & X \\ X^T & 0_{m \times m} \end{bmatrix} \]

where \(n\) and \(m\) are the number of rows and columns of X, respectively.

Usage:

vec s = svds(X, k)
vec s = svds(X, k, tol)

svds(vec s, X, k)
svds(vec s, X, k, tol)

svds(mat U, vec s, mat V, sp_mat X, k)
svds(mat U, vec s, mat V, sp_mat X, k, tol)

svds(cx_mat U, vec s, cx_mat V, sp_cx_mat X, k)
svds(cx_mat U, vec s, cx_mat V, sp_cx_mat X, k, tol)

k specifies the number of singular values and singular vectors.

The singular values are in descending order.

The argument tol is optional; it specifies the tolerance for convergence, and it is passed as tol / sqrt(2) to eigs_sym.

If the decomposition fails:

Caveats:

Examples

[[cpp11::register]] list svds1_(const doubles_matrix<>& x, const int& k) {
  sp_mat X = as_SpMat(x);

  // convert all values below 0.1 to zero
  X.transform([](double val) { return (std::abs(val) < 0.1) ? 0 : val; });

  mat U;
  vec s;
  mat V;

  bool ok = svds(U, s, V, X, k);

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

  return out;
}

Solve a system of linear equations

Solve a sparse system of linear equations, \(A \cdot X = B\), where \(A\) is a sparse matrix, \(B\) is a dense matrix or vector, and \(X\) is unknown.

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

Usage:

// A = matrix, b = vector
vec x = spsolve(A, b)
vec x = spsolve(A, b, solver)
vec x = spsolve(A, b, solver, opts)

// A = matrix, B = matrix
mat X = spsolve(A, B)

spsolve(x, A, b)
spsolve(x, A, b, solver)
spsolve(x, A, b, solver, opts)

The solver argument is optional; solver is either "superlu" (default) or "lapack".

The opts argument is optional and applicable to the SuperLU solver, and is an instance of the superlu_opts structure:

struct superlu_opts
{
  bool             allow_ugly;   // default: false
  bool             equilibrate;  // default: false
  bool             symmetric;    // default: false
  double           pivot_thresh; // default: 1.0
  permutation_type permutation;  // default: superlu_opts::COLAMD
  refine_type      refine;       // default: superlu_opts::REF_NONE
};

If no solution is found:

The SuperLU solver is mainly useful for very large and/or highly sparse matrices.

To reuse the SuperLU factorisation of \(A\) for finding solutions where \(B\) is iteratively changed, see the spsolve_factoriser class.

If there is sufficient memory to store a dense version of matrix \(A\), the LAPACK solver can be faster.

Examples

[[cpp11::register]] doubles spsolve1_(const doubles_matrix<>& a,
                                      const doubles& b,
                                      const char* method) {
  sp_mat A = as_SpMat(a);
  vec B = as_Col(b);

  vec X = spsolve(A, B, method);

  return as_doubles(X);
}

Factorise a sparse matrix for solving linear systems

Class for factorisation of sparse matrix \(A\) for solving systems of linear equations in the form \(AX = B\).

Allows the SuperLU factorisation of \(A\) to be reused for finding solutions in cases where \(B\) is iteratively changed.

For an instance of spsolve_factoriser named as SF, the member functions are:

Caveats:

Examples

[[cpp11::register]] list spsolve_factoriser1_(const doubles_matrix<>& a,
                                              const list& b) {
  sp_mat A = as_SpMat(a);

  bool status = SF.factorise(A);

  if (status == false) {
    stop("factorisation failed");
  }
  
  double rcond_value = SF.rcond();

  vec B1 = as_Col(b[0]);
  vec B2 = as_Col(b[1]);

  vec X1, X2;

  bool ok1 = SF.solve(X1, B1);
  bool ok2 = SF.solve(X2, B2);

  if (ok1 == false) {
    stop("couldn't find X1");
  }

  if (ok2 == false) {
    stop("couldn't find X2");
  }

  writable::list out(3);
  out[0] = writable::logicals({status && ok1 && ok2});
  out[1] = as_doubles(X1);
  out[2] = as_doubles(X2);

  return out;
}