Work with eigenvalue/eigenvector solver?

[duetosymmetry]

Since this provides a variety of factorizations, I thought it might be possible to use this package to solve an infinite-dimensional eigenvalue/eigenvector problem. For example, can ArnoldiMethod.jl be made to use an infinite-dimensional factorization to find the closest eigenvalue/vector pair to some initial guess (via shift-and-invert)? I was not able to make it work out of the box, but I'm also not fluent in Julia.

[dlfivefifty]

What you mentioned is exactly the type of application we have in mind for the package. If you give a MWE I can try to get it working

[duetosymmetry]

The problem I'm interested in is the computation of quasinormal mode frequencies of spinning black holes... part of that problem involves computing so-called "spin-weighted spheroidal harmonics" as a spectral decomposition in the basis of "spin-weighted spherical harmonics". These are a generalization of both the spherical harmonics and the spheroidal harmonics (https://dlmf.nist.gov/30). I implemented this in python here, simply truncating the complex symmetric pentadiagonal matrix and using a canned eigenvalue/vector solve. The method is described in Sec. IID of this paper. In this gist I write a Julia function to generate the infinite matrix — note that this is parameterized over the arguments s, m (integers, though the problem could be defined for half-integers as well) and c (complex, could be Complex{BigFloat} for arbitrary precision studies).

In the paper linked above, the rows/columns are labeled by l which starts from l_min(s, m) = max(abs(s), abs(m)), taking values of l = l_min(s,m) + [0, 1, 2, 3...]. Then the eigenvalues should be near l(l+1)-s(s+1) for allowed values of l, and the corresponding eigenvectors near [0, 0, ... 0, 1, 0, ...].

[dlfivefifty]

There’s no way a MWE requires knowing about spin-weighted SH.

Usually for SH the operators decouple across Fourier modes so there is probably a simple example of a banded matrix. I’m thinking like a 2 line example, the first line makes the operator and the 2nd line plugs it in to arnoldi

[duetosymmetry]

Sorry, I didn't see that you wanted a simpler example (in my pre-coffee state, I thought it was just any example). The gist is pretty short, but I'll try to think of a simpler problem.

[duetosymmetry]

@dlfivefifty I created a more minimal example. Sorry that it's still not just 2 lines. It's just that the problems I'm familiar with where you might want an unknown eigenvalue of an infinite-dimensional matrix problem are this same type of problem. In this same gist, I added a jupyter notebook constructing the infinite matrix for a simpler problem, and getting a specific eigenvalue in 3 different ways: using builtin LinearAlgebra, using KrylovKit, and using ArnoldiMethod. But none of these can work with a ℵ₀×ℵ₀ matrix.

I think I can re-implement the Arnoldi method so it doesn't attempt to pre-allocate an infinite amount of memory... but (i) can any currently-written code be made to work? (ii) I assume there are details I don't know yet, like if the (lazy) basis vectors of the Krylov subspace need to be extended at each iteration, etc.

[dlfivefifty]

There is nothing in the gist that I can obviously copy-and-paste so its not very helpful

[duetosymmetry]

I think it just takes a few seconds to render the ipynb. I just converted it to a plain .jl file here. The fundamental issue is that the existing packages try to preallocate storage and fill it with 0's, but this just makes the runtime try to grab ever-increasing amounts of memory, filling it with 0's.

[duetosymmetry]

@dlfivefifty any suggestions for how I can contribute to make InfiniteLinearAlgebra.jl able to solve eigen-problems? What would be your strategy to make this happen?

[dlfivefifty]

You just need low order discrete eigenvalues? If so it would be easy to port over the ApproxFun code:

https://github.com/JuliaApproximation/ApproxFunBase.jl/blob/master/src/eigen.jl

[duetosymmetry]

I have no problem truncating the infinite matrix and then using other routines... but I thought the spirit for this package would be that the eigenvectors come back as an infinite array, with more components computed lazily when needed

[dlfivefifty]

That could be doable but assuming you only have discrete eigenvalues

Eigenvectors should be an ∞ matrix btw.

For continuous spectrum it would have to return a QuasiArray. Actually I had some work on this in SpectralMeasures.jl

[duetosymmetry]

Yes, I am thinking of the simpler case of a discrete spectrum. I agree that it ought to be an infinite matrix. I guess this means the number of eigenvectors should also be lazy, and more are computed as needed? This seems very tricky to me — do you have to be able to tell how many of your (Arnoldi) vectors are numerically converged to your desired tolerance?

[dlfivefifty]

As long as we don’t claim to guarantee an ordering it isn’t that hard I think: just do simealtanous inverse iteration with Rayleigh shifts. Since eigenvectors have decay we would compute the whole vector a la \