From f33c15d16166a2b4017a4b012bc4987b3b5c73f9 Mon Sep 17 00:00:00 2001 From: Simeon David Schaub Date: Thu, 5 Feb 2026 18:34:50 +0100 Subject: [PATCH 01/13] [RFC] Symplectic Arnoldi See the discussion in #150. On my test case, this already gives me significantly better performance, and improved accuracy, if using re-skew-orthogonalization. An initial sketch was written by Copilot, but I revised and reviewed the result extensively. closes #150 --- src/KrylovKit.jl | 7 +- src/algorithms.jl | 86 ++++++++++- src/factorizations/arnoldi.jl | 106 +++++++++----- src/innerproductvec.jl | 99 +++++++------ src/orthonormal.jl | 266 ++++++++++++++++++++++++++++++++-- test/runtests.jl | 3 + test/symplectic_arnoldi.jl | 53 +++++++ 7 files changed, 522 insertions(+), 98 deletions(-) create mode 100644 test/symplectic_arnoldi.jl diff --git a/src/KrylovKit.jl b/src/KrylovKit.jl index 27d8e615..74596171 100644 --- a/src/KrylovKit.jl +++ b/src/KrylovKit.jl @@ -32,14 +32,19 @@ export linsolve, reallinsolve, lssolve, reallssolve export eigsolve, geneigsolve, realeigsolve, schursolve, svdsolve, bieigsolve export exponentiate, expintegrator export orthogonalize, orthogonalize!!, orthonormalize, orthonormalize!! +export skeworthogonalize, skeworthogonalize!!, skeworthonormalize, skeworthonormalize!! export basis, rayleighquotient, residual, normres, rayleighextension export initialize, initialize!, expand!, shrink! export ClassicalGramSchmidt, ClassicalGramSchmidt2, ClassicalGramSchmidtIR export ModifiedGramSchmidt, ModifiedGramSchmidt2, ModifiedGramSchmidtIR +export ClassicalSymplecticGramSchmidt, ClassicalSymplecticGramSchmidt2 +export ClassicalSymplecticGramSchmidtIR +export ModifiedSymplecticGramSchmidt, ModifiedSymplecticGramSchmidt2 +export ModifiedSymplecticGramSchmidtIR export LanczosIterator, BlockLanczosIterator, ArnoldiIterator, GKLIterator, BiArnoldiIterator export CG, GMRES, BiCGStab, Lanczos, BlockLanczos, Arnoldi, GKL, GolubYe, LSMR, BiArnoldi export KrylovDefaults, EigSorter -export RecursiveVec, InnerProductVec, Block +export RecursiveVec, InnerProductVec, Block, SymplecticBasis, numpairs # Multithreading const _NTHREADS = Ref(1) diff --git a/src/algorithms.jl b/src/algorithms.jl index 6aa2cdcc..d22bfb1c 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -79,6 +79,82 @@ struct ModifiedGramSchmidtIR{S <: Real} <: Reorthogonalizer end ModifiedGramSchmidtIR() = ModifiedGramSchmidtIR(1 / sqrt(2)) # Daniel-Gragg-Kaufman-Stewart +# Skew-orthogonalization for symplectic bases +""" + abstract type SkewOrthogonalizer + +Supertype for representing different skew-orthogonalization strategies or algorithms that +produce symplectic (Darboux) bases. A symplectic basis satisfies the relations +`⟨u_{2m-1}, u_{2n}⟩ = δ_{mn}`, `⟨u_{2m}, u_{2n}⟩ = 0`, and `⟨u_{2m-1}, u_{2n-1}⟩ = 0`. + +See also: [`ClassicalSymplecticGramSchmidt`](@ref), [`ModifiedSymplecticGramSchmidt`](@ref), +[`ClassicalSymplecticGramSchmidt2`](@ref), [`ModifiedSymplecticGramSchmidt2`](@ref), +[`ClassicalSymplecticGramSchmidtIR`](@ref), [`ModifiedSymplecticGramSchmidtIR`](@ref). +""" +abstract type SkewOrthogonalizer end +abstract type SkewReorthogonalizer <: SkewOrthogonalizer end + +# Simple +""" + ClassicalSymplecticGramSchmidt() + +Represents the classical symplectic Gram Schmidt algorithm for skew-orthogonalizing vectors +to produce a symplectic (Darboux) basis, typically not an optimal choice. +""" +struct ClassicalSymplecticGramSchmidt <: SkewOrthogonalizer end + +""" + ModifiedSymplecticGramSchmidt() + +Represents the modified symplectic Gram Schmidt algorithm for skew-orthogonalizing vectors +to produce a symplectic (Darboux) basis. +""" +struct ModifiedSymplecticGramSchmidt <: SkewOrthogonalizer end + +# A single reorthogonalization always +""" + ClassicalSymplecticGramSchmidt2() + +Represents the classical symplectic Gram Schmidt algorithm with a second reskew- +orthogonalization step always taking place. +""" +struct ClassicalSymplecticGramSchmidt2 <: SkewReorthogonalizer end + +""" + ModifiedSymplecticGramSchmidt2() + +Represents the modified symplectic Gram Schmidt algorithm with a second reskew- +orthogonalization step always taking place. +""" +struct ModifiedSymplecticGramSchmidt2 <: SkewReorthogonalizer end + +# Iterative reorthogonalization +""" + ClassicalSymplecticGramSchmidtIR(η::Real = 1/sqrt(2)) + +Represents the classical symplectic Gram Schmidt algorithm with iterative (i.e. zero or +more) reskew-orthogonalization until the norm of the vector after a skew-orthogonalization +step has not decreased by a factor smaller than `η` with respect to the norm before the +step. The default value corresponds to the Daniel-Gragg-Kaufman-Stewart condition. +""" +struct ClassicalSymplecticGramSchmidtIR{S <: Real} <: SkewReorthogonalizer + η::S +end +ClassicalSymplecticGramSchmidtIR() = ClassicalSymplecticGramSchmidtIR(1 / sqrt(2)) + +""" + ModifiedSymplecticGramSchmidtIR(η::Real = 1/sqrt(2)) + +Represents the modified symplectic Gram Schmidt algorithm with iterative (i.e. zero or +more) reskew-orthogonalization until the norm of the vector after a skew-orthogonalization +step has not decreased by a factor smaller than `η` with respect to the norm before the +step. The default value corresponds to the Daniel-Gragg-Kaufman-Stewart condition. +""" +struct ModifiedSymplecticGramSchmidtIR{S <: Real} <: SkewReorthogonalizer + η::S +end +ModifiedSymplecticGramSchmidtIR() = ModifiedSymplecticGramSchmidtIR(1 / sqrt(2)) + # Solving eigenvalue problems abstract type KrylovAlgorithm end @@ -145,7 +221,7 @@ The initial block size determines the maximum multiplicity of the target eigenva The iteration stops when either the norm of the residual is below `tol` or a sufficient number of eigenvectors have converged. [Reference](https://www.netlib.org/utk/people/JackDongarra/etemplates/node250.html) -Use `Arnoldi` for non-symmetric or non-Hermitian linear operators. +Use `Arnoldi` for non-symmetric or non-Hermitian linear operators. See also: [Factorization types](@ref), [`eigsolve`](@ref), [`Arnoldi`](@ref), [`Orthogonalizer`](@ref) """ @@ -357,7 +433,7 @@ end """ GMRES(; krylovdim=KrylovDefaults.krylovdim[], maxiter=KrylovDefaults.maxiter[], - tol=KrylovDefaults.tol[], + tol=KrylovDefaults.tol[], orth::Orthogonalizer=KrylovDefaults.orth, verbosity=KrylovDefaults.verbosity[]) @@ -460,7 +536,7 @@ end Construct an instance of the Biconjugate gradient algorithm with specified parameters, which can be passed to `linsolve` in order to iteratively solve a linear system general linear map. The `BiCGStab` method will search for the optimal `x` in a Krylov subspace - of maximal size `maxiter`, or stop when `norm(A*x - b) < tol`. The default verbosity level + of maximal size `maxiter`, or stop when `norm(A*x - b) < tol`. The default verbosity level `verbosity` amounts to printing warnings upon lack of convergence. See also: [`linsolve`](@ref), [`GMRES`](@ref), [`CG`](@ref), [`BiCG`](@ref), [`LSMR`](@ref), @@ -484,7 +560,7 @@ abstract type LeastSquaresSolver <: KrylovAlgorithm end """ LSMR(; krylovdim=1, maxiter=KrylovDefaults.maxiter[], - tol=KrylovDefaults.tol[], + tol=KrylovDefaults.tol[], orth::Orthogonalizer=ModifiedGramSchmidt(), verbosity=KrylovDefaults.verbosity[]) @@ -495,7 +571,7 @@ algebraically equivalent to applying MINRES to the normal equations ``(A^*A + λ^2I)x = A^*b``, but has better numerical properties, especially if ``A`` is ill-conditioned. -The `LSMR` method will search for the optimal ``x`` in a Krylov subspace of maximal size +The `LSMR` method will search for the optimal ``x`` in a Krylov subspace of maximal size `maxiter`, or stop when ``norm(A'*(A*x - b) + λ^2 * x) < tol``. The parameter `krylovdim` does in this case not indicate that a subspace of that size will be built, but represents the number of most recent vectors that will be kept to which the next vector will be reorthogonalized. diff --git a/src/factorizations/arnoldi.jl b/src/factorizations/arnoldi.jl index 5b1ead15..8456056d 100644 --- a/src/factorizations/arnoldi.jl +++ b/src/factorizations/arnoldi.jl @@ -1,5 +1,5 @@ """ - mutable struct ArnoldiFactorization{T,S} <: KrylovFactorization{T,S} + mutable struct ArnoldiFactorization{T,S,B} <: KrylovFactorization{T,S} Structure to store an Arnoldi factorization of a linear map `A` of the form @@ -9,15 +9,15 @@ A * V = V * B + r * b' For a given Arnoldi factorization `fact` of length `k = length(fact)`, the basis `V` is obtained via [`basis(fact)`](@ref basis) and is an instance of [`OrthonormalBasis{T}`](@ref -Basis), with also `length(V) == k` and where `T` denotes the type of vector like objects -used in the problem. The Rayleigh quotient `B` is obtained as -[`rayleighquotient(fact)`](@ref) and is of type [`B::PackedHessenberg{S<:Number}`](@ref -PackedHessenberg) with `size(B) == (k,k)`. The residual `r` is obtained as -[`residual(fact)`](@ref) and is of type `T`. One can also query [`normres(fact)`](@ref) to -obtain `norm(r)`, the norm of the residual. The vector `b` has no dedicated name but can be -obtained via [`rayleighextension(fact)`](@ref). It takes the default value ``e_k``, i.e. the -unit vector of all zeros and a one in the last entry, which is represented using -[`SimpleBasisVector`](@ref). +Basis) or [`SymplecticBasis{T}`](@ref) (when using a [`SkewOrthogonalizer`](@ref)), with +also `length(V) == k` and where `T` denotes the type of vector like objects used in the +problem. The Rayleigh quotient `B` is obtained as [`rayleighquotient(fact)`](@ref) and is +of type [`B::PackedHessenberg{S<:Number}`](@ref PackedHessenberg) with `size(B) == (k,k)`. +The residual `r` is obtained as [`residual(fact)`](@ref) and is of type `T`. One can also +query [`normres(fact)`](@ref) to obtain `norm(r)`, the norm of the residual. The vector `b` +has no dedicated name but can be obtained via [`rayleighextension(fact)`](@ref). It takes +the default value ``e_k``, i.e. the unit vector of all zeros and a one in the last entry, +which is represented using [`SimpleBasisVector`](@ref). An Arnoldi factorization `fact` can be destructured as `V, B, r, nr, b = fact` with `nr = norm(r)`. @@ -28,9 +28,9 @@ Arnoldi factorizations of a given linear map and a starting vector. See [`LanczosFactorization`](@ref) and [`LanczosIterator`](@ref) for a Krylov factorization that is optimized for real symmetric or complex hermitian linear maps. """ -mutable struct ArnoldiFactorization{T, S} <: KrylovFactorization{T, S} +mutable struct ArnoldiFactorization{T, S, B <: Basis{T}} <: KrylovFactorization{T, S} k::Int # current Krylov dimension - V::OrthonormalBasis{T} # basis of length k + V::B # basis of length k (OrthonormalBasis or SymplecticBasis) H::Vector{S} # stores the Hessenberg matrix in packed form r::T # residual end @@ -52,8 +52,8 @@ rayleighextension(F::ArnoldiFactorization) = SimpleBasisVector(F.k, F.k) # Arnoldi iteration for constructing the orthonormal basis of a Krylov subspace. """ - struct ArnoldiIterator{F,T,O<:Orthogonalizer} <: KrylovIterator{F,T} - ArnoldiIterator(f, v₀, [orth::Orthogonalizer = KrylovDefaults.orth]) + struct ArnoldiIterator{F,T,O<:Union{Orthogonalizer, SkewOrthogonalizer}} <: KrylovIterator{F,T} + ArnoldiIterator(f, v₀, [orth::Union{Orthogonalizer, SkewOrthogonalizer} = KrylovDefaults.orth]) Iterator that takes a general linear map `f::F` and an initial vector `v₀::T` and generates an expanding `ArnoldiFactorization` thereof. In particular, `ArnoldiIterator` iterates over @@ -63,9 +63,9 @@ progressively expanding Arnoldi factorizations using the The argument `f` can be a matrix, or a function accepting a single argument `v`, so that `f(v)` implements the action of the linear map on the vector `v`. -The optional argument `orth` specifies which [`Orthogonalizer`](@ref) to be used. The -default value in [`KrylovDefaults`](@ref) is to use [`ModifiedGramSchmidtIR`](@ref), which -possibly uses reorthogonalization steps. +The optional argument `orth` specifies which [`Orthogonalizer`](@ref) or +[`SkewOrthogonalizer`](@ref) to be used. The default value in [`KrylovDefaults`](@ref) is +to use [`ModifiedGramSchmidtIR`](@ref), which possibly uses reorthogonalization steps. When iterating over an instance of `ArnoldiIterator`, the values being generated are instances of [`ArnoldiFactorization`](@ref), which can be immediately destructured into a @@ -108,7 +108,7 @@ factorization in place. See also [`initialize!(::KrylovIterator, information will be discarded) and [`shrink!(::KrylovFactorization, k)`](@ref) to shrink an existing factorization down to length `k`. """ -struct ArnoldiIterator{F, T, O <: Orthogonalizer} <: KrylovIterator{F, T} +struct ArnoldiIterator{F, T, O <: Union{Orthogonalizer, SkewOrthogonalizer}} <: KrylovIterator{F, T} operator::F x₀::T orth::O @@ -148,33 +148,40 @@ function initialize(iter::ArnoldiIterator; verbosity::Int = KrylovDefaults.verbo else r = scale!!(Ax₀, 1 / β₀) end - βold = norm(r) - r = add!!(r, v, -α) - β = norm(r) - # possibly reorthogonalize - if iter.orth isa Union{ClassicalGramSchmidt2, ModifiedGramSchmidt2} - dα = inner(v, r) - α += dα - r = add!!(r, v, -dα) + if iter.orth isa Orthogonalizer + βold = norm(r) + r = add!!(r, v, -α) β = norm(r) - elseif iter.orth isa Union{ClassicalGramSchmidtIR, ModifiedGramSchmidtIR} - while eps(one(β)) < β < iter.orth.η * βold - βold = β + # possibly reorthogonalize + if iter.orth isa Union{ClassicalGramSchmidt2, ModifiedGramSchmidt2} dα = inner(v, r) α += dα r = add!!(r, v, -dα) β = norm(r) + elseif iter.orth isa Union{ClassicalGramSchmidtIR, ModifiedGramSchmidtIR} + while eps(one(β)) < β < alg.η * βold + βold = β + dα = inner(v, r) + α += dα + r = add!!(r, v, -dα) + β = norm(r) + end end + V = OrthonormalBasis([v]) + H = T[α, β] + else + iszero(α) && throw(ArgumentError("initial vector and its image are symplectically orthogonal")) + V = SymplecticBasis([v]) + H = T[zero(α), α] end - V = OrthonormalBasis([v]) - H = T[α, β] if verbosity > EACHITERATION_LEVEL @info "Arnoldi initiation at dimension 1: subspace normres = $(normres2string(β))" end - return state = ArnoldiFactorization(1, V, H, r) + return ArnoldiFactorization(1, V, H, r) end + function initialize!( - iter::ArnoldiIterator, state::ArnoldiFactorization; + iter::ArnoldiIterator{<:Any, <:Any, <:Orthogonalizer}, state::ArnoldiFactorization; verbosity::Int = KrylovDefaults.verbosity[] ) x₀ = iter.x₀ @@ -186,10 +193,17 @@ function initialize!( V[1] = scale!!(V[1], x₀, 1 / norm(x₀)) w = apply(iter.operator, V[1]) - r, α = orthogonalize!!(w, V[1], iter.orth) - β = norm(r) + if iter.orth isa Orthogonalizer + r, α = orthogonalize!!(w, V[1], iter.orth) + β = norm(r) + push!(H, α, β) + else + α = inner(V[1], w) + iszero(α) && throw(ArgumentError("initial vector and its image are symplectically orthogonal")) + r = w + push!(H, zero(α), α) + end state.k = 1 - push!(H, α, β) state.r = r if verbosity > EACHITERATION_LEVEL @info "Arnoldi initiation at dimension 1: subspace normres = $(normres2string(β))" @@ -205,7 +219,11 @@ function expand!( V = state.V H = state.H r = state.r - β = normres(state) + if iter.orth isa Orthogonalizer + β = normres(state) + else + β = iseven(k) ? normres(state) : norm(r) + end push!(V, scale(r, 1 / β)) m = length(H) resize!(H, m + k + 1) @@ -227,7 +245,11 @@ function shrink!(state::ArnoldiFactorization, k; verbosity::Int = KrylovDefaults r = pop!(V) resize!(H, (k * k + 3 * k) >> 1) state.k = k - β = normres(state) + if iter.orth isa Orthogonalizer + β = normres(state) + else + β = iseven(k) ? normres(state) : norm(r) + end if verbosity > EACHITERATION_LEVEL @info "Arnoldi reduction to dimension $k: subspace normres = $(normres2string(β))" end @@ -243,3 +265,11 @@ function arnoldirecurrence!!( r, h = orthogonalize!!(w, V, h, orth) return r, norm(r) end + +function arnoldirecurrence!!( + operator, V::SymplecticBasis, h::AbstractVector, orth::SkewOrthogonalizer + ) + w = apply(operator, last(V)) + r, h = skeworthogonalize!!(w, V, h, orth) + return r, iseven(length(V)) ? norm(r) : inner(last(V), r) +end diff --git a/src/innerproductvec.jl b/src/innerproductvec.jl index d04e9df1..dea411c8 100644 --- a/src/innerproductvec.jl +++ b/src/innerproductvec.jl @@ -1,5 +1,5 @@ """ - v = InnerProductVec(vec, dotf) + v = InnerProductVec(vec, dotf, [normf]) Create a new vector `v` from an existing vector `dotf` with a modified inner product given by `inner`. The vector `vec`, which can be any type (not necessarily `Vector`) that supports @@ -9,50 +9,58 @@ with scalars (both out of place and in place using `mul!`, `rmul!`, `axpy!` and applied to `v` is simply forwarded to `v.vec`. The inner product between vectors `v1 = InnerProductVec(vec1, dotf)` and `v2 = InnerProductVec(vec2, dotf)` is computed as `dot(v1, v2) = dotf(v1.vec, v2.vec) = dotf(vec1, vec2)`. The inner product between vectors -with different `dotf` functions is not defined. Similarly, The norm of `v::InnerProductVec` -is defined as `v = sqrt(real(dot(v, v))) = sqrt(real(dotf(vec, vec)))`. +with different `dotf` functions is not defined. + +By default, the norm of `v::InnerProductVec` is defined as +`norm(v) = sqrt(real(dot(v, v))) = sqrt(real(dotf(vec, vec)))`. However, an optional third +argument `normf` can be provided to define a different norm function, which is useful for +cases like symplectic forms where the inner product is skew-symmetric and cannot define a +norm. In this case, `norm(v) = normf(v.vec)`. In a (linear) map applied to `v`, the original vector can be obtained as `v.vec` or simply as `v[]`. """ -struct InnerProductVec{F, T} +struct InnerProductVec{F, N, T} vec::T dotf::F + normf::N end -Base.:-(v::InnerProductVec) = InnerProductVec(-v.vec, v.dotf) -function Base.:+(v::InnerProductVec{F}, w::InnerProductVec{F}) where {F} - return InnerProductVec(v.vec + w.vec, v.dotf) +InnerProductVec(vec, dotf) = InnerProductVec(vec, dotf, nothing) + +Base.:-(v::InnerProductVec) = InnerProductVec(-v.vec, v.dotf, v.normf) +function Base.:+(v::InnerProductVec{F, N}, w::InnerProductVec{F, N}) where {F, N} + return InnerProductVec(v.vec + w.vec, v.dotf, v.normf) end -function Base.:-(v::InnerProductVec{F}, w::InnerProductVec{F}) where {F} - return InnerProductVec(v.vec - w.vec, v.dotf) +function Base.:-(v::InnerProductVec{F, N}, w::InnerProductVec{F, N}) where {F, N} + return InnerProductVec(v.vec - w.vec, v.dotf, v.normf) end -Base.:*(v::InnerProductVec, a::Number) = InnerProductVec(v.vec * a, v.dotf) -Base.:*(a::Number, v::InnerProductVec) = InnerProductVec(a * v.vec, v.dotf) -Base.:/(v::InnerProductVec, a::Number) = InnerProductVec(v.vec / a, v.dotf) -Base.:\(a::Number, v::InnerProductVec) = InnerProductVec(a \ v.vec, v.dotf) +Base.:*(v::InnerProductVec, a::Number) = InnerProductVec(v.vec * a, v.dotf, v.normf) +Base.:*(a::Number, v::InnerProductVec) = InnerProductVec(a * v.vec, v.dotf, v.normf) +Base.:/(v::InnerProductVec, a::Number) = InnerProductVec(v.vec / a, v.dotf, v.normf) +Base.:\(a::Number, v::InnerProductVec) = InnerProductVec(a \ v.vec, v.dotf, v.normf) function Base.similar(v::InnerProductVec, (::Type{T}) = scalartype(v)) where {T} - return InnerProductVec(similar(v.vec), v.dotf) + return InnerProductVec(similar(v.vec), v.dotf, v.normf) end Base.getindex(v::InnerProductVec) = v.vec -function Base.copy!(w::InnerProductVec{F}, v::InnerProductVec{F}) where {F} +function Base.copy!(w::InnerProductVec{F, N}, v::InnerProductVec{F, N}) where {F, N} copy!(w.vec, v.vec) return w end function LinearAlgebra.mul!( - w::InnerProductVec{F}, a::Number, v::InnerProductVec{F} - ) where {F} + w::InnerProductVec{F, N}, a::Number, v::InnerProductVec{F, N} + ) where {F, N} mul!(w.vec, a, v.vec) return w end function LinearAlgebra.mul!( - w::InnerProductVec{F}, v::InnerProductVec{F}, a::Number - ) where {F} + w::InnerProductVec{F, N}, v::InnerProductVec{F, N}, a::Number + ) where {F, N} mul!(w.vec, v.vec, a) return w end @@ -63,75 +71,76 @@ function LinearAlgebra.rmul!(v::InnerProductVec, a::Number) end function LinearAlgebra.axpy!( - a::Number, v::InnerProductVec{F}, w::InnerProductVec{F} - ) where {F} + a::Number, v::InnerProductVec{F, N}, w::InnerProductVec{F, N} + ) where {F, N} axpy!(a, v.vec, w.vec) return w end function LinearAlgebra.axpby!( - a::Number, v::InnerProductVec{F}, b, w::InnerProductVec{F} - ) where {F} + a::Number, v::InnerProductVec{F, N}, b, w::InnerProductVec{F, N} + ) where {F, N} axpby!(a, v.vec, b, w.vec) return w end -function LinearAlgebra.dot(v::InnerProductVec{F}, w::InnerProductVec{F}) where {F} +function LinearAlgebra.dot(v::InnerProductVec{F, N}, w::InnerProductVec{F, N}) where {F, N} return v.dotf(v.vec, w.vec) end -VectorInterface.scalartype(::Type{<:InnerProductVec{F, T}}) where {F, T} = scalartype(T) +VectorInterface.scalartype(::Type{<:InnerProductVec{F, N, T}}) where {F, N, T} = scalartype(T) function VectorInterface.zerovector(v::InnerProductVec, T::Type{<:Number}) - return InnerProductVec(zerovector(v.vec, T), v.dotf) + return InnerProductVec(zerovector(v.vec, T), v.dotf, v.normf) end function VectorInterface.zerovector!(v::InnerProductVec) - return InnerProductVec(zerovector!(v.vec), v.dotf) + return InnerProductVec(zerovector!(v.vec), v.dotf, v.normf) end function VectorInterface.zerovector!!(v::InnerProductVec) - return InnerProductVec(zerovector!!(v.vec), v.dotf) + return InnerProductVec(zerovector!!(v.vec), v.dotf, v.normf) end function VectorInterface.scale(v::InnerProductVec, a::Number) - return InnerProductVec(scale(v.vec, a), v.dotf) + return InnerProductVec(scale(v.vec, a), v.dotf, v.normf) end function VectorInterface.scale!!(v::InnerProductVec, a::Number) - return InnerProductVec(scale!!(v.vec, a), v.dotf) + return InnerProductVec(scale!!(v.vec, a), v.dotf, v.normf) end function VectorInterface.scale!(v::InnerProductVec, a::Number) scale!(v.vec, a) return v end function VectorInterface.scale!!( - w::InnerProductVec{F}, v::InnerProductVec{F}, a::Number - ) where {F} - return InnerProductVec(scale!!(w.vec, v.vec, a), w.dotf) + w::InnerProductVec{F, N}, v::InnerProductVec{F, N}, a::Number + ) where {F, N} + return InnerProductVec(scale!!(w.vec, v.vec, a), w.dotf, w.normf) end function VectorInterface.scale!( - w::InnerProductVec{F}, v::InnerProductVec{F}, a::Number - ) where {F} + w::InnerProductVec{F, N}, v::InnerProductVec{F, N}, a::Number + ) where {F, N} scale!(w.vec, v.vec, a) return w end function VectorInterface.add( - v::InnerProductVec{F}, w::InnerProductVec{F}, a::Number, b::Number - ) where {F} - return InnerProductVec(add(v.vec, w.vec, a, b), v.dotf) + v::InnerProductVec{F, N}, w::InnerProductVec{F, N}, a::Number, b::Number + ) where {F, N} + return InnerProductVec(add(v.vec, w.vec, a, b), v.dotf, v.normf) end function VectorInterface.add!!( - v::InnerProductVec{F}, w::InnerProductVec{F}, a::Number, b::Number - ) where {F} - return InnerProductVec(add!!(v.vec, w.vec, a, b), v.dotf) + v::InnerProductVec{F, N}, w::InnerProductVec{F, N}, a::Number, b::Number + ) where {F, N} + return InnerProductVec(add!!(v.vec, w.vec, a, b), v.dotf, v.normf) end function VectorInterface.add!( - v::InnerProductVec{F}, w::InnerProductVec{F}, a::Number, b::Number - ) where {F} + v::InnerProductVec{F, N}, w::InnerProductVec{F, N}, a::Number, b::Number + ) where {F, N} add!(v.vec, w.vec, a, b) return v end -function VectorInterface.inner(v::InnerProductVec{F}, w::InnerProductVec{F}) where {F} +function VectorInterface.inner(v::InnerProductVec{F, N}, w::InnerProductVec{F, N}) where {F, N} return v.dotf(v.vec, w.vec) end -VectorInterface.norm(v::InnerProductVec) = sqrt(real(inner(v, v))) +VectorInterface.norm(v::InnerProductVec{F, Nothing}) where {F} = sqrt(real(inner(v, v))) +VectorInterface.norm(v::InnerProductVec) = v.normf(v.vec) diff --git a/src/orthonormal.jl b/src/orthonormal.jl index 9089b8f3..0eec9474 100644 --- a/src/orthonormal.jl +++ b/src/orthonormal.jl @@ -68,15 +68,15 @@ _use_multithreaded_array_kernel(::Type) = false function _use_multithreaded_array_kernel(::Type{<:Array{T}}) where {T <: Number} return isbitstype(T) && get_num_threads() > 1 end -function _use_multithreaded_array_kernel(::Type{<:OrthonormalBasis{T}}) where {T} +function _use_multithreaded_array_kernel(::Type{<:Basis{T}}) where {T} return _use_multithreaded_array_kernel(T) end """ - project!!(y::AbstractVector, b::OrthonormalBasis, x, + project!!(y::AbstractVector, b::Basis, x, [α::Number = 1, β::Number = 0, r = Base.OneTo(length(b))]) -For a given orthonormal basis `b`, compute the expansion coefficients `y` resulting from +For a given basis `b`, compute the expansion coefficients `y` resulting from projecting the vector `x` onto the subspace spanned by `b`; more specifically this computes ``` @@ -86,7 +86,7 @@ projecting the vector `x` onto the subspace spanned by `b`; more specifically th for all ``j ∈ r``. """ function project!!( - y::AbstractVector, b::OrthonormalBasis, x, + y::AbstractVector, b::Basis, x, α::Number = true, β::Number = false, r = Base.OneTo(length(b)) ) # no specialized routine for IndexLinear x because reduction dimension is large dimension @@ -118,10 +118,10 @@ function project!!( end """ - unproject!!(y, b::OrthonormalBasis, x::AbstractVector, + unproject!!(y, b::Basis, x::AbstractVector, [α::Number = 1, β::Number = 0, r = Base.OneTo(length(b))]) -For a given orthonormal basis `b`, reconstruct the vector-like object `y` that is defined by +For a given basis `b`, reconstruct the vector-like object `y` that is defined by expansion coefficients with respect to the basis vectors in `b` in `x`; more specifically this computes @@ -130,7 +130,7 @@ this computes ``` """ function unproject!!( - y, b::OrthonormalBasis, x::AbstractVector, + y, b::Basis, x::AbstractVector, α::Number = true, β::Number = false, r = Base.OneTo(length(b)) ) if _use_multithreaded_array_kernel(y) @@ -149,7 +149,7 @@ function unproject!!( return y end function unproject_linear_multithreaded!( - y::AbstractArray, b::OrthonormalBasis{<:AbstractArray}, x::AbstractVector, + y::AbstractArray, b::Basis{<:AbstractArray}, x::AbstractVector, α::Number = true, β::Number = false, r = Base.OneTo(length(b)) ) # multi-threaded implementation, similar to BLAS level 2 matrix vector multiplication @@ -172,7 +172,7 @@ function unproject_linear_multithreaded!( return y end function unproject_linear_kernel!( - y::AbstractArray, b::OrthonormalBasis{<:AbstractArray}, x::AbstractVector, + y::AbstractArray, b::Basis{<:AbstractArray}, x::AbstractVector, I, α::Number, β::Number, r ) return @inbounds begin @@ -554,3 +554,251 @@ and its concrete subtypes [`ClassicalGramSchmidt`](@ref), [`ModifiedGramSchmidt` [`ClassicalGramSchmidtIR`](@ref) and [`ModifiedGramSchmidtIR`](@ref). """ orthonormalize, orthonormalize!! + +# Definition of a symplectic (Darboux) basis +""" + SymplecticBasis{T} <: Basis{T} + +A list of vector like objects of type `T` that form a symplectic (Darboux) basis with +respect to a skew-symmetric bilinear form `ω`. A symplectic basis satisfies the relations +`ω(u_{2m-1}, u_{2n}) = δ_{mn}`, `ω(u_{2m}, u_{2n}) = 0`, and `ω(u_{2m-1}, u_{2n-1}) = 0`. +See also [`Basis`](@ref). + +Skew-orthonormality of the vectors contained in an instance `b` of `SymplecticBasis` is not +checked when elements are added; it is up to the algorithm that constructs `b` to guarantee +skew-orthonormality. + +Vectors are added in pairs: odd-indexed vectors `u_{2m-1}` and their symplectic partners +`u_{2m}`. The function [`skeworthogonalize`](@ref) or [`skeworthonormalize`](@ref) can be +used to skew-orthogonalize a new vector with respect to the existing basis. These functions +require a skew-symmetric form `ω` to be provided. Note that in-place versions +[`skeworthogonalize!!`](@ref) or [`skeworthonormalize!!`](@ref) are also available. +""" +struct SymplecticBasis{T} <: Basis{T} + basis::Vector{T} +end +SymplecticBasis{T}() where {T} = SymplecticBasis{T}(Vector{T}(undef, 0)) + +# Iterator methods for SymplecticBasis +Base.IteratorSize(::Type{<:SymplecticBasis}) = Base.HasLength() +Base.IteratorEltype(::Type{<:SymplecticBasis}) = Base.HasEltype() + +Base.length(b::SymplecticBasis) = length(b.basis) +Base.eltype(b::SymplecticBasis{T}) where {T} = T + +Base.iterate(b::SymplecticBasis) = Base.iterate(b.basis) +Base.iterate(b::SymplecticBasis, state) = Base.iterate(b.basis, state) + +Base.getindex(b::SymplecticBasis, i) = getindex(b.basis, i) +Base.setindex!(b::SymplecticBasis, i, q) = setindex!(b.basis, i, q) +Base.firstindex(b::SymplecticBasis) = firstindex(b.basis) +Base.lastindex(b::SymplecticBasis) = lastindex(b.basis) + +Base.first(b::SymplecticBasis) = first(b.basis) +Base.last(b::SymplecticBasis) = last(b.basis) + +Base.popfirst!(b::SymplecticBasis) = popfirst!(b.basis) +Base.pop!(b::SymplecticBasis) = pop!(b.basis) +Base.push!(b::SymplecticBasis{T}, q::T) where {T} = (push!(b.basis, q); return b) +Base.empty!(b::SymplecticBasis) = (empty!(b.basis); return b) +Base.sizehint!(b::SymplecticBasis, k::Int) = (sizehint!(b.basis, k); return b) +Base.resize!(b::SymplecticBasis, k::Int) = (resize!(b.basis, k); return b) + +""" + numpairs(b::SymplecticBasis) -> Int + +Return the number of symplectic pairs in the basis `b`. This equals `div(length(b), 2)`. +""" +numpairs(b::SymplecticBasis) = div(length(b), 2) + +# Skew-orthogonalization of a vector against a given SymplecticBasis +skeworthogonalize(v, args...) = skeworthogonalize!!(scale(v, true), args...) + +function skeworthogonalize!!( + v::T, b::SymplecticBasis{T}, alg::SkewOrthogonalizer + ) where {T} + S = promote_type(scalartype(v), scalartype(T)) + c = Vector{S}(undef, length(b)) + return skeworthogonalize!!(v, b, c, alg) +end + +# See pages 3-4 of https://people.math.ethz.ch/%7Eacannas/Papers/lsg.pdf +function skeworthogonalize!!( + v::T, b::SymplecticBasis{T}, x::AbstractVector, ::ClassicalSymplecticGramSchmidt + ) where {T} + np = numpairs(b) + idx_odd = 1:2:(2np - 1) + idx_even = 2:2:2np + project!!(view(x, idx_odd), b, v, -1, 0, idx_even) + project!!(view(x, idx_even), b, v, 1, 0, idx_odd) + v = unproject!!(v, b, x, -1, 1) + return (v, x) +end + +function reskeworthogonalize!!( + v::T, b::SymplecticBasis{T}, x::AbstractVector, ::ClassicalSymplecticGramSchmidt + ) where {T} + s = similar(x) ## EXTRA ALLOCATION + (v, x) = skeworthogonalize!!(v, b, s, ClassicalSymplecticGramSchmidt()) + x .+= s + return (v, x) +end + +function skeworthogonalize!!( + v::T, b::SymplecticBasis{T}, x::AbstractVector, ::ClassicalSymplecticGramSchmidt2 + ) where {T} + (v, x) = skeworthogonalize!!(v, b, x, ClassicalSymplecticGramSchmidt()) + return reskeworthogonalize!!(v, b, x, ClassicalSymplecticGramSchmidt()) +end + +function skeworthogonalize!!( + v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ClassicalSymplecticGramSchmidtIR + ) where {T} + nold = norm(v) + (v, x) = skeworthogonalize!!(v, b, x, ClassicalSymplecticGramSchmidt()) + nnew = norm(v) + while eps(one(nnew)) < nnew < alg.η * nold + nold = nnew + (v, x) = reskeworthogonalize!!(v, b, x, ClassicalSymplecticGramSchmidt()) + nnew = norm(v) + end + return (v, x) +end + +function skeworthogonalize!!( + v::T, b::SymplecticBasis{T}, x::AbstractVector, ::ModifiedSymplecticGramSchmidt + ) where {T} + np = numpairs(b) + for m in 1:np + i_odd = 2m - 1 + i_even = 2m + h_e = inner(b[i_odd], v) + h_f = inner(b[i_even], v) + x[i_odd] = -h_f + x[i_even] = h_e + v = add!!(v, b[i_odd], h_f) + v = add!!(v, b[i_even], -h_e) + end + return (v, x) +end + +function reskeworthogonalize!!( + v::T, b::SymplecticBasis{T}, x::AbstractVector, ::ModifiedSymplecticGramSchmidt + ) where {T} + np = numpairs(b) + for m in 1:np + i_odd = 2m - 1 + i_even = 2m + h_e = inner(b[i_odd], v) + h_f = inner(b[i_even], v) + s_odd = -h_f + s_even = h_e + x[i_odd] += s_odd + x[i_even] += s_even + v = add!!(v, b[i_odd], h_f) + v = add!!(v, b[i_even], -h_e) + end + return (v, x) +end + +function skeworthogonalize!!( + v::T, b::SymplecticBasis{T}, x::AbstractVector, ::ModifiedSymplecticGramSchmidt2 + ) where {T} + (v, x) = skeworthogonalize!!(v, b, x, ModifiedSymplecticGramSchmidt()) + return reskeworthogonalize!!(v, b, x, ModifiedSymplecticGramSchmidt()) +end + +function skeworthogonalize!!( + v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ModifiedSymplecticGramSchmidtIR + ) where {T} + nold = norm(v) + (v, x) = skeworthogonalize!!(v, b, x, ModifiedSymplecticGramSchmidt()) + nnew = norm(v) + while eps(one(nnew)) < nnew < alg.η * nold + nold = nnew + (v, x) = reskeworthogonalize!!(v, b, x, ModifiedSymplecticGramSchmidt()) + nnew = norm(v) + end + return (v, x) +end + +# Skew-orthonormalization: skew-orthogonalization and normalization +# For odd vectors: normalize with standard norm +# For even vectors: scale so that ω(partner, v) = 1 +skeworthonormalize(v, args...) = skeworthonormalize!!(scale(v, VectorInterface.One()), args...) + +function skeworthonormalize!!(v, b::SymplecticBasis, args...) + out = skeworthogonalize!!(v, b, args...) # out[1] === v + if iseven(length(b)) + # Adding odd vector: normalize with standard norm + β = norm(v) + v = scale!!(v, inv(β)) + else + # Adding even vector: scale so that ω(partner, v) = 1 + # The partner is the last vector in b (the odd vector of the current pair) + β = inner(last(b), v) + v = scale!!(v, inv(β)) + end + return tuple(v, β, Base.tail(out)...) +end + +""" + skeworthogonalize(v, b::SymplecticBasis, [x::AbstractVector,] alg::SkewOrthogonalizer) -> w, x + skeworthogonalize!!(v, b::SymplecticBasis, [x::AbstractVector,] alg::SkewOrthogonalizer) -> w, x + +Skew-orthogonalize vector `v` against all the vectors in the symplectic basis `b` using the +skew-orthogonalization algorithm `alg` of type [`SkewOrthogonalizer`](@ref), and return the +resulting vector `w` and the overlap coefficients `x` of `v` with the basis vectors in `b`. + +The skew-orthogonalization uses the skew-symmetric form `ω = inner`, which is expected to +satisfy `ω(u, v) = -ω(v, u)` when vectors are appropriately wrapped (e.g., using +`InnerProductVec`). For a symplectic basis with pairs `(u_{2m-1}, u_{2m})` where +`ω(u_{2m-1}, u_{2m}) = 1`, the skew-orthogonalization ensures: +- `ω(u_{2m-1}, w) = 0` for all `m` +- `ω(u_{2m}, w) = 0` for all `m` + +In case of `skeworthogonalize!!`, the vector `v` is mutated in place. In both functions, +storage for the overlap coefficients `x` can be provided as optional argument +`x::AbstractVector` with `length(x) >= length(b)`. + +Note that `w` is not normalized, see also [`skeworthonormalize`](@ref). + +For more information on possible skew-orthogonalization algorithms, see +[`SkewOrthogonalizer`](@ref) and its concrete subtypes +[`ClassicalSymplecticGramSchmidt`](@ref), [`ModifiedSymplecticGramSchmidt`](@ref), +[`ClassicalSymplecticGramSchmidt2`](@ref), [`ModifiedSymplecticGramSchmidt2`](@ref), +[`ClassicalSymplecticGramSchmidtIR`](@ref) and [`ModifiedSymplecticGramSchmidtIR`](@ref). +""" +skeworthogonalize, skeworthogonalize!! + +""" + skeworthonormalize(v, b::SymplecticBasis, [x::AbstractVector,] alg::SkewOrthogonalizer) -> w, β, x + skeworthonormalize!!(v, b::SymplecticBasis, [x::AbstractVector,] alg::SkewOrthogonalizer) -> w, β, x + +Skew-orthonormalize vector `v` against all the vectors in the symplectic basis `b` using +the skew-orthogonalization algorithm `alg` of type [`SkewOrthogonalizer`](@ref). + +The normalization depends on the current length of the basis: + +**When `length(b)` is even** (adding an odd vector): returns the resulting vector `w` +normalized to unit norm (`‖w‖ = 1`), the norm `β = ‖v‖` after skew-orthogonalizing, and +the overlap coefficients `x`. + +**When `length(b)` is odd** (adding an even vector): returns the resulting vector `w` +scaled such that `ω(last(b), w) = 1`, where `last(b)` is the odd partner of the pair. +Returns the scaling factor `β = ω(last(b), v)` after skew-orthogonalizing, and the overlap +coefficients `x`. + +In case of `skeworthonormalize!!`, the vector `v` is mutated in place. In both functions, +storage for the overlap coefficients `x` can be provided as optional argument +`x::AbstractVector` with `length(x) >= length(b)`. + +See [`skeworthogonalize`](@ref) if `w` does not need to be normalized. + +For more information on possible skew-orthogonalization algorithms, see +[`SkewOrthogonalizer`](@ref) and its concrete subtypes +[`ClassicalSymplecticGramSchmidt`](@ref), [`ModifiedSymplecticGramSchmidt`](@ref), +[`ClassicalSymplecticGramSchmidt2`](@ref), [`ModifiedSymplecticGramSchmidt2`](@ref), +[`ClassicalSymplecticGramSchmidtIR`](@ref) and [`ModifiedSymplecticGramSchmidtIR`](@ref). +""" +skeworthonormalize, skeworthonormalize!! diff --git a/test/runtests.jl b/test/runtests.jl index 8268ef79..b952a653 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -29,6 +29,9 @@ t = time() @testset "Krylov factorisations" verbose = true begin include("factorize.jl") end +@testset "Symplectic Arnoldi" verbose = true begin + include("symplectic_arnoldi.jl") +end @testset "Linear problems with linsolve" verbose = true begin include("linsolve.jl") end diff --git a/test/symplectic_arnoldi.jl b/test/symplectic_arnoldi.jl new file mode 100644 index 00000000..112e705d --- /dev/null +++ b/test/symplectic_arnoldi.jl @@ -0,0 +1,53 @@ +@testset "Symplectic Arnoldi" begin + ε(x, y) = sign(y - x) / 2 + + n = 8 + domain = 0:(n - 1) + w(x) = 1 + + skew_dot(u, v) = sum(ε(x, y) * w(x) * w(y) * u[x + 1] * v[y + 1] + for x in domain, y in domain) + + function run_symplectic_arnoldi(orth) + v₀ = InnerProductVec(fill(1 / sqrt(sum(w, domain)), n), skew_dot, norm) + itr = ArnoldiIterator( + u -> InnerProductVec(domain .* u[], u.dotf, u.normf), + v₀, + orth, + ) + fact = initialize(itr) + for _ in 1:(n - 1) + expand!(itr, fact) + end + return stack(getindex, basis(fact).basis) + end + + function max_symplectic_error(W) + max_err = 0.0 + for i in axes(W, 2), j in axes(W, 2) + val = skew_dot(W[:, i], W[:, j]) + if isodd(i) && j == i + 1 + max_err = max(max_err, abs(val - 1)) + elseif isodd(j) && i == j + 1 + max_err = max(max_err, abs(val + 1)) + else + max_err = max(max_err, abs(val)) + end + end + return max_err + end + + algs = ( + ClassicalSymplecticGramSchmidt(), + ModifiedSymplecticGramSchmidt(), + ClassicalSymplecticGramSchmidt2(), + ModifiedSymplecticGramSchmidt2(), + ClassicalSymplecticGramSchmidtIR(0.75), + ModifiedSymplecticGramSchmidtIR(0.75), + ) + + for alg in algs + W = run_symplectic_arnoldi(alg) + @test max_symplectic_error(W) < 1e-8 + end +end From 434743b8d0f69ef85f3b844fbb166cee618d2252 Mon Sep 17 00:00:00 2001 From: Simeon David Schaub Date: Fri, 6 Feb 2026 13:52:55 +0100 Subject: [PATCH 02/13] fix formatting --- test/symplectic_arnoldi.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/test/symplectic_arnoldi.jl b/test/symplectic_arnoldi.jl index 112e705d..b4505b41 100644 --- a/test/symplectic_arnoldi.jl +++ b/test/symplectic_arnoldi.jl @@ -5,8 +5,10 @@ domain = 0:(n - 1) w(x) = 1 - skew_dot(u, v) = sum(ε(x, y) * w(x) * w(y) * u[x + 1] * v[y + 1] - for x in domain, y in domain) + skew_dot(u, v) = sum( + ε(x, y) * w(x) * w(y) * u[x + 1] * v[y + 1] + for x in domain, y in domain + ) function run_symplectic_arnoldi(orth) v₀ = InnerProductVec(fill(1 / sqrt(sum(w, domain)), n), skew_dot, norm) @@ -48,6 +50,6 @@ for alg in algs W = run_symplectic_arnoldi(alg) - @test max_symplectic_error(W) < 1e-8 + @test max_symplectic_error(W) < 1.0e-8 end end From 05c0d7619aaa39fb8c946b99cb7c35373c12ceb2 Mon Sep 17 00:00:00 2001 From: Simeon David Schaub Date: Fri, 6 Feb 2026 18:47:05 +0100 Subject: [PATCH 03/13] fix typos --- src/factorizations/arnoldi.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/factorizations/arnoldi.jl b/src/factorizations/arnoldi.jl index 8456056d..d0988e48 100644 --- a/src/factorizations/arnoldi.jl +++ b/src/factorizations/arnoldi.jl @@ -159,7 +159,7 @@ function initialize(iter::ArnoldiIterator; verbosity::Int = KrylovDefaults.verbo r = add!!(r, v, -dα) β = norm(r) elseif iter.orth isa Union{ClassicalGramSchmidtIR, ModifiedGramSchmidtIR} - while eps(one(β)) < β < alg.η * βold + while eps(one(β)) < β < iter.orth.η * βold βold = β dα = inner(v, r) α += dα @@ -245,7 +245,7 @@ function shrink!(state::ArnoldiFactorization, k; verbosity::Int = KrylovDefaults r = pop!(V) resize!(H, (k * k + 3 * k) >> 1) state.k = k - if iter.orth isa Orthogonalizer + if V isa OrthonormalBasis β = normres(state) else β = iseven(k) ? normres(state) : norm(r) From 3f2939eb1e15547ae17a2d8bb21905df6cf7dd0e Mon Sep 17 00:00:00 2001 From: Simeon David Schaub Date: Fri, 6 Feb 2026 19:04:01 +0100 Subject: [PATCH 04/13] fix inference issue --- src/factorizations/biarnoldi.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/factorizations/biarnoldi.jl b/src/factorizations/biarnoldi.jl index ce259830..ac42ff66 100644 --- a/src/factorizations/biarnoldi.jl +++ b/src/factorizations/biarnoldi.jl @@ -1,6 +1,6 @@ -mutable struct BiArnoldiFactorization{T, S} <: KrylovFactorization{T, S} - VH::ArnoldiFactorization{T, S} - WK::ArnoldiFactorization{T, S} +mutable struct BiArnoldiFactorization{T, S, B} <: KrylovFactorization{T, S} + VH::ArnoldiFactorization{T, S, B} + WK::ArnoldiFactorization{T, S, B} end Base.length(F::BiArnoldiFactorization) = length(F.VH) From 39ea1221dd41adb04550e5553dd0f342b05247a9 Mon Sep 17 00:00:00 2001 From: Simeon David Schaub Date: Fri, 6 Feb 2026 19:13:54 +0100 Subject: [PATCH 05/13] remove outdated type restriction again --- src/factorizations/arnoldi.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/factorizations/arnoldi.jl b/src/factorizations/arnoldi.jl index d0988e48..7f7725ac 100644 --- a/src/factorizations/arnoldi.jl +++ b/src/factorizations/arnoldi.jl @@ -181,7 +181,7 @@ function initialize(iter::ArnoldiIterator; verbosity::Int = KrylovDefaults.verbo end function initialize!( - iter::ArnoldiIterator{<:Any, <:Any, <:Orthogonalizer}, state::ArnoldiFactorization; + iter::ArnoldiIterator, state::ArnoldiFactorization; verbosity::Int = KrylovDefaults.verbosity[] ) x₀ = iter.x₀ From 8b86171672cdc859a04d6336722c1965b6b3f767 Mon Sep 17 00:00:00 2001 From: Simeon David Schaub Date: Fri, 6 Feb 2026 22:05:44 +0100 Subject: [PATCH 06/13] fix docs build --- docs/src/man/implementation.md | 25 +++++++++++++++++++++++++ src/orthonormal.jl | 5 ----- 2 files changed, 25 insertions(+), 5 deletions(-) diff --git a/docs/src/man/implementation.md b/docs/src/man/implementation.md index c631bf22..99d4bf55 100644 --- a/docs/src/man/implementation.md +++ b/docs/src/man/implementation.md @@ -13,6 +13,11 @@ that case, the specific implementation `OrthonormalBasis{T}` can be used: KrylovKit.OrthonormalBasis ``` +For symplectic (Darboux) bases, the specific implementation `SymplecticBasis{T}` can be used: +```@docs +KrylovKit.SymplecticBasis +``` + We can orthogonalize or orthonormalize a given vector to another vector (assumed normalized) or to a given [`KrylovKit.OrthonormalBasis`](@ref) using ```@docs @@ -25,6 +30,26 @@ KrylovKit.orthogonalize!! KrylovKit.orthonormalize!! ``` +For symplectic bases, the corresponding skew-orthogonalization and skew-orthonormalization +interfaces are +```@docs +KrylovKit.skeworthogonalize +KrylovKit.skeworthonormalize +KrylovKit.skeworthogonalize!! +KrylovKit.skeworthonormalize!! +``` + +The available skew-orthogonalization algorithms are +```@docs +KrylovKit.SkewOrthogonalizer +KrylovKit.ClassicalSymplecticGramSchmidt +KrylovKit.ModifiedSymplecticGramSchmidt +KrylovKit.ClassicalSymplecticGramSchmidt2 +KrylovKit.ModifiedSymplecticGramSchmidt2 +KrylovKit.ClassicalSymplecticGramSchmidtIR +KrylovKit.ModifiedSymplecticGramSchmidtIR +``` + The expansion coefficients of a general vector in terms of a given orthonormal basis can be obtained as ```@docs KrylovKit.project!! diff --git a/src/orthonormal.jl b/src/orthonormal.jl index 0eec9474..dcb7790f 100644 --- a/src/orthonormal.jl +++ b/src/orthonormal.jl @@ -604,11 +604,6 @@ Base.empty!(b::SymplecticBasis) = (empty!(b.basis); return b) Base.sizehint!(b::SymplecticBasis, k::Int) = (sizehint!(b.basis, k); return b) Base.resize!(b::SymplecticBasis, k::Int) = (resize!(b.basis, k); return b) -""" - numpairs(b::SymplecticBasis) -> Int - -Return the number of symplectic pairs in the basis `b`. This equals `div(length(b), 2)`. -""" numpairs(b::SymplecticBasis) = div(length(b), 2) # Skew-orthogonalization of a vector against a given SymplecticBasis From 8525cbdc3974b01cfa15d2c285a8384a31b37e16 Mon Sep 17 00:00:00 2001 From: Simeon David Schaub Date: Sat, 7 Feb 2026 16:52:46 +0100 Subject: [PATCH 07/13] implement more stable ESR2, make it the default --- docs/Project.toml | 1 + src/KrylovKit.jl | 1 + src/algorithms.jl | 69 ++++++++++++++++++++++++++--------- src/factorizations/arnoldi.jl | 21 ++++++++--- src/innerproductvec.jl | 3 ++ src/orthonormal.jl | 52 +++++++++++++++++++------- test/symplectic_arnoldi.jl | 49 +++++++++++++++++-------- 7 files changed, 144 insertions(+), 52 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index 1814eb33..5a4b3ba2 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,5 +1,6 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" [compat] Documenter = "1" diff --git a/src/KrylovKit.jl b/src/KrylovKit.jl index 74596171..2e0d3875 100644 --- a/src/KrylovKit.jl +++ b/src/KrylovKit.jl @@ -37,6 +37,7 @@ export basis, rayleighquotient, residual, normres, rayleighextension export initialize, initialize!, expand!, shrink! export ClassicalGramSchmidt, ClassicalGramSchmidt2, ClassicalGramSchmidtIR export ModifiedGramSchmidt, ModifiedGramSchmidt2, ModifiedGramSchmidtIR +export ESR, ESR1, ESR2, ESR3 export ClassicalSymplecticGramSchmidt, ClassicalSymplecticGramSchmidt2 export ClassicalSymplecticGramSchmidtIR export ModifiedSymplecticGramSchmidt, ModifiedSymplecticGramSchmidt2 diff --git a/src/algorithms.jl b/src/algorithms.jl index d22bfb1c..cd45b256 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -80,6 +80,19 @@ end ModifiedGramSchmidtIR() = ModifiedGramSchmidtIR(1 / sqrt(2)) # Daniel-Gragg-Kaufman-Stewart # Skew-orthogonalization for symplectic bases +# TODO: implement ESR3 +""" + @enum ESR ESR1 ESR2 ESR3 + +Enum for selecting Elementary Symplectic factorization (ESR) variant: +- `ESR1`: r11 = ||x₁||, r12 = 0 +- `ESR2`: r11 = ||x₁||, r12 = s₁ᵀx₂ (most stable, s₁ and s₂ orthogonal) +- `ESR3`: r11 = ||⟨x₁, x₂⟩||, r12 = 0 (TODO) + +See https://journal.austms.org.au/ojs/index.php/ANZIAMJ/article/view/9380/1920, page 2 for details. +""" +@enum ESR ESR1 ESR2# ESR3 + """ abstract type SkewOrthogonalizer @@ -96,64 +109,84 @@ abstract type SkewReorthogonalizer <: SkewOrthogonalizer end # Simple """ - ClassicalSymplecticGramSchmidt() + ClassicalSymplecticGramSchmidt(esr::ESR = ESR2) Represents the classical symplectic Gram Schmidt algorithm for skew-orthogonalizing vectors -to produce a symplectic (Darboux) basis, typically not an optimal choice. +to produce a symplectic (Darboux) basis, typically not an optimal choice. The `esr` parameter +selects the Elementary Symplectic factorization variant (default: ESR2, most stable). """ -struct ClassicalSymplecticGramSchmidt <: SkewOrthogonalizer end +struct ClassicalSymplecticGramSchmidt <: SkewOrthogonalizer + esr::ESR +end +ClassicalSymplecticGramSchmidt() = ClassicalSymplecticGramSchmidt(ESR2) """ - ModifiedSymplecticGramSchmidt() + ModifiedSymplecticGramSchmidt(esr::ESR = ESR2) Represents the modified symplectic Gram Schmidt algorithm for skew-orthogonalizing vectors -to produce a symplectic (Darboux) basis. +to produce a symplectic (Darboux) basis. The `esr` parameter selects the Elementary +Symplectic factorization variant (default: ESR2, most stable). """ -struct ModifiedSymplecticGramSchmidt <: SkewOrthogonalizer end +struct ModifiedSymplecticGramSchmidt <: SkewOrthogonalizer + esr::ESR +end +ModifiedSymplecticGramSchmidt() = ModifiedSymplecticGramSchmidt(ESR2) # A single reorthogonalization always """ - ClassicalSymplecticGramSchmidt2() + ClassicalSymplecticGramSchmidt2(esr::ESR = ESR2) Represents the classical symplectic Gram Schmidt algorithm with a second reskew- -orthogonalization step always taking place. +orthogonalization step always taking place. The `esr` parameter selects the Elementary +Symplectic factorization variant (default: ESR2, most stable). """ -struct ClassicalSymplecticGramSchmidt2 <: SkewReorthogonalizer end +struct ClassicalSymplecticGramSchmidt2 <: SkewReorthogonalizer + esr::ESR +end +ClassicalSymplecticGramSchmidt2() = ClassicalSymplecticGramSchmidt2(ESR2) """ - ModifiedSymplecticGramSchmidt2() + ModifiedSymplecticGramSchmidt2(esr::ESR = ESR2) Represents the modified symplectic Gram Schmidt algorithm with a second reskew- -orthogonalization step always taking place. +orthogonalization step always taking place. The `esr` parameter selects the Elementary +Symplectic factorization variant (default: ESR2, most stable). """ -struct ModifiedSymplecticGramSchmidt2 <: SkewReorthogonalizer end +struct ModifiedSymplecticGramSchmidt2 <: SkewReorthogonalizer + esr::ESR +end +ModifiedSymplecticGramSchmidt2() = ModifiedSymplecticGramSchmidt2(ESR2) # Iterative reorthogonalization """ - ClassicalSymplecticGramSchmidtIR(η::Real = 1/sqrt(2)) + ClassicalSymplecticGramSchmidtIR(η::Real = 1/sqrt(2), esr::ESR = ESR2) Represents the classical symplectic Gram Schmidt algorithm with iterative (i.e. zero or more) reskew-orthogonalization until the norm of the vector after a skew-orthogonalization step has not decreased by a factor smaller than `η` with respect to the norm before the -step. The default value corresponds to the Daniel-Gragg-Kaufman-Stewart condition. +step. The default value corresponds to the Daniel-Gragg-Kaufman-Stewart condition. The `esr` +parameter selects the Elementary Symplectic factorization variant (default: ESR2, most stable). """ struct ClassicalSymplecticGramSchmidtIR{S <: Real} <: SkewReorthogonalizer η::S + esr::ESR end -ClassicalSymplecticGramSchmidtIR() = ClassicalSymplecticGramSchmidtIR(1 / sqrt(2)) +ClassicalSymplecticGramSchmidtIR(η::Real = 1 / sqrt(2)) = ClassicalSymplecticGramSchmidtIR(η, ESR2) """ - ModifiedSymplecticGramSchmidtIR(η::Real = 1/sqrt(2)) + ModifiedSymplecticGramSchmidtIR(η::Real = 1/sqrt(2), esr::ESR = ESR2) Represents the modified symplectic Gram Schmidt algorithm with iterative (i.e. zero or more) reskew-orthogonalization until the norm of the vector after a skew-orthogonalization step has not decreased by a factor smaller than `η` with respect to the norm before the -step. The default value corresponds to the Daniel-Gragg-Kaufman-Stewart condition. +step. The default value corresponds to the Daniel-Gragg-Kaufman-Stewart condition. The `esr` +parameter selects the Elementary Symplectic factorization variant (default: ESR2, most stable). """ struct ModifiedSymplecticGramSchmidtIR{S <: Real} <: SkewReorthogonalizer η::S + esr::ESR end -ModifiedSymplecticGramSchmidtIR() = ModifiedSymplecticGramSchmidtIR(1 / sqrt(2)) +ModifiedSymplecticGramSchmidtIR(η::Real = 1 / sqrt(2)) = ModifiedSymplecticGramSchmidtIR(η, ESR2) # Solving eigenvalue problems abstract type KrylovAlgorithm end diff --git a/src/factorizations/arnoldi.jl b/src/factorizations/arnoldi.jl index 7f7725ac..1d1870c2 100644 --- a/src/factorizations/arnoldi.jl +++ b/src/factorizations/arnoldi.jl @@ -172,7 +172,13 @@ function initialize(iter::ArnoldiIterator; verbosity::Int = KrylovDefaults.verbo else iszero(α) && throw(ArgumentError("initial vector and its image are symplectically orthogonal")) V = SymplecticBasis([v]) - H = T[zero(α), α] + if iter.orth.esr == ESR2 + r12 = standard_dot(v, r) + r = add!!(r, v, -r12) + H = T[r12, α] + else + H = T[zero(α), α] + end end if verbosity > EACHITERATION_LEVEL @info "Arnoldi initiation at dimension 1: subspace normres = $(normres2string(β))" @@ -201,7 +207,13 @@ function initialize!( α = inner(V[1], w) iszero(α) && throw(ArgumentError("initial vector and its image are symplectically orthogonal")) r = w - push!(H, zero(α), α) + if iter.orth.esr == ESR2 + r12 = standard_dot(V[1], r) + r = add!!(r, V[1], -r12) + push!(H, r12, α) + else + push!(H, zero(α), α) + end end state.k = 1 state.r = r @@ -253,13 +265,12 @@ function shrink!(state::ArnoldiFactorization, k; verbosity::Int = KrylovDefaults if verbosity > EACHITERATION_LEVEL @info "Arnoldi reduction to dimension $k: subspace normres = $(normres2string(β))" end - state.r = scale!!(r, β) return state end # Arnoldi recurrence: simply use provided orthonormalization routines function arnoldirecurrence!!( - operator, V::OrthonormalBasis, h::AbstractVector, orth::Orthogonalizer + operator, V::OrthonormalBasis, h::AbstractVector, orth::Orthogonalizer, ) w = apply(operator, last(V)) r, h = orthogonalize!!(w, V, h, orth) @@ -267,7 +278,7 @@ function arnoldirecurrence!!( end function arnoldirecurrence!!( - operator, V::SymplecticBasis, h::AbstractVector, orth::SkewOrthogonalizer + operator, V::SymplecticBasis, h::AbstractVector, orth::SkewOrthogonalizer, ) w = apply(operator, last(V)) r, h = skeworthogonalize!!(w, V, h, orth) diff --git a/src/innerproductvec.jl b/src/innerproductvec.jl index dea411c8..8e0a4176 100644 --- a/src/innerproductvec.jl +++ b/src/innerproductvec.jl @@ -144,3 +144,6 @@ end VectorInterface.norm(v::InnerProductVec{F, Nothing}) where {F} = sqrt(real(inner(v, v))) VectorInterface.norm(v::InnerProductVec) = v.normf(v.vec) + +standard_dot(v::InnerProductVec, w::InnerProductVec) = dot(v.vec, w.vec) +standard_dot(v, w) = dot(v, w) diff --git a/src/orthonormal.jl b/src/orthonormal.jl index dcb7790f..9504e804 100644 --- a/src/orthonormal.jl +++ b/src/orthonormal.jl @@ -619,19 +619,26 @@ end # See pages 3-4 of https://people.math.ethz.ch/%7Eacannas/Papers/lsg.pdf function skeworthogonalize!!( - v::T, b::SymplecticBasis{T}, x::AbstractVector, ::ClassicalSymplecticGramSchmidt + v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ClassicalSymplecticGramSchmidt ) where {T} np = numpairs(b) idx_odd = 1:2:(2np - 1) idx_even = 2:2:2np project!!(view(x, idx_odd), b, v, -1, 0, idx_even) project!!(view(x, idx_even), b, v, 1, 0, idx_odd) + if isodd(length(b)) + if alg.esr == ESR2 + x[2np + 1] = standard_dot(last(b), v) + else + x[2np + 1] = zero(eltype(x)) + end + end v = unproject!!(v, b, x, -1, 1) return (v, x) end function reskeworthogonalize!!( - v::T, b::SymplecticBasis{T}, x::AbstractVector, ::ClassicalSymplecticGramSchmidt + v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ClassicalSymplecticGramSchmidt ) where {T} s = similar(x) ## EXTRA ALLOCATION (v, x) = skeworthogonalize!!(v, b, s, ClassicalSymplecticGramSchmidt()) @@ -640,28 +647,30 @@ function reskeworthogonalize!!( end function skeworthogonalize!!( - v::T, b::SymplecticBasis{T}, x::AbstractVector, ::ClassicalSymplecticGramSchmidt2 + v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ClassicalSymplecticGramSchmidt2 ) where {T} - (v, x) = skeworthogonalize!!(v, b, x, ClassicalSymplecticGramSchmidt()) - return reskeworthogonalize!!(v, b, x, ClassicalSymplecticGramSchmidt()) + csgs = ClassicalSymplecticGramSchmidt(alg.esr) + (v, x) = skeworthogonalize!!(v, b, x, csgs) + return reskeworthogonalize!!(v, b, x, csgs) end function skeworthogonalize!!( v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ClassicalSymplecticGramSchmidtIR ) where {T} + csgs = ClassicalSymplecticGramSchmidt(alg.esr) nold = norm(v) - (v, x) = skeworthogonalize!!(v, b, x, ClassicalSymplecticGramSchmidt()) + (v, x) = skeworthogonalize!!(v, b, x, csgs) nnew = norm(v) while eps(one(nnew)) < nnew < alg.η * nold nold = nnew - (v, x) = reskeworthogonalize!!(v, b, x, ClassicalSymplecticGramSchmidt()) + (v, x) = reskeworthogonalize!!(v, b, x, csgs) nnew = norm(v) end return (v, x) end function skeworthogonalize!!( - v::T, b::SymplecticBasis{T}, x::AbstractVector, ::ModifiedSymplecticGramSchmidt + v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ModifiedSymplecticGramSchmidt ) where {T} np = numpairs(b) for m in 1:np @@ -674,11 +683,19 @@ function skeworthogonalize!!( v = add!!(v, b[i_odd], h_f) v = add!!(v, b[i_even], -h_e) end + if isodd(length(b)) + if alg.esr == ESR2 + x[2np + 1] = standard_dot(last(b), v) + v = add!!(v, last(b), -x[2np + 1]) + else + x[2np + 1] = zero(eltype(x)) + end + end return (v, x) end function reskeworthogonalize!!( - v::T, b::SymplecticBasis{T}, x::AbstractVector, ::ModifiedSymplecticGramSchmidt + v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ModifiedSymplecticGramSchmidt ) where {T} np = numpairs(b) for m in 1:np @@ -693,25 +710,32 @@ function reskeworthogonalize!!( v = add!!(v, b[i_odd], h_f) v = add!!(v, b[i_even], -h_e) end + if isodd(length(b)) && alg.esr == ESR2 + r11 = standard_dot(last(b), v) + x[2np + 1] += r11 + v = add!!(v, last(b), -r11) + end return (v, x) end function skeworthogonalize!!( - v::T, b::SymplecticBasis{T}, x::AbstractVector, ::ModifiedSymplecticGramSchmidt2 + v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ModifiedSymplecticGramSchmidt2 ) where {T} - (v, x) = skeworthogonalize!!(v, b, x, ModifiedSymplecticGramSchmidt()) - return reskeworthogonalize!!(v, b, x, ModifiedSymplecticGramSchmidt()) + msgs = ModifiedSymplecticGramSchmidt(alg.esr) + (v, x) = skeworthogonalize!!(v, b, x, msgs) + return reskeworthogonalize!!(v, b, x, msgs) end function skeworthogonalize!!( v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ModifiedSymplecticGramSchmidtIR ) where {T} + msgs = ModifiedSymplecticGramSchmidt(alg.esr) nold = norm(v) - (v, x) = skeworthogonalize!!(v, b, x, ModifiedSymplecticGramSchmidt()) + (v, x) = skeworthogonalize!!(v, b, x, msgs) nnew = norm(v) while eps(one(nnew)) < nnew < alg.η * nold nold = nnew - (v, x) = reskeworthogonalize!!(v, b, x, ModifiedSymplecticGramSchmidt()) + (v, x) = reskeworthogonalize!!(v, b, x, msgs) nnew = norm(v) end return (v, x) diff --git a/test/symplectic_arnoldi.jl b/test/symplectic_arnoldi.jl index b4505b41..5c9ba781 100644 --- a/test/symplectic_arnoldi.jl +++ b/test/symplectic_arnoldi.jl @@ -1,7 +1,6 @@ @testset "Symplectic Arnoldi" begin ε(x, y) = sign(y - x) / 2 - n = 8 domain = 0:(n - 1) w(x) = 1 @@ -10,7 +9,7 @@ for x in domain, y in domain ) - function run_symplectic_arnoldi(orth) + function run_symplectic_arnoldi(orth, inplace_init) v₀ = InnerProductVec(fill(1 / sqrt(sum(w, domain)), n), skew_dot, norm) itr = ArnoldiIterator( u -> InnerProductVec(domain .* u[], u.dotf, u.normf), @@ -18,10 +17,13 @@ orth, ) fact = initialize(itr) - for _ in 1:(n - 1) + if inplace_init + fact = initialize!(itr, fact) + end + for i in 1:(n - 1) expand!(itr, fact) end - return stack(getindex, basis(fact).basis) + return stack(getindex, basis(fact).basis), rayleighquotient(fact), residual(fact), rayleighextension(fact) end function max_symplectic_error(W) @@ -39,17 +41,34 @@ return max_err end - algs = ( - ClassicalSymplecticGramSchmidt(), - ModifiedSymplecticGramSchmidt(), - ClassicalSymplecticGramSchmidt2(), - ModifiedSymplecticGramSchmidt2(), - ClassicalSymplecticGramSchmidtIR(0.75), - ModifiedSymplecticGramSchmidtIR(0.75), - ) + for esr in (ESR1, ESR2) + algs = ( + ClassicalSymplecticGramSchmidt(esr), + ModifiedSymplecticGramSchmidt(esr), + ClassicalSymplecticGramSchmidt2(esr), + ModifiedSymplecticGramSchmidt2(esr), + ClassicalSymplecticGramSchmidtIR(0.75, esr), + ModifiedSymplecticGramSchmidtIR(0.75, esr), + ) + @testset "$alg" for alg in algs + W1, H1, r1, b1 = run_symplectic_arnoldi(alg, false) + @test max_symplectic_error(W1) < 1.0e-8 + if alg == ClassicalSymplecticGramSchmidt2(ESR1) + @test_broken Diagonal(domain) * W1 ≈ W1 * H1 + r1[] * b1' atol = 1.0e-8 + else + @test Diagonal(domain) * W1 ≈ W1 * H1 + r1[] * b1' atol = 1.0e-8 + end - for alg in algs - W = run_symplectic_arnoldi(alg) - @test max_symplectic_error(W) < 1.0e-8 + W2, H2, r2, b2 = run_symplectic_arnoldi(alg, true) + @test W1 ≈ W2 + @test H1 ≈ H2 + @test r1[] ≈ r2[] atol = 1.0e-8 + @test max_symplectic_error(W2) < 1.0e-8 + if alg == ClassicalSymplecticGramSchmidt2(ESR1) + @test_broken Diagonal(domain) * W2 ≈ W2 * H2 + r2[] * b2' atol = 1.0e-8 + else + @test Diagonal(domain) * W2 ≈ W2 * H2 + r2[] * b2' atol = 1.0e-8 + end + end end end From ff9ae2ec8ce830e83e9c3188ed4940a685bbff75 Mon Sep 17 00:00:00 2001 From: Simeon David Schaub Date: Sat, 7 Feb 2026 16:52:46 +0100 Subject: [PATCH 08/13] introduce ESR3m, fix `reskeworthogonalize!!` bug --- src/KrylovKit.jl | 2 +- src/algorithms.jl | 5 ++--- src/factorizations/arnoldi.jl | 19 ++++++++++++++----- src/orthonormal.jl | 8 +++----- test/symplectic_arnoldi.jl | 14 +++----------- 5 files changed, 23 insertions(+), 25 deletions(-) diff --git a/src/KrylovKit.jl b/src/KrylovKit.jl index 2e0d3875..be11da31 100644 --- a/src/KrylovKit.jl +++ b/src/KrylovKit.jl @@ -37,7 +37,7 @@ export basis, rayleighquotient, residual, normres, rayleighextension export initialize, initialize!, expand!, shrink! export ClassicalGramSchmidt, ClassicalGramSchmidt2, ClassicalGramSchmidtIR export ModifiedGramSchmidt, ModifiedGramSchmidt2, ModifiedGramSchmidtIR -export ESR, ESR1, ESR2, ESR3 +export ESR, ESR1, ESR2, ESR3m export ClassicalSymplecticGramSchmidt, ClassicalSymplecticGramSchmidt2 export ClassicalSymplecticGramSchmidtIR export ModifiedSymplecticGramSchmidt, ModifiedSymplecticGramSchmidt2 diff --git a/src/algorithms.jl b/src/algorithms.jl index cd45b256..d9f918c9 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -80,18 +80,17 @@ end ModifiedGramSchmidtIR() = ModifiedGramSchmidtIR(1 / sqrt(2)) # Daniel-Gragg-Kaufman-Stewart # Skew-orthogonalization for symplectic bases -# TODO: implement ESR3 """ @enum ESR ESR1 ESR2 ESR3 Enum for selecting Elementary Symplectic factorization (ESR) variant: - `ESR1`: r11 = ||x₁||, r12 = 0 - `ESR2`: r11 = ||x₁||, r12 = s₁ᵀx₂ (most stable, s₁ and s₂ orthogonal) -- `ESR3`: r11 = ||⟨x₁, x₂⟩||, r12 = 0 (TODO) +- `ESR3m`: r11 = 1, r12 = 0, r22 = ⟨x₁, x₂⟩ (original ESR3: r11 = ||⟨x₁, x₂⟩||, r12 = 0, r22 = ±1) See https://journal.austms.org.au/ojs/index.php/ANZIAMJ/article/view/9380/1920, page 2 for details. """ -@enum ESR ESR1 ESR2# ESR3 +@enum ESR ESR1 ESR2 ESR3m """ abstract type SkewOrthogonalizer diff --git a/src/factorizations/arnoldi.jl b/src/factorizations/arnoldi.jl index 1d1870c2..1b4da90a 100644 --- a/src/factorizations/arnoldi.jl +++ b/src/factorizations/arnoldi.jl @@ -47,7 +47,7 @@ Base.eltype(::Type{<:ArnoldiFactorization{<:Any, S}}) where {S} = S basis(F::ArnoldiFactorization) = F.V rayleighquotient(F::ArnoldiFactorization) = PackedHessenberg(F.H, F.k) residual(F::ArnoldiFactorization) = F.r -@inbounds normres(F::ArnoldiFactorization) = abs(F.H[end]) +@inbounds normres(F::ArnoldiFactorization) = F.H[end] rayleighextension(F::ArnoldiFactorization) = SimpleBasisVector(F.k, F.k) # Arnoldi iteration for constructing the orthonormal basis of a Krylov subspace. @@ -135,7 +135,11 @@ end function initialize(iter::ArnoldiIterator; verbosity::Int = KrylovDefaults.verbosity[]) # initialize without using eltype x₀ = iter.x₀ - β₀ = norm(x₀) + if iter.orth isa Orthogonalizer || iter.orth.esr != ESR3 + β₀ = norm(x₀) + else + β₀ = one(scalartype(x₀)) + end iszero(β₀) && throw(ArgumentError("initial vector should not have norm zero")) Ax₀ = apply(iter.operator, x₀) α = inner(x₀, Ax₀) / (β₀ * β₀) @@ -197,7 +201,12 @@ function initialize!( end H = empty!(state.H) - V[1] = scale!!(V[1], x₀, 1 / norm(x₀)) + if iter.orth isa Orthogonalizer || iter.orth.esr != ESR3 + β₀ = norm(x₀) + else + β₀ = one(scalartype(x₀)) + end + V[1] = scale!!(V[1], x₀, 1 / β₀) w = apply(iter.operator, V[1]) if iter.orth isa Orthogonalizer r, α = orthogonalize!!(w, V[1], iter.orth) @@ -231,7 +240,7 @@ function expand!( V = state.V H = state.H r = state.r - if iter.orth isa Orthogonalizer + if iter.orth isa Orthogonalizer || iter.orth.esr == ESR3 β = normres(state) else β = iseven(k) ? normres(state) : norm(r) @@ -282,5 +291,5 @@ function arnoldirecurrence!!( ) w = apply(operator, last(V)) r, h = skeworthogonalize!!(w, V, h, orth) - return r, iseven(length(V)) ? norm(r) : inner(last(V), r) + return r, iseven(length(V)) ? (orth.esr == ESR3 ? one(scalartype(r)) : norm(r)) : inner(last(V), r) end diff --git a/src/orthonormal.jl b/src/orthonormal.jl index 9504e804..b986ebaa 100644 --- a/src/orthonormal.jl +++ b/src/orthonormal.jl @@ -641,7 +641,7 @@ function reskeworthogonalize!!( v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ClassicalSymplecticGramSchmidt ) where {T} s = similar(x) ## EXTRA ALLOCATION - (v, x) = skeworthogonalize!!(v, b, s, ClassicalSymplecticGramSchmidt()) + (v, s) = skeworthogonalize!!(v, b, s, ClassicalSymplecticGramSchmidt()) x .+= s return (v, x) end @@ -703,10 +703,8 @@ function reskeworthogonalize!!( i_even = 2m h_e = inner(b[i_odd], v) h_f = inner(b[i_even], v) - s_odd = -h_f - s_even = h_e - x[i_odd] += s_odd - x[i_even] += s_even + x[i_odd] -= h_f + x[i_even] += h_e v = add!!(v, b[i_odd], h_f) v = add!!(v, b[i_even], -h_e) end diff --git a/test/symplectic_arnoldi.jl b/test/symplectic_arnoldi.jl index 5c9ba781..13265bfe 100644 --- a/test/symplectic_arnoldi.jl +++ b/test/symplectic_arnoldi.jl @@ -41,7 +41,7 @@ return max_err end - for esr in (ESR1, ESR2) + for esr in (ESR1, ESR2, ESR3m) algs = ( ClassicalSymplecticGramSchmidt(esr), ModifiedSymplecticGramSchmidt(esr), @@ -53,22 +53,14 @@ @testset "$alg" for alg in algs W1, H1, r1, b1 = run_symplectic_arnoldi(alg, false) @test max_symplectic_error(W1) < 1.0e-8 - if alg == ClassicalSymplecticGramSchmidt2(ESR1) - @test_broken Diagonal(domain) * W1 ≈ W1 * H1 + r1[] * b1' atol = 1.0e-8 - else - @test Diagonal(domain) * W1 ≈ W1 * H1 + r1[] * b1' atol = 1.0e-8 - end + @test Diagonal(domain) * W1 ≈ W1 * H1 + r1[] * b1' atol = 1.0e-8 W2, H2, r2, b2 = run_symplectic_arnoldi(alg, true) @test W1 ≈ W2 @test H1 ≈ H2 @test r1[] ≈ r2[] atol = 1.0e-8 @test max_symplectic_error(W2) < 1.0e-8 - if alg == ClassicalSymplecticGramSchmidt2(ESR1) - @test_broken Diagonal(domain) * W2 ≈ W2 * H2 + r2[] * b2' atol = 1.0e-8 - else - @test Diagonal(domain) * W2 ≈ W2 * H2 + r2[] * b2' atol = 1.0e-8 - end + @test Diagonal(domain) * W2 ≈ W2 * H2 + r2[] * b2' atol = 1.0e-8 end end end From d1956eb52fb5b841fc1e5c8a844c49a10bd876b7 Mon Sep 17 00:00:00 2001 From: Simeon David Schaub Date: Wed, 11 Feb 2026 17:00:36 +0100 Subject: [PATCH 09/13] fix some typos --- src/algorithms.jl | 2 +- src/factorizations/arnoldi.jl | 8 ++++---- src/orthonormal.jl | 8 ++++---- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/algorithms.jl b/src/algorithms.jl index d9f918c9..771a6771 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -81,7 +81,7 @@ ModifiedGramSchmidtIR() = ModifiedGramSchmidtIR(1 / sqrt(2)) # Daniel-Gragg-Kauf # Skew-orthogonalization for symplectic bases """ - @enum ESR ESR1 ESR2 ESR3 + @enum ESR ESR1 ESR2 ESR3m Enum for selecting Elementary Symplectic factorization (ESR) variant: - `ESR1`: r11 = ||x₁||, r12 = 0 diff --git a/src/factorizations/arnoldi.jl b/src/factorizations/arnoldi.jl index 1b4da90a..010fda7f 100644 --- a/src/factorizations/arnoldi.jl +++ b/src/factorizations/arnoldi.jl @@ -135,7 +135,7 @@ end function initialize(iter::ArnoldiIterator; verbosity::Int = KrylovDefaults.verbosity[]) # initialize without using eltype x₀ = iter.x₀ - if iter.orth isa Orthogonalizer || iter.orth.esr != ESR3 + if iter.orth isa Orthogonalizer || iter.orth.esr != ESR3m β₀ = norm(x₀) else β₀ = one(scalartype(x₀)) @@ -201,7 +201,7 @@ function initialize!( end H = empty!(state.H) - if iter.orth isa Orthogonalizer || iter.orth.esr != ESR3 + if iter.orth isa Orthogonalizer || iter.orth.esr != ESR3m β₀ = norm(x₀) else β₀ = one(scalartype(x₀)) @@ -240,7 +240,7 @@ function expand!( V = state.V H = state.H r = state.r - if iter.orth isa Orthogonalizer || iter.orth.esr == ESR3 + if iter.orth isa Orthogonalizer || iter.orth.esr == ESR3m β = normres(state) else β = iseven(k) ? normres(state) : norm(r) @@ -291,5 +291,5 @@ function arnoldirecurrence!!( ) w = apply(operator, last(V)) r, h = skeworthogonalize!!(w, V, h, orth) - return r, iseven(length(V)) ? (orth.esr == ESR3 ? one(scalartype(r)) : norm(r)) : inner(last(V), r) + return r, iseven(length(V)) ? (orth.esr == ESR3m ? one(scalartype(r)) : norm(r)) : inner(last(V), r) end diff --git a/src/orthonormal.jl b/src/orthonormal.jl index b986ebaa..5210f41d 100644 --- a/src/orthonormal.jl +++ b/src/orthonormal.jl @@ -744,11 +744,11 @@ end # For even vectors: scale so that ω(partner, v) = 1 skeworthonormalize(v, args...) = skeworthonormalize!!(scale(v, VectorInterface.One()), args...) -function skeworthonormalize!!(v, b::SymplecticBasis, args...) - out = skeworthogonalize!!(v, b, args...) # out[1] === v +function skeworthonormalize!!(v, b::SymplecticBasis, x::AbstractVector, alg::SkewOrthogonalizer) + out = skeworthogonalize!!(v, b, x, alg) if iseven(length(b)) # Adding odd vector: normalize with standard norm - β = norm(v) + β = alg.esr == ESR3m ? one(scalartype(v)) : norm(v) v = scale!!(v, inv(β)) else # Adding even vector: scale so that ω(partner, v) = 1 @@ -756,7 +756,7 @@ function skeworthonormalize!!(v, b::SymplecticBasis, args...) β = inner(last(b), v) v = scale!!(v, inv(β)) end - return tuple(v, β, Base.tail(out)...) + return (v, β, Base.tail(out)...) end """ From 7e59b3ca32a182ee0eb17aa9b42b859266dc69a9 Mon Sep 17 00:00:00 2001 From: Simeon David Schaub Date: Thu, 19 Feb 2026 10:58:17 +0100 Subject: [PATCH 10/13] introduce `SymplecticFormVec` instead of extending `InnerProductVec` --- src/KrylovKit.jl | 3 +- src/factorizations/arnoldi.jl | 14 +- src/innerproductvec.jl | 258 ++++++++++++++++++++++++++-------- src/orthonormal.jl | 47 ++++--- test/symplectic_arnoldi.jl | 4 +- 5 files changed, 240 insertions(+), 86 deletions(-) diff --git a/src/KrylovKit.jl b/src/KrylovKit.jl index be11da31..f0015553 100644 --- a/src/KrylovKit.jl +++ b/src/KrylovKit.jl @@ -45,7 +45,8 @@ export ModifiedSymplecticGramSchmidtIR export LanczosIterator, BlockLanczosIterator, ArnoldiIterator, GKLIterator, BiArnoldiIterator export CG, GMRES, BiCGStab, Lanczos, BlockLanczos, Arnoldi, GKL, GolubYe, LSMR, BiArnoldi export KrylovDefaults, EigSorter -export RecursiveVec, InnerProductVec, Block, SymplecticBasis, numpairs +export RecursiveVec, InnerProductVec, SymplecticFormVec, Block, SymplecticBasis, numpairs +export symplecticform # Multithreading const _NTHREADS = Ref(1) diff --git a/src/factorizations/arnoldi.jl b/src/factorizations/arnoldi.jl index 010fda7f..e93a1cba 100644 --- a/src/factorizations/arnoldi.jl +++ b/src/factorizations/arnoldi.jl @@ -142,7 +142,11 @@ function initialize(iter::ArnoldiIterator; verbosity::Int = KrylovDefaults.verbo end iszero(β₀) && throw(ArgumentError("initial vector should not have norm zero")) Ax₀ = apply(iter.operator, x₀) - α = inner(x₀, Ax₀) / (β₀ * β₀) + if iter.orth isa Orthogonalizer + α = inner(x₀, Ax₀) / (β₀ * β₀) + else + α = symplecticform(x₀, Ax₀) / (β₀ * β₀) + end T = typeof(α) # scalar type of the Rayleigh quotient # this line determines the vector type that we will henceforth use # vector scalar type can be different from `T`, e.g. for real inner products @@ -177,7 +181,7 @@ function initialize(iter::ArnoldiIterator; verbosity::Int = KrylovDefaults.verbo iszero(α) && throw(ArgumentError("initial vector and its image are symplectically orthogonal")) V = SymplecticBasis([v]) if iter.orth.esr == ESR2 - r12 = standard_dot(v, r) + r12 = inner(v, r) r = add!!(r, v, -r12) H = T[r12, α] else @@ -213,11 +217,11 @@ function initialize!( β = norm(r) push!(H, α, β) else - α = inner(V[1], w) + α = symplecticform(V[1], w) iszero(α) && throw(ArgumentError("initial vector and its image are symplectically orthogonal")) r = w if iter.orth.esr == ESR2 - r12 = standard_dot(V[1], r) + r12 = inner(V[1], r) r = add!!(r, V[1], -r12) push!(H, r12, α) else @@ -291,5 +295,5 @@ function arnoldirecurrence!!( ) w = apply(operator, last(V)) r, h = skeworthogonalize!!(w, V, h, orth) - return r, iseven(length(V)) ? (orth.esr == ESR3m ? one(scalartype(r)) : norm(r)) : inner(last(V), r) + return r, iseven(length(V)) ? (orth.esr == ESR3m ? one(scalartype(r)) : norm(r)) : symplecticform(last(V), r) end diff --git a/src/innerproductvec.jl b/src/innerproductvec.jl index 8e0a4176..9fb15669 100644 --- a/src/innerproductvec.jl +++ b/src/innerproductvec.jl @@ -1,5 +1,5 @@ """ - v = InnerProductVec(vec, dotf, [normf]) + v = InnerProductVec(vec, dotf) Create a new vector `v` from an existing vector `dotf` with a modified inner product given by `inner`. The vector `vec`, which can be any type (not necessarily `Vector`) that supports @@ -9,58 +9,50 @@ with scalars (both out of place and in place using `mul!`, `rmul!`, `axpy!` and applied to `v` is simply forwarded to `v.vec`. The inner product between vectors `v1 = InnerProductVec(vec1, dotf)` and `v2 = InnerProductVec(vec2, dotf)` is computed as `dot(v1, v2) = dotf(v1.vec, v2.vec) = dotf(vec1, vec2)`. The inner product between vectors -with different `dotf` functions is not defined. - -By default, the norm of `v::InnerProductVec` is defined as -`norm(v) = sqrt(real(dot(v, v))) = sqrt(real(dotf(vec, vec)))`. However, an optional third -argument `normf` can be provided to define a different norm function, which is useful for -cases like symplectic forms where the inner product is skew-symmetric and cannot define a -norm. In this case, `norm(v) = normf(v.vec)`. +with different `dotf` functions is not defined. Similarly, The norm of `v::InnerProductVec` +is defined as `v = sqrt(real(dot(v, v))) = sqrt(real(dotf(vec, vec)))`. In a (linear) map applied to `v`, the original vector can be obtained as `v.vec` or simply as `v[]`. """ -struct InnerProductVec{F, N, T} +struct InnerProductVec{F, T} vec::T dotf::F - normf::N end -InnerProductVec(vec, dotf) = InnerProductVec(vec, dotf, nothing) - -Base.:-(v::InnerProductVec) = InnerProductVec(-v.vec, v.dotf, v.normf) -function Base.:+(v::InnerProductVec{F, N}, w::InnerProductVec{F, N}) where {F, N} - return InnerProductVec(v.vec + w.vec, v.dotf, v.normf) +Base.:-(v::InnerProductVec) = InnerProductVec(-v.vec, v.dotf) +function Base.:+(v::InnerProductVec{F}, w::InnerProductVec{F}) where {F} + return InnerProductVec(v.vec + w.vec, v.dotf) end -function Base.:-(v::InnerProductVec{F, N}, w::InnerProductVec{F, N}) where {F, N} - return InnerProductVec(v.vec - w.vec, v.dotf, v.normf) +function Base.:-(v::InnerProductVec{F}, w::InnerProductVec{F}) where {F} + return InnerProductVec(v.vec - w.vec, v.dotf) end -Base.:*(v::InnerProductVec, a::Number) = InnerProductVec(v.vec * a, v.dotf, v.normf) -Base.:*(a::Number, v::InnerProductVec) = InnerProductVec(a * v.vec, v.dotf, v.normf) -Base.:/(v::InnerProductVec, a::Number) = InnerProductVec(v.vec / a, v.dotf, v.normf) -Base.:\(a::Number, v::InnerProductVec) = InnerProductVec(a \ v.vec, v.dotf, v.normf) +Base.:*(v::InnerProductVec, a::Number) = InnerProductVec(v.vec * a, v.dotf) +Base.:*(a::Number, v::InnerProductVec) = InnerProductVec(a * v.vec, v.dotf) +Base.:/(v::InnerProductVec, a::Number) = InnerProductVec(v.vec / a, v.dotf) +Base.:\(a::Number, v::InnerProductVec) = InnerProductVec(a \ v.vec, v.dotf) function Base.similar(v::InnerProductVec, (::Type{T}) = scalartype(v)) where {T} - return InnerProductVec(similar(v.vec), v.dotf, v.normf) + return InnerProductVec(similar(v.vec), v.dotf) end Base.getindex(v::InnerProductVec) = v.vec -function Base.copy!(w::InnerProductVec{F, N}, v::InnerProductVec{F, N}) where {F, N} +function Base.copy!(w::InnerProductVec{F}, v::InnerProductVec{F}) where {F} copy!(w.vec, v.vec) return w end function LinearAlgebra.mul!( - w::InnerProductVec{F, N}, a::Number, v::InnerProductVec{F, N} - ) where {F, N} + w::InnerProductVec{F}, a::Number, v::InnerProductVec{F} + ) where {F} mul!(w.vec, a, v.vec) return w end function LinearAlgebra.mul!( - w::InnerProductVec{F, N}, v::InnerProductVec{F, N}, a::Number - ) where {F, N} + w::InnerProductVec{F}, v::InnerProductVec{F}, a::Number + ) where {F} mul!(w.vec, v.vec, a) return w end @@ -71,79 +63,233 @@ function LinearAlgebra.rmul!(v::InnerProductVec, a::Number) end function LinearAlgebra.axpy!( - a::Number, v::InnerProductVec{F, N}, w::InnerProductVec{F, N} - ) where {F, N} + a::Number, v::InnerProductVec{F}, w::InnerProductVec{F} + ) where {F} axpy!(a, v.vec, w.vec) return w end function LinearAlgebra.axpby!( - a::Number, v::InnerProductVec{F, N}, b, w::InnerProductVec{F, N} - ) where {F, N} + a::Number, v::InnerProductVec{F}, b, w::InnerProductVec{F} + ) where {F} axpby!(a, v.vec, b, w.vec) return w end -function LinearAlgebra.dot(v::InnerProductVec{F, N}, w::InnerProductVec{F, N}) where {F, N} +function LinearAlgebra.dot(v::InnerProductVec{F}, w::InnerProductVec{F}) where {F} return v.dotf(v.vec, w.vec) end -VectorInterface.scalartype(::Type{<:InnerProductVec{F, N, T}}) where {F, N, T} = scalartype(T) +VectorInterface.scalartype(::Type{<:InnerProductVec{F, T}}) where {F, T} = scalartype(T) function VectorInterface.zerovector(v::InnerProductVec, T::Type{<:Number}) - return InnerProductVec(zerovector(v.vec, T), v.dotf, v.normf) + return InnerProductVec(zerovector(v.vec, T), v.dotf) end function VectorInterface.zerovector!(v::InnerProductVec) - return InnerProductVec(zerovector!(v.vec), v.dotf, v.normf) + return InnerProductVec(zerovector!(v.vec), v.dotf) end function VectorInterface.zerovector!!(v::InnerProductVec) - return InnerProductVec(zerovector!!(v.vec), v.dotf, v.normf) + return InnerProductVec(zerovector!!(v.vec), v.dotf) end function VectorInterface.scale(v::InnerProductVec, a::Number) - return InnerProductVec(scale(v.vec, a), v.dotf, v.normf) + return InnerProductVec(scale(v.vec, a), v.dotf) end function VectorInterface.scale!!(v::InnerProductVec, a::Number) - return InnerProductVec(scale!!(v.vec, a), v.dotf, v.normf) + return InnerProductVec(scale!!(v.vec, a), v.dotf) end function VectorInterface.scale!(v::InnerProductVec, a::Number) scale!(v.vec, a) return v end function VectorInterface.scale!!( - w::InnerProductVec{F, N}, v::InnerProductVec{F, N}, a::Number - ) where {F, N} - return InnerProductVec(scale!!(w.vec, v.vec, a), w.dotf, w.normf) + w::InnerProductVec{F}, v::InnerProductVec{F}, a::Number + ) where {F} + return InnerProductVec(scale!!(w.vec, v.vec, a), w.dotf) end function VectorInterface.scale!( - w::InnerProductVec{F, N}, v::InnerProductVec{F, N}, a::Number - ) where {F, N} + w::InnerProductVec{F}, v::InnerProductVec{F}, a::Number + ) where {F} scale!(w.vec, v.vec, a) return w end function VectorInterface.add( - v::InnerProductVec{F, N}, w::InnerProductVec{F, N}, a::Number, b::Number - ) where {F, N} - return InnerProductVec(add(v.vec, w.vec, a, b), v.dotf, v.normf) + v::InnerProductVec{F}, w::InnerProductVec{F}, a::Number, b::Number + ) where {F} + return InnerProductVec(add(v.vec, w.vec, a, b), v.dotf) end function VectorInterface.add!!( - v::InnerProductVec{F, N}, w::InnerProductVec{F, N}, a::Number, b::Number - ) where {F, N} - return InnerProductVec(add!!(v.vec, w.vec, a, b), v.dotf, v.normf) + v::InnerProductVec{F}, w::InnerProductVec{F}, a::Number, b::Number + ) where {F} + return InnerProductVec(add!!(v.vec, w.vec, a, b), v.dotf) end function VectorInterface.add!( - v::InnerProductVec{F, N}, w::InnerProductVec{F, N}, a::Number, b::Number - ) where {F, N} + v::InnerProductVec{F}, w::InnerProductVec{F}, a::Number, b::Number + ) where {F} add!(v.vec, w.vec, a, b) return v end -function VectorInterface.inner(v::InnerProductVec{F, N}, w::InnerProductVec{F, N}) where {F, N} +function VectorInterface.inner(v::InnerProductVec{F}, w::InnerProductVec{F}) where {F} return v.dotf(v.vec, w.vec) end -VectorInterface.norm(v::InnerProductVec{F, Nothing}) where {F} = sqrt(real(inner(v, v))) -VectorInterface.norm(v::InnerProductVec) = v.normf(v.vec) +VectorInterface.norm(v::InnerProductVec) = sqrt(real(inner(v, v))) + +""" + v = SymplecticFormVec(vec, skewf) + +Create a new vector `v` from an existing vector `vec` with a custom symplectic form given +by `skewf`. The vector `vec`, which can be any type that supports the basic vector interface +required by KrylovKit, is wrapped in a custom struct `v::SymplecticFormVec`. All vector +space functionality is forwarded to `v.vec`, and the inner product `inner(v1, v2)` is +computed as the standard inner product `inner(v1.vec, v2.vec)`. The symplectic form between +vectors `v1 = SymplecticFormVec(vec1, skewf)` and `v2 = SymplecticFormVec(vec2, skewf)` is +computed as `symplecticform(v1, v2) = skewf(vec1, vec2)`. + +In a (linear) map applied to `v`, the original vector can be obtained as `v.vec` or simply +as `v[]`. +""" +struct SymplecticFormVec{F, T} + vec::T + skewf::F +end + +Base.:-(v::SymplecticFormVec) = SymplecticFormVec(-v.vec, v.skewf) +function Base.:+(v::SymplecticFormVec{F}, w::SymplecticFormVec{F}) where {F} + return SymplecticFormVec(v.vec + w.vec, v.skewf) +end +function Base.:-(v::SymplecticFormVec{F}, w::SymplecticFormVec{F}) where {F} + return SymplecticFormVec(v.vec - w.vec, v.skewf) +end +Base.:*(v::SymplecticFormVec, a::Number) = SymplecticFormVec(v.vec * a, v.skewf) +Base.:*(a::Number, v::SymplecticFormVec) = SymplecticFormVec(a * v.vec, v.skewf) +Base.:/(v::SymplecticFormVec, a::Number) = SymplecticFormVec(v.vec / a, v.skewf) +Base.:\(a::Number, v::SymplecticFormVec) = SymplecticFormVec(a \ v.vec, v.skewf) + +function Base.similar(v::SymplecticFormVec, (::Type{T}) = scalartype(v)) where {T} + return SymplecticFormVec(similar(v.vec), v.skewf) +end + +Base.getindex(v::SymplecticFormVec) = v.vec + +function Base.copy!(w::SymplecticFormVec{F}, v::SymplecticFormVec{F}) where {F} + copy!(w.vec, v.vec) + return w +end -standard_dot(v::InnerProductVec, w::InnerProductVec) = dot(v.vec, w.vec) -standard_dot(v, w) = dot(v, w) +function LinearAlgebra.mul!( + w::SymplecticFormVec{F}, a::Number, v::SymplecticFormVec{F} + ) where {F} + mul!(w.vec, a, v.vec) + return w +end + +function LinearAlgebra.mul!( + w::SymplecticFormVec{F}, v::SymplecticFormVec{F}, a::Number + ) where {F} + mul!(w.vec, v.vec, a) + return w +end + +function LinearAlgebra.rmul!(v::SymplecticFormVec, a::Number) + rmul!(v.vec, a) + return v +end + +function LinearAlgebra.axpy!( + a::Number, v::SymplecticFormVec{F}, w::SymplecticFormVec{F} + ) where {F} + axpy!(a, v.vec, w.vec) + return w +end +function LinearAlgebra.axpby!( + a::Number, v::SymplecticFormVec{F}, b, w::SymplecticFormVec{F} + ) where {F} + axpby!(a, v.vec, b, w.vec) + return w +end + +function LinearAlgebra.dot(v::SymplecticFormVec{F}, w::SymplecticFormVec{F}) where {F} + return dot(v.vec, w.vec) +end + +VectorInterface.scalartype(::Type{<:SymplecticFormVec{F, T}}) where {F, T} = scalartype(T) + +function VectorInterface.zerovector(v::SymplecticFormVec, T::Type{<:Number}) + return SymplecticFormVec(zerovector(v.vec, T), v.skewf) +end +function VectorInterface.zerovector!(v::SymplecticFormVec) + return SymplecticFormVec(zerovector!(v.vec), v.skewf) +end +function VectorInterface.zerovector!!(v::SymplecticFormVec) + return SymplecticFormVec(zerovector!!(v.vec), v.skewf) +end + +function VectorInterface.scale(v::SymplecticFormVec, a::Number) + return SymplecticFormVec(scale(v.vec, a), v.skewf) +end +function VectorInterface.scale!!(v::SymplecticFormVec, a::Number) + return SymplecticFormVec(scale!!(v.vec, a), v.skewf) +end +function VectorInterface.scale!(v::SymplecticFormVec, a::Number) + scale!(v.vec, a) + return v +end +function VectorInterface.scale!!( + w::SymplecticFormVec{F}, v::SymplecticFormVec{F}, a::Number + ) where {F} + return SymplecticFormVec(scale!!(w.vec, v.vec, a), w.skewf) +end +function VectorInterface.scale!( + w::SymplecticFormVec{F}, v::SymplecticFormVec{F}, a::Number + ) where {F} + scale!(w.vec, v.vec, a) + return w +end + +function VectorInterface.add( + v::SymplecticFormVec{F}, w::SymplecticFormVec{F}, a::Number, b::Number + ) where {F} + return SymplecticFormVec(add(v.vec, w.vec, a, b), v.skewf) +end +function VectorInterface.add!!( + v::SymplecticFormVec{F}, w::SymplecticFormVec{F}, a::Number, b::Number + ) where {F} + return SymplecticFormVec(add!!(v.vec, w.vec, a, b), v.skewf) +end +function VectorInterface.add!( + v::SymplecticFormVec{F}, w::SymplecticFormVec{F}, a::Number, b::Number + ) where {F} + add!(v.vec, w.vec, a, b) + return v +end + +function VectorInterface.inner(v::SymplecticFormVec{F}, w::SymplecticFormVec{F}) where {F} + return inner(v.vec, w.vec) +end + +VectorInterface.norm(v::SymplecticFormVec) = sqrt(real(inner(v, v))) + +# Symplectic form +""" + symplecticform(v, w) + +Compute the symplectic form `ω(v, w)` between two vectors. For `SymplecticFormVec` vectors, +this uses the custom skew-symmetric form provided at construction. For `AbstractVector`, +this computes the standard symplectic form `ω(v, w) = Σ_i (v[2i-1]*w[2i] - v[2i]*w[2i-1])`. +""" +function symplecticform(v::SymplecticFormVec{F}, w::SymplecticFormVec{F}) where {F} + return v.skewf(v.vec, w.vec) +end + +function symplecticform(v::AbstractVector, w::AbstractVector) + length(v) == length(w) || throw(DimensionMismatch()) + n = length(v) + T = promote_type(eltype(v), eltype(w)) + result = zero(T) + @inbounds for i in 1:2:n + result += v[i] * w[i + 1] - v[i + 1] * w[i] + end + return result +end diff --git a/src/orthonormal.jl b/src/orthonormal.jl index 5210f41d..455e0a01 100644 --- a/src/orthonormal.jl +++ b/src/orthonormal.jl @@ -74,20 +74,23 @@ end """ project!!(y::AbstractVector, b::Basis, x, - [α::Number = 1, β::Number = 0, r = Base.OneTo(length(b))]) + [α::Number = 1, β::Number = 0, r = Base.OneTo(length(b))]; + innerfun = inner) For a given basis `b`, compute the expansion coefficients `y` resulting from projecting the vector `x` onto the subspace spanned by `b`; more specifically this computes ``` - y[j] = β*y[j] + α * inner(b[r[j]], x) + y[j] = β*y[j] + α * innerfun(b[r[j]], x) ``` -for all ``j ∈ r``. +for all ``j ∈ r``. The keyword `innerfun` allows using a custom bilinear form instead of +the default `inner`. """ function project!!( y::AbstractVector, b::Basis, x, - α::Number = true, β::Number = false, r = Base.OneTo(length(b)) + α::Number = true, β::Number = false, r = Base.OneTo(length(b)); + innerfun = inner ) # no specialized routine for IndexLinear x because reduction dimension is large dimension length(y) == length(r) || throw(DimensionMismatch()) @@ -96,9 +99,9 @@ function project!!( Threads.@spawn for j in $J @inbounds begin if β == 0 - y[j] = α * inner(b[r[j]], x) + y[j] = α * innerfun(b[r[j]], x) else - y[j] = β * y[j] + α * inner(b[r[j]], x) + y[j] = β * y[j] + α * innerfun(b[r[j]], x) end end end @@ -107,9 +110,9 @@ function project!!( for j in 1:length(r) @inbounds begin if β == 0 - y[j] = α * inner(b[r[j]], x) + y[j] = α * innerfun(b[r[j]], x) else - y[j] = β * y[j] + α * inner(b[r[j]], x) + y[j] = β * y[j] + α * innerfun(b[r[j]], x) end end end @@ -624,11 +627,11 @@ function skeworthogonalize!!( np = numpairs(b) idx_odd = 1:2:(2np - 1) idx_even = 2:2:2np - project!!(view(x, idx_odd), b, v, -1, 0, idx_even) - project!!(view(x, idx_even), b, v, 1, 0, idx_odd) + project!!(view(x, idx_odd), b, v, -1, 0, idx_even; innerfun = symplecticform) + project!!(view(x, idx_even), b, v, 1, 0, idx_odd; innerfun = symplecticform) if isodd(length(b)) if alg.esr == ESR2 - x[2np + 1] = standard_dot(last(b), v) + x[2np + 1] = inner(last(b), v) else x[2np + 1] = zero(eltype(x)) end @@ -676,8 +679,8 @@ function skeworthogonalize!!( for m in 1:np i_odd = 2m - 1 i_even = 2m - h_e = inner(b[i_odd], v) - h_f = inner(b[i_even], v) + h_e = symplecticform(b[i_odd], v) + h_f = symplecticform(b[i_even], v) x[i_odd] = -h_f x[i_even] = h_e v = add!!(v, b[i_odd], h_f) @@ -685,7 +688,7 @@ function skeworthogonalize!!( end if isodd(length(b)) if alg.esr == ESR2 - x[2np + 1] = standard_dot(last(b), v) + x[2np + 1] = inner(last(b), v) v = add!!(v, last(b), -x[2np + 1]) else x[2np + 1] = zero(eltype(x)) @@ -701,15 +704,15 @@ function reskeworthogonalize!!( for m in 1:np i_odd = 2m - 1 i_even = 2m - h_e = inner(b[i_odd], v) - h_f = inner(b[i_even], v) + h_e = symplecticform(b[i_odd], v) + h_f = symplecticform(b[i_even], v) x[i_odd] -= h_f x[i_even] += h_e v = add!!(v, b[i_odd], h_f) v = add!!(v, b[i_even], -h_e) end if isodd(length(b)) && alg.esr == ESR2 - r11 = standard_dot(last(b), v) + r11 = inner(last(b), v) x[2np + 1] += r11 v = add!!(v, last(b), -r11) end @@ -753,7 +756,7 @@ function skeworthonormalize!!(v, b::SymplecticBasis, x::AbstractVector, alg::Ske else # Adding even vector: scale so that ω(partner, v) = 1 # The partner is the last vector in b (the odd vector of the current pair) - β = inner(last(b), v) + β = symplecticform(last(b), v) v = scale!!(v, inv(β)) end return (v, β, Base.tail(out)...) @@ -767,10 +770,10 @@ Skew-orthogonalize vector `v` against all the vectors in the symplectic basis `b skew-orthogonalization algorithm `alg` of type [`SkewOrthogonalizer`](@ref), and return the resulting vector `w` and the overlap coefficients `x` of `v` with the basis vectors in `b`. -The skew-orthogonalization uses the skew-symmetric form `ω = inner`, which is expected to -satisfy `ω(u, v) = -ω(v, u)` when vectors are appropriately wrapped (e.g., using -`InnerProductVec`). For a symplectic basis with pairs `(u_{2m-1}, u_{2m})` where -`ω(u_{2m-1}, u_{2m}) = 1`, the skew-orthogonalization ensures: +The skew-orthogonalization uses the symplectic form [`symplecticform`](@ref), which is +expected to satisfy `ω(u, v) = -ω(v, u)`. For vectors wrapped in `SymplecticFormVec`, a +custom form function can be provided. For a symplectic basis with pairs `(u_{2m-1}, u_{2m})` +where `ω(u_{2m-1}, u_{2m}) = 1`, the skew-orthogonalization ensures: - `ω(u_{2m-1}, w) = 0` for all `m` - `ω(u_{2m}, w) = 0` for all `m` diff --git a/test/symplectic_arnoldi.jl b/test/symplectic_arnoldi.jl index 13265bfe..f36cf3ea 100644 --- a/test/symplectic_arnoldi.jl +++ b/test/symplectic_arnoldi.jl @@ -10,9 +10,9 @@ ) function run_symplectic_arnoldi(orth, inplace_init) - v₀ = InnerProductVec(fill(1 / sqrt(sum(w, domain)), n), skew_dot, norm) + v₀ = SymplecticFormVec(fill(1 / sqrt(sum(w, domain)), n), skew_dot) itr = ArnoldiIterator( - u -> InnerProductVec(domain .* u[], u.dotf, u.normf), + u -> SymplecticFormVec(domain .* u[], u.skewf), v₀, orth, ) From 24ccb87e2a918aae48320c722eb20d9edd8ae740 Mon Sep 17 00:00:00 2001 From: Simeon David Schaub Date: Sat, 21 Feb 2026 11:32:22 +0100 Subject: [PATCH 11/13] implement review suggestions --- docs/src/man/implementation.md | 22 +- src/KrylovKit.jl | 6 + src/innerproductvec.jl | 158 --------------- src/orthonormal.jl | 293 ++------------------------- src/skeworthonormal.jl | 359 +++++++++++++++++++++++++++++++++ src/symplecticformvec.jl | 191 ++++++++++++++++++ 6 files changed, 591 insertions(+), 438 deletions(-) create mode 100644 src/skeworthonormal.jl create mode 100644 src/symplecticformvec.jl diff --git a/docs/src/man/implementation.md b/docs/src/man/implementation.md index 99d4bf55..2753a50a 100644 --- a/docs/src/man/implementation.md +++ b/docs/src/man/implementation.md @@ -50,11 +50,31 @@ KrylovKit.ClassicalSymplecticGramSchmidtIR KrylovKit.ModifiedSymplecticGramSchmidtIR ``` +The skew-orthogonalization algorithms require a symplectic form to be defined for the +vector type. This can be done by defining `KrylovKit.symplecticform(v, w)` for your vector +type, or by wrapping your vectors in [`SymplecticFormVec`](@ref): +```@docs +KrylovKit.symplecticform +KrylovKit.SymplecticFormVec +``` + +The normalization step in skew-orthonormalization is controlled by the elementary symplectic +reflection (ESR) variant: +```@docs +KrylovKit.ESR +``` + The expansion coefficients of a general vector in terms of a given orthonormal basis can be obtained as ```@docs KrylovKit.project!! ``` -whereas the inverse calculation is obtained as + +For symplectic bases, the analogous projection using the symplectic form is +```@docs +KrylovKit.skewproject!! +``` + +The inverse calculation, reconstructing a vector from expansion coefficients, is obtained as ```@docs KrylovKit.unproject!! ``` diff --git a/src/KrylovKit.jl b/src/KrylovKit.jl index f0015553..6876c885 100644 --- a/src/KrylovKit.jl +++ b/src/KrylovKit.jl @@ -174,6 +174,9 @@ include("algorithms.jl") # OrthonormalBasis, orthogonalization and orthonormalization methods include("orthonormal.jl") +# SymplecticBasis, skew-orthogonalization and skew-orthonormalization methods +include("skeworthonormal.jl") + # Dense linear algebra structures and functions used in the algorithms below include("dense/givens.jl") include("dense/linalg.jl") @@ -250,6 +253,9 @@ end # vectors with modified inner product include("innerproductvec.jl") +# vectors with custom symplectic form +include("symplecticformvec.jl") + # support for real _realinner(v, w) = real(inner(v, w)) const RealVec{V} = InnerProductVec{typeof(_realinner), V} diff --git a/src/innerproductvec.jl b/src/innerproductvec.jl index 9fb15669..d04e9df1 100644 --- a/src/innerproductvec.jl +++ b/src/innerproductvec.jl @@ -135,161 +135,3 @@ function VectorInterface.inner(v::InnerProductVec{F}, w::InnerProductVec{F}) whe end VectorInterface.norm(v::InnerProductVec) = sqrt(real(inner(v, v))) - -""" - v = SymplecticFormVec(vec, skewf) - -Create a new vector `v` from an existing vector `vec` with a custom symplectic form given -by `skewf`. The vector `vec`, which can be any type that supports the basic vector interface -required by KrylovKit, is wrapped in a custom struct `v::SymplecticFormVec`. All vector -space functionality is forwarded to `v.vec`, and the inner product `inner(v1, v2)` is -computed as the standard inner product `inner(v1.vec, v2.vec)`. The symplectic form between -vectors `v1 = SymplecticFormVec(vec1, skewf)` and `v2 = SymplecticFormVec(vec2, skewf)` is -computed as `symplecticform(v1, v2) = skewf(vec1, vec2)`. - -In a (linear) map applied to `v`, the original vector can be obtained as `v.vec` or simply -as `v[]`. -""" -struct SymplecticFormVec{F, T} - vec::T - skewf::F -end - -Base.:-(v::SymplecticFormVec) = SymplecticFormVec(-v.vec, v.skewf) -function Base.:+(v::SymplecticFormVec{F}, w::SymplecticFormVec{F}) where {F} - return SymplecticFormVec(v.vec + w.vec, v.skewf) -end -function Base.:-(v::SymplecticFormVec{F}, w::SymplecticFormVec{F}) where {F} - return SymplecticFormVec(v.vec - w.vec, v.skewf) -end -Base.:*(v::SymplecticFormVec, a::Number) = SymplecticFormVec(v.vec * a, v.skewf) -Base.:*(a::Number, v::SymplecticFormVec) = SymplecticFormVec(a * v.vec, v.skewf) -Base.:/(v::SymplecticFormVec, a::Number) = SymplecticFormVec(v.vec / a, v.skewf) -Base.:\(a::Number, v::SymplecticFormVec) = SymplecticFormVec(a \ v.vec, v.skewf) - -function Base.similar(v::SymplecticFormVec, (::Type{T}) = scalartype(v)) where {T} - return SymplecticFormVec(similar(v.vec), v.skewf) -end - -Base.getindex(v::SymplecticFormVec) = v.vec - -function Base.copy!(w::SymplecticFormVec{F}, v::SymplecticFormVec{F}) where {F} - copy!(w.vec, v.vec) - return w -end - -function LinearAlgebra.mul!( - w::SymplecticFormVec{F}, a::Number, v::SymplecticFormVec{F} - ) where {F} - mul!(w.vec, a, v.vec) - return w -end - -function LinearAlgebra.mul!( - w::SymplecticFormVec{F}, v::SymplecticFormVec{F}, a::Number - ) where {F} - mul!(w.vec, v.vec, a) - return w -end - -function LinearAlgebra.rmul!(v::SymplecticFormVec, a::Number) - rmul!(v.vec, a) - return v -end - -function LinearAlgebra.axpy!( - a::Number, v::SymplecticFormVec{F}, w::SymplecticFormVec{F} - ) where {F} - axpy!(a, v.vec, w.vec) - return w -end -function LinearAlgebra.axpby!( - a::Number, v::SymplecticFormVec{F}, b, w::SymplecticFormVec{F} - ) where {F} - axpby!(a, v.vec, b, w.vec) - return w -end - -function LinearAlgebra.dot(v::SymplecticFormVec{F}, w::SymplecticFormVec{F}) where {F} - return dot(v.vec, w.vec) -end - -VectorInterface.scalartype(::Type{<:SymplecticFormVec{F, T}}) where {F, T} = scalartype(T) - -function VectorInterface.zerovector(v::SymplecticFormVec, T::Type{<:Number}) - return SymplecticFormVec(zerovector(v.vec, T), v.skewf) -end -function VectorInterface.zerovector!(v::SymplecticFormVec) - return SymplecticFormVec(zerovector!(v.vec), v.skewf) -end -function VectorInterface.zerovector!!(v::SymplecticFormVec) - return SymplecticFormVec(zerovector!!(v.vec), v.skewf) -end - -function VectorInterface.scale(v::SymplecticFormVec, a::Number) - return SymplecticFormVec(scale(v.vec, a), v.skewf) -end -function VectorInterface.scale!!(v::SymplecticFormVec, a::Number) - return SymplecticFormVec(scale!!(v.vec, a), v.skewf) -end -function VectorInterface.scale!(v::SymplecticFormVec, a::Number) - scale!(v.vec, a) - return v -end -function VectorInterface.scale!!( - w::SymplecticFormVec{F}, v::SymplecticFormVec{F}, a::Number - ) where {F} - return SymplecticFormVec(scale!!(w.vec, v.vec, a), w.skewf) -end -function VectorInterface.scale!( - w::SymplecticFormVec{F}, v::SymplecticFormVec{F}, a::Number - ) where {F} - scale!(w.vec, v.vec, a) - return w -end - -function VectorInterface.add( - v::SymplecticFormVec{F}, w::SymplecticFormVec{F}, a::Number, b::Number - ) where {F} - return SymplecticFormVec(add(v.vec, w.vec, a, b), v.skewf) -end -function VectorInterface.add!!( - v::SymplecticFormVec{F}, w::SymplecticFormVec{F}, a::Number, b::Number - ) where {F} - return SymplecticFormVec(add!!(v.vec, w.vec, a, b), v.skewf) -end -function VectorInterface.add!( - v::SymplecticFormVec{F}, w::SymplecticFormVec{F}, a::Number, b::Number - ) where {F} - add!(v.vec, w.vec, a, b) - return v -end - -function VectorInterface.inner(v::SymplecticFormVec{F}, w::SymplecticFormVec{F}) where {F} - return inner(v.vec, w.vec) -end - -VectorInterface.norm(v::SymplecticFormVec) = sqrt(real(inner(v, v))) - -# Symplectic form -""" - symplecticform(v, w) - -Compute the symplectic form `ω(v, w)` between two vectors. For `SymplecticFormVec` vectors, -this uses the custom skew-symmetric form provided at construction. For `AbstractVector`, -this computes the standard symplectic form `ω(v, w) = Σ_i (v[2i-1]*w[2i] - v[2i]*w[2i-1])`. -""" -function symplecticform(v::SymplecticFormVec{F}, w::SymplecticFormVec{F}) where {F} - return v.skewf(v.vec, w.vec) -end - -function symplecticform(v::AbstractVector, w::AbstractVector) - length(v) == length(w) || throw(DimensionMismatch()) - n = length(v) - T = promote_type(eltype(v), eltype(w)) - result = zero(T) - @inbounds for i in 1:2:n - result += v[i] * w[i + 1] - v[i + 1] * w[i] - end - return result -end diff --git a/src/orthonormal.jl b/src/orthonormal.jl index 455e0a01..a79417fc 100644 --- a/src/orthonormal.jl +++ b/src/orthonormal.jl @@ -73,24 +73,24 @@ function _use_multithreaded_array_kernel(::Type{<:Basis{T}}) where {T} end """ - project!!(y::AbstractVector, b::Basis, x, - [α::Number = 1, β::Number = 0, r = Base.OneTo(length(b))]; - innerfun = inner) + project!!(y::AbstractVector, b::OrthonormalBasis, x, + [α::Number = 1, β::Number = 0, r = Base.OneTo(length(b))]) -For a given basis `b`, compute the expansion coefficients `y` resulting from +For a given orthonormal basis `b`, compute the expansion coefficients `y` resulting from projecting the vector `x` onto the subspace spanned by `b`; more specifically this computes ``` - y[j] = β*y[j] + α * innerfun(b[r[j]], x) + y[j] = β*y[j] + α * inner(b[r[j]], x) ``` -for all ``j ∈ r``. The keyword `innerfun` allows using a custom bilinear form instead of -the default `inner`. +for all ``j ∈ r``. + +See also the method for [`SymplecticBasis`](@ref) which uses [`symplecticform`](@ref) +instead of `inner`. """ function project!!( - y::AbstractVector, b::Basis, x, - α::Number = true, β::Number = false, r = Base.OneTo(length(b)); - innerfun = inner + y::AbstractVector, b::OrthonormalBasis, x, + α::Number = true, β::Number = false, r = Base.OneTo(length(b)) ) # no specialized routine for IndexLinear x because reduction dimension is large dimension length(y) == length(r) || throw(DimensionMismatch()) @@ -99,9 +99,9 @@ function project!!( Threads.@spawn for j in $J @inbounds begin if β == 0 - y[j] = α * innerfun(b[r[j]], x) + y[j] = α * inner(b[r[j]], x) else - y[j] = β * y[j] + α * innerfun(b[r[j]], x) + y[j] = β * y[j] + α * inner(b[r[j]], x) end end end @@ -110,9 +110,9 @@ function project!!( for j in 1:length(r) @inbounds begin if β == 0 - y[j] = α * innerfun(b[r[j]], x) + y[j] = α * inner(b[r[j]], x) else - y[j] = β * y[j] + α * innerfun(b[r[j]], x) + y[j] = β * y[j] + α * inner(b[r[j]], x) end end end @@ -557,268 +557,3 @@ and its concrete subtypes [`ClassicalGramSchmidt`](@ref), [`ModifiedGramSchmidt` [`ClassicalGramSchmidtIR`](@ref) and [`ModifiedGramSchmidtIR`](@ref). """ orthonormalize, orthonormalize!! - -# Definition of a symplectic (Darboux) basis -""" - SymplecticBasis{T} <: Basis{T} - -A list of vector like objects of type `T` that form a symplectic (Darboux) basis with -respect to a skew-symmetric bilinear form `ω`. A symplectic basis satisfies the relations -`ω(u_{2m-1}, u_{2n}) = δ_{mn}`, `ω(u_{2m}, u_{2n}) = 0`, and `ω(u_{2m-1}, u_{2n-1}) = 0`. -See also [`Basis`](@ref). - -Skew-orthonormality of the vectors contained in an instance `b` of `SymplecticBasis` is not -checked when elements are added; it is up to the algorithm that constructs `b` to guarantee -skew-orthonormality. - -Vectors are added in pairs: odd-indexed vectors `u_{2m-1}` and their symplectic partners -`u_{2m}`. The function [`skeworthogonalize`](@ref) or [`skeworthonormalize`](@ref) can be -used to skew-orthogonalize a new vector with respect to the existing basis. These functions -require a skew-symmetric form `ω` to be provided. Note that in-place versions -[`skeworthogonalize!!`](@ref) or [`skeworthonormalize!!`](@ref) are also available. -""" -struct SymplecticBasis{T} <: Basis{T} - basis::Vector{T} -end -SymplecticBasis{T}() where {T} = SymplecticBasis{T}(Vector{T}(undef, 0)) - -# Iterator methods for SymplecticBasis -Base.IteratorSize(::Type{<:SymplecticBasis}) = Base.HasLength() -Base.IteratorEltype(::Type{<:SymplecticBasis}) = Base.HasEltype() - -Base.length(b::SymplecticBasis) = length(b.basis) -Base.eltype(b::SymplecticBasis{T}) where {T} = T - -Base.iterate(b::SymplecticBasis) = Base.iterate(b.basis) -Base.iterate(b::SymplecticBasis, state) = Base.iterate(b.basis, state) - -Base.getindex(b::SymplecticBasis, i) = getindex(b.basis, i) -Base.setindex!(b::SymplecticBasis, i, q) = setindex!(b.basis, i, q) -Base.firstindex(b::SymplecticBasis) = firstindex(b.basis) -Base.lastindex(b::SymplecticBasis) = lastindex(b.basis) - -Base.first(b::SymplecticBasis) = first(b.basis) -Base.last(b::SymplecticBasis) = last(b.basis) - -Base.popfirst!(b::SymplecticBasis) = popfirst!(b.basis) -Base.pop!(b::SymplecticBasis) = pop!(b.basis) -Base.push!(b::SymplecticBasis{T}, q::T) where {T} = (push!(b.basis, q); return b) -Base.empty!(b::SymplecticBasis) = (empty!(b.basis); return b) -Base.sizehint!(b::SymplecticBasis, k::Int) = (sizehint!(b.basis, k); return b) -Base.resize!(b::SymplecticBasis, k::Int) = (resize!(b.basis, k); return b) - -numpairs(b::SymplecticBasis) = div(length(b), 2) - -# Skew-orthogonalization of a vector against a given SymplecticBasis -skeworthogonalize(v, args...) = skeworthogonalize!!(scale(v, true), args...) - -function skeworthogonalize!!( - v::T, b::SymplecticBasis{T}, alg::SkewOrthogonalizer - ) where {T} - S = promote_type(scalartype(v), scalartype(T)) - c = Vector{S}(undef, length(b)) - return skeworthogonalize!!(v, b, c, alg) -end - -# See pages 3-4 of https://people.math.ethz.ch/%7Eacannas/Papers/lsg.pdf -function skeworthogonalize!!( - v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ClassicalSymplecticGramSchmidt - ) where {T} - np = numpairs(b) - idx_odd = 1:2:(2np - 1) - idx_even = 2:2:2np - project!!(view(x, idx_odd), b, v, -1, 0, idx_even; innerfun = symplecticform) - project!!(view(x, idx_even), b, v, 1, 0, idx_odd; innerfun = symplecticform) - if isodd(length(b)) - if alg.esr == ESR2 - x[2np + 1] = inner(last(b), v) - else - x[2np + 1] = zero(eltype(x)) - end - end - v = unproject!!(v, b, x, -1, 1) - return (v, x) -end - -function reskeworthogonalize!!( - v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ClassicalSymplecticGramSchmidt - ) where {T} - s = similar(x) ## EXTRA ALLOCATION - (v, s) = skeworthogonalize!!(v, b, s, ClassicalSymplecticGramSchmidt()) - x .+= s - return (v, x) -end - -function skeworthogonalize!!( - v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ClassicalSymplecticGramSchmidt2 - ) where {T} - csgs = ClassicalSymplecticGramSchmidt(alg.esr) - (v, x) = skeworthogonalize!!(v, b, x, csgs) - return reskeworthogonalize!!(v, b, x, csgs) -end - -function skeworthogonalize!!( - v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ClassicalSymplecticGramSchmidtIR - ) where {T} - csgs = ClassicalSymplecticGramSchmidt(alg.esr) - nold = norm(v) - (v, x) = skeworthogonalize!!(v, b, x, csgs) - nnew = norm(v) - while eps(one(nnew)) < nnew < alg.η * nold - nold = nnew - (v, x) = reskeworthogonalize!!(v, b, x, csgs) - nnew = norm(v) - end - return (v, x) -end - -function skeworthogonalize!!( - v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ModifiedSymplecticGramSchmidt - ) where {T} - np = numpairs(b) - for m in 1:np - i_odd = 2m - 1 - i_even = 2m - h_e = symplecticform(b[i_odd], v) - h_f = symplecticform(b[i_even], v) - x[i_odd] = -h_f - x[i_even] = h_e - v = add!!(v, b[i_odd], h_f) - v = add!!(v, b[i_even], -h_e) - end - if isodd(length(b)) - if alg.esr == ESR2 - x[2np + 1] = inner(last(b), v) - v = add!!(v, last(b), -x[2np + 1]) - else - x[2np + 1] = zero(eltype(x)) - end - end - return (v, x) -end - -function reskeworthogonalize!!( - v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ModifiedSymplecticGramSchmidt - ) where {T} - np = numpairs(b) - for m in 1:np - i_odd = 2m - 1 - i_even = 2m - h_e = symplecticform(b[i_odd], v) - h_f = symplecticform(b[i_even], v) - x[i_odd] -= h_f - x[i_even] += h_e - v = add!!(v, b[i_odd], h_f) - v = add!!(v, b[i_even], -h_e) - end - if isodd(length(b)) && alg.esr == ESR2 - r11 = inner(last(b), v) - x[2np + 1] += r11 - v = add!!(v, last(b), -r11) - end - return (v, x) -end - -function skeworthogonalize!!( - v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ModifiedSymplecticGramSchmidt2 - ) where {T} - msgs = ModifiedSymplecticGramSchmidt(alg.esr) - (v, x) = skeworthogonalize!!(v, b, x, msgs) - return reskeworthogonalize!!(v, b, x, msgs) -end - -function skeworthogonalize!!( - v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ModifiedSymplecticGramSchmidtIR - ) where {T} - msgs = ModifiedSymplecticGramSchmidt(alg.esr) - nold = norm(v) - (v, x) = skeworthogonalize!!(v, b, x, msgs) - nnew = norm(v) - while eps(one(nnew)) < nnew < alg.η * nold - nold = nnew - (v, x) = reskeworthogonalize!!(v, b, x, msgs) - nnew = norm(v) - end - return (v, x) -end - -# Skew-orthonormalization: skew-orthogonalization and normalization -# For odd vectors: normalize with standard norm -# For even vectors: scale so that ω(partner, v) = 1 -skeworthonormalize(v, args...) = skeworthonormalize!!(scale(v, VectorInterface.One()), args...) - -function skeworthonormalize!!(v, b::SymplecticBasis, x::AbstractVector, alg::SkewOrthogonalizer) - out = skeworthogonalize!!(v, b, x, alg) - if iseven(length(b)) - # Adding odd vector: normalize with standard norm - β = alg.esr == ESR3m ? one(scalartype(v)) : norm(v) - v = scale!!(v, inv(β)) - else - # Adding even vector: scale so that ω(partner, v) = 1 - # The partner is the last vector in b (the odd vector of the current pair) - β = symplecticform(last(b), v) - v = scale!!(v, inv(β)) - end - return (v, β, Base.tail(out)...) -end - -""" - skeworthogonalize(v, b::SymplecticBasis, [x::AbstractVector,] alg::SkewOrthogonalizer) -> w, x - skeworthogonalize!!(v, b::SymplecticBasis, [x::AbstractVector,] alg::SkewOrthogonalizer) -> w, x - -Skew-orthogonalize vector `v` against all the vectors in the symplectic basis `b` using the -skew-orthogonalization algorithm `alg` of type [`SkewOrthogonalizer`](@ref), and return the -resulting vector `w` and the overlap coefficients `x` of `v` with the basis vectors in `b`. - -The skew-orthogonalization uses the symplectic form [`symplecticform`](@ref), which is -expected to satisfy `ω(u, v) = -ω(v, u)`. For vectors wrapped in `SymplecticFormVec`, a -custom form function can be provided. For a symplectic basis with pairs `(u_{2m-1}, u_{2m})` -where `ω(u_{2m-1}, u_{2m}) = 1`, the skew-orthogonalization ensures: -- `ω(u_{2m-1}, w) = 0` for all `m` -- `ω(u_{2m}, w) = 0` for all `m` - -In case of `skeworthogonalize!!`, the vector `v` is mutated in place. In both functions, -storage for the overlap coefficients `x` can be provided as optional argument -`x::AbstractVector` with `length(x) >= length(b)`. - -Note that `w` is not normalized, see also [`skeworthonormalize`](@ref). - -For more information on possible skew-orthogonalization algorithms, see -[`SkewOrthogonalizer`](@ref) and its concrete subtypes -[`ClassicalSymplecticGramSchmidt`](@ref), [`ModifiedSymplecticGramSchmidt`](@ref), -[`ClassicalSymplecticGramSchmidt2`](@ref), [`ModifiedSymplecticGramSchmidt2`](@ref), -[`ClassicalSymplecticGramSchmidtIR`](@ref) and [`ModifiedSymplecticGramSchmidtIR`](@ref). -""" -skeworthogonalize, skeworthogonalize!! - -""" - skeworthonormalize(v, b::SymplecticBasis, [x::AbstractVector,] alg::SkewOrthogonalizer) -> w, β, x - skeworthonormalize!!(v, b::SymplecticBasis, [x::AbstractVector,] alg::SkewOrthogonalizer) -> w, β, x - -Skew-orthonormalize vector `v` against all the vectors in the symplectic basis `b` using -the skew-orthogonalization algorithm `alg` of type [`SkewOrthogonalizer`](@ref). - -The normalization depends on the current length of the basis: - -**When `length(b)` is even** (adding an odd vector): returns the resulting vector `w` -normalized to unit norm (`‖w‖ = 1`), the norm `β = ‖v‖` after skew-orthogonalizing, and -the overlap coefficients `x`. - -**When `length(b)` is odd** (adding an even vector): returns the resulting vector `w` -scaled such that `ω(last(b), w) = 1`, where `last(b)` is the odd partner of the pair. -Returns the scaling factor `β = ω(last(b), v)` after skew-orthogonalizing, and the overlap -coefficients `x`. - -In case of `skeworthonormalize!!`, the vector `v` is mutated in place. In both functions, -storage for the overlap coefficients `x` can be provided as optional argument -`x::AbstractVector` with `length(x) >= length(b)`. - -See [`skeworthogonalize`](@ref) if `w` does not need to be normalized. - -For more information on possible skew-orthogonalization algorithms, see -[`SkewOrthogonalizer`](@ref) and its concrete subtypes -[`ClassicalSymplecticGramSchmidt`](@ref), [`ModifiedSymplecticGramSchmidt`](@ref), -[`ClassicalSymplecticGramSchmidt2`](@ref), [`ModifiedSymplecticGramSchmidt2`](@ref), -[`ClassicalSymplecticGramSchmidtIR`](@ref) and [`ModifiedSymplecticGramSchmidtIR`](@ref). -""" -skeworthonormalize, skeworthonormalize!! diff --git a/src/skeworthonormal.jl b/src/skeworthonormal.jl new file mode 100644 index 00000000..e83fcc16 --- /dev/null +++ b/src/skeworthonormal.jl @@ -0,0 +1,359 @@ +# Definition of a symplectic (Darboux) basis +""" + SymplecticBasis{T} <: Basis{T} + +A list of vector like objects of type `T` that form a symplectic (Darboux) basis with +respect to a skew-symmetric bilinear form `ω`. A symplectic basis satisfies the relations +`ω(u_{2m-1}, u_{2n}) = δ_{mn}`, `ω(u_{2m}, u_{2n}) = 0`, and `ω(u_{2m-1}, u_{2n-1}) = 0`. +See also [`Basis`](@ref). + +Skew-orthonormality of the vectors contained in an instance `b` of `SymplecticBasis` is not +checked when elements are added; it is up to the algorithm that constructs `b` to guarantee +skew-orthonormality. + +Vectors are added in pairs: odd-indexed vectors `u_{2m-1}` and their symplectic partners +`u_{2m}`. The function [`skeworthogonalize`](@ref) or [`skeworthonormalize`](@ref) can be +used to skew-orthogonalize a new vector with respect to the existing basis. These functions +require a skew-symmetric form `ω` to be defined for the vector type via +[`symplecticform`](@ref). + +# Symplectic form requirements + +To use a `SymplecticBasis` and the associated [`SkewOrthogonalizer`](@ref) algorithms with +your vector type, you must ensure that [`symplecticform(v, w)`](@ref symplecticform) is +defined. There are two approaches: + +1. **Define `KrylovKit.symplecticform(v, w)`** for your vector type. A default + implementation is provided for `AbstractVector` types, computing the canonical symplectic + form `ω(v, w) = vᵀ J w = Σᵢ (v[2i-1] w[2i] - v[2i] w[2i-1])`. + +2. **Wrap your vectors** in [`SymplecticFormVec(v, skewf)`](@ref SymplecticFormVec), which + carries a custom symplectic form function `skewf`. + +See also [`skeworthogonalize`](@ref), [`skeworthonormalize`](@ref), +[`SkewOrthogonalizer`](@ref), [`symplecticform`](@ref). +""" +struct SymplecticBasis{T} <: Basis{T} + basis::Vector{T} +end +SymplecticBasis{T}() where {T} = SymplecticBasis{T}(Vector{T}(undef, 0)) + +# Iterator methods for SymplecticBasis +Base.IteratorSize(::Type{<:SymplecticBasis}) = Base.HasLength() +Base.IteratorEltype(::Type{<:SymplecticBasis}) = Base.HasEltype() + +Base.length(b::SymplecticBasis) = length(b.basis) +Base.eltype(b::SymplecticBasis{T}) where {T} = T + +Base.iterate(b::SymplecticBasis) = Base.iterate(b.basis) +Base.iterate(b::SymplecticBasis, state) = Base.iterate(b.basis, state) + +Base.getindex(b::SymplecticBasis, i) = getindex(b.basis, i) +Base.setindex!(b::SymplecticBasis, i, q) = setindex!(b.basis, i, q) +Base.firstindex(b::SymplecticBasis) = firstindex(b.basis) +Base.lastindex(b::SymplecticBasis) = lastindex(b.basis) + +Base.first(b::SymplecticBasis) = first(b.basis) +Base.last(b::SymplecticBasis) = last(b.basis) + +Base.popfirst!(b::SymplecticBasis) = popfirst!(b.basis) +Base.pop!(b::SymplecticBasis) = pop!(b.basis) +Base.push!(b::SymplecticBasis{T}, q::T) where {T} = (push!(b.basis, q); return b) +Base.empty!(b::SymplecticBasis) = (empty!(b.basis); return b) +Base.sizehint!(b::SymplecticBasis, k::Int) = (sizehint!(b.basis, k); return b) +Base.resize!(b::SymplecticBasis, k::Int) = (resize!(b.basis, k); return b) + +numpairs(b::SymplecticBasis) = div(length(b), 2) + +# Projection using the symplectic form for SymplecticBasis +""" + skewproject!!(y::AbstractVector, b::SymplecticBasis, x, + [α::Number = 1, β::Number = 0, r = Base.OneTo(length(b))]) + +Project the vector `x` onto the symplectic basis `b` using the symplectic form +[`symplecticform`](@ref). The projection coefficients are computed using the symplectic +partner of each output index, with appropriate signs dictated by the Darboux basis +structure: + +- for odd `j`: `y[j] = β*y[j] - α * symplecticform(b[r[j+1]], x)` +- for even `j`: `y[j] = β*y[j] + α * symplecticform(b[r[j-1]], x)` + +That is, each coefficient is computed from the symplectic form with the *partner* basis +vector (even partner for odd indices, odd partner for even indices), with a sign flip for +odd indices. This cross-partner convention ensures correct skew-orthogonal projection in a +symplectic basis where `ω(u_{2m-1}, u_{2m}) = 1`. +""" +function skewproject!!( + y::AbstractVector, b::SymplecticBasis, x, + α::Number = true, β::Number = false, r = Base.OneTo(length(b)) + ) + length(y) == length(r) || throw(DimensionMismatch()) + if get_num_threads() > 1 + @sync for J in splitrange(1:length(r), get_num_threads()) + Threads.@spawn for j in $J + @inbounds begin + i = isodd(j) ? r[j + 1] : r[j - 1] + σ = isodd(j) ? -1 : 1 + if β == 0 + y[j] = α * σ * symplecticform(b[i], x) + else + y[j] = β * y[j] + α * σ * symplecticform(b[i], x) + end + end + end + end + else + for j in 1:length(r) + @inbounds begin + i = isodd(j) ? r[j + 1] : r[j - 1] + σ = isodd(j) ? -1 : 1 + if β == 0 + y[j] = α * σ * symplecticform(b[i], x) + else + y[j] = β * y[j] + α * σ * symplecticform(b[i], x) + end + end + end + end + return y +end + +# Skew-orthogonalization of a vector against a given SymplecticBasis +skeworthogonalize(v, args...) = skeworthogonalize!!(scale(v, true), args...) + +function skeworthogonalize!!( + v::T, b::SymplecticBasis{T}, alg::SkewOrthogonalizer + ) where {T} + S = promote_type(scalartype(v), scalartype(T)) + c = Vector{S}(undef, length(b)) + return skeworthogonalize!!(v, b, c, alg) +end + +# See pages 3-4 of https://people.math.ethz.ch/%7Eacannas/Papers/lsg.pdf +function skeworthogonalize!!( + v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ClassicalSymplecticGramSchmidt + ) where {T} + np = numpairs(b) + idx_pairs = 1:2np + skewproject!!(view(x, idx_pairs), b, v, 1, 0, idx_pairs) + if isodd(length(b)) + if alg.esr == ESR2 + x[2np + 1] = inner(last(b), v) + else + x[2np + 1] = zero(eltype(x)) + end + end + v = unproject!!(v, b, x, -1, 1) + return (v, x) +end + +function reskeworthogonalize!!( + v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ClassicalSymplecticGramSchmidt + ) where {T} + s = similar(x) ## EXTRA ALLOCATION + (v, s) = skeworthogonalize!!(v, b, s, alg) + x .+= s + return (v, x) +end + +function skeworthogonalize!!( + v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ClassicalSymplecticGramSchmidt2 + ) where {T} + csgs = ClassicalSymplecticGramSchmidt(alg.esr) + (v, x) = skeworthogonalize!!(v, b, x, csgs) + return reskeworthogonalize!!(v, b, x, csgs) +end + +function skeworthogonalize!!( + v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ClassicalSymplecticGramSchmidtIR + ) where {T} + csgs = ClassicalSymplecticGramSchmidt(alg.esr) + nold = norm(v) + (v, x) = skeworthogonalize!!(v, b, x, csgs) + nnew = norm(v) + while eps(one(nnew)) < nnew < alg.η * nold + nold = nnew + (v, x) = reskeworthogonalize!!(v, b, x, csgs) + nnew = norm(v) + end + return (v, x) +end + +function skeworthogonalize!!( + v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ModifiedSymplecticGramSchmidt + ) where {T} + np = numpairs(b) + for m in 1:np + i_odd = 2m - 1 + i_even = 2m + h_e = symplecticform(b[i_odd], v) + h_f = symplecticform(b[i_even], v) + x[i_odd] = -h_f + x[i_even] = h_e + v = add!!(v, b[i_odd], h_f) + v = add!!(v, b[i_even], -h_e) + end + if isodd(length(b)) + if alg.esr == ESR2 + x[2np + 1] = inner(last(b), v) + v = add!!(v, last(b), -x[2np + 1]) + else + x[2np + 1] = zero(eltype(x)) + end + end + return (v, x) +end + +function reskeworthogonalize!!( + v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ModifiedSymplecticGramSchmidt + ) where {T} + np = numpairs(b) + for m in 1:np + i_odd = 2m - 1 + i_even = 2m + h_e = symplecticform(b[i_odd], v) + h_f = symplecticform(b[i_even], v) + x[i_odd] -= h_f + x[i_even] += h_e + v = add!!(v, b[i_odd], h_f) + v = add!!(v, b[i_even], -h_e) + end + if isodd(length(b)) && alg.esr == ESR2 + r11 = inner(last(b), v) + x[2np + 1] += r11 + v = add!!(v, last(b), -r11) + end + return (v, x) +end + +function skeworthogonalize!!( + v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ModifiedSymplecticGramSchmidt2 + ) where {T} + msgs = ModifiedSymplecticGramSchmidt(alg.esr) + (v, x) = skeworthogonalize!!(v, b, x, msgs) + return reskeworthogonalize!!(v, b, x, msgs) +end + +function skeworthogonalize!!( + v::T, b::SymplecticBasis{T}, x::AbstractVector, alg::ModifiedSymplecticGramSchmidtIR + ) where {T} + msgs = ModifiedSymplecticGramSchmidt(alg.esr) + nold = norm(v) + (v, x) = skeworthogonalize!!(v, b, x, msgs) + nnew = norm(v) + while eps(one(nnew)) < nnew < alg.η * nold + nold = nnew + (v, x) = reskeworthogonalize!!(v, b, x, msgs) + nnew = norm(v) + end + return (v, x) +end + +# Skew-orthonormalization: skew-orthogonalization and normalization +# For odd vectors: normalize with standard norm +# For even vectors: scale so that ω(partner, v) = 1 +skeworthonormalize(v, args...) = skeworthonormalize!!(scale(v, VectorInterface.One()), args...) + +function skeworthonormalize!!(v, b::SymplecticBasis, x::AbstractVector, alg::SkewOrthogonalizer) + out = skeworthogonalize!!(v, b, x, alg) + if iseven(length(b)) + # Adding odd vector: normalize with standard norm / don't normalize if ESR3m + β = alg.esr == ESR3m ? one(scalartype(v)) : norm(v) + v = scale!!(v, inv(β)) + else + # Adding even vector: scale so that ω(partner, v) = 1 + # The partner is the last vector in b (the odd vector of the current pair) + β = symplecticform(last(b), v) + v = scale!!(v, inv(β)) + end + return (v, β, Base.tail(out)...) +end + +""" + skeworthogonalize(v, b::SymplecticBasis, [x::AbstractVector,] alg::SkewOrthogonalizer) -> w, x + skeworthogonalize!!(v, b::SymplecticBasis, [x::AbstractVector,] alg::SkewOrthogonalizer) -> w, x + +Skew-orthogonalize vector `v` against all the vectors in the symplectic basis `b` using the +skew-orthogonalization algorithm `alg` of type [`SkewOrthogonalizer`](@ref), and return the +resulting vector `w` and the overlap coefficients `x` of `v` with the basis vectors in `b`. + +The skew-orthogonalization uses the symplectic form [`symplecticform`](@ref), which is +expected to satisfy `ω(u, v) = -ω(v, u)`. For vectors wrapped in +[`SymplecticFormVec`](@ref), a custom form function can be provided. For a symplectic basis +with pairs `(u_{2m-1}, u_{2m})` where `ω(u_{2m-1}, u_{2m}) = 1`, the +skew-orthogonalization ensures: +- `ω(u_{2m-1}, w) = 0` for all `m` +- `ω(u_{2m}, w) = 0` for all `m` + +When the basis has odd length (i.e. after an odd vector has been added but before its even +partner), the behaviour depends on the [`ESR`](@ref) variant of the algorithm: +- `ESR1`: no additional projection is performed against the unpaired odd vector. +- `ESR2`: additionally projects `v` against the unpaired odd vector using the standard + inner product, i.e. `x[end] = inner(last(b), v)` and `w = v - x[end] * last(b)`. This + makes `s₁` and `s₂` orthogonal and tends to improve numerical stability. +- `ESR3m`: same as `ESR1` (no additional projection). + +See [Salam (2014)](https://journal.austms.org.au/ojs/index.php/ANZIAMJ/article/view/9380/1920) +for details on the elementary symplectic reflections and their effect on stability. + +In case of `skeworthogonalize!!`, the vector `v` is mutated in place. In both functions, +storage for the overlap coefficients `x` can be provided as optional argument +`x::AbstractVector` with `length(x) >= length(b)`. + +Note that `w` is not normalized, see also [`skeworthonormalize`](@ref). + +For more information on possible skew-orthogonalization algorithms, see +[`SkewOrthogonalizer`](@ref) and its concrete subtypes +[`ClassicalSymplecticGramSchmidt`](@ref), [`ModifiedSymplecticGramSchmidt`](@ref), +[`ClassicalSymplecticGramSchmidt2`](@ref), [`ModifiedSymplecticGramSchmidt2`](@ref), +[`ClassicalSymplecticGramSchmidtIR`](@ref) and [`ModifiedSymplecticGramSchmidtIR`](@ref). +""" +skeworthogonalize, skeworthogonalize!! + +""" + skeworthonormalize(v, b::SymplecticBasis, [x::AbstractVector,] alg::SkewOrthogonalizer) -> w, β, x + skeworthonormalize!!(v, b::SymplecticBasis, [x::AbstractVector,] alg::SkewOrthogonalizer) -> w, β, x + +Skew-orthonormalize vector `v` against all the vectors in the symplectic basis `b` using +the skew-orthogonalization algorithm `alg` of type [`SkewOrthogonalizer`](@ref). + +The normalization depends on the current length of the basis and the [`ESR`](@ref) variant +of the algorithm. Vectors are added in pairs `(s₁, s₂)` via an elementary symplectic +reflection (ESR) factorization of `(x₁, x₂)` (the input vectors before and after +skew-orthogonalization). The scaling factors `r₁₁` and `r₂₂` in the ESR depend on the +variant: + +- **`ESR1`**: `r₁₁ = ‖x₁‖`, `r₁₂ = 0`, `r₂₂ = ω(s₁, x₂)`. The odd vector `s₁` is + normalized to unit norm. +- **`ESR2`**: `r₁₁ = ‖x₁‖`, `r₁₂ = s₁ᵀ x₂`, `r₂₂ = ω(s₁, y)` where + `y = x₂ - r₁₂ s₁`. The odd vector `s₁` is normalized to unit norm, and the even + vector is additionally projected against `s₁` using the standard inner product before + symplectic scaling, improving stability. +- **`ESR3m`**: `r₁₁ = 1` (no norm scaling), `r₁₂ = 0`, `r₂₂ = ω(s₁, x₂)`. The odd + vector `s₁ = x₁` is used without normalization; `r₁₁` is recorded as + `one(scalartype(v))` rather than `‖x₁‖`. + +See [Salam (2014)](https://journal.austms.org.au/ojs/index.php/ANZIAMJ/article/view/9380/1920) +for details on the elementary symplectic reflections. + +**When `length(b)` is even** (adding an odd vector `s₁`): returns the resulting vector `w` +normalized to unit norm (`‖w‖ = 1`) for `ESR1`/`ESR2`, or unnormalized for `ESR3m`. The +scaling factor `β = r₁₁` is returned. + +**When `length(b)` is odd** (adding an even vector `s₂`): returns the resulting vector `w` +scaled such that `ω(last(b), w) = 1`, where `last(b)` is the odd partner of the pair. +The scaling factor `β = r₂₂ = ω(last(b), v)` after skew-orthogonalizing is returned. + +In case of `skeworthonormalize!!`, the vector `v` is mutated in place. In both functions, +storage for the overlap coefficients `x` can be provided as optional argument +`x::AbstractVector` with `length(x) >= length(b)`. + +See [`skeworthogonalize`](@ref) if `w` does not need to be normalized. + +For more information on possible skew-orthogonalization algorithms, see +[`SkewOrthogonalizer`](@ref) and its concrete subtypes +[`ClassicalSymplecticGramSchmidt`](@ref), [`ModifiedSymplecticGramSchmidt`](@ref), +[`ClassicalSymplecticGramSchmidt2`](@ref), [`ModifiedSymplecticGramSchmidt2`](@ref), +[`ClassicalSymplecticGramSchmidtIR`](@ref) and [`ModifiedSymplecticGramSchmidtIR`](@ref). +""" +skeworthonormalize, skeworthonormalize!! diff --git a/src/symplecticformvec.jl b/src/symplecticformvec.jl new file mode 100644 index 00000000..904d573d --- /dev/null +++ b/src/symplecticformvec.jl @@ -0,0 +1,191 @@ +""" + v = SymplecticFormVec(vec, skewf) + +Create a new vector `v` from an existing vector `vec` with a custom symplectic form given +by `skewf`. The vector `vec`, which can be any type that supports the basic vector interface +required by KrylovKit, is wrapped in a custom struct `v::SymplecticFormVec`. All vector +space functionality is forwarded to `v.vec`, and the inner product `inner(v1, v2)` is +computed as the standard inner product `inner(v1.vec, v2.vec)`. The symplectic form between +vectors `v1 = SymplecticFormVec(vec1, skewf)` and `v2 = SymplecticFormVec(vec2, skewf)` is +computed as `symplecticform(v1, v2) = skewf(vec1, vec2)`. + +In a (linear) map applied to `v`, the original vector can be obtained as `v.vec` or simply +as `v[]`. + +# Usage with SkewOrthogonalizer / Symplectic Arnoldi + +In order to use a [`SkewOrthogonalizer`](@ref) algorithm in an +[`ArnoldiFactorization`](@ref), the vector type used must support the +[`symplecticform`](@ref) function. There are two ways to achieve this: + +1. **Define `KrylovKit.symplecticform(v, w)`** directly for your vector type. For + `AbstractVector` types, a default implementation is already provided that computes the + standard symplectic form `ω(v, w) = Σᵢ (v[2i-1] w[2i] - v[2i] w[2i-1])`. + +2. **Wrap your vectors** in `SymplecticFormVec(v, skewf)`, where `skewf(v, w)` computes the + desired symplectic form between the unwrapped vectors `v` and `w`. This is analogous to + [`InnerProductVec`](@ref) for custom inner products. + +See also [`symplecticform`](@ref), [`SkewOrthogonalizer`](@ref). +""" +struct SymplecticFormVec{F, T} + vec::T + skewf::F +end + +Base.:-(v::SymplecticFormVec) = SymplecticFormVec(-v.vec, v.skewf) +function Base.:+(v::SymplecticFormVec{F}, w::SymplecticFormVec{F}) where {F} + return SymplecticFormVec(v.vec + w.vec, v.skewf) +end +function Base.:-(v::SymplecticFormVec{F}, w::SymplecticFormVec{F}) where {F} + return SymplecticFormVec(v.vec - w.vec, v.skewf) +end +Base.:*(v::SymplecticFormVec, a::Number) = SymplecticFormVec(v.vec * a, v.skewf) +Base.:*(a::Number, v::SymplecticFormVec) = SymplecticFormVec(a * v.vec, v.skewf) +Base.:/(v::SymplecticFormVec, a::Number) = SymplecticFormVec(v.vec / a, v.skewf) +Base.:\(a::Number, v::SymplecticFormVec) = SymplecticFormVec(a \ v.vec, v.skewf) + +function Base.similar(v::SymplecticFormVec, (::Type{T}) = scalartype(v)) where {T} + return SymplecticFormVec(similar(v.vec), v.skewf) +end + +Base.getindex(v::SymplecticFormVec) = v.vec + +function Base.copy!(w::SymplecticFormVec{F}, v::SymplecticFormVec{F}) where {F} + copy!(w.vec, v.vec) + return w +end + +function LinearAlgebra.mul!( + w::SymplecticFormVec{F}, a::Number, v::SymplecticFormVec{F} + ) where {F} + mul!(w.vec, a, v.vec) + return w +end + +function LinearAlgebra.mul!( + w::SymplecticFormVec{F}, v::SymplecticFormVec{F}, a::Number + ) where {F} + mul!(w.vec, v.vec, a) + return w +end + +function LinearAlgebra.rmul!(v::SymplecticFormVec, a::Number) + rmul!(v.vec, a) + return v +end + +function LinearAlgebra.axpy!( + a::Number, v::SymplecticFormVec{F}, w::SymplecticFormVec{F} + ) where {F} + axpy!(a, v.vec, w.vec) + return w +end +function LinearAlgebra.axpby!( + a::Number, v::SymplecticFormVec{F}, b, w::SymplecticFormVec{F} + ) where {F} + axpby!(a, v.vec, b, w.vec) + return w +end + +function LinearAlgebra.dot(v::SymplecticFormVec{F}, w::SymplecticFormVec{F}) where {F} + return dot(v.vec, w.vec) +end + +VectorInterface.scalartype(::Type{<:SymplecticFormVec{F, T}}) where {F, T} = scalartype(T) + +function VectorInterface.zerovector(v::SymplecticFormVec, T::Type{<:Number}) + return SymplecticFormVec(zerovector(v.vec, T), v.skewf) +end +function VectorInterface.zerovector!(v::SymplecticFormVec) + return SymplecticFormVec(zerovector!(v.vec), v.skewf) +end +function VectorInterface.zerovector!!(v::SymplecticFormVec) + return SymplecticFormVec(zerovector!!(v.vec), v.skewf) +end + +function VectorInterface.scale(v::SymplecticFormVec, a::Number) + return SymplecticFormVec(scale(v.vec, a), v.skewf) +end +function VectorInterface.scale!!(v::SymplecticFormVec, a::Number) + return SymplecticFormVec(scale!!(v.vec, a), v.skewf) +end +function VectorInterface.scale!(v::SymplecticFormVec, a::Number) + scale!(v.vec, a) + return v +end +function VectorInterface.scale!!( + w::SymplecticFormVec{F}, v::SymplecticFormVec{F}, a::Number + ) where {F} + return SymplecticFormVec(scale!!(w.vec, v.vec, a), w.skewf) +end +function VectorInterface.scale!( + w::SymplecticFormVec{F}, v::SymplecticFormVec{F}, a::Number + ) where {F} + scale!(w.vec, v.vec, a) + return w +end + +function VectorInterface.add( + v::SymplecticFormVec{F}, w::SymplecticFormVec{F}, a::Number, b::Number + ) where {F} + return SymplecticFormVec(add(v.vec, w.vec, a, b), v.skewf) +end +function VectorInterface.add!!( + v::SymplecticFormVec{F}, w::SymplecticFormVec{F}, a::Number, b::Number + ) where {F} + return SymplecticFormVec(add!!(v.vec, w.vec, a, b), v.skewf) +end +function VectorInterface.add!( + v::SymplecticFormVec{F}, w::SymplecticFormVec{F}, a::Number, b::Number + ) where {F} + add!(v.vec, w.vec, a, b) + return v +end + +function VectorInterface.inner(v::SymplecticFormVec{F}, w::SymplecticFormVec{F}) where {F} + return inner(v.vec, w.vec) +end + +VectorInterface.norm(v::SymplecticFormVec) = sqrt(real(inner(v, v))) + +# Symplectic form +""" + symplecticform(v, w) + +Compute the symplectic form `ω(v, w)` between two vectors `v` and `w`. + +The symplectic form is a skew-symmetric bilinear form, i.e. `ω(v, w) = -ω(w, v)`. + +# Methods + +For [`SymplecticFormVec`](@ref) vectors `v = SymplecticFormVec(vec, skewf)`, the symplectic +form is computed using the custom function: `symplecticform(v, w) = skewf(v.vec, w.vec)`. + +For `AbstractVector` types, a default implementation computes the standard (canonical) +symplectic form: + +```math +ω(v, w) = v^T J w = \\sum_i \\bigl(v_{2i-1}\\, w_{2i} - v_{2i}\\, w_{2i-1}\\bigr) +``` + +For other custom vector types, define a method `KrylovKit.symplecticform(v::MyVec, w::MyVec)` +to enable the use of [`SkewOrthogonalizer`](@ref) algorithms. + +See also [`SymplecticFormVec`](@ref), [`SkewOrthogonalizer`](@ref), +[`SymplecticBasis`](@ref). +""" +function symplecticform(v::SymplecticFormVec{F}, w::SymplecticFormVec{F}) where {F} + return v.skewf(v.vec, w.vec) +end + +function symplecticform(v::AbstractVector, w::AbstractVector) + length(v) == length(w) || throw(DimensionMismatch()) + n = length(v) + T = promote_type(eltype(v), eltype(w)) + result = zero(T) + @inbounds for i in 1:2:n + result += v[i] * w[i + 1] - v[i + 1] * w[i] + end + return result +end From 6e770efb048619cbee141ea682069ae6d7e0ed94 Mon Sep 17 00:00:00 2001 From: Simeon David Schaub Date: Sat, 21 Feb 2026 11:50:03 +0100 Subject: [PATCH 12/13] use `abs(F.H[end])` for `normres` with orthonormal basis --- src/factorizations/arnoldi.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/factorizations/arnoldi.jl b/src/factorizations/arnoldi.jl index e93a1cba..d63ee1f9 100644 --- a/src/factorizations/arnoldi.jl +++ b/src/factorizations/arnoldi.jl @@ -47,7 +47,8 @@ Base.eltype(::Type{<:ArnoldiFactorization{<:Any, S}}) where {S} = S basis(F::ArnoldiFactorization) = F.V rayleighquotient(F::ArnoldiFactorization) = PackedHessenberg(F.H, F.k) residual(F::ArnoldiFactorization) = F.r -@inbounds normres(F::ArnoldiFactorization) = F.H[end] +@inbounds normres(F::ArnoldiFactorization{<:Any, <:Any, <:OrthonormalBasis}) = abs(F.H[end]) +@inbounds normres(F::ArnoldiFactorization{<:Any, <:Any, <:SymplecticBasis}) = F.H[end] rayleighextension(F::ArnoldiFactorization) = SimpleBasisVector(F.k, F.k) # Arnoldi iteration for constructing the orthonormal basis of a Krylov subspace. From 362331172110bffbf2befb27d7b39cc14f709d06 Mon Sep 17 00:00:00 2001 From: Simeon David Schaub Date: Sat, 21 Feb 2026 12:00:39 +0100 Subject: [PATCH 13/13] revert accidentally committed change to `shrink!` --- src/factorizations/arnoldi.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/factorizations/arnoldi.jl b/src/factorizations/arnoldi.jl index d63ee1f9..580b9822 100644 --- a/src/factorizations/arnoldi.jl +++ b/src/factorizations/arnoldi.jl @@ -279,6 +279,7 @@ function shrink!(state::ArnoldiFactorization, k; verbosity::Int = KrylovDefaults if verbosity > EACHITERATION_LEVEL @info "Arnoldi reduction to dimension $k: subspace normres = $(normres2string(β))" end + state.r = scale!!(r, β) return state end