ArnoldiMethod.jl

https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable
Admin
ArnoldiMethod.jl provides an iterative, matrix-free method to compute approximate solutions to...
Tokens:5,084
Snippets:48
Trust Score:9.7
Update:1 month ago
Show doc for...
Tokens:
Raw
### Install ArnoldiMethod.jl

Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable

Use the Julia package manager to add the ArnoldiMethod package.

```julia
(v1.6) pkg> add ArnoldiMethod
```

--------------------------------

### Allocate ArnoldiWorkspace examples

Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable

Demonstrates various ways to initialize the Arnoldi workspace, including automatic and manual memory allocation.

```julia
# allocates workspace for 20-dimensional Krylov subspace
arnoldi = ArnoldiWorkspace(100, 20)

# allocate workspace for 20-dimensional Krylov subspace, with initial vector ones(100) copied into
# the first column of V
arnoldi = ArnoldiWorkspace(ones(100), 20)

# manually allocate workspace with V, H
V = Matrix{Float64}(undef, 100, 21)
H = Matrix{Float64}(undef, 21, 20)
arnoldi = ArnoldiWorkspace(V, H)

# manually allocate all arrays, including temporaries
V, tmp = Matrix{Float64}(undef, 100, 21), Matrix{Float64}(undef, 100, 21)
H, Q = Matrix{Float64}(undef, 21, 20), Matrix{Float64}(undef, 20, 20)
arnoldi = ArnoldiWorkspace(V, H, V_tmp = tmp, Q = Q)
```

--------------------------------

### Load ArnoldiMethod.jl

Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable

Import the package into your Julia session.

```julia
using ArnoldiMethod
```

--------------------------------

### Resuming partial Schur decomposition

Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable

Demonstrates initializing an ArnoldiWorkspace and performing iterative partial Schur decomposition updates.

```julia
A = rand(100, 100)

# Pre-allocate the relevant Krylov subspace arrays
V, H = rand(100, 21), rand(21, 20)
arnoldi = ArnoldiWorkspace(V, H)

# Find a few eigenvalues
_, info_1 = partialschur!(A, arnoldi, nev = 3, tol = 1e-12)

# Then continue to find a couple more. Notice: 5 in total, so 2 more. Allow larger errors by
# changin `tol`.
F, info_2 = partialschur!(A, arnoldi, nev = 5, start_from = info_1.nconverged + 1 , tol = 1e-8)
@show norm(A * F.Q - F.Q * F.R)
```

--------------------------------

### ArnoldiWorkspace

Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable

Constructor for pre-allocating memory for the Arnoldi relation matrices.

```APIDOC
## ArnoldiWorkspace(n, k) / ArnoldiWorkspace(v1, k) / ArnoldiWorkspace(V, H; V_tmp, Q)

### Description
Holds the large arrays for the Arnoldi relation: Vₖ₊₁ and Hₖ. Used for pre-allocating memory or working with custom matrix types.

### Parameters
- **n** (Int) - Required - Dimension of the Krylov subspace.
- **k** (Int) - Required - Size of the subspace.
- **v1** (Vector) - Optional - Initial vector.
- **V, H, V_tmp, Q** (Matrix) - Optional - Manually allocated matrices for the Arnoldi relation.
```

--------------------------------

### partialschur(A; v1, workspace, nev, which, tol, mindim, maxdim, restarts)

Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable

Computes an approximate, partial Schur decomposition of a matrix A. It runs iteratively until eigenpairs are approximated to the prescribed tolerance or restarts are exhausted.

```APIDOC
## partialschur(A; v1, workspace, nev, which, tol, mindim, maxdim, restarts)

### Description
Finds `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`.

### Parameters
- **A** (Any) - Required - The linear map/matrix to decompose.
- **nev** (Int) - Optional - Number of eigenvalues (Default: min(6, size(A, 1))).
- **which** (Symbol or Target) - Optional - Target eigenvalue selection (Default: :LM).
- **tol** (Real) - Optional - Tolerance for convergence (Default: √eps).
- **v1** (AbstractVector) - Optional - Starting vector for the Krylov subspace.
- **mindim** (Int) - Optional - Minimum Krylov dimension.
- **maxdim** (Int) - Optional - Maximum Krylov dimension.
- **restarts** (Int) - Optional - Maximum number of restarts (Default: 200).

### Returns
- **decomp** (PartialSchur) - The partial Schur decomposition struct.
- **history** (History) - Information about the convergence of the method.
```

