Getting started

Installing

In Julia open the package manager in the REPL via ] and run:

(v1.0) pkg> add ArnoldiMethod

Then use the package.

using ArnoldiMethod

Construct a partial Schur decomposition

ArnoldiMethod.jl exports the partialschur function which can be used to obtain a partial Schur decomposition of any matrix A.

ArnoldiMethod.partialschur — Function.

partialschur(A; nev, which, tol, mindim, maxdim, restarts) → PartialSchur, History

Find nev approximate eigenpairs of A with eigenvalues near a specified target.

The matrix A can be any linear map that implements mul!(y, A, x), eltype and size.

The method will run iteratively until the eigenpairs are approximated to the prescribed tolerance or until restarts restarts have passed.

Arguments

The most important keyword arguments:

Keyword	Type	Default	Description
`nev`	`Int`	`min(6, size(A, 1))`	Number of eigenvalues
`which`	`Target`	`LM()`	One of `LM()`, `LR()`, `SR()`, `LI()`, `SI()`, see below.
`tol`	`Real`	`√eps`	Tolerance for convergence: ‖Ax - xλ‖₂ < tol * ‖λ‖

The target which can be any of subtypes(ArnoldiMethod.Target):

Target	Description
`LM()`	Largest magnitude: `abs(λ)` is largest
`LR()`	Largest real part: `real(λ)` is largest
`SR()`	Smallest real part: `real(λ)` is smallest
`LI()`	Largest imaginary part: `imag(λ)` is largest
`SI()`	Smallest imaginary part: `imag(λ)` is smallest

Note

The targets LI() and SI() only make sense in complex arithmetic. In real arithmetic λ is an eigenvalue iff conj(λ) is an eigenvalue and this conjugate pair converges simultaneously.

Return values

The function returns a tuple

decomp, history = partialschur(A, ...)

where decomp is a PartialSchur struct which forms a partial Schur decomposition of A to a prescribed tolerance:

> norm(A * decomp.Q - decomp.Q * decomp.R)

history is a History struct that holds some basic information about convergence of the method:

> history.converged
true
> @show history
Converged after 359 matrix-vector products

Advanced usage

Further there are advanced keyword arguments for tuning the algorithm:

Keyword	Type	Default	Description
`mindim`	`Int`	`min(max(10, nev), size(A,1))`	Minimum Krylov dimension (≥ nev)
`maxdim`	`Int`	`min(max(20, 2nev), size(A,1))`	Maximum Krylov dimension (≥ min)
`restarts`	`Int`	`200`	Maximum number of restarts

When the algorithm does not converge, one can increase restarts. When the algorithm converges too slowly, one can play with mindim and maxdim. It is suggested to keep mindim equal to or slightly larger than nev, and maxdim is usually about two times mindim.

source

From a Schur decomposition to an eigendecomposition

The eigenvalues and eigenvectors are obtained from the Schur form with the partialeigen function that is exported by ArnoldiMethod.jl:

ArnoldiMethod.partialeigen — Function.

partialeigen(P::PartialSchur) → (Vector{<:Union{Real,Complex}}, Matrix{<:Union{Real,Complex}})

Transforms a partial Schur decomposition into an eigendecomposition.

Note

For real-symmetric and Hermitian matrices the Schur vectors coincide with the eigenvectors, and hence it is not necessary to call this function in that case.

The method still relies on LAPACK to compute the eigenvectors of the (quasi) upper triangular matrix R from the partial Schur decomposition.

Note

This method is currently type unstable for real matrices, since we have not yet decided how to deal with complex conjugate pairs of eigenvalues. E.g. if almost all eigenvalues are real, but there are just a few conjugate pairs, should all eigenvectors be complex-valued?

source

Example

Here we compute the first ten eigenvalues and eigenvectors of a tridiagonal sparse matrix.

julia> using ArnoldiMethod, LinearAlgebra, SparseArrays
julia> A = spdiagm(
           -1 => fill(-1.0, 99),
            0 => fill(2.0, 100), 
            1 => fill(-1.0, 99)
       );
julia> decomp, history = partialschur(A, nev=10, tol=1e-6, which=SR());
julia> decomp
PartialSchur decomposition (Float64) of dimension 10
eigenvalues:
10-element Array{Complex{Float64},1}:
 0.0009674354160236865 + 0.0im
  0.003868805732811139 + 0.0im
  0.008701304061962657 + 0.0im
   0.01546025527344699 + 0.0im
  0.024139120518486677 + 0.0im
    0.0347295035554728 + 0.0im
   0.04722115887278571 + 0.0im
   0.06160200160067088 + 0.0im
    0.0778581192025522 + 0.0im
   0.09597378493453936 + 0.0im
julia> history
Converged: 10 of 10 eigenvalues in 174 matrix-vector products
julia> norm(A * decomp.Q - decomp.Q * decomp.R)
6.39386920955869e-8
julia> λs, X = partialeigen(decomp);
julia> norm(A * X - X * Diagonal(λs))
6.393869211477937e-8

The PartialSchur and History structs

For completeness, the return values of the partialschur function:

ArnoldiMethod.PartialSchur — Type.

PartialSchur(Q, R, eigenvalues)

Holds an orthonormal basis Q and a (quasi) upper triangular matrix R.

For convenience the eigenvalues that appear on the diagonal of R are also listed as eigenvalues, which is in particular useful in the case of real matrices with complex eigenvalues. Note that the eigenvalues are always a complex, even when the matrix R is real.

source

ArnoldiMethod.History — Type.

History(mvproducts, nconverged, converged, nev)

History shows whether the method has converged (when nconverged ≥ nev) and how many matrix-vector products were necessary to do so.

source