[RFC] Symplectic Arnoldi by simeonschaub · Pull Request #151 · Jutho/KrylovKit.jl

simeonschaub · 2026-02-06T12:48:18Z

See the discussion in #150. On my test case, this already gives me
significantly better performance, and improved accuracy, if using
re-skew-orthogonalization, than my hand-rolled Arnoldi implementation. An initial sketch was written by Copilot, but
I revised and reviewed the result extensively.

closes #150

github-actions · 2026-02-06T12:48:49Z

Your PR no longer requires formatting changes. Thank you for your contribution!

codecov · 2026-02-06T20:43:24Z

Codecov Report

❌ Patch coverage is 70.92652% with 91 lines in your changes missing coverage. Please review.
✅ Project coverage is 85.37%. Comparing base (6c0cc8b) to head (d8267cd).

Files with missing lines	Patch %	Lines
src/symplecticformvec.jl	21.91%	57 Missing ⚠️
src/skeworthonormal.jl	84.09%	21 Missing ⚠️
src/factorizations/arnoldi.jl	92.63%	7 Missing ⚠️
src/algorithms.jl	50.00%	6 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #151      +/-   ##
==========================================
- Coverage   88.42%   85.37%   -3.05%     
==========================================
  Files          36       38       +2     
  Lines        3964     4219     +255     
==========================================
+ Hits         3505     3602      +97     
- Misses        459      617     +158

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

simeonschaub · 2026-02-10T13:48:24Z

Ok, I think I am happy with the implementation now from my side. This should enable some initial experimentation with symplectic Gram-Schmidt methods. Future work might include blocked methods such as described in https://journal.austms.org.au/ojs/index.php/ANZIAMJ/article/view/9380/1920, but I don't want to make this PR much bigger than it already is.

Jutho · 2026-02-18T11:47:23Z

Ok, I will have to study this a bit further and let it sink in, but my initial impression is that I am not a big fan of "recycling" inner and InnerProductVec for computing the symplectic form, which than results in hacky solutions such as having to introduce a separate norm and a standard_dot function. While I agree that could duplication is typically harder to maintain, in this case I wouldn't mind if the symplectic functionality is provided through a separate SymplecticFormVec that copies a lot of the functionality of InnerProductVec, and uses a new dedicated sympleticform (or other shorter name) primitive for applying the symplectic form. I prefer the mathematical correctness of such an approach (inner and symplecticform have fundamentally different behavior under interchanging the two arguments), and think that will lead to cleaner and easier to comprehend code.

simeonschaub · 2026-02-18T18:01:37Z

Sure, I can change the design! I wasn't too sure about extending InnerProductVec anyways, probably makes sense to introduce a new type for symplectic forms

Jutho · 2026-02-20T21:50:33Z

 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]


What is the rationale behind this change? Can normres be negative in the skew-orthogonal case? In the orthogonal case, I agree that F.H[end] is positive by construction, but in the complex case, it is a positive number encoded in a complex data type, and I would like normres to return a real number. So is real(F.H[end]) a solution?

Yes, for symplectic bases normres is allowed to be negative, since all the even vectors need to be scaled in a way such that ω(last(b), v) = +1. So if we don't allow normres to be negative ω(last(b), v) might come out to be -1 instead of +1, which is not what we want.

I haven't really played around with complex vectors yet, but I imagine the phase would be relevant there as well. Would a good compromise be to have two methods:

@inbounds normres(F::ArnoldiFactorization{<:Any, <:Any, <:OrthonormalBasis}) = abs(F.H[end]) @inbounds normres(F::ArnoldiFactorization{<:Any, <:Any, <:SymplecticBasis}) = F.H[end]

?

Jutho · 2026-02-20T23:04:55Z

Thanks for the changes; this is starting to look really nice. I have a few remaining comments and suggestions:

