Support starting vec (#141)

Adds a keyword arg `partialschur(...; v1)` which the vector that induces the Krylov subspace.

Co-authored-by: Samuel Powell <samuelpowell@users.noreply.github.com>

Author: haampie

src/ArnoldiMethod.jl

@@ -8,24 +8,45 @@ using Base: RefValue, OneTo
 export partialschur, LM, SR, LR, SI, LI, partialeigen
 
 """
-    Arnoldi(n, k) → Arnoldi
+ArnoldiWorkspace(n, k) → ArnoldiWorkspace
 
-Pre-allocated Arnoldi relation of the Vₖ₊₁ and Hₖ matrices that satisfy
-A * Vₖ = Vₖ₊₁ * Hₖ, where Vₖ₊₁ is orthonormal of size n × (k+1) and Hₖ upper 
-Hessenberg of size (k+1) × k. The constructor will just allocate sufficient
-space, but will *not* initialize the first vector of `v₁`. For the latter see
-`reinitialize!`.
+Holds the large arrays for the Arnoldi relation: Vₖ₊₁ and Hₖ are matrices that
+satisfy A * Vₖ = Vₖ₊₁ * Hₖ, where Vₖ₊₁ is orthonormal of size n × (k+1) and Hₖ upper 
+Hessenberg of size (k+1) × k. The matrices V_tmp and Q are used for restarts, and
+have similar size as Vₖ₊₁ and Hₖ (but Q is k × k, not k+1 × k).
 """
-struct Arnoldi{T,TV<:StridedMatrix{T},TH<:StridedMatrix{T}}
+struct ArnoldiWorkspace{T,TV<:AbstractMatrix{T},TH<:AbstractMatrix{T}}
+    # Orthonormal basis of the Krylov subspace.
     V::TV
+
+    # The Hessenberg matrix of the Arnoldi relation.
     H::TH
 
-    function Arnoldi{T}(matrix_order::Int, krylov_dimension::Int) where {T}
+    # Temporary matrix similar to V, used to restart.
+    V_tmp::TV
+
+    # Unitary matrix size of (square) H to do a change of basis.
+    Q::TH
+
+    function ArnoldiWorkspace(::Type{T}, matrix_order::Int, krylov_dimension::Int) where {T}
+        # Without an initial vector.
         krylov_dimension <= matrix_order ||
             throw(ArgumentError("Krylov dimension should be less than matrix order."))
         V = Matrix{T}(undef, matrix_order, krylov_dimension + 1)
+        V_tmp = similar(V)
         H = zeros(T, krylov_dimension + 1, krylov_dimension)
-        return new{T,typeof(V),typeof(H)}(V, H)
+        Q = similar(H, krylov_dimension, krylov_dimension)
+        return new{T,typeof(V),typeof(H)}(V, H, V_tmp, Q)
+    end
+
+    function ArnoldiWorkspace(v1::AbstractVector{T}, krylov_dimension::Int) where {T}
+        # From an initial vector v1.
+        V = similar(v1, length(v1), krylov_dimension + 1)
+        V_tmp = similar(V)
+        H = similar(V, krylov_dimension + 1, krylov_dimension)
+        fill!(H, zero(eltype(H)))
+        Q = similar(H, krylov_dimension, krylov_dimension)
+        return new{T,typeof(V),typeof(H)}(V, H, V_tmp, Q)
     end
 end
 

src/expansion.jl

@@ -9,12 +9,16 @@ Generate a random `j+1`th column orthonormal against V[:,1:j]
 Returns true if the column is a valid new basis vector.
 Returns false if the column is numerically in the span of the previous vectors.
 """
-function reinitialize!(arnoldi::Arnoldi{T}, j::Int = 0) where {T}
+function reinitialize!(
+    arnoldi::ArnoldiWorkspace{T},
+    j::Int = 0,
+    populate! = rand!,
+) where {T}
     V = arnoldi.V
     v = view(V, :, j + 1)
 
     # Generate a new random column
-    rand!(v)
+    populate!(v)
 
     # Norm before orthogonalization
     rnorm = norm(v)
@@ -62,7 +66,7 @@ Orthogonalize arnoldi.V[:, j+1] against arnoldi.V[:, 1:j].
 Returns true if the column is a valid new basis vector.
 Returns false if the column is numerically in the span of the previous vectors.
 """
-function orthogonalize!(arnoldi::Arnoldi{T}, j::Integer) where {T}
+function orthogonalize!(arnoldi::ArnoldiWorkspace{T}, j::Integer) where {T}
     V = arnoldi.V
     H = arnoldi.H
 
@@ -109,7 +113,7 @@ end
 
 Perform Arnoldi from `from` to `to`.
 """
-function iterate_arnoldi!(A, arnoldi::Arnoldi{T}, range::UnitRange{Int}) where {T}
+function iterate_arnoldi!(A, arnoldi::ArnoldiWorkspace{T}, range::UnitRange{Int}) where {T}
     V, H = arnoldi.V, arnoldi.H
 
     for j in range

src/run.jl

@@ -13,7 +13,7 @@ end
 
 """
 ```julia
-partialschur(A; nev, which, tol, mindim, maxdim, restarts) → PartialSchur, History
+partialschur(A; v1, workspace, nev, which, tol, mindim, maxdim, restarts) → PartialSchur, History
 ```
 
 Find `nev` approximate eigenpairs of `A` with eigenvalues near a specified target.
@@ -33,6 +33,7 @@ The most important keyword arguments:
 | `nev` | `Int` | `min(6, size(A, 1))` |Number of eigenvalues |
 | `which` | `Symbol` or `Target` | `:LM` | One of `:LM`, `:LR`, `:SR`, `:LI`, `:SI`, see below. |
 | `tol` | `Real` | `√eps` | Tolerance for convergence: ‖Ax - xλ‖₂ < tol * ‖λ‖ |
+| `v1` | `AbstractVector` | `nothing` | Optional starting vector for the Krylov subspace |
 
 The target `which` can be any of:
 
@@ -93,6 +94,7 @@ is usually about two times `mindim`.
 """
 function partialschur(
     A;
+    v1::Union{AbstractVector,Nothing} = nothing,
     nev::Int = min(6, size(A, 1)),
     which::Union{Target,Symbol} = LM(),
     tol::Real = sqrt(eps(real(vtype(A)))),
@@ -101,14 +103,24 @@ function partialschur(
     restarts::Int = 200,
 )
     s = checksquare(A)
-    if nev < 1
-        throw(ArgumentError("nev cannot be less than 1"))
-    end
+    nev < 1 && throw(ArgumentError("nev cannot be less than 1"))
     nev ≤ mindim ≤ maxdim ≤ s || throw(
-        ArgumentError("nev ≤ mindim ≤ maxdim does not hold, got $nev ≤ $mindim ≤ $maxdim"),
+        ArgumentError(
+            "nev ≤ mindim ≤ maxdim ≤ size(A, 1) does not hold, got $nev ≤ $mindim ≤ $maxdim ≤ $s",
+        ),
     )
     _which = which isa Target ? which : _symbol_to_target(which)
-    _partialschur(A, vtype(A), mindim, maxdim, nev, tol, restarts, _which)
+
+    if v1 === nothing
+        arnoldi = ArnoldiWorkspace(vtype(A), size(A, 1), maxdim)
+        reinitialize!(arnoldi, 0, v -> rand!(v))
+    else
+        length(v1) == size(A, 1) ||
+            throw(ArgumentError("v1 should have the same dimension as A"))
+        arnoldi = ArnoldiWorkspace(v1, maxdim)
+        reinitialize!(arnoldi, 0, v -> copyto!(v, v1))
+    end
+    _partialschur(A, arnoldi, mindim, maxdim, nev, tol, restarts, _which)
 end
 
 _symbol_to_target(sym::Symbol) =
@@ -156,33 +168,25 @@ end
 
 function _partialschur(
     A,
-    ::Type{T},
+    arnoldi::ArnoldiWorkspace{T},
     mindim::Int,
     maxdim::Int,
     nev::Int,
     tol::Ttol,
     restarts::Int,
     which::Target,
 ) where {T,Ttol<:Real}
-    n = size(A, 1)
-
-    # Pre-allocated Arnoldi decomp
-    arnoldi = Arnoldi{T}(n, maxdim)
-
     # Unpack for convenience
     H = arnoldi.H
     V = arnoldi.V
+    Vtmp = arnoldi.V_tmp
+    Q = arnoldi.Q
 
-    # For a change of basis we have Vtmp as working space
-    Vtmp = Matrix{T}(undef, n, maxdim)
-
-    # Unitary matrix used for change of basis of V.
-    Q = Matrix{T}(undef, maxdim, maxdim)
-
-    # We only need to store one eignvector of the Hessenberg matrix
+    # We only need to store one eigenvector of the Hessenberg matrix.
     x = zeros(complex(T), maxdim)
 
     # And we store the reflector to transform H back to Hessenberg separately
+    # TODO: can be in-place in H now.
     G = Reflector{T}(maxdim)
 
     # Approximate residual norms for all Ritz values, and Ritz values
@@ -205,7 +209,6 @@ function _partialschur(
     prods = mindim
 
     # Initialize an Arnoldi relation of size `mindim`
-    reinitialize!(arnoldi)
     iterate_arnoldi!(A, arnoldi, 1:mindim)
 
     for iter = 1:restarts

test/expansion.jl

@@ -1,10 +1,10 @@
 # Tests the Arnoldi relation AV = VH when expanding the search subspace
 
 using Test, LinearAlgebra, SparseArrays
-using ArnoldiMethod: reinitialize!, Arnoldi, iterate_arnoldi!
+using ArnoldiMethod: reinitialize!, ArnoldiWorkspace, iterate_arnoldi!
 
 @testset "Initialization" begin
-    arnoldi = Arnoldi{Float64}(5, 3)
+    arnoldi = ArnoldiWorkspace(Float64, 5, 3)
     reinitialize!(arnoldi)
     @test norm(arnoldi.V[:, 1]) ≈ 1
 end
@@ -15,7 +15,7 @@ end
     for T in (Float64, BigFloat)
         A = sprand(T, n, n, 0.1) + I
 
-        arnoldi = Arnoldi{T}(n, max)
+        arnoldi = ArnoldiWorkspace(T, n, max)
         reinitialize!(arnoldi)
         V, H = arnoldi.V, arnoldi.H
 
@@ -40,7 +40,7 @@ end
         ]
 
         # and an initial vector [1; 0; ... 0]
-        vh = Arnoldi{T}(8, 5)
+        vh = ArnoldiWorkspace(T, 8, 5)
         V, H = vh.V, vh.H
         V[:, 1] .= zero(T)
         V[1, 1] = one(T)

test/householder.jl

@@ -70,7 +70,7 @@ end
     Q = qr(randn(T, k, k)).Q * Matrix(1.0I, k, k)
 
     A = rand(T, n, n)
-    arn = Arnoldi{T}(n, k)
+    arn = ArnoldiWorkspace(T, n, k)
     reinitialize!(arn)
     iterate_arnoldi!(A, arn, 1:k)
 

test/partial_schur.jl

@@ -61,6 +61,21 @@ end
     @test_throws ArgumentError partialschur(A, nev = 10)
 end
 
+
+@testset "Eigenvector as initial vector is not problematic" begin
+    A = rand(30, 30)
+    A += A'
+    λs, X = eigen(Symmetric(A))  # ensure real eigenvectors
+
+    λ, x = λs[end], X[:, end]
+    decomp, history = partialschur(A, v1 = x, nev = 2, tol = 1e-8)
+
+    @test history.converged
+    @test norm(A * decomp.Q - decomp.Q * decomp.R) < 1e-7
+    @test abs(maximum(real(decomp.eigenvalues)) - λ) < 1e-7
+end
+
+
 @testset "Target non-dominant eigenvalues" begin
     # Dominant eigenvalues 50, 51, 52, 53, but we target the smallest real part
     A = Diagonal([1:0.1:10; 50:53])