[RFC] Symplectic Arnoldi

docs/Project.toml

@@ -1,5 +1,6 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 
 [compat]
 Documenter = "1"

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,11 +30,51 @@ 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 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!!
 ```

src/KrylovKit.jl

@@ -32,14 +32,21 @@ 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 ESR, ESR1, ESR2, ESR3m
+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, SymplecticFormVec, Block, SymplecticBasis, numpairs
+export symplecticform
 
 # Multithreading
 const _NTHREADS = Ref(1)
@@ -167,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")
@@ -243,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}

src/algorithms.jl

@@ -79,6 +79,114 @@ struct ModifiedGramSchmidtIR{S <: Real} <: Reorthogonalizer
 end
 ModifiedGramSchmidtIR() = ModifiedGramSchmidtIR(1 / sqrt(2)) # Daniel-Gragg-Kaufman-Stewart
 
+# Skew-orthogonalization for symplectic bases
+"""
+    @enum ESR ESR1 ESR2 ESR3m
+
+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)
+- `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 ESR3m
+
+"""
+    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(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. The `esr` parameter
+selects the Elementary Symplectic factorization variant (default: ESR2, most stable).
+"""
+struct ClassicalSymplecticGramSchmidt <: SkewOrthogonalizer
+    esr::ESR
+end
+ClassicalSymplecticGramSchmidt() = ClassicalSymplecticGramSchmidt(ESR2)
+
+"""
+    ModifiedSymplecticGramSchmidt(esr::ESR = ESR2)
+
+Represents the modified symplectic Gram Schmidt algorithm for skew-orthogonalizing vectors
+to produce a symplectic (Darboux) basis. The `esr` parameter selects the Elementary
+Symplectic factorization variant (default: ESR2, most stable).
+"""
+struct ModifiedSymplecticGramSchmidt <: SkewOrthogonalizer
+    esr::ESR
+end
+ModifiedSymplecticGramSchmidt() = ModifiedSymplecticGramSchmidt(ESR2)
+
+# A single reorthogonalization always
+"""
+    ClassicalSymplecticGramSchmidt2(esr::ESR = ESR2)
+
+Represents the classical symplectic Gram Schmidt algorithm with a second reskew-
+orthogonalization step always taking place. The `esr` parameter selects the Elementary
+Symplectic factorization variant (default: ESR2, most stable).
+"""
+struct ClassicalSymplecticGramSchmidt2 <: SkewReorthogonalizer
+    esr::ESR
+end
+ClassicalSymplecticGramSchmidt2() = ClassicalSymplecticGramSchmidt2(ESR2)
+
+"""
+    ModifiedSymplecticGramSchmidt2(esr::ESR = ESR2)
+
+Represents the modified symplectic Gram Schmidt algorithm with a second reskew-
+orthogonalization step always taking place. The `esr` parameter selects the Elementary
+Symplectic factorization variant (default: ESR2, most stable).
+"""
+struct ModifiedSymplecticGramSchmidt2 <: SkewReorthogonalizer
+    esr::ESR
+end
+ModifiedSymplecticGramSchmidt2() = ModifiedSymplecticGramSchmidt2(ESR2)
+
+# Iterative reorthogonalization
+"""
+    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. The `esr`
+parameter selects the Elementary Symplectic factorization variant (default: ESR2, most stable).
+"""
+struct ClassicalSymplecticGramSchmidtIR{S <: Real} <: SkewReorthogonalizer
+    η::S
+    esr::ESR
+end
+ClassicalSymplecticGramSchmidtIR(η::Real = 1 / sqrt(2)) = ClassicalSymplecticGramSchmidtIR(η, ESR2)
+
+"""
+    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. The `esr`
+parameter selects the Elementary Symplectic factorization variant (default: ESR2, most stable).
+"""
+struct ModifiedSymplecticGramSchmidtIR{S <: Real} <: SkewReorthogonalizer
+    η::S
+    esr::ESR
+end
+ModifiedSymplecticGramSchmidtIR(η::Real = 1 / sqrt(2)) = ModifiedSymplecticGramSchmidtIR(η, ESR2)
+
 # Solving eigenvalue problems
 abstract type KrylovAlgorithm end
 
@@ -145,7 +253,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 +465,7 @@ end
 """
     GMRES(; krylovdim=KrylovDefaults.krylovdim[],
             maxiter=KrylovDefaults.maxiter[],
-            tol=KrylovDefaults.tol[], 
+            tol=KrylovDefaults.tol[],
             orth::Orthogonalizer=KrylovDefaults.orth,
             verbosity=KrylovDefaults.verbosity[])
 
@@ -460,7 +568,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 +592,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 +603,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.