I think it would be easier if SymplecticFormVec just lives in its own source file, similar to InnerProductVec.
Similarly, I would move everything related to SymplecticBasis and skew-orthogonalization in a separate file "skeworthonormal.jl", that lives next to "orthonormal.jl".
I don't really like the innerfun keyword to project, since it is really completely determined by the type of b. I think this can more cleanly be solved by renaming project to orthproject and similarly have a skewproject, or simply by keeping the name project but have a separate method definition for OrthonormalBasis and SymplecticBasis. Again, I wouldn't mind that there is some code duplication here.
It would be great if you could at a bit more documentation, explaining that in order to use a SkewOrthogonalizer algorithm in ArnoldiFactorization, your vector type either needs to define KrylovKit.symplecticform, or declare which method is representing the symplectic form via SimplecticFormVec(v, form), and also explaining the default implementation for AbstractVector types.

simeonschaub · 2026-02-21T11:16:46Z

It would be great if you could at a bit more documentation, explaining that in order to use a SkewOrthogonalizer algorithm in ArnoldiFactorization, your vector type either needs to define KrylovKit.symplecticform, or declare which method is representing the symplectic form via SimplecticFormVec(v, form), and also explaining the default implementation for AbstractVector types.

Did you want me to add some more docs regarding this to the docstrings of ArnoldiFactorization and ArnoldiIteration as well? I just thought since InnerProductVec isn't mentioned there either, it might be more appropriate to add that explanation to the docstrings of SymplecticBasis and the (re)skeortho(gon|norm)alize(!!) functions instead and link to that in the implementation details docs. I can add it there as well though if you think it makes sense.

simeonschaub · 2026-03-26T11:47:35Z

Did you have a chance to take another look?

Jutho · 2026-04-24T16:19:50Z

My apologies for the long silence. Overall I think this is a great PR with very nice and consistent coding patterns.

The only part of which I am not convinced is including both the orthogonal and symplectic Arnoldi process in the same ArnoldiFactorization object (by introducing the type of basis as a new type parameter), i.e. the changes in the "factorizations/arnoldi.jl" file, for two reasons. Conceptually, I think most people think of Arnoldi as having an orthogonal basis for the Krylov subspace (so the symplectic part is quite hidden in this type parameter). Codewise, I don't like the resulting if ... else ... structures.

I would propose to have an explicit SymplecticArnoldiFactorization structure, and for those methods that would be identical for both structures, to use a Union{ArnoldiFactorization,SymplecticArnoldiFactorization} type specification for the factorization argument.

See the discussion in Jutho#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 Jutho#150

Instead of parameterizing `ArnoldiFactorization` on the basis type, the symplectic Arnoldi process now uses a dedicated `SymplecticArnoldiFactorization`. Methods that are identical for both structures use a `Union{ArnoldiFactorization, SymplecticArnoldiFactorization}` argument, while `initialize(!)`, `expand!` and `shrink!` now have separate method definitions for the orthogonal and symplectic case instead of `if ... else ...` branches on the orthogonalizer type. Assisted-By: Claude Fable 5 <noreply@anthropic.com>

- `shrink!` for `SymplecticArnoldiFactorization` rescaled the residual by `norm(r)` when shrinking to odd length, but the popped even-index basis vector is ω-normalized, so the stored coefficient `H[end]` is the correct scale. Since that coefficient is also what `norm(r)` evaluates to in the odd-index case, both `shrink!` and `expand!` now simply use `normres(state)` unconditionally, like their orthonormal counterparts. - the `Base.iterate` termination check now uses `abs(normres(state))`, as `normres` of a `SymplecticArnoldiFactorization` is a signed (or complex) scaling factor; a negative value previously terminated iteration immediately. - `expand!` now throws on symplectic breakdown (zero `normres`) instead of silently producing Inf/NaN basis vectors, matching the guard in `initialize`. - `skeworthonormalize!!` normalized the input binding instead of the vector returned by `skeworthogonalize!!`, which is wrong for vector types whose `!!` methods return new objects. - `skewproject!!` read `r[j + 1]` out of bounds under `@inbounds` for the default index range on an odd-length basis; the default now covers complete pairs and uneven ranges throw an `ArgumentError`. - `reskeworthogonalize!!` no longer allocates a fresh coefficient vector on every reorthogonalization pass of `ClassicalSymplecticGramSchmidtIR`. - add tests for `shrink!`, the `Base.iterate` protocol and `skeworthonormalize`. Assisted-By: Claude Fable 5 <noreply@anthropic.com>

