--- title: "Decompositions, factorisations, inverses and equation solvers (sparse matrices)" output: rmarkdown::html_vignette bibliography: "references.bib" vignette: > %\VignetteIndexEntry{Decompositions, factorisations, inverses and equation solvers (sparse 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). # Eigen decomposition of a symmetric matrix {#eig_sym} Obtain a limited number of eigenvalues and eigenvectors of sparse symmetric real matrix `X`. Usage: ```cpp 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: ```cpp 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: - The number of obtained eigenvalues/eigenvectors may be lower than requested, depending on the given data. - If the decomposition fails, try first increasing `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). - For an alternative to the `"sm"` form, use the shift-invert mode with `sigma` set to `0.0`. - The implementation in Armadillo 12.6 is considerably faster than earlier versions; further speedups can be obtained by enabling OpenMP in your compiler (e.g., `-fopenmp` in GCC and clang). ## Examples ```cpp [[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 {#eig_gen} Obtain a limited number of eigenvalues and eigenvectors of sparse general (non-symmetric/non-hermitian) square matrix `X`. Usage: ```cpp 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: ```cpp 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: - The number of obtained eigenvalues/eigenvectors may be lower than requested, depending on the given data. - If the decomposition fails, try first increasing `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). - For an alternative to the `"sm"` form, use the shift-invert mode with `sigma` set to `0.0`. ## Examples ```cpp [[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 {#svds} 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: ```cpp 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. - Depending on the given matrix, `svds` may find fewer singular values than specified. - The implementation in Armadillo 12.6 is considerably faster than earlier versions; further speedups can be obtained by enabling OpenMP in your compiler (e.g., `-fopenmp` in GCC and clang). ## Examples ```cpp [[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 {#spsolve} 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: ```cpp // 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"`. - For `"superlu"`, `ARMA_USE_SUPERLU` must be enabled in `armadillo/config.hpp`. - For `"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: ```cpp 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. ## Examples ```cpp [[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 {#spsolve_factoriser} 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 the `opts` argument as per the `spsolve()` 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: - If the factorisation of $A$ does not need to be reused, use `spsolve()` instead. - This class internally uses the SuperLU solver; `ARMA_USE_SUPERLU` must be enabled in `config.hpp`. ## Examples ```cpp [[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; } ```