General real eigenproblems

This page is a walkthrough showing how to use ezarpack::arpack_solver<Asymmetric, Backend> in your C++ code to compute a few eigenpairs \((\lambda,\mathbf{x})\) of

\[\hat A \mathbf{x} = \lambda \hat M \mathbf{x}\]

with a real matrix \(\hat A\) and a real symmetric positive semi-definite matrix \(\hat M\). The asymmetric solver class supports a few computational modes, where the original eigenproblem is recast into

\[\hat O \mathbf{x} = \mu \mathbf{x}.\]

Operator \(\hat O\) acts in a vector space equipped with an inner product defined by the matrix \(\hat B = \hat M\),

\[\langle \mathbf{x}, \mathbf{y} \rangle = \mathbf{x}^\dagger \hat B \mathbf{y}.\]

There are explicit relations between the original eigenvalues \(\lambda\) and their transformed counterparts \(\mu\). \(\lambda\) are either real or come in complex-conjugate pairs. Eigenvectors \(\mathbf{x}\) corresponding to the real eigenvalues \(\lambda\) are also real, and those corresponding to the conjugate pairs come as complex conjugate pairs themselves.

Note

If the transformed matrix \(\hat O\) is symmetric w.r.t. the \(\hat B\)-weighted inner product,

\[\langle \mathbf{x}, \hat O \mathbf{y} \rangle = \langle \hat O \mathbf{x}, \mathbf{y} \rangle,\]

or, equivalently,

\[\hat B \hat O = \hat O^T \hat B,\]

then the faster ezarpack::arpack_solver<Symmetric, Backend> should be used instead.

Typical steps needed to compute the eigenpairs are as follows.

  1. Decide what storage backend you want to use or whether it is appropriate to implement a new one. In the following, we will assume that the Eigen backend has been selected.

  2. Include <ezarpack/arpack_solver.hpp> and the relevant backend header.

#include <ezarpack/arpack_solver.hpp>
#include <ezarpack/storages/eigen.hpp>

Note

<ezarpack/arpack_solver.hpp> includes all three specializations of ezarpack::arpack_solver at once. If you want to speed up compilation a little bit, you can include <ezarpack/solver_base.hpp> and <ezarpack/solver_asymmetric.hpp> instead.

  1. Create a solver object.

const int N = 1000; // Size of matrices A and M

using namespace ezarpack;

// Shorthand for solver's type.
using solver_t = arpack_solver<Asymmetric, eigen_storage>;
// Solver object.
solver_t solver(N);

  1. Fill a params_t structure with calculation parameters.

// params_t is a structure holding parameters of
// the Implicitly Restarted Arnoldi iteration.
using params_t = solver_t::params_t;

// Requested number of eigenvalues to compute.
const int nev = 10;

// Compute the eigenvalues with the largest magnitudes.
auto eigenvalues_select = params_t::LargestMagnitude;

// Compute Ritz vectors (eigenvectors).
params_t::compute_vectors_t compute_vectors = params_t::Ritz;

// params_t's constructor takes three arguments -- mandatory parameters
// that need be set explicitly.
params_t params(nev, eigenvalues_select, compute_vectors);

The following table contains an annotated list of all supported parameters.

Parameter name

Type

Default value

Description

n_eigenvalues

unsigned int

n/a

Number of eigenvalues to compute.

eigenvalues_select

params_t::eigenvalues_select_t (enumeration)

n/a

Part of the spectrum to target. Acceptable values are LargestMagnitude (largest eigenvalues in magnitude), SmallestMagnitude (smallest eigenvalues in magnitude), LargestReal (eigenvalues of largest real part), SmallestReal (eigenvalues of smallest real part), LargestImag (eigenvalues of largest imaginary part) and SmallestImag (eigenvalues of smallest imaginary part).

ncv

int

min(2 * n_eigenvalues + 2, N)

How many Arnoldi vectors to generate at each iteration.

compute_vectors

compute_vectors_t (enumeration)

n/a

Schur – compute only Schur vectors (orthogonal basis vectors of the n_eigenvalues-dimensional subspace), Ritz – compute Ritz vectors (eigenvectors) in addition to the Schur vectors, None – compute neither Schur nor Ritz vectors.

random_residual_vector

bool

true

Use a randomly generated initial residual vector?

sigma

std::complex<double>

0

Complex eigenvalue shift \(\sigma\) for spectral transformation modes.

tolerance

double

Machine precision

Relative tolerance for Ritz value (eigenvalue) convergence.

max_iter

unsigned int

INT_MAX

Maximum number of Arnoldi update iterations allowed.

Note

