Function arpack_rssolve 

pub fn arpack_rssolve<F: FnMut(&[f64], &mut [f64])>(
    n: usize,
    matvec: F,
    options: &ArpackOptions,
    storage: Option<&mut ArpackStorage>,
) -> Result<ArpackSymmetricResult>

Eigenvalues and eigenvectors of a symmetric linear operator given as a Rust closure, with ARPACK’s implicitly restarted Lanczos method (igraph_arpack_rssolve).

matvec(x, y) must store the product A x into y (which is zeroed before each call, so accumulating into it is fine); both slices have length n. The matrix is never formed: this is the method of choice for large sparse or structured matrices. In ArpackMode::ShiftInvert the closure must compute (A - sigma I)^-1 x instead. A panic in the closure stops ARPACK and is resumed by this function; a product containing NaN or infinite values stops it with an ErrorKind::InvalidValue error.

Time complexity: depends on the matrix-vector product; usually a few iterations suffice, so with an O(n) product the eigenvalues are found in about O(n) time.

In regular mode, 1 × 1 and 2 × 2 problems are solved in closed form by igraph without running ARPACK. For 2 × 2 problems igraph 1.0.0 and 1.0.1 confuse the magnitude-based choices with the algebraic ones and return unnormalized eigenvectors; this wrapper computes both eigenpairs and selects and normalizes the requested ones itself.

Singular operators: the start vector is multiplied first

ARPACK does not use the start vector v0 itself as the first Krylov vector: it first applies the operator and starts from A v0 (to force the start vector into the range of the operator, see dgetv0). For a singular symmetric operator the range is orthogonal to the null space, so the component of v0 along the null vectors is wiped out: when A v0 has no null component exactly (e.g. integer matrices and start vectors, as for graph Laplacians), the eigenvalue 0 can be invisible to the iteration, and ARPACK then silently returns other eigenvalues instead. A start vector lying in the null space fails outright with ArpackError::ZeroStart. The random default start vector usually recovers the null space through rounding errors, but there is no guarantee.

The robust fix is a shift: solve for A + s I with some s > 0 that makes the operator non-singular (for a positive semidefinite A such as a Laplacian any s > 0 works, and “smallest magnitude” becomes “smallest algebraic”), then subtract s from the eigenvalues; the eigenvectors are unchanged. See the second example below.

Nesting

igraph’s ARPACK keeps the state of the running iteration in thread-local statics, so it is not re-entrant: matvec (and any interruption or progress handler, which igraph calls from inside the ARPACK loop) must not start another ARPACK computation on the same thread. The ARPACK-based functions of this module (this function, arpack_rnsolve, the SparseMat ARPACK solvers, the ARPACK paths of eigen_matrix_symmetric and eigen_symmetric_fn, Graph::eigen_adjacency and the spectral embeddings) check this: started while another one of them runs on the thread, they fail with an ErrorKind::Failure error instead of corrupting the running solver. So do the ARPACK-based functions of other modules: Graph::eigenvector_centrality, Graph::hub_and_authority_scores, Graph::centralization_eigenvector_centrality, PageRank with PageRankAlgo::Arpack and Graph::community_leading_eigenvector. Calling other igraph functions from matvec is fine, even failing ones: the closure runs in its own level of igraph’s cleanup stack.

Binds igraph_arpack_rssolve.

See also SparseMat::arpack_rssolve for a stored sparse matrix, eigen_symmetric_fn for a front-end that can also use LAPACK, and the graph-level spectral routines built on ARPACK, e.g. Graph::eigenvector_centrality and Graph::eigen_adjacency.

Errors

Invalid options (see ArpackOptions) give ErrorKind::InvalidValue; ARPACK failures give ErrorKind::Arpack, see ArpackError::from_error for the ARPACK error condition.

Examples

The cycle graph C_n has adjacency eigenvalues 2 cos(2 pi k / n): the largest is 2, with the constant eigenvector.

use igraph::linalg::{arpack_rssolve, ArpackOptions, ArpackWhich};
let n = 12;
let cycle = |x: &[f64], y: &mut [f64]| {
    for i in 0..n {
        y[i] = x[(i + 1) % n] + x[(i + n - 1) % n];
    }
};
let opts = ArpackOptions::default().with_nev(1).with_which(ArpackWhich::LargestAlgebraic);
let res = arpack_rssolve(n, cycle, &opts, None).unwrap();
assert!((res.values[0] - 2.0).abs() < 1e-10);
let v = res.vectors.column(0);
assert!(v.iter().all(|x| (x.abs() - 1.0 / (n as f64).sqrt()).abs() < 1e-8));

The Laplacian of the path on 10 vertices is singular (eigenvalue 0, with the constant eigenvector). The ramp v0 = (0, 1, ..., 9) has a constant component, but ARPACK iterates from L v0 = e_9 - e_0, which has exactly none (and, by the mirror symmetry of the path, neither have the following Krylov vectors): the eigenvalue 0 is missed altogether. Shifting to L + I keeps the constant component of v0 and finds it:

use igraph::linalg::{arpack_rssolve, ArpackOptions, ArpackWhich};
use std::f64::consts::PI;
let n = 10;
// y = (L + shift I) x for the path Laplacian L.
let laplacian = move |shift: f64| {
    move |x: &[f64], y: &mut [f64]| {
        for i in 0..n {
            let degree = if i == 0 || i == n - 1 { 1.0 } else { 2.0 };
            y[i] = (degree + shift) * x[i];
            if i > 0 {
                y[i] -= x[i - 1];
            }
            if i + 1 < n {
                y[i] -= x[i + 1];
            }
        }
    }
};
let ramp: Vec<f64> = (0..n).map(|i| i as f64).collect();
let lambda1 = 2.0 - 2.0 * (PI / n as f64).cos(); // the smallest non-zero eigenvalue

let opts = ArpackOptions::default()
    .with_nev(2)
    .with_which(ArpackWhich::SmallestMagnitude)
    .with_start(ramp);
let plain = arpack_rssolve(n, laplacian(0.0), &opts, None).unwrap();
assert!(plain.values.iter().all(|v| v.abs() > 0.09)); // 0 is missed
assert!(plain.values.iter().any(|v| (v - lambda1).abs() < 1e-8));

let shifted_opts = opts.with_which(ArpackWhich::SmallestAlgebraic);
let shifted = arpack_rssolve(n, laplacian(1.0), &shifted_opts, None).unwrap();
let values: Vec<f64> = shifted.values.iter().map(|v| v - 1.0).collect();
assert!(values[0].abs() < 1e-10);
assert!((values[1] - lambda1).abs() < 1e-8);