--------------------------------

### Define History struct

Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable

Tracks convergence status and computational effort of the Arnoldi method.

```julia
History(mvproducts, nconverged, converged, nev)
```

--------------------------------

### Compute Eigenvalues and Eigenvectors

Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable

Compute the first ten eigenvalues and eigenvectors of a tridiagonal sparse matrix using partialschur and partialeigen.

```julia
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
```

--------------------------------

### Check Convergence History

Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable

Access the convergence status and information from the history object.

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

--------------------------------

### partialschur!(A, arnoldi; start_from, initialize, nev, which, tol, mindim, maxdim, restarts)

Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable

A variant of partialschur that operates on a pre-allocated ArnoldiWorkspace instance to avoid allocations or continue from an existing decomposition.

```APIDOC
## partialschur!(A, arnoldi; start_from, initialize, nev, which, tol, mindim, maxdim, restarts)

### Description
Operates on a pre-allocated `ArnoldiWorkspace` instance. Useful for continuing from an existing decomposition, using custom array types, or avoiding repeated allocations.

### Parameters
- **A** (Any) - Required - The linear map/matrix.
- **arnoldi** (ArnoldiWorkspace) - Required - Pre-allocated workspace.
- **start_from** (Int) - Optional - Number of Schur vectors to continue from.
- **initialize** (Bool) - Optional - Whether to initialize the Krylov subspace (Default: true).
- **nev** (Int) - Optional - Number of eigenvalues.
- **which** (Symbol or Target) - Optional - Target eigenvalue selection.
- **tol** (Real) - Optional - Tolerance for convergence.
- **mindim** (Int) - Optional - Minimum Krylov dimension.
- **maxdim** (Int) - Optional - Maximum Krylov dimension.
- **restarts** (Int) - Optional - Maximum number of restarts.

### Returns
- **decomp** (PartialSchur) - The partial Schur decomposition struct containing views of the workspace.
- **history** (History) - Information about the convergence of the method.
```

--------------------------------

### Shift-and-invert for smallest eigenvalues

Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable

Uses LinearMaps.jl to invert matrix A, allowing the solver to find the largest eigenvalues of A⁻¹, which correspond to the smallest eigenvalues of A.

```julia
using ArnoldiMethod, LinearAlgebra, LinearMaps

# Define a matrix whose eigenvalues you want
A = rand(100,100)

# Factorizes A and builds a linear map that applies inv(A) to a vector.
function construct_linear_map(A)
    F = factorize(A)
    LinearMap{eltype(A)}((y, x) -> ldiv!(y, F, x), size(A,1), ismutating=true)
end

# Target the largest eigenvalues of the inverted problem
decomp, = partialschur(construct_linear_map(A), nev=4, tol=1e-5, restarts=100, which=:LM)
λs_inv, X = partialeigen(decomp)

# Eigenvalues have to be inverted to find the smallest eigenvalues of the non-inverted problem.
λs = 1 ./ λs_inv
 
# Show that Ax = xλ
@show norm(A * X - X * Diagonal(λs)) # 7.38473677258669e-6
```

--------------------------------

### Shift-and-invert for generalized eigenvalue problems

Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable

Recasts Ax=Bxλ into A⁻¹Bx=xθ to target eigenvalues closest to the origin.

```julia
using ArnoldiMethod, LinearAlgebra, LinearMaps

# Define the matrices of the generalized eigenvalue problem
A, B = rand(100, 100), rand(100, 100)

struct ShiftAndInvert{TA,TB,TT}
    A_lu::TA
    B::TB
    temp::TT
end

function (M::ShiftAndInvert)(y, x)
    mul!(M.temp, M.B, x)
    ldiv!(y, M.A_lu, M.temp)
end

function construct_linear_map(A, B)
    a = ShiftAndInvert(factorize(A), B, Vector{eltype(A)}(undef, size(A, 1)))
    LinearMap{eltype(A)}(a, size(A, 1), ismutating = true)
end

# Target the largest eigenvalues of the inverted problem
decomp, = partialschur(
    construct_linear_map(A, B),
    nev = 4,
    tol = 1e-5,
    restarts = 100,
    which = :LM,
)
λs_inv, X = partialeigen(decomp)

# Eigenvalues have to be inverted to find the smallest eigenvalues of the non-inverted problem.
λs = 1 ./ λs_inv

# Show that Ax = Bxλ
@show norm(A * X - B * X * Diagonal(λs)) # 2.8043149027575927e-6
```

