Iterative Eigensolvers

Iterative Eigensolvers

Julia provides bindings to ARPACK, which can be used to perform iterative solutions for eigensystems (using eigs) or singular value decompositions (using svds).

eigs calculates the eigenvalues and, optionally, eigenvectors of its input(s) using implicitly restarted Lanczos or Arnoldi iterations for real symmetric or general nonsymmetric matrices respectively.

For the single matrix version,

eigs(A; nev=6, ncv=max(20,2*nev+1), which=:LM, tol=0.0, maxiter=300, sigma=nothing, ritzvec=true, v0=zeros((0,))) -> (d,[v,],nconv,niter,nmult,resid)

the following keyword arguments are supported:

whichtype of eigenvalues
:LMeigenvalues of largest magnitude (default)
:SMeigenvalues of smallest magnitude
:LReigenvalues of largest real part
:SReigenvalues of smallest real part
:LIeigenvalues of largest imaginary part (nonsymmetric or complex A only)
:SIeigenvalues of smallest imaginary part (nonsymmetric or complex A only)
:BEcompute half of the eigenvalues from each end of the spectrum, biased in favor of the high end. (real symmetric A only)

We can see the various keywords in action in the following examples:

julia> using IterativeEigensolvers

julia> A = Diagonal(1:4);

julia> λ, ϕ = eigs(A, nev = 2, which=:SM);

julia> λ
2-element Array{Float64,1}:
 1.0000000000000002
 2.0

julia> B = Diagonal([1., 2., -3im, 4im])
4×4 Diagonal{Complex{Float64},Array{Complex{Float64},1}}:
 1.0+0.0im      ⋅          ⋅          ⋅
     ⋅      2.0+0.0im      ⋅          ⋅
     ⋅          ⋅      0.0-3.0im      ⋅
     ⋅          ⋅          ⋅      0.0+4.0im

julia> λ, ϕ = eigs(B, nev=1, which=:LI);

julia> λ
1-element Array{Complex{Float64},1}:
 -4.440892098500626e-16 + 3.999999999999998im

julia> λ, ϕ = eigs(B, nev=1, which=:SI);

julia> λ
1-element Array{Complex{Float64},1}:
 1.3877787807814457e-16 - 2.999999999999999im

julia> λ, ϕ = eigs(B, nev=1, which=:LR);

julia> λ
1-element Array{Complex{Float64},1}:
 2.0 + 4.242754940683747e-17im

julia> λ, ϕ = eigs(B, nev=1, which=:SR);

julia> λ
1-element Array{Complex{Float64},1}:
 4.440892098500626e-16 + 4.0000000000000036im

julia> λ, ϕ = eigs(B, nev=1, sigma=1.5);

julia> λ
1-element Array{Complex{Float64},1}:
 1.9999999999999996 + 2.4290457684137336e-17im
Note

The sigma and which keywords interact: the description of eigenvalues searched for by which do not necessarily refer to the eigenvalues of A, but rather the linear operator constructed by the specification of the iteration mode implied by sigma.

sigmaiteration modewhich refers to eigenvalues of
nothingordinary (forward)$A$
real or complexinverse with level shift sigma$(A - \\sigma I )^{-1}$
Note

Although tol has a default value, the best choice depends strongly on the matrix A. We recommend that users always specify a value for tol which suits their specific needs.

For details of how the errors in the computed eigenvalues are estimated, see:

  • B. N. Parlett, "The Symmetric Eigenvalue Problem", SIAM: Philadelphia, 2/e (1998), Ch. 13.2, "Accessing Accuracy in Lanczos Problems", pp. 290-292 ff.

  • R. B. Lehoucq and D. C. Sorensen, "Deflation Techniques for an Implicitly Restarted Arnoldi Iteration", SIAM Journal on Matrix Analysis and Applications (1996), 17(4), 789–821. doi:10.1137/S0895479895281484

For the two-input generalized eigensolution version,

eigs(A, B; nev=6, ncv=max(20,2*nev+1), which=:LM, tol=0.0, maxiter=300, sigma=nothing, ritzvec=true, v0=zeros((0,))) -> (d,[v,],nconv,niter,nmult,resid)

the following keyword arguments are supported:

whichtype of eigenvalues
:LMeigenvalues of largest magnitude (default)
:SMeigenvalues of smallest magnitude
:LReigenvalues of largest real part
:SReigenvalues of smallest real part
:LIeigenvalues of largest imaginary part (nonsymmetric or complex A only)
:SIeigenvalues of smallest imaginary part (nonsymmetric or complex A only)
:BEcompute half of the eigenvalues from each end of the spectrum, biased in favor of the high end. (real symmetric A only)