In the spectral transformation modes, values of eigenvalues_select refer to the spectrum of the transformed problem, not the original one. For instance, LargestMagnitude used in the shift-and-invert mode with a real shift \(\sigma\) will pick eigenvalues \(\lambda\) closest to \(\sigma\), because they correspond to the eigenvalues \(\mu = 1/(\lambda - \sigma)\) that have the largest magnitude.

  1. Optionally set the initial vector for Arnoldi iteration if a better choice than a random vector is known. random_residual_vector parameter must be set to false for the changes made to the initial vector to take effect.

    A view of the residual vector is accessible via method residual_vector() of the solver.

    // Set all components of the initial vector to 1.
auto rv = solver.residual_vector();
for(int i = 0; i < N; ++i) rv[i] = 1.0;

    One may also call residual_vector() later, after a diagonalization run has started, to retrieve the current residual vector.

  2. Choose one of supported computational modes and perform diagonalization. In this part, user is supposed to call the solver object and pass the parameter structure as well as callable objects (e.g. lambda-functions) that represent action of operators \(\hat O\) and \(\hat B\) on a given vector. The supplied objects will be called to generate Arnoldi vectors. Syntax and semantics of the C++ code vary between the computational modes and will be explained individually for each of them.

    • Standard mode (for standard eigenproblems, \(\hat M = \hat I\)).

      using vector_view_t = solver_t::vector_view_t;
using vector_const_view_t = solver_t::vector_const_view_t;

auto Aop = [](vector_const_view_t in, vector_view_t out) {
  // Code implementing action of matrix A on vector 'in':
  // out = A * in
};

solver(Aop, params);

    • Regular inverse mode (for positive-definite \(\hat M\)).

      In this mode, the transformed eigenproblem is defined by \(\hat O = \hat M^{-1} \hat A\), \(\hat B = \hat M\) and \(\lambda = \mu\).

      using vector_view_t = solver_t::vector_view_t;
using vector_const_view_t = solver_t::vector_const_view_t;

auto op = [](vector_const_view_t in, vector_view_t out) {
  // Code implementing action of matrix M^{-1} A on vector 'in':
  // out = M^{-1} * A * in
};
auto Bop = [](vector_const_view_t in, vector_view_t out) {
  // Code implementing action of matrix M on vector 'in':
  // out = M * in
};

solver(op, Bop, solver_t::Inverse, params);

      Inverting a sparse matrix \(\hat M\) will likely make it dense, which is usually undesirable from the storage standpoint. A more practical solution is to compute the sparse LU or Cholesky factorization of \(\hat M\) once (outside of the lambda-function’s body), and write the lambda-function so that it computes out as the solution of the linear system M * out = A * in using the precomputed factorization.

    • Complex Shift-and-Invert mode in real arithmetic I.

      In this mode, the transformed eigenproblem is defined by \(\hat O = \Re[(\hat A - \sigma\hat M)^{-1} \hat M]\), \(\hat B = \hat M\) and \(\mu = \frac{1}{2}\left[\frac{1}{\lambda-\sigma} + \frac{1}{\lambda-\sigma^*}\right]\). The complex spectral shift \(\sigma\) must be set in the parameters structure, see the list of parameters.

      using vector_view_t = solver_t::vector_view_t;
using vector_const_view_t = solver_t::vector_const_view_t;

auto op = [](vector_const_view_t in, vector_view_t out) {
  // Code implementing action of matrix
  // Re[(A - sigma*M)^{-1} * M] on 'in':
  // out = Re[(A - sigma*M)^{-1} * M] * in
};
auto Bop = [](vector_const_view_t in, vector_view_t out) {
  // Code implementing action of matrix M on vector 'in':
  // out = M * in
};

solver(op, Bop, solver_t::ShiftAndInvertReal, params);

    • Complex Shift-and-Invert mode in real arithmetic II.

      In this mode, the transformed eigenproblem is defined by \(\hat O = \Im[(\hat A - \sigma\hat M)^{-1} \hat M]\), \(\hat B = \hat M\) and \(\mu = \frac{1}{2i}\left[\frac{1}{\lambda-\sigma} - \frac{1}{\lambda-\sigma^*}\right]\). The complex spectral shift \(\sigma\) must be set in the parameters structure, see the list of parameters.

      using vector_view_t = solver_t::vector_view_t;
using vector_const_view_t = solver_t::vector_const_view_t;

auto op = [](vector_const_view_t in, vector_view_t out) {
  // Code implementing action of matrix
  // Im[(A - sigma*M)^{-1} * M] on 'in':
  // out = Im[(A - sigma*M)^{-1} * M] * in
};
auto Bop = [](vector_const_view_t in, vector_view_t out) {
  // Code implementing action of matrix M on vector 'in':
  // out = M * in
};

