This vignette is adapted from the official Armadillo documentation.
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:
"lm": obtain eigenvalues with largest magnitude (default operation)."sm": obtain eigenvalues with smallest magnitude (see the caveats below)."la": obtain eigenvalues with largest algebraic value.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)
};tol specifies the tolerance for convergence.maxiter specifies the maximum number of Arnoldi iterations.subdim specifies the dimension of the Krylov subspace, with the constraint k < subdim <= X.n_rows; the recommended value is subdim >= 2*k.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:
eigval = eigs_sym(X,k) resets eigval and throws an error.eigs_sym(eigval,X,k) resets eigval and returns a bool set to false (an error is not thrown).eigs_sym(eigval,eigvec,X,k) resets eigval and eigvec and returns a bool set to false (an error is not thrown).Caveats:
opts.subdim (Krylov subspace dimension), and, as secondary options, try increasing opts.maxiter (maximum number of iterations), and/or opts.tol (tolerance for convergence), and/or k (number of eigenvalues)."sm" form, use the shift-invert mode with sigma set to 0.0.-fopenmp in GCC and clang).[[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;
}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:
"lm": obtain eigenvalues with largest magnitude (default operation)."sm": obtain eigenvalues with smallest magnitude (see the caveats below)."lr": obtain eigenvalues with largest real part."sr": obtain eigenvalues with smallest real part."li": obtain eigenvalues with largest imaginary part."si": obtain eigenvalues with smallest imaginary part.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)
};tol specifies the tolerance for convergence.maxiter specifies the maximum number of Arnoldi iterations.subdim specifies the dimension of the Krylov subspace, with the constraint k < subdim <= X.n_rows; the recommended value is subdim >= 2*k.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:
eigval = eigs_gen(X, k) resets eigval and throws an error.eigs_gen(eigval, X, k) resets eigval and returns a bool set to false (an error is not thrown).eigs_gen(eigval,eigvec, X, k) resets eigval and eigvec and returns a bool set to false (an error is not thrown).Caveats:
opts.subdim (Krylov subspace dimension), and, as secondary options, try increasing opts.maxiter (maximum number of iterations), and/or opts.tol (tolerance for convergence), and/or k (number of eigenvalues)."sm" form, use the shift-invert mode with sigma set to 0.0.[[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;
}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:
s = svds(X, k) resets s and throws an error.svds(s, X, k) resets s and returns a bool set to false (an error is not thrown).svds(U, s, V, X, k) resets U, s, and V and returns a bool set to false (an error is not thrown).Caveats:
svds is intended only for finding a few singular values from a large sparse matrix; to find all singular values, use svd instead.svds may find fewer singular values than specified.-fopenmp in GCC and clang).[[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 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".
"superlu", ARMA_USE_SUPERLU must be enabled in armadillo/config.hpp."lapack", the sparse matrix \(A\) is converted to a dense matrix before using the LAPACK solver. This considerably increases memory usage.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
};allow_ugly is either true or false; it indicates whether to keep solutions of systems singular to working precision.equilibrate is either true or false; it indicates whether to equilibrate the system (scale the rows and columns of \(A\) to have unit norm).symmetric is either true or false; it indicates whether to use SuperLU symmetric mode, which gives preference to diagonal pivots.pivot_thresh is in the range [0.0, 1.0], used for determining whether a diagonal entry is an acceptable pivot (details in SuperLU documentation).permutation specifies the type of column permutation; it is one of:
superlu_opts::NATURAL: natural ordering.superlu_opts::MMD_ATA: minimum degree ordering on structure of \(A^T \cdot A\).superlu_opts::MMD_AT_PLUS_A: minimum degree ordering on structure of \(A^T + A\).superlu_opts::COLAMD: approximate minimum degree column ordering.refine specifies the type of iterative refinement; it is one of:
superlu_opts::REF_NONE: no refinement.superlu_opts::REF_SINGLE: iterative refinement in single precision.superlu_opts::REF_DOUBLE: iterative refinement in double precision.superlu_opts::REF_EXTRA: iterative refinement in extra precision.If no solution is found:
x = spsolve(A, b) resets x and throws an error.spsolve(x, A, b) resets x and returns a bool set to false (an error is not thrown).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.
[[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);
}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:
SF.factorise(A. opts): factorise square-sized sparse matrix $A. Optional settings are given in theoptsargument as per thespsolve()` function. If the factorisation fails, a bool set to false is returned.SF.solve(X, B): using the given dense matrix \(B\) and the computed factorisation, store in \(X\) the solution to $AX = B`. If computing the solution fails, \(X\) is reset and a bool set to false is returned.SF.rcond(): return the 1-norm estimate of the reciprocal condition number computed during the factorisation. Values close to 1 suggest that the factorised matrix is well-conditioned. Values close to 0 suggest that the factorised matrix is badly conditioned.SF.reset(): reset the instance and release all memory related to the stored factorisation; this is automatically done when the instance goes out of scope.Caveats:
spsolve() instead.ARMA_USE_SUPERLU must be enabled in config.hpp.[[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;
}