# Finding \*all\* the eigenvalues of a sparse matrix?

**URL:** <https://discourse.julialang.org/t/finding-all-the-eigenvalues-of-a-sparse-matrix/116895>\
**Category:** Numerics\
**Tags:** linearalgebra, sparse, eigenvalues, matrix\
**Created:** [July 11, 2024, 9:50am UTC](https://discourse.julialang.org/t/finding-all-the-eigenvalues-of-a-sparse-matrix/116895 "2024-07-11T09:50:48Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![rafaelbailo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafaelbailo/32/215573_2.png) [@rafaelbailo](https://discourse.julialang.org/u/rafaelbailo)\
**Post date:** [July 11, 2024, 9:50am UTC](https://discourse.julialang.org/t/finding-all-the-eigenvalues-of-a-sparse-matrix/116895/1 "2024-07-11T09:50:48Z")

</div>

Seems that LinearAlgebra.jl does not work with sparse matrices, and the libraries that do (e.g. Arpack.jl, ArnoldiMethod.jl, KrylovKit.jl) all use iterative methods, which are not efficient to find _all_ the eigenvalues. Short of using LinearAlgebra.jl with a regular matrix, do I have an alternative?

---

<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:** [July 11, 2024, 10:24am UTC](https://discourse.julialang.org/t/finding-all-the-eigenvalues-of-a-sparse-matrix/116895/2 "2024-07-11T10:24:51Z")

</div>

in general, there do not exist efficient algorithms for all the eigenvalues of a spare matrix.

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [July 11, 2024, 3:03pm UTC](https://discourse.julialang.org/t/finding-all-the-eigenvalues-of-a-sparse-matrix/116895/3 "2024-07-11T15:03:21Z")

</div>

To elaborate a bit on @Oscar_Smith’s answer: The main reason why there is in general in no efficient algorithm is because the matrix of eigenvectors is generally not sparse even if your original matrix is. Thus if you compute all eigenvectors you are computing a dense matrix anyways and there is no benefit of starting with a sparse one.

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [July 11, 2024, 5:37pm UTC](https://discourse.julialang.org/t/finding-all-the-eigenvalues-of-a-sparse-matrix/116895/4 "2024-07-11T17:37:33Z")

</div>

For a small sparse matrix (1000x1000 or so), convert to dense. Otherwise, ask yourself “why do I think I need ALL eigenvalues?”.

---

<div class="post-metadata">

**Author:** ![Thomas\_Rey](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/thomas_rey/32/210591_2.png) [@Thomas\_Rey](https://discourse.julialang.org/u/Thomas_Rey)\
**Post date:** [July 11, 2024, 8:17pm UTC](https://discourse.julialang.org/t/finding-all-the-eigenvalues-of-a-sparse-matrix/116895/5 "2024-07-11T20:17:36Z")

</div>

This post was temporarily hidden by the community for possibly being off-topic, unfocused, inappropriate, or spammy.

---

<div class="post-metadata">

**Author:** ![rafaelbailo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafaelbailo/32/215573_2.png) [@rafaelbailo](https://discourse.julialang.org/u/rafaelbailo)\
**Post date:** [July 11, 2024, 8:50pm UTC](https://discourse.julialang.org/t/finding-all-the-eigenvalues-of-a-sparse-matrix/116895/6 "2024-07-11T20:50:29Z")

</div>

Definitely. @PetrKryslUCSD, I could do without rhetorical questions about my thought process. There are problems in mathematics (kinetic equations, asymptotics, multi-scale problems among others) where knowing only _some_ eigenvalues of a very large matrix is not good enough. Hence the very specific title to my very specific question!

---

<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:** [July 11, 2024, 9:18pm UTC](https://discourse.julialang.org/t/finding-all-the-eigenvalues-of-a-sparse-matrix/116895/8 "2024-07-11T21:18:04Z")

</div>

I don’t think it’s necessary to get rancorous. As others have said, the general answer to your question is that it currently seems impossible (or at least impractical) to avoid dense methods for general sparse matrices if you want precise values for all the eigenvalues. (Unless your sparsity pattern is very special?)

For many problems in applied mathematics that at first seem to require the whole spectrum, people have worked hard to find sparse-compatible workarounds that compute a final result in a different way while avoiding computing the explicit spectrum as an intermediate step. Maybe you can explore this route? If you want to share more about your ultimate goal —what are you using the spectrum to accomplish?—then maybe someone here will have an idea.

---

<div class="post-metadata">

**Author:** ![kbarros](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kbarros/32/7723_2.png) [@kbarros](https://discourse.julialang.org/u/kbarros)\
**Post date:** [July 11, 2024, 9:51pm UTC](https://discourse.julialang.org/t/finding-all-the-eigenvalues-of-a-sparse-matrix/116895/9 "2024-07-11T21:51:42Z")

</div>

> For many problems in applied mathematics that at first seem to require the whole spectrum, people have worked hard to find sparse-compatible workarounds that compute a final result in a different way while avoiding computing the explicit spectrum as an intermediate step.

This is a good point. For some problems, it is enough to approximate the _density_ of the eigenvalue spectrum. Or, perhaps, to approximate the trace of some function of the given sparse matrix. For such problems, there are powerful stochastic approximation techniques that avoid explicit computation of the individual eigenvectors. From the physics perspective, see, for example, [the kernel polynomial method](https://arxiv.org/abs/cond-mat/0504627) and subsequent [“probing” techniques](https://onlinelibrary.wiley.com/doi/full/10.1002/nla.779) to accelerate stochastic convergence.

---

<div class="post-metadata">

**Author:** ![GregVernon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gregvernon/32/5611_2.png) [@GregVernon](https://discourse.julialang.org/u/GregVernon)\
**Post date:** [July 11, 2024, 9:53pm UTC](https://discourse.julialang.org/t/finding-all-the-eigenvalues-of-a-sparse-matrix/116895/10 "2024-07-11T21:53:27Z")

</div>

In a grad-school class based on [Trefethen’s text](https://epubs.siam.org/doi/10.1137/1.9781611977165) I also came across a need to compute all eigenvalues/eigenvectors of sparse matrices. I ended up just using my own hand-written Arnoldi method, but in my literature review I came across a figure:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/1/f/1fec647f253f1ea37e6f0bc0626ed609955ce248.png)

(Google Image search finds [this presentation](https://blogs.fau.de/essex/files/2012/11/ESSEX-SPPEXA-EuroPar2014.pdf))

Which would suggest Kernel Polynomial Methods based (unsurprisingly) on Chebyshev polynomials. I never investigated this further as I only needed a one-off approach and my own naive code was good enough.

@kbarros – jinx!

---

<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:** [July 11, 2024, 10:17pm UTC](https://discourse.julialang.org/t/finding-all-the-eigenvalues-of-a-sparse-matrix/116895/11 "2024-07-11T22:17:58Z")

</div>

> [@GregVernon](#):
>
> I also came across a need to compute all eigenvalues/eigenvectors of sparse matrices. I ended up just using my own hand-written Arnoldi method

If you run Arnoldi long enough to compute all of the eigenvalues (i.e. m steps without restarting for an m \times m matrix), isn’t it essentially equivalent in cost (or worse) to just using a dense eigensolver from the beginning?

Or did you somehow partition your spectrum into chunks and run Arnoldi with shift-and-invert to separately compute eigenvalues a few at a time in each chunk?

> [@GregVernon](#):
>
> Which would suggest Kernel Polynomial Methods based (unsurprisingly) on Chebyshev polynomials.

Don’t those compute continuous approximations for the spectral density / DOS (as well as other functions of the spectrum), rather than individual eigenvalues _per se_?

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [July 11, 2024, 11:05pm UTC](https://discourse.julialang.org/t/finding-all-the-eigenvalues-of-a-sparse-matrix/116895/12 "2024-07-11T23:05:54Z")

</div>

Do tell: if some eigenvalues of a VERY large matrix (what is very large in your neck of the woods?) are not enough for some people, what do they do?

---

<div class="post-metadata">

**Author:** ![GregVernon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gregvernon/32/5611_2.png) [@GregVernon](https://discourse.julialang.org/u/GregVernon)\
**Post date:** [July 11, 2024, 11:46pm UTC](https://discourse.julialang.org/t/finding-all-the-eigenvalues-of-a-sparse-matrix/116895/13 "2024-07-11T23:46:47Z")

</div>

IIRC at the time there was some shaking going on in the Julia ecosystem regarding eigenvalue solvers? In any case, the class was numerical linear algebra and I’d just written an Arnoldi solver in Matlab for an assignment (which I then used to make [this animation on the Arnoldi Iteration](https://en.wikipedia.org/wiki/Arnoldi_iteration#:~:text=is%20more%20complete.-,Arnoldi%20iteration,-demonstrating%20convergence%20of) page in Wikipedia) I thought it would be fun to rewrite in Julia.

![undefined](https://global.discourse-cdn.com/julialang/original/3X/0/8/08f61380e73a1f6a97261a5f728263464cdcfca9.gif)

---

<div class="post-metadata">

**Author:** ![GregVernon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gregvernon/32/5611_2.png) [@GregVernon](https://discourse.julialang.org/u/GregVernon)\
**Post date:** [July 11, 2024, 11:54pm UTC](https://discourse.julialang.org/t/finding-all-the-eigenvalues-of-a-sparse-matrix/116895/14 "2024-07-11T23:54:19Z")

</div>

> [@stevengj](#):
>
> Don’t those compute continuous approximations for the spectral density / DOS (as well as other functions of the spectrum), rather than individual eigenvalues _per se_?

Admittedly I don’t know much about the method(s), but found the graphic at least interesting on the surface - I’d never before needed to calculate all eigenvectors of a (large) matrix before and didn’t even know about [FEAST](https://www.feast-solver.org/index.htm) (or [PFEAST](https://vkalantzi.github.io/papers/pfeast2.pdf)) until I saw the graphic.

---

<div class="post-metadata">

**Author:** ![Ralph\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ralph_smith/32/10344_2.png) [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)\
**Post date:** [July 12, 2024, 2:38am UTC](https://discourse.julialang.org/t/finding-all-the-eigenvalues-of-a-sparse-matrix/116895/15 "2024-07-12T02:38:33Z")

</div>

FEAST and others build efficient solvers for large numbers of eigenpairs by separating the spectrum into pieces and distributing the pieces over compute nodes. Although they take more operations than the standard (also iterative!) dense-matrix methods, they would win out for the full spectrum of big enough problems because of the cleaner parallelism and memory use. I’m not aware of any implementations or wrappers of these in Julia as of now.

Some such approaches use Chebychev series of matvec products or solves to emphasize segments of the real line (for Hermitian matrices), a significant cousin of the shift-and-invert strategy. This is not the same as KPM, which “sets aside the requirement for a complete and exact knowledge of the spectrum” in favor of an approximate density (which is a preferable goal in some applications).

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [July 15, 2024, 6:53pm UTC](https://discourse.julialang.org/t/finding-all-the-eigenvalues-of-a-sparse-matrix/116895/16 "2024-07-15T18:53:49Z")

</div>

6 posts were split to a new topic: [Meta discussion about _all_ the eigenvalues](https://discourse.julialang.org/t/meta-discussion-about-all-the-eigenvalues/117064)

---

<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:** [July 12, 2024, 11:55am UTC](https://discourse.julialang.org/t/finding-all-the-eigenvalues-of-a-sparse-matrix/116895/19 "2024-07-12T11:55:57Z")

</div>

> [@Ralph\_Smith](#):
>
> Although they take more operations than the standard (also iterative!) dense-matrix methods, they would win out for the full spectrum of big enough problems because of the cleaner parallelism and memory use. I’m not aware of any implementations or wrappers of these in Julia as of now

I’m skeptical of this claim, do you have any references for it? It’s usually possible to get very good parallel scaling for dense \Theta(m^3) methods on m \times m matrices, because doing \Theta(m^3) work on \Theta(m^2) data makes it possible to hide the communications latency behind computation. (This is the same reason as why BLAS3 algorithms on a _single_ CPU, properly implemented, can [hide the memory latency](https://github.com/mitmath/18335/blob/spring21/notes/Memory-and-Matrices-latest.ipynb) and hit nearly the peak flop rate even for matrices that exceed the cache size.)

For example, here is a graph from [Petschow et al. (2013)](https://epubs.siam.org/doi/abs/10.1137/110848803) ([arXiv](https://ar5iv.labs.arxiv.org/html/1205.2107)) showing that [Elemental’s](https://github.com/LLNL/Elemental) dense Hermitian eigensolver (accessible in Julia via [Elemental.jl](https://github.com/JuliaParallel/Elemental.jl?tab=readme-ov-file)) attains nearly linear scaling up to 2048 cores on a fixed 20,000 \times 20,000 matrix (which is only 3GB in double precision). An even bigger matrix should scale even better (because the computational work scales faster than the communication).  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/1/9/19bde5e19ae69623cb4190d873d94496b3f1d622.png)  
See also [Davidovic (2021)](https://doi.org/10.23919/MIPRO52101.2021.9596900) for a review of more parallel dense eigensolvers.

So, if using FEAST for _all_ eigenvalues takes more computation (more arithmetic operations) than a dense eigensolver, which seems plausible(?), and both FEAST and the dense eigensolver get linear speedup from parallelism, then it seems like FEAST should still lose to a dense eigensolver for computing all eigenvalues.

The downside of good parallel scaling of \Theta(m^3) algorithms is that, even if you give me 1000x more processors, I can only solve a 10x larger matrix in the same amount of time. And if this comes from a 3d physical problem where m \sim \mathrm{diameter}^3, I can therefore only increase the diameter of my computational volume by \sqrt[3]{10} \approx 2.15 times with 1000x more CPUs. So dense-matrix parallelism doesn’t accomplish what you usually _really_ want, which is to enable you to attack much bigger problems. That’s why it’s still so important to find clever approaches that avoid computing the entire spectrum explicitly if you can.

---

<div class="post-metadata">

**Author:** ![Bruno\_Amorim](https://avatars.discourse-cdn.com/v4/letter/b/f6c823/32.png) [@Bruno\_Amorim](https://discourse.julialang.org/u/Bruno_Amorim)\
**Post date:** [July 12, 2024, 6:42pm UTC](https://discourse.julialang.org/t/finding-all-the-eigenvalues-of-a-sparse-matrix/116895/20 "2024-07-12T18:42:59Z")

</div>

Another possible option is the Eigenvalues Slicing Library: [[1802.05215] The Eigenvalues Slicing Library (EVSL): Algorithms, Implementation, and Software](https://arxiv.org/abs/1802.05215)  
(No wrapper in Julia that I am aware)

---

<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:** [July 12, 2024, 7:02pm UTC](https://discourse.julialang.org/t/finding-all-the-eigenvalues-of-a-sparse-matrix/116895/21 "2024-07-12T19:02:29Z")

</div>

> [@Bruno\_Amorim](#):
>
> Another possible option is the Eigenvalues Slicing Library: [[1802.05215] The Eigenvalues Slicing Library (EVSL): Algorithms, Implementation, and Software](https://arxiv.org/abs/1802.05215)

That sounds like a nice approach, but as they say in the concluding remarks, it is still apparently a research problem to make a scalable parallel implementation of this kind of algorithm (“A fully parallel version of EVSL using MPI is currently being developed.”). (The examples in the paper where they computed _all_ the eigenvalues appear to have been only for relatively small matrices, e.g. 2000x2000, where probably a dense method would be faster; in the larger examples of the paper they computed at most a few hundred eigenvalues.)

Unfortunately, if you look at the [github EVSL repository](https://github.com/eigs/EVSL), it doesn’t seem like it’s seen major development since the paper came out (2019), and in particular there doesn’t seem to have been any progress on a parallel version.

---

<div class="post-metadata">

**Author:** ![Ralph\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ralph_smith/32/10344_2.png) [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)\
**Post date:** [July 12, 2024, 9:19pm UTC](https://discourse.julialang.org/t/finding-all-the-eigenvalues-of-a-sparse-matrix/116895/22 "2024-07-12T21:19:13Z")

</div>

Petschow’s paper concerns the tridiagonal problem, not the full dense case. The reduction to tridiagonal (or Hessenberg) form was traditionally dominated by level 2 BLAS and thus memory/communication bound, so was the bottleneck that inspired my claim. TIL that the [ELPA](https://elpa.mpcdf.mpg.de/index.html) and [SLATE](https://icl.utk.edu/slate/) teams have built blocking schemes which radically improve this for the Hermitian case, so you are correct there. I haven’t found similar results for the nonsymmetric (Hessenberg) case, but the SLATE team seems to be working on it.

---

<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:** [July 14, 2024, 12:03pm UTC](https://discourse.julialang.org/t/finding-all-the-eigenvalues-of-a-sparse-matrix/116895/23 "2024-07-14T12:03:51Z")

</div>

> [@Ralph\_Smith](#):
>
> The reduction to tridiagonal (or Hessenberg) form was traditionally dominated by level 2 BLAS and thus memory/communication bound, so was the bottleneck that inspired my claim.

Petschow’s paper explicitly looked at this for ScaLAPACK in section 2.1 of their paper, and contradicts your claim — they found that the reduction to tridiagonal form was less than 10% of the compute time and that the clear parallelization bottleneck was the tridiagonal eigensolver. That’s why they focused in the paper on improving this stage of the computation.

[Next page](https://discourse.julialang.org/t/finding-all-the-eigenvalues-of-a-sparse-matrix/116895.md?page=2)