solver(op, Bop, solver_t::ShiftAndInvertImag, params);

    Shift-and-Invert modes in real arithmetic (ShiftAndInvertReal and ShiftAndInvertImag) make extraction of eigenvalues more involved (see below) and should only be used when the amount of storage used by complex arithmetic is prohibitive. Otherwise, the ShiftAndInvert mode of ezarpack::arpack_solver<Complex, Backend> is preferable. Both Shift-and-Invert modes work well close to the complex shift \(\sigma\). However, for large \(\lambda\) operator \(\hat O\) in the mode II dampens the eigenvalues more strongly than that from the mode I. If \(\sigma\) is purely real, \(\hat O = 0\) in the mode II is ill-defined and ShiftAndInvertReal is the only choice.

    Note

    In most computational modes above, it is seemingly necessary to apply operator \(\hat B\) to the same vector twice per generated Arnoldi vector, once in functor op and once in Bop. It is actually possible to spare one of the applications. Calling solver.Bx_available() inside op will tell whether Bop has already been called at the current iteration, and solver.Bx_vector() will return a constant view of the application result \(\hat B \mathbf{x}\).

    The in and out views passed to the callable objects always expose one of three length-\(N\) vectors stored inside the solver object. There is another, indirect way to access them.

    // Get index (0-2) of the current 'in' vector and request a view of it
auto in_view = solver.workspace_vector(solver.in_vector_n());
// Similar for the 'out' vector
auto out_view = solver.workspace_vector(solver.out_vector_n());

    In advanced usage scenarios, the implicit restarting procedure can be customized via an extra argument of solver’s call operator. See ‘Advanced: Customization of Lanczos/Arnoldi implicit restarting’ for more details.

    auto shifts_f = [](solver_t::real_vector_const_view_t ritz_values_re,
                   solver_t::real_vector_const_view_t ritz_values_im,
                   solver_t::real_vector_const_view_t ritz_bounds,
                   solver_t::real_vector_view_t shifts_re,
                   solver_t::real_vector_view_t shifts_im) {
                     // Compute shifts for the implicit restarting
                   };

// Standard mode
solver(op, params, shifts_f);
// Other modes, e.g. Inverse
solver(op, Bop, solver_t::Inverse, params, shifts_f);

    solver_t::operator() can throw two special exception types.

    • ezarpack::maxiter_reached - Maximum number of implicitly restarted Arnoldi iterations has been reached.

    • ezarpack::ncv_insufficient - No shifts could be applied during a cycle of the Implicitly restarted Arnoldi iteration. Consider increasing the number of Arnoldi vectors generated at each iteration (ncv parameter).

    The rest of possible problems reported by ARPACK-NG result in generic std::runtime_error exceptions.

  3. Request computed eigenvalues. In the standard and Inverse computational modes this is done by simply calling

    auto lambda = solver.eigenvalues();

    In the Shift-and-Invert modes, however, the situation is more tricky. Original eigenvalues \(\lambda\) have to be extracted as solutions of a quadratic equation with coefficients given in terms of computed eigenvalues \(\mu\) and of the shift \(\sigma\). ARPACK-NG does not solve that quadratic equation, because figuring out which of the two roots is the actual eigenvalue is not always trivial. Instead, ezARPACK allows to derive \(\lambda\) from computed eigenvectors \(\mathbf{x}\) (compute_vectors parameter must be set to params_t::Ritz) as Rayleigh quotients

    \[\lambda = \frac{\langle \mathbf{x}, \hat A \mathbf{x} \rangle} {\langle \mathbf{x}, \mathbf{x} \rangle} = \frac{\mathbf{x}^\dagger \hat A \mathbf{x}} {\mathbf{x}^\dagger \hat B \mathbf{x}}.\]

    The corresponding C++ code reads

    auto Aop = [](vector_const_view_t in, vector_view_t out) {
  // Code implementing action of matrix A on vector 'in':
  // out = A * in
};

auto lambda = solver.eigenvalues(Aop);

    If the diagonalization run has ended prematurely (for example, when the maximum number of iterations has been reached), then it is still possible to extract solver.nconv() converged eigenpairs.

  4. Optionally request computed eigenvectors (provided the compute_vectors parameter has been set to params_t::Ritz):

    auto vecs = solver.eigenvectors();

    The eigenvectors are solver.nconv() columns of the complex matrix vecs.

  5. Optionally request the Schur vectors, i.e. \(\hat B\)-orthogonal basis vectors of the relevant vector subspace (compute_vectors must be either params_t::Schur or params_t::Ritz).

    auto basis = solver.schur_vectors();

    The basis vectors are solver.nconv() columns of the real matrix view basis.

  6. Optionally request statistics about the completed run.

// Print some computation statistics
auto stats = solver.stats();

std::cout << "Number of Arnoldi update iterations: " << stats.n_iter
          << std::endl;
std::cout << "Total number of O*x operations: " << stats.n_op_x_operations
          << std::endl;
std::cout << "Total number of B*x operations: " << stats.n_b_x_operations
          << std::endl;
std::cout << "Total number of steps of re-orthogonalization: "
          << stats.n_reorth_steps << std::endl;