--------------------------------

### partialschur

Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable

Computes the partial Schur decomposition of a matrix A. Users can provide an initial guess vector v1 to speed up convergence.

```APIDOC
## partialschur(A, v1=my_initial_vector)

### Description
Computes the partial Schur decomposition of a matrix A. An optional initial guess vector v1 can be provided to build the Krylov subspace.

### Parameters
- **A** (Matrix) - Required - The matrix to decompose.
- **v1** (Vector) - Optional - Initial guess for the target eigenvector.
```

--------------------------------

### partialeigen(P::PartialSchur)

Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable

Transforms a partial Schur decomposition into an eigendecomposition by computing the eigendecomposition of the upper triangular matrix R.

```APIDOC
## partialeigen(P::PartialSchur)

### Description
Transforms a partial Schur decomposition into an eigendecomposition. This is a post-processing step that computes the eigendecomposition of the upper triangular matrix R from the partial Schur decomposition.

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

### Parameters
- **P** (PartialSchur) - Required - The partial Schur decomposition object to be transformed.

### Returns
- **Vector{<:Union{Real,Complex}}** - The eigenvalues.
- **Matrix{<:Union{Real,Complex}}** - The eigenvectors.
```

--------------------------------

### Verify Decomposition Norm

Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable

Check the residual norm of the computed partial Schur decomposition.

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

--------------------------------

### Partial Schur Decomposition Function Signature

Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable

The primary interface for computing a partial Schur decomposition of a linear map A.

```julia
partialschur(A; v1, workspace, nev, which, tol, mindim, maxdim, restarts) → PartialSchur, History
```

--------------------------------

### Partial Schur Return Values

Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable

The function returns a tuple containing the decomposition structure and the convergence history.

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

--------------------------------

### Define PartialSchur struct

Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable

Represents the orthonormal basis and quasi-upper triangular matrix resulting from the partial Schur decomposition.

```julia
PartialSchur(Q, R, eigenvalues)
```

--------------------------------

### Generalized eigenvalue problem with symmetric positive-definite B

Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable

Uses Cholesky decomposition of B to transform the generalized problem into a standard Schur decomposition.

```julia
using ArnoldiMethod, LinearAlgebra, LinearMaps
struct WithBInnerProduct{TA,TL}
    A::TA
    L::TL
end

function (M::WithBInnerProduct)(y, x)
    # Julia unfortunately does not have in-place CHOLMOD solve, so this is far from optimal.
    tmp = M.L \ x
    mul!(y, M.A, tmp)
    y .= M.L' \ y
    return y
end

# Define the matrices of the generalized eigenvalue problem
A = rand(100, 100)
B = Diagonal(range(1.0, 2.0, length = 100))

# Reformulate the problem as standard Schur decomposition
F = cholesky(B)
linear_map = LinearMap{eltype(A)}(WithBInnerProduct(A, F.L), size(A, 1), ismutating = true)
decomp, info = partialschur(linear_map, nev = 4, which = :LM, tol = 1e-10)

# Translate back to the original problem
Q = F.L' \ decomp.Q

@show norm(Q' * A * Q - decomp.R)  # 3.883933945390996e-14
@show norm(Q' * B * Q - I)  # 3.1672155003480583e-15
```

--------------------------------

### Compute partial eigendecomposition

Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable

Transforms a PartialSchur object into a tuple containing eigenvalues and eigenvectors. Note that this is unnecessary for real-symmetric or complex-Hermitian matrices.

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

--------------------------------

### In-place Partial Schur Decomposition

Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable

A variant of partialschur that uses a pre-allocated ArnoldiWorkspace to avoid allocations and support custom array types.

```julia
partialschur!(A, arnoldi; start_from, initialize, nev, which, tol, mindim, maxdim, restarts) → PartialSchur, History
```

--------------------------------

### Define ArnoldiWorkspace struct

Source: https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable

Constructs the workspace for storing large arrays used in the Arnoldi relation.

```julia
ArnoldiWorkspace(n, k) → ArnoldiWorkspace
ArnoldiWorkspace(v1, k) → ArnoldiWorkspace
ArnoldiWorkspace(V, H; V_tmp, Q) → ArnoldiWorkspace
```
Tokens:
Raw