simeonschaub · 2026-08-04T08:31:06Z

Sorry, I was a little busy with other stuff these last few months, but finally had time to come back to this! I addressed your review comments. I also had claude review the existing code and it found some correctness issues I had it fix and add some additional tests. I reviewed the fixes and they all seemed sensible to me

simeonschaub mentioned this pull request Feb 6, 2026

support for skew-orthogonalization / symplectic Gram-Schmidt #150

Open

simeonschaub force-pushed the sds/symplectic branch from c68189d to 8525cbd Compare February 9, 2026 22:11

Jutho reviewed Feb 20, 2026

View reviewed changes

simeonschaub force-pushed the sds/symplectic branch from 3623311 to c7d0666 Compare August 4, 2026 08:22

simeonschaub added 15 commits August 4, 2026 10:23

fix formatting

88b30bc

fix typos

692d319

fix inference issue

8cbf663

remove outdated type restriction again

8e2a789

fix docs build

d741b60

implement more stable ESR2, make it the default

81a1769

introduce ESR3m, fix reskeworthogonalize!! bug

a136773

fix some typos

a0eaaa1

introduce SymplecticFormVec instead of extending InnerProductVec

877bb0f

implement review suggestions

6d871de

use abs(F.H[end]) for normres with orthonormal basis

7cd8277

revert accidentally committed change to shrink!

94979b9

simeonschaub force-pushed the sds/symplectic branch from c7d0666 to d8267cd Compare August 4, 2026 08:24

Navigation Menu

[RFC] Symplectic Arnoldi - #151

[RFC] Symplectic Arnoldi#151
simeonschaub wants to merge 15 commits into
Jutho:masterfrom
simeonschaub:sds/symplectic

simeonschaub commented Feb 6, 2026 •

edited

Loading

Uh oh!

github-actions Bot commented Feb 6, 2026 •

edited

Loading

Uh oh!

codecov Bot commented Feb 6, 2026 •

edited

Loading

Uh oh!

simeonschaub commented Feb 10, 2026

Uh oh!

Jutho commented Feb 18, 2026

Uh oh!

simeonschaub commented Feb 18, 2026

Uh oh!

Jutho Feb 20, 2026

Uh oh!

simeonschaub Feb 21, 2026

Uh oh!

Jutho commented Feb 20, 2026

Uh oh!

simeonschaub commented Feb 21, 2026

Uh oh!

simeonschaub commented Mar 26, 2026

Uh oh!

Jutho commented Apr 24, 2026

Uh oh!

simeonschaub commented Aug 4, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Conversation

simeonschaub commented Feb 6, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

github-actions Bot commented Feb 6, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

codecov Bot commented Feb 6, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

simeonschaub commented Feb 10, 2026

Uh oh!

Jutho commented Feb 18, 2026

Uh oh!

simeonschaub commented Feb 18, 2026

Uh oh!

Jutho Feb 20, 2026

Choose a reason for hiding this comment

Uh oh!

simeonschaub Feb 21, 2026

Choose a reason for hiding this comment

Uh oh!

Jutho commented Feb 20, 2026

Uh oh!

simeonschaub commented Feb 21, 2026

Uh oh!

simeonschaub commented Mar 26, 2026

Uh oh!

Jutho commented Apr 24, 2026

Uh oh!

simeonschaub commented Aug 4, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

simeonschaub commented Feb 6, 2026 •

edited

Loading

github-actions Bot commented Feb 6, 2026 •

edited

Loading

codecov Bot commented Feb 6, 2026 •

edited

Loading