eigs returns the nev requested eigenvalues in d, the corresponding Ritz vectors v (only if ritzvec=true), the number of converged eigenvalues nconv, the number of iterations niter and the number of matrix vector multiplications nmult, as well as the final residual vector resid.

We can see the various keywords in action in the following examples:

julia> using IterativeEigensolvers

julia> A = sparse(1.0I, 4, 4); B = Diagonal(1:4);

julia> λ, ϕ = eigs(A, B, nev = 2);

julia> λ
2-element Array{Float64,1}:
 1.0
 0.4999999999999999

julia> A = sparse(1.0I, 4, 4); B = Diagonal([1, -2im, 3, 4im]);

julia> λ, ϕ = eigs(A, B, nev=1, which=:SI);

julia> λ
1-element Array{Complex{Float64},1}:
 0.03291282838780993 - 2.0627621271174514im

julia> λ, ϕ = eigs(A, B, nev=1, which=:LI);

julia> λ
1-element Array{Complex{Float64},1}:
 -0.6428551411711136 + 2.1820633510068994im
Note

The sigma and which keywords interact: the description of eigenvalues searched for by which do not necessarily refer to the eigenvalue problem $Av = Bv\\lambda$, but rather the linear operator constructed by the specification of the iteration mode implied by sigma.

sigmaiteration modewhich refers to the problem
nothingordinary (forward)$Av = Bv\\lambda$
real or complexinverse with level shift sigma$(A - \\sigma B )^{-1}B = v\\nu$
eigs(A; nev=6, ncv=max(20,2*nev+1), which=:LM, tol=0.0, maxiter=300, sigma=nothing, ritzvec=true, v0=zeros((0,))) -> (d,[v,],nconv,niter,nmult,resid)

Computes eigenvalues d of A using implicitly restarted Lanczos or Arnoldi iterations for real symmetric or general nonsymmetric matrices respectively. See the manual for more information.

eigs returns the nev requested eigenvalues in d, the corresponding Ritz vectors v (only if ritzvec=true), the number of converged eigenvalues nconv, the number of iterations niter and the number of matrix vector multiplications nmult, as well as the final residual vector resid.

Examples

julia> using IterativeEigensolvers

julia> A = Diagonal(1:4);

julia> λ, ϕ = eigs(A, nev = 2);

julia> λ
2-element Array{Float64,1}:
 4.0
 3.0
source
eigs(A, B; nev=6, ncv=max(20,2*nev+1), which=:LM, tol=0.0, maxiter=300, sigma=nothing, ritzvec=true, v0=zeros((0,))) -> (d,[v,],nconv,niter,nmult,resid)

Computes generalized eigenvalues d of A and B using implicitly restarted Lanczos or Arnoldi iterations for real symmetric or general nonsymmetric matrices respectively. See the manual for more information.

source
svds(A; nsv=6, ritzvec=true, tol=0.0, maxiter=1000, ncv=2*nsv, v0=zeros((0,))) -> (SVD([left_sv,] s, [right_sv,]), nconv, niter, nmult, resid)

Computes the largest singular values s of A using implicitly restarted Lanczos iterations derived from eigs.

Inputs

  • A: Linear operator whose singular values are desired. A may be represented as a subtype of AbstractArray, e.g., a sparse matrix, or any other type supporting the four methods size(A), eltype(A), A * vector, and A' * vector.

  • nsv: Number of singular values. Default: 6.

  • ritzvec: If true, return the left and right singular vectors left_sv and right_sv. If false, omit the singular vectors. Default: true.

  • tol: tolerance, see eigs.

  • maxiter: Maximum number of iterations, see eigs. Default: 1000.

  • ncv: Maximum size of the Krylov subspace, see eigs (there called nev). Default: 2*nsv.

  • v0: Initial guess for the first Krylov vector. It may have length min(size(A)...), or 0.

Outputs

  • svd: An SVD object containing the left singular vectors, the requested values, and the right singular vectors. If ritzvec = false, the left and right singular vectors will be empty.

  • nconv: Number of converged singular values.

  • niter: Number of iterations.

  • nmult: Number of matrix–vector products used.

  • resid: Final residual vector.

Examples

julia> A = Diagonal(1:4);

julia> s = svds(A, nsv = 2)[1];

julia> s.S
2-element Array{Float64,1}:
 4.0
 2.9999999999999996
Implementation

svds(A) is formally equivalent to calling eigs to perform implicitly restarted Lanczos tridiagonalization on the Hermitian matrix $A^\prime A$ or $AA^\prime$ such that the size is smallest.

source