# 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:** 12\
**Page:** 2

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [June 5, 2024, 2:25pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/22 "2024-06-05T14:25:28Z")

</div>

> [@bremez](#):
>
> @juthohaegeman Any thoughts on a way to pass on the set of eigenvectors of M\_iMiM\_i as the initial guesses to `eigensolve` of M\_{i+1}Mi+1M\_{i+1}?

Haven’t read the whole thread, but LOBPCG does this. There’s an implementation in iterativesolvers.jl, and we also have one in [DFTK.jl/src/eigen/lobpcg\_hyper\_impl.jl at master · JuliaMolSim/DFTK.jl · GitHub](https://github.com/JuliaMolSim/DFTK.jl/blob/master/src/eigen/lobpcg_hyper_impl.jl)

---

<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, 2:59pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/23 "2024-06-05T14:59:35Z")

</div>

I’ll have a look! Perhaps unsurprisingly, my use case is also about computing bandstructures 😉

---

<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 11, 2024, 12:57pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/24 "2024-06-11T12:57:33Z")

</div>

@antoine-levitt I’ve had a chance to try `IterativeSolvers.lobpcg!`. Without reusing vectors, it is 2 to 3x slower than `KrylovKit` and `ArnoldiMethod`. With reusing the vectors from the previous k-point, it roughly halves the number of iterations needed to converge (~15 instead of ~30), making it comparable, but typically somewhat slower.

I thought that perhaps instead of just reusing vectors (feeding V\_{n-1}, the matrix of eigenvectors of the previous k-point) I instead try to extrapolate to the new k-point, i.e. feed into `lobpcg` as initial guess 2 V\_{n-1} - V\_{n-2}, would lead to better convergence. However, this actually hurt convergence and returned me to the ~30 iterations as with not reusing vectors at all. Any other ideas how to accelerate that convergence?

---

<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 11, 2024, 1:37pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/25 "2024-06-11T13:37:20Z")

</div>

> [@bremez](#):
>
> @juthohaegeman Any thoughts on a way to pass on the set of eigenvectors of M\_iMiM\_i as the initial guesses to `eigensolve` of M\_{i+1}Mi+1M\_{i+1}?

Since KrylovKit.jl does not have block implementations, you can only pass one starting vector. You could add a linear combination of the previous eigenvectors. If you would do this without changing the matrix, the individual eigenvectors would be reconstructed from the linear combination in a number of function applications equal to that of the number of vectors involved.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [June 11, 2024, 2:43pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/26 "2024-06-11T14:43:26Z")

</div>

2Vn-1 - Vn-2 is not great as it would not be able to accomodate band crossings and the like. If you want to use the information of Vn-2, a better solution is to do a rayleigh-ritz method (diagonalize the operator in the basis formed of the union of Vn-1 and Vn-2) if you can afford it, although I’m not sure how much iterations it will buy you. Or even pass to LOBPCG more vectors than you actually need, and setting the stopping criterion to only look at the ones you actually need (I don’t know about IterativeSolvers but our solver has an option for this)

---

<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 11, 2024, 5:59pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/27 "2024-06-11T17:59:53Z")

</div>

> [@bremez](#):
>
> I thought that perhaps instead of just reusing vectors (feeding V\_{n-1}, the matrix of eigenvectors of the previous k-point) I instead try to extrapolate to the new k-point

This is called a [numerical continuation](https://en.wikipedia.org/wiki/Numerical_continuation) problem, and in related problems (tracking [nonlinear eigenvalues for lasing modes](http://doi.org/10.1103/PhysRevA.90.023816)) we’ve had good results with [BifurcationKit.jl](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/) (which can do higher-order extrapolation and handle certain kinds of singularities); I haven’t tried it for band structures, though.

(But you still have to be careful of band crossings, i.e. where you are computing N bands but the N+1-th band crosses the N-th band, as @antoine-levitt discusses above.)

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [June 11, 2024, 6:21pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/28 "2024-06-11T18:21:53Z")

</div>

Presumably for band structures it’s a complex operator with arbitrary phase on the eigenvectors, so if you want to extrapolate you’d need at least to align the vectors before, ie find the best unitary R such that Vn-2 R is closest to Vn-1 (which you do by an svd of the overlap between Vn-2 and Vn-1)

---

<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 11, 2024, 6:31pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/29 "2024-06-11T18:31:29Z")

</div>

> [@antoine-levitt](#):
>
> if you want to extrapolate you’d need at least to align the vectors before, ie find the best unitary R such that Vn-2 R is closest to Vn-1

Note that for LOBPCG you should only need to extrapolate the subspace, i.e. you should only need to extrapolate to a set of vectors that nearly spans the new eigenvectors.

But for extrapolating from multiple k points I agree that you need to be careful about choosing the bases. One possibility is described in [https://onlinelibrary.wiley.com/doi/abs/10.1002/pssb.202000260](https://onlinelibrary.wiley.com/doi/abs/10.1002/pssb.202000260)

Also [https://doi.org/10.1021/acs.jctc.0c00684](https://doi.org/10.1021/acs.jctc.0c00684)

---

<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 12, 2024, 1:28pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/30 "2024-06-12T13:28:28Z")

</div>

@antoine-levitt @stevengj Fair, I did not consider that the phase might be an issue. For the low order extrapolation I was proposing, would simply fixing the phase of each basis vector, i.e. enforce `Vn[1,:]` are real, solve most of the problem? (Admittedly, this does not deal with the arbitrary superposition of degenerate eigenvectors. I am dealing with a low-symmetry case, so I expect at most doublets. Perhaps @stevengj 's comment on needing only extrapolated superposition applies.) Otherwise, could you be a bit more explicit about the SVD approach you outlined?

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [June 13, 2024, 1:49pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/31 "2024-06-13T13:49:00Z")

</div>

> [@bremez](#):
>
> For the low order extrapolation I was proposing, would simply fixing the phase of each basis vector, i.e. enforce `Vn[1,:]` are real, solve most of the problem?

Probably, as long as the first component doesn’t vanish. A better way is to align them, ie find the phase phi such that vn-1 phi is closest to vn-2. This comes out as vn-1 becomes vn-2 \<vn-2, vn-1\>/|\<vn-2, vn-1\>| (or something like this).

In the general multivector case, phi becomes a unitary matrix which can be obtained from an SVD of the overlap matrix between Vn-1 and Vn-2. I’m sure it’s very standard but I don’t know offhand of a place where it’s written down properly…

---

<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 14, 2024, 1:31pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/32 "2024-06-14T13:31:21Z")

</div>

[Would this be it](https://en.wikipedia.org/wiki/Singular_value_decomposition#Nearest_orthogonal_matrix)?

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [June 14, 2024, 2:08pm UTC](https://discourse.julialang.org/t/accelerating-repeated-eigenvalue-searches/114989/33 "2024-06-14T14:08:23Z")

</div>

Yep!

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