# Accelerating repeated eigenvalue searches

**URL:** <https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989>\
**Category:** Numerics\
**Tags:** linear-algebra, eigenvalues\
**Created:** [May 30, 2024, 10:47pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989 "2024-05-30T22:47:44Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![bremez](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bremez/32/38777_2.png) [@bremez](https://discourse.julialang.org/u/bremez)\
**Post date:** [May 30, 2024, 10:47pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/1 "2024-05-30T22:47:44Z")

</div>

I need to compute the eigenvalues of 2N\times2N (`N` ~ 1000) complex, hermitian, dense, matrices, the number of matrices being ~100. They have the following properties:

1. All matrices M\_i take the form M\_i = A + D\_i, where A is common to all, and D\_i are different tridiagonal matrices.
2. The matrices are positive-semidefinite, and I am interested only in eigenvalues in the range, say, [0, 2]. I know that the number of such eigenvalues is of order 20.

Is there a way to exploit the common structure of the matrices (e.g. an advantageous factorization of A), or the desired interval, or both, to speed up the calculation of these multiple eigenspectra?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 31, 2024, 5:35pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/2 "2024-05-31T17:35:47Z")

</div>

> [@bremez](#):
>
> The matrices are positive-semidefinite, and I am interested only in eigenvalues in the range, say, [0, 2]. I know that the number of such eigenvalues is of order 20

You can call `eigvals(Hermitian(Mᵢ), 0, 2)` to only compute eigenvalues in [0,2], for example.

I don’t know how to exploit an A + D\_i structure, however.

---

<div class="post-metadata">

**Author:** ![bremez](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bremez/32/38777_2.png) [@bremez](https://discourse.julialang.org/u/bremez)\
**Post date:** [June 1, 2024, 6:40pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/4 "2024-06-01T18:40:27Z")

</div>

Does this work substantially differently to `eigvals!(Hermitian(M), 1:20)`, which is what I have been doing so far? Asking for the ~20 lowest eigenvalues, I then just filter for those \le 2. However, my impression was that requesting a small number had the same runtime irrespective of the number (so `1:5`, etc., would take about as long), and comparable to requesting the full ~1000 eigenspectrum.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 1, 2024, 6:47pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/5 "2024-06-01T18:47:54Z")

</div>

> [@bremez](#):
>
> However, my impression was that requesting a small number had the same runtime irrespective of the number (so `1:5`, etc., would take about as long), and comparable to requesting the full ~1000 eigenspectrum.

Unfortunately, whether you ask for a few eigenvalues or all the eigenvalues, the first step in the dense eigensolver is to compute the tridiagonal/Hessenberg reduction of the full matrix, which is what takes the majority of the time.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [June 1, 2024, 7:16pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/6 "2024-06-01T19:16:31Z")

</div>

> [@bremez](#):
>
> Asking for the ~20 lowest eigenvalues,

This is the piece of information I was looking for. You can try iterative methods

1. shift invert should give you the smallest eigenvalues, you can even ask for the first 20. See for example `eigsolve` in `KrylovKit` or `ArnoldiMethod`
2. you can also try [Home · JacobiDavidson.jl](https://haampie.github.io/JacobiDavidson.jl/stable/) which is kindof a newton algo to find the eigenelements, I did not try that much in the past.

For this to work, you have to solve `Mi sol = rhs` fast. Now, I would say that solving `(I+A\D)sol = A\rhs` with iterative method could be fast depending on the properties of `A\D` (do not form this!). Even better, I would consider that giant matrix `diag(M1, ..., MN)` of size 200N~200000. Indeed, ( I\_{100}\otimes A)^{-1} is simple to find 😉 and so you could try shift invert on this big matrix

---

<div class="post-metadata">

**Author:** ![bremez](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bremez/32/38777_2.png) [@bremez](https://discourse.julialang.org/u/bremez)\
**Post date:** [June 3, 2024, 3:07pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/7 "2024-06-03T15:07:55Z")

</div>

> Now, I would say that solving `(I+A\D)sol = A\rhs` with iterative method could be fast depending on the properties of `A\D`

By “solving with iterative methods”, you mean you would expect that e.g. `eigsolve(M)` would converge rapidly? Otherwise, I’m not sure what in practice you are suggesting (if not the formation of `A\D`, which you warn against).

> [@rveltz](#):
>
> Even better, I would consider that giant matrix `diag(M1, ..., MN)` of size 200N~200000. Indeed, ( I\_{100}\otimes A)^{-1} is simple to find.

What do I benefit from considering this large matrix? I suppose it can be realized sparsely, but otherwise seems like extra compute work. I also want to know which eivenvalue to attribute to which M\_i. The inverse is as hard as inverting A, which I’m not sure is not singular.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [June 3, 2024, 6:37pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/8 "2024-06-03T18:37:36Z")

</div>

> you mean you would expect that e.g. `eigsolve(M)` would converge rapidly?

Yes but I think `eigsolve(x->M\x)` would converge faster. Also, for computing `x->M\x` I would use `linsolve` with `A` serving as a preconditioner. A basic example is in this [tutorial](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/tutorialsSH3d/#sh3d).

> What do I benefit from considering this large matrix?

That is the title of your post! To try use threading. Being block diagonal, the eigenvalues are easy to attribute to the relevant block

---

<div class="post-metadata">

**Author:** ![bremez](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bremez/32/38777_2.png) [@bremez](https://discourse.julialang.org/u/bremez)\
**Post date:** [June 3, 2024, 6:59pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/9 "2024-06-03T18:59:25Z")

</div>

Ah! Though I imagine ideally a better threading performance would be to just parallelize over the different M\_i, kept separate?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [June 3, 2024, 7:42pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/10 "2024-06-03T19:42:56Z")

</div>

yeah I would expect that to work better.

---

<div class="post-metadata">

**Author:** ![juthohaegeman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juthohaegeman/32/8620_2.png) [@juthohaegeman](https://discourse.julialang.org/u/juthohaegeman)\
**Post date:** [June 3, 2024, 8:10pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/11 "2024-06-03T20:10:46Z")

</div>

> [@rveltz](#):
>
> Yes but I think `eigsolve(x->M\x)` would converge faster. Also, for computing `x->M\x` I would use `linsolve` with `A` serving as a preconditioner. A basic example is in this [tutorial](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/tutorialsSH3d/#sh3d).

Do you think a shift and invert (without the shift, so just invert) strategy would lead to faster convergence here? If it is guaranteed that the matrix is positive definite, then all its eigenvalues are on the positive real axis, and the range [0,2] contains extremal eigenvalues (use `:SR` in `KrylovKit.eigsolve`).

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [June 4, 2024, 6:49am UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/12 "2024-06-04T06:49:27Z")

</div>

it is worth a try

---

<div class="post-metadata">

**Author:** ![bremez](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bremez/32/38777_2.png) [@bremez](https://discourse.julialang.org/u/bremez)\
**Post date:** [June 4, 2024, 4:47pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/13 "2024-06-04T16:47:54Z")

</div>

After some several tinkerings, some comments:

1. Perhaps I’m miss-applying the preconditioner, but from what I see, say, `cg(Mi, v; Pl = lu(A), log = true).isconverged` returns false (for a random vector `v`). Same for changing LU to QR factorization, and `cg` to `gmres`. Here I used `IterativeSolvers.jl`. Perhaps I should have mentioned this before, but I expect the spectrum of M\_i to be dominated by D\_i whereas A is a “small perturbation” - does this make it a bad preconditioner? I am not too familiar with how to asses/pick preconditioners. (Note also A is NOT positive definite, only the sums A + D\_i)
2. One alternative is to use, say, M\_1 as the preconditioner. This seems to at least converge. However, this is not much faster than `M\v`.
3. Using `KrylovKit.eigsolve(Mi,...)` is much much faster than `KrylovKit.eigsolve(v -> Mi\v, ...)`, and in light of the only modest gains of the linear solvers in the preceding point, I’m not sure this approach has much to gain. What was the original motivation to think this might be quicker?
4. `KrylovKit.eigsolve` does beats `LinearAlgebra.eigvals!` for ~20ish eigenvalues that I need. However, `KrylovKit` allocates a lot of memory as my matrix size increases (I might need to go to larger `N` than expected). Any way to reduce that?
5. Perhaps related to the previous item: Using [this example](https://docs.julialang.org/en/v1/manual/multi-threading/#Using-@threads-without-data-races) to multithread over the set M\_i, `nthreads=4`, (i) `LinearAlgebra.eigvals!` speeds up roughly 3x, whereas (ii) `KrylovKit.eigsolve` **is 2x slower** than single-threaded, coincidentally wiping out the advantage and giving comparable runtimes (but with the extra memory usage). Perhaps memory bandwidth issue? So from best to worse, single-thread Krylov → (multithreaded Krylov, multithreaded `eigvals`) → single-thread `eigvals`. Thoughts?
6. I should have been more precise: My D\_i are not just tridiagonal, they are in fact 2\times2 block-diagonal. Is that of any use?

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [June 4, 2024, 5:59pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/14 "2024-06-04T17:59:48Z")

</div>

> [@bremez](#):
>
> but I expect the spectrum of M\_i to be dominated by D\_i whereas A is a “small perturbation” - does this make it a bad preconditioner?

Basically, you want P^{-1}M\_i=I+K\_i where K is a compact operator.

Just to be sure, in your notations, M = diag(M\_1,\cdots, M\_n)?

---

<div class="post-metadata">

**Author:** ![bremez](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bremez/32/38777_2.png) [@bremez](https://discourse.julialang.org/u/bremez)\
**Post date:** [June 4, 2024, 6:46pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/15 "2024-06-04T18:46:04Z")

</div>

No, I had meant any individual M\_i, now corrected. Compact K just in terms of “bandwidth”/sparsity for fast multiplication? Any conditions on its spectrum, etc.?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 4, 2024, 8:10pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/16 "2024-06-04T20:10:43Z")

</div>

> [@rveltz](#):
>
> Yes but I think `eigsolve(x->M\x)` would converge faster. Also, for computing `x->M\x` I would use `linsolve` with `A` serving as a preconditioner.

It’s certainly worth trying, but I’m not too optimistic that a nested iterative method like this will be faster than the dense eigensolver for a 1000x1000 dense matrix, unless the `linsolve` (preconditioned CG) converges extremely quickly (which requires your tridiagonal matrices D\_i to not shift the spectrum of A + D\_i much compared to A, i.e. it could work well if \Vert D\_i \Vert \ll \Vert A \Vert).

> [@bremez](#):
>
> Perhaps I should have mentioned this before, but I expect the spectrum of M\_i to be dominated by D\_i whereas A is a “small perturbation” - does this make it a bad preconditioner?

Yes, in this case you want to do something completely different.

I guess you could try an iterative linear solve for M\_i \setminus x preconditioned by D\_i \setminus x, since tridiagonal linear solve is linear time, and if that converges sufficiently quickly try plugging it into `eigsolve`.

---

<div class="post-metadata">

**Author:** ![juthohaegeman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juthohaegeman/32/8620_2.png) [@juthohaegeman](https://discourse.julialang.org/u/juthohaegeman)\
**Post date:** [June 4, 2024, 10:18pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/17 "2024-06-04T22:18:44Z")

</div>

> [@bremez](#):
>
> Perhaps related to the previous item: Using [this example](https://docs.julialang.org/en/v1/manual/multi-threading/#Using-@threads-without-data-races) to multithread over the set M\_iMiM\_i, `nthreads=4`, (i) `LinearAlgebra.eigvals!` speeds up roughly 3x, whereas (ii) `KrylovKit.eigsolve` **is 2x slower** than single-threaded, coincidentally wiping out the advantage and giving comparable runtimes (but with the extra memory usage). Perhaps memory bandwidth issue? So from best to worse, single-thread Krylov → (multithreaded Krylov, multithreaded `eigvals`) → single-thread `eigvals`. Thoughts?

Not sure why this happens, but you can disable the multithreading that is within KrylovKit using `KrylovKit.disable_threads()`. I assume there is some conflict between the threading of BLAS in the matrix vector multiplication, and the threading in KrylovKit, although the latter is purely restricted to the “native KrylovKit” parts of the code where there are no BLAS calls or function calls to the linear operator going on.

---

<div class="post-metadata">

**Author:** ![bremez](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bremez/32/38777_2.png) [@bremez](https://discourse.julialang.org/u/bremez)\
**Post date:** [June 4, 2024, 11:23pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/18 "2024-06-04T23:23:47Z")

</div>

It seems like you suggested I can use (LU of) D\_i to precondiction CG. However, turns out this is not as rapidly convergent as preconditioning GMRES with M\_1 (surprising?) [GMRES fails to converge without preconditioning or with D\_i,A as preconditioners]. This is the fastest iterative realization I found, about 0.5 the time of `Mi \ v`. (Digression: for some reason `Hermitian(Mi) \ v` takes longer.)

However, I think this direction might be moot. I have noticed that using `eigsolve(inv(Mi), 20, :LR,...)` reports 6~7 iterations, and ~70 “linear map applications” of the matrix multiplication at my typical use cases. Now, in terms of runtime, `inv(Mi)` itself takes a while - but turns out, less than 70 runs of preconditioned GMRES. So I might as well just compute M\_i^{-1} explicitly. The subsequent runtime of `eigsolve` is about half as much as the inversion. The total is approximately 33% the time of the explicit eigensolver `LinearAlgebra.eigvals(Hermitian(Mi), 1:20)`.

Finally: is the inversion trick worth it? Yes! `eigsolve(Mi, 20, :SR,...)` takes ~50 iterations and 430 linear map applications. Since this is more than 3x70, that is more than the combined cost of explicit inversion and `eigsolve`ing the inverse. Thanks for all the comments, I sure learned a lot!

One final point: Another thing I forgot to mention is that the ~100 different D\_i are essentially a sampling of a continuous matrix D(x). Therefore, I expect that between M\_i and M\_{i+1} the eigenspectrum and eigenvectors to only change by a little. Is there a way to feed **all 20** eigenvectors as initial guesses to some iterative solver? It seems to me that at least `KrylovKit.eigsolve` only takes one initial guess vector. Perhaps looping over 20 calls to `ArnoldiMethod.partialschur!`? If this will dramatically reduce the number of iterations/linear maps for M\_{2,3,...}, perhaps the cost of the inversion can be avoided.

---

<div class="post-metadata">

**Author:** ![bremez](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bremez/32/38777_2.png) [@bremez](https://discourse.julialang.org/u/bremez)\
**Post date:** [June 4, 2024, 11:52pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/19 "2024-06-04T23:52:04Z")

</div>

Thanks for the tip. Even when I disable `KrylovKit`’s threads, sometimes single threaded is faster. However, it seems I regain the performance advantage of multithreading if the matrices are small enough, say below 400x400 (dense) instead of 1000x1000. As the matrix size increases this benefit vanishes and sometimes there is slight regression. Once I also `BLAS.set_num_threads(1)`, I finally see the performance separation between my own higher-level multithreading (its mostly my “single-thread” version that slows down). So perhaps the issue is how BLAS decide to multithread over matrices of size ~1000. Curiously, the benefit I get for multithreading (still with single-threaded BLAS) is closer to 2.0x~2.5x, instead of the 4x I’d expect with `Threads.nthreads==4`.

---

<div class="post-metadata">

**Author:** ![juthohaegeman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juthohaegeman/32/8620_2.png) [@juthohaegeman](https://discourse.julialang.org/u/juthohaegeman)\
**Post date:** [June 5, 2024, 8:51am UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/20 "2024-06-05T08:51:19Z")

</div>

What is “your own multithreading” in this case? You mean an outer threaded for loop over the ~100 matrices that you need to find the eigenvalues of? In that case I would indeed expect a speedup proportional to the number of threads, as the overhead of the threading should be small compared to the cost of the computation. However, it is indeed important then to disable threading at all lower levels, i.e. both in KrylovKit and in BLAS.

---

<div class="post-metadata">

**Author:** ![bremez](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bremez/32/38777_2.png) [@bremez](https://discourse.julialang.org/u/bremez)\
**Post date:** [June 5, 2024, 1:34pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/21 "2024-06-05T13:34:18Z")

</div>

My solution follows the Julia manual. Here `system` is a struct that defines the physical system I am solving, and contains the precalculated matrix A. `xvals` is the value of the parameter with which I sample the block-diagonals D\_i \in \lbrace D(x) | x \in \text{xvals}\rbrace (see my previous post). These diagonals are computed (very cheaply) on the fly inside `_computeSpectrum_singlecore`. `args` passes on my desired number of eigenvalues per value of x, as well as a flag for which solver to use -`LinearAlgebra` (direct), `KrylovKit`, or `ArnoldiMethod` (shift-invert and `:LR`).

Right now I disabled `KrylovKit`’s threading globally to simply comparisons, so treat it as off.

After manually inverting, it seems `KrylovKit` and `ArnoldiMethod` give similar runtime performance, with the latter using fewer allocations, so I am leaning towards making that the default.

```julia
function computeSpectrum(system, xvals; multicore = true, args...)
    # TODO: consider a heuristic when to multicore and when not to?
    if multicore
        blas_threads = BLAS.get_num_threads()
        BLAS.set_num_threads(1); # each instance should be single-threaded
        kchunks = Iterators.partition(xvals, length(xvals) ÷ Threads.nthreads())
        tasks = map(kchunks) do chunk
            Threads.@spawn _computeSpectrum_singlecore(system, chunk; args...)
        end
        spectra = fetch.(tasks)
        # TODO catch possible exceptions to guarantee this happens
        BLAS.set_num_threads(blas_threads)
        return vcat(spectra...)
    else
        return _computeSpectrum_singlecore(system, xval; args...)
    end
end

```

@juthohaegeman Any thoughts on a way to pass on the set of eigenvectors of M\_i as the initial guesses to `eigensolve` of M\_{i+1}?

[Next page](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989.md?page=2)