Option to diagonalize after every new Krylov vector in eigsolve

[antoine-levitt]

eigsolve(epsfun, length(scfres.ρ.real), 1, :SR, verbosity=3, tol=1e-4)

[ Info: Arnoldi iteration step 1: normres = 0.016207911009266315
[ Info: Arnoldi iteration step 2: normres = 0.5590671934901157
[ Info: Arnoldi iteration step 3: normres = 0.905377226214006
[ Info: Arnoldi iteration step 4: normres = 0.23298405811493555
[ Info: Arnoldi iteration step 5: normres = 0.15201652536133126
[ Info: Arnoldi iteration step 6: normres = 0.592603021837977
[ Info: Arnoldi iteration step 7: normres = 0.12176731843540282
[ Info: Arnoldi iteration step 8: normres = 0.09097000303823394
[ Info: Arnoldi iteration step 9: normres = 0.04323381192254804
[ Info: Arnoldi iteration step 10: normres = 0.07849715070693251
[ Info: Arnoldi iteration step 11: normres = 0.129007401971026
[ Info: Arnoldi iteration step 12: normres = 0.7133902458129476
[ Info: Arnoldi iteration step 13: normres = 0.05367711291850232
[ Info: Arnoldi iteration step 14: normres = 0.03533800454486114
[ Info: Arnoldi iteration step 15: normres = 0.7191209788317857
[ Info: Arnoldi iteration step 16: normres = 0.08411095972078231
[ Info: Arnoldi iteration step 17: normres = 0.06905579200587413
[ Info: Arnoldi iteration step 18: normres = 0.03099370046380324
[ Info: Arnoldi iteration step 19: normres = 0.17225062546432904
[ Info: Arnoldi iteration step 20: normres = 0.28845620299403657
[ Info: Arnoldi iteration step 21: normres = 0.15783635449059297
[ Info: Arnoldi iteration step 22: normres = 0.03996345360229037
[ Info: Arnoldi iteration step 23: normres = 0.025087829332985617
[ Info: Arnoldi iteration step 24: normres = 0.020227416987651632
[ Info: Arnoldi iteration step 25: normres = 0.12157529694425907
[ Info: Arnoldi iteration step 26: normres = 0.02732767047408143
[ Info: Arnoldi iteration step 27: normres = 0.08297860675872647
[ Info: Arnoldi iteration step 28: normres = 0.04088602000015171
[ Info: Arnoldi iteration step 29: normres = 0.015330135139814571
[ Info: Arnoldi iteration step 30: normres = 0.029899360771800297
[ Info: Arnoldi schursolve in iter 1: 1 values converged, normres = (1.34e-05)
┌ Info: Arnoldi eigsolve finished after 1 iterations:
│ *  1 eigenvalues converged
│ *  norm of residuals = (1.3359144756588437e-5,)
└ *  number of operations = 30

The normres printouts here are strange: they can't be residuals, right?

[Jutho]

normres is the residual of the final vector of the Krylov subspace, i.e. in the Krylov relation

A V_m = V_m H_m + r em' = V_m Hm + beta * v{m+1} * em'

where A is the linear map, V_m the basis of the Krylov subspace of length m, H_m the reduced matrix V_m' A Vm, r = beta v{m+1} the residual, which has norm beta and its direction becomes the next Krylov vector (and em the vector with all zeros except a 1).

The value being printed is beta. However, it's not because beta is not going to zero, that the subspace itself cannot contain an approximate eigenvector whose residual is going to zero. I agree that it would probably be more useful to print that information. I will have to look into it.

[antoine-levitt]

Alright, that makes sense, thanks!

[Jutho]

Now I see, you have verbosity 3, so it also prints the beta after every time A is applied, i.e. in every expansion step of the Krylov subspace. At that point, the information about the residual of the eigenvector is not available, since that computation is only done after the Krylov subspace has been built up to the specified krylovdim. At that time, this information is printed:

[ Info: Arnoldi schursolve in iter 1: 1 values converged, normres = (1.34e-05)

Here, however, this is only printed once, because it is immediately converged.

[Jutho]

With lower verbosity, but in a calculation that needs several restarts before convergence, you would only see the information message about the residual of the eigenvector.

The extra information that you get by verbosity > 2 is completely controlled by the routine that builds the Krylov subspace (independent of whether it is an eigenvalue problem, linear problem, ...) and so that can only print the value of beta at every step.

[Jutho]

I will close this; I hope you agree that the printing is then fine?

[antoine-levitt]

I think it's fine for the intended use of krylovkit. My use case is probably a bit different than most, in that each application of A is extremely expensive (seconds to hours) but the linalg is cheap. It's then useful to get a printout as often as possible. In fact often the eigensolver converges in less than krylovdim iterations, but it doesn't realize this because it waits until the subspace is built to diagobalize the matrix. Is there any way to accommodate this use case in krylovkit? Otherwise I'll probably just build my own Arnoldi routine, diagonalizing at each step.

[Jutho]

This is currently not supported, but also in the application domain in which I am active this is sometimes useful. So maybe this should become an option yes.

[antoine-levitt]

Also while we're at wishlist stuff, it would be great to get the low-rank approximation generated by the Arnoldi process.

[Jutho]

Does schursolve do what you want? Or really just the arnoldi proces itself?

[antoine-levitt]

Huh, I looked into the docs and found I can just iterate on the arnoldi directly, which is the easiest way to do what I want (and more). That's completely awesome, thanks for making that possible! Btw, the doc has a typo over at https://jutho.github.io/KrylovKit.jl/latest/man/implementation/: for V,B,r,nr,b in ArnoldiIterator(f, v₀) should be for (V,B,r,nr,b) in ArnoldiIterator(f, v₀)

[Jutho]

Yes the KrylovIterators was going to be my next suggestion. Thanks.

[antoine-levitt]

if anybody's interested, here's what I cooked up for my use case (getting the extremal eigenvalues):

using KrylovKit

N = 100
A = randn(N, N)
# A = A+A'
x0 = randn(N)

for (V, B, r, nr, b) in ArnoldiIterator(A, x0)
    # A * V = V * B + r * b'.
    V = hcat(V...)
    AV = V*B + r*b'

    ew, ev = eigen(B)
    Vr = V*ev
    AVr = AV*ev
    R = AVr - Vr * Diagonal(ew)

    N = size(V, 2)
    inds = 1:min(N, 4)
    inds = [inds..., (N .+ 1 .- inds)...]
    display([ew[inds] [norm(r) for r in eachcol(R[:, inds])]])
end

[Jutho]

I have added an eager keyword to try to use the Krylov subspace for the relevant task (eigsolve, svdsolve, expintegrator) after every expansion step. Available in v0.5.0 which will hopefully be registered later today.

[antoine-levitt]

You rock, thanks!