support for skew-orthogonalization / symplectic Gram-Schmidt

[simeonschaub]

I have a use case for an Arnoldi iteration, but where instead of being orthonormal, my U in A * U = U * H is actually symplectic, such that Uᵀ * Jₘ * U = Jₙ where J₂ₙ = Iₙ ⊗ [0 1; -1 0]. This can be achieved with a skew-symmetric version of the Gram-Schmidt process, such as described on pages 3-4 of https://people.math.ethz.ch/%7Eacannas/Papers/lsg.pdf. Most of the literature I could find on this only discusses symplectic Arnoldi methods in the context of Hamiltonian systems, but I am interested in calculating skew-orthogonal polynomials w.r.t. a skew-symmetric inner product, so my operator A is not Hamiltonian.

My question is, would something like this be in scope for KrylovKit.jl? Looking at the design in https://github.com/Jutho/KrylovKit.jl/blob/master/src/orthonormal.jl, I think it might lend itself quite nicely to extending it to different notions of "orthogonality" and there could then be different SkewOrthogonalizers, which do different versions of symplectic Gram-Schmidt. Any thoughts regarding such a design? If that sounds interesting, I can try putting up a proof of concept.

[Jutho]

That sounds very interesting. While I am aware of symplectic diagonalisation, I was not aware about Krylov methods for this (I assume you will need a dense symlectic diagonalization as well for the subproblem?).

I am definitely very supportive of having this functionality available in KrylovKit, but with very little bandwidth to do anything myself in the next few months. However, if you want to put a PR together, then I am definitely willing to make some time for feedback and comments.

[simeonschaub]

Thanks for the feedback, that sounds great! I will see what it takes to incorporate this into KrylovKit then.

(I assume you will need a dense symplectic diagonalization as well for the subproblem?)

In my case, I actually only care about the Arnoldi factorization itself, since that already gives me the evaluated first n skew-orthogonal polynomials as U and the recurrence relation, from which the polynomials can be construced explicitly, as H. But I could also see applications where you would want (skew) Eigen/singular values of such an operator.

[Jutho]

Ok, let's start with the Arnoldi factorization itself. That will yield a smaller PR which is easier to review anyway 😄 . Maybe I can use this as a building block of a symplectic diagonalisation at some later stage.

[simeonschaub]

Ok, so one issue I am currently encountering is that IIUC, v₀ is not part of the basis V, so it is not included in the orthogonalization procedure. I am assuming this was deliberate? For my use case at least, it would be much cleaner to not have a manual initialization step and just treat [v₀] as the first basis. Then the procedure can just start from k = 0. Unfortunately this doesn't really fit in with how ArnoldiFactorization is defined.

Any thoughts on how to overcome this? I don't really want to define a new ArnoldiFactorizationIncludingV₀ factorization that handles this case appropriately, but maybe that is the best option to maintain backwards compatibility?

[Jutho]

It is, but the initialization is somewhat involved to account for changing element types:

function initialize(iter::ArnoldiIterator; verbosity::Int = KrylovDefaults.verbosity[])
    # initialize without using eltype
    x₀ = iter.x₀
    β₀ = norm(x₀)
    iszero(β₀) && throw(ArgumentError("initial vector should not have norm zero"))
    Ax₀ = apply(iter.operator, x₀)
    α = inner(x₀, Ax₀) / (β₀ * β₀)
    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
    v = add!!(scale(Ax₀, zero(α)), x₀, 1 / β₀)

Note that the last line is really saying: v is the result of adding x₀ (which is v₀ in some parts of the code, that is a bit inconsistent) with a factor 1 / β₀ (where β₀ the norm of x₀) to a vector corresponding to A * x₀ scaled by a factor zero, i.e. adding x₀ / norm(x₀) to a zero vector. Hence, v really is just the normalized version of v₀, but with potentially a different type.

The reason for this is the following: suppose a user enters an initial vector of integers. Then normalization will have to change this to floating point numbers. But then maybe the linear map A may for example contain complex numbers, and all the other vectors in the Krylov subspace might have to be complex. Since I want all vectors to have a stable type, this code tries to figure all of this out in the initialization, by already applying A to x₀ once, and then further multiplying this with a zero scalar of the same type as the result of taking an inner product and dividing it by a norm.

So as a worked out example: if x₀ is a Vector{Int}, and A amounts to just multiplying with im, then Ax₀ would be of type Vector{Complex{Int}}, inner(x₀, Ax₀) would of type Complex{Int}, but β₀ would be of type Float64, α would be of type Complex{Float64}, and scale(Ax₀, zero(α)) would be of type Vector{Complex{Float64}}, which is what v ends up being. But ultimately, the value of v is equivalent to a normalized version of v₀, and that is the first vector that is added to the Krylov subspace.

The result of Ax₀ is furthermore not discarded, but used to compute the first residual (which will end up being the second Krylov vector, after normalization, at least in the Arnoldi routine).

[simeonschaub]

Ah, I misunderstood what the code was doing. Thanks for the detailed explanation! I have put up #151 as a starting point for discussion, curious to hear what you think