# Special method for finding eigenvalues of a large sparse symmetric matrix?

**URL:** <https://discourse.julialang.org/t/special-method-for-finding-eigenvalues-of-a-large-sparse-symmetric-matrix/9621>\
**Category:** Numerics\
**Created:** [March 9, 2018, 8:09pm UTC](https://discourse.julialang.org/t/special-method-for-finding-eigenvalues-of-a-large-sparse-symmetric-matrix/9621 "2018-03-09T20:09:05Z")\
**Posts on this page:** 18\
**Page:** 1

<div class="post-metadata">

**Author:** ![406700](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/406700/32/3705_2.png) [@406700](https://discourse.julialang.org/u/406700)\
**Post date:** [March 9, 2018, 8:09pm UTC](https://discourse.julialang.org/t/special-method-for-finding-eigenvalues-of-a-large-sparse-symmetric-matrix/9621/1 "2018-03-09T20:09:05Z")

</div>

Hi,  
I am looking to solve the eigenvalues of sparse symmetric matrix of ~ 14000x14000.  
Eigfact is too slow, and didn’t give results after an hour running. (8gb ram i7 processor)  
Are there special iteritive methods?  
One post mentions julia pro linking to MKL but I can’t imagine it would give the type of speedup I need.  
I’m new to Julia and programming in general so let me know if I’m missing something obvious.  
Thanks!

---

<div class="post-metadata">

**Author:** ![simonbyrne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simonbyrne/32/19_2.png) [@simonbyrne](https://discourse.julialang.org/u/simonbyrne)\
**Post date:** [March 9, 2018, 9:12pm UTC](https://discourse.julialang.org/t/special-method-for-finding-eigenvalues-of-a-large-sparse-symmetric-matrix/9621/2 "2018-03-09T21:12:14Z")

</div>

Try [`eigs`](https://docs.julialang.org/en/stable/stdlib/linalg/#Base.LinAlg.eigs-Tuple%7BAny%7D).

---

<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:** [March 9, 2018, 10:57pm UTC](https://discourse.julialang.org/t/special-method-for-finding-eigenvalues-of-a-large-sparse-symmetric-matrix/9621/3 "2018-03-09T22:57:28Z")

</div>

Do you really need _all_ the eigenvalues? If so, why? (There might be a better approach to your ultimate goal.)

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [March 9, 2018, 11:10pm UTC](https://discourse.julialang.org/t/special-method-for-finding-eigenvalues-of-a-large-sparse-symmetric-matrix/9621/4 "2018-03-09T23:10:14Z")

</div>

Running on N cores gives roughly an N-fold speedup with MKL, I believe.

---

<div class="post-metadata">

**Author:** ![406700](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/406700/32/3705_2.png) [@406700](https://discourse.julialang.org/u/406700)\
**Post date:** [March 10, 2018, 5:20pm UTC](https://discourse.julialang.org/t/special-method-for-finding-eigenvalues-of-a-large-sparse-symmetric-matrix/9621/5 "2018-03-10T17:20:48Z")

</div>

@ simonbyrne eigs seems to actually be a little slower.  
@ stevengj I don’t necessarily need all the eigenvalue but they have to be at leat the n smallest values. Is there a better function for this?  
@ dpsanders, I’ll try with MKL and report if there is a significant speedup.

Thank for the answers!

---

<div class="post-metadata">

**Author:** ![406700](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/406700/32/3705_2.png) [@406700](https://discourse.julialang.org/u/406700)\
**Post date:** [March 10, 2018, 5:20pm UTC](https://discourse.julialang.org/t/special-method-for-finding-eigenvalues-of-a-large-sparse-symmetric-matrix/9621/6 "2018-03-10T17:20:53Z")

</div>

I just noticed that running the MKL version htop is showing only 100% on 4 cores with low values on the others while without hkl I see 100% on all 8. The core i7 uses hyper-threading which is why I see 8 cores.  
Any inight on this?

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [March 10, 2018, 7:07pm UTC](https://discourse.julialang.org/t/special-method-for-finding-eigenvalues-of-a-large-sparse-symmetric-matrix/9621/7 "2018-03-10T19:07:41Z")

</div>

You could try simeaultanious inverse iteration with Rayleigh shifts.

---

<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:** [March 10, 2018, 7:44pm UTC](https://discourse.julialang.org/t/special-method-for-finding-eigenvalues-of-a-large-sparse-symmetric-matrix/9621/8 "2018-03-10T19:44:20Z")

</div>

> [@406700](#):
>
> I just noticed that running the MKL version htop is showing only 100% on 4 cores with low values on the others while without hkl I see 100% on all 8. The core i7 uses hyper-threading which is why I see 8 cores.
> 
> Any inight on this?

Linear algebra like this usually doesn’t benefit from hyperthreading; it’s better to run full out on whatever fits in cache. MKL knows this but for some reason OpenBLAS defaults to using all logical cores. Top shows 100% even if some cores are waiting for food.

For small eigenvalues, `eigs` has to factor a matrix at every iteration. If you are using a sparse matrix type, that is done on one core, and is expensive. 14000 is not too large, so you are better off converting to a full matrix, which is factored fast with BLAS.

Make sure your matrix is really symmetric (`A=(A+A')/2`, if need be). The general algorithms are much more expensive. Even `eigfact` should finish in less than 30 min for a dense symmetric matrix of rank 14000.

---

<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:** [March 10, 2018, 8:28pm UTC](https://discourse.julialang.org/t/special-method-for-finding-eigenvalues-of-a-large-sparse-symmetric-matrix/9621/9 "2018-03-10T20:28:49Z")

</div>

> [@406700](#):
>
> I don’t necessarily need all the eigenvalue but they have to be at leat the n smallest values. Is there a better function for this?

Yes, use `eigs`. If `n` is much smaller than the size of the matrix, as is typical in many applications, this can be computed _much_ more efficiently than computing all the eigenvalues.

---

<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:** [March 10, 2018, 8:47pm UTC](https://discourse.julialang.org/t/special-method-for-finding-eigenvalues-of-a-large-sparse-symmetric-matrix/9621/10 "2018-03-10T20:47:06Z")

</div>

> [@Ralph\_Smith](#):
>
> For small eigenvalues, `eigs` has to factor a matrix at every iteration. If you are using a sparse matrix type, that is done on one core, and is expensive. 14000 is not too large, so you are better off converting to a full matrix, which is factored fast with BLAS.

If you only want a few eigenvalues, I’m skeptical that it is advantageous to convert a matrix of that size to dense.

For the smallest-magnitude eigenvalues, it should be able to re-use the same (sparse) factorization over and over again—`eigs` does _not_ need to re-factorize on every iteration—in which case solving for a few eigenvalues should be vastly faster than a dense solve.

Indeed, I just checked that this is what `eigs` seems to do, and it appears to be reasonably efficient for sparse problems. I can solve for the 20 smallest-magnitude eigenvalues of a 100000x100000 real-symmetric matrix (coming from a 2d discretized Laplacian) in about 2 seconds on my laptop, orders of magnitude faster than a dense eigensolver would be even if I had enough memory (which I don’t — it would require more than 80GB of RAM). (Try the [second example from this notebook](http://nbviewer.jupyter.org/url/math.mit.edu/~stevenj/18.303/min-max-examples.ipynb), changing to `N=500` to increase the matrix size.)

---

<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:** [March 10, 2018, 9:28pm UTC](https://discourse.julialang.org/t/special-method-for-finding-eigenvalues-of-a-large-sparse-symmetric-matrix/9621/11 "2018-03-10T21:28:26Z")

</div>

My mistake, my recollection of the algorithm was confused.

What I meant to advocate was converting to dense and then calling `eigs`, since the Krylov subspace schemes will indeed win for small `nev`. I got 700s for sparse and 33s for dense 14000x14000 with `eigs` (random symmetric structure, about 10% fill), vs. about 800s for a dense `eigfact`. For very large and very sparse cases like your example, of course things will be as you show, especially with simple structure. My claim is just that when level 2 and 3 BLAS **are** practical, they often beat the sparse versions.

---

<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:** [March 10, 2018, 9:57pm UTC](https://discourse.julialang.org/t/special-method-for-finding-eigenvalues-of-a-large-sparse-symmetric-matrix/9621/12 "2018-03-10T21:57:15Z")

</div>

For the record, it seems this is particularly good performance because it’s a 2D sparsity pattern. With a modified version of your code for the first eigenvalue of the Laplacian, with about 21k dof in both cases, I get 0.09s in 2D, 1.29s in 3D (which is of course better than a dense eigensolver anyway). I seem to remember there’s an asymptotic difference in performance from the dimensionality of the problem.

---

<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:** [March 11, 2018, 12:05am UTC](https://discourse.julialang.org/t/special-method-for-finding-eigenvalues-of-a-large-sparse-symmetric-matrix/9621/13 "2018-03-11T00:05:51Z")

</div>

10% fill is huge. There is usually no point in using a sparse solver if you have that many nonzero entries. (On the order of 10 nonzero entries per row is common in real sparse-matrix applications, which corresponds to \< 0.1% fill for a 14000x14000 matrix.)

And sparse-direct solvers will generally perform poorly for random fill. The algorithms work best for sparsity that comes from local interactions on a mesh or grid.

1d grids are best, and 2d grids/meshes are still very good. For 3d grids, the scaling of sparse-direct methods is worse (this is a consequence of the difficulty of partitioning meshes in higher dimensions, essentially…see e.g. Davis’ book on sparse-direct methods). For very large 3d problems you usually have to switch to purely iterative methods (GMRES, CG, etc) and hope you can find a good preconditioner.

---

<div class="post-metadata">

**Author:** ![406700](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/406700/32/3705_2.png) [@406700](https://discourse.julialang.org/u/406700)\
**Post date:** [March 13, 2018, 8:03pm UTC](https://discourse.julialang.org/t/special-method-for-finding-eigenvalues-of-a-large-sparse-symmetric-matrix/9621/14 "2018-03-13T20:03:21Z")

</div>

Wow thanks for all the suggestions. I have to say I wasn’t expecting such a good response on my first post here.

I actually found out that I had a dumb mistake in the code. I converted to a dense matrix to solve using eigfact,  
and forgot to convert back to sparse to sparse for eigs. I have ~ .0003% filling so sparse is the way to go.  
I can now solve the code with eigs:  
@timev: 0.901941 seconds (10.80 k allocations: 4.219 MiB)

---

<div class="post-metadata">

**Author:** ![406700](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/406700/32/3705_2.png) [@406700](https://discourse.julialang.org/u/406700)\
**Post date:** [March 13, 2018, 8:13pm UTC](https://discourse.julialang.org/t/special-method-for-finding-eigenvalues-of-a-large-sparse-symmetric-matrix/9621/15 "2018-03-13T20:13:30Z")

</div>

> [@406700](#):
>
> 03

sorry correction more like .03% fill

---

<div class="post-metadata">

**Author:** ![406700](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/406700/32/3705_2.png) [@406700](https://discourse.julialang.org/u/406700)\
**Post date:** [March 13, 2018, 8:27pm UTC](https://discourse.julialang.org/t/special-method-for-finding-eigenvalues-of-a-large-sparse-symmetric-matrix/9621/16 "2018-03-13T20:27:17Z")

</div>

for anyone curious here is the full comparison of @timev for eigs with sparse and dense matrices  
sparse:  
0.901941 seconds (10.80 k allocations: 4.219 MiB)  
elapsed time (ns): 901940886  
bytes allocated: 4423584  
pool allocs: 10789  
non-pool GC allocs:1  
malloc() calls: 6

dense:  
194.259476 seconds (149.69 k allocations: 1.541 GiB, 0.00% gc time)  
elapsed time (ns): 194259475636  
gc time (ns): 1994472  
bytes allocated: 1655020183  
pool allocs: 149651  
non-pool GC allocs:30  
malloc() calls: 5  
GC pauses: 1

---

<div class="post-metadata">

**Author:** ![Carol](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carol/32/8140_2.png) [@Carol](https://discourse.julialang.org/u/Carol)\
**Post date:** [May 29, 2019, 5:23pm UTC](https://discourse.julialang.org/t/special-method-for-finding-eigenvalues-of-a-large-sparse-symmetric-matrix/9621/17 "2019-05-29T17:23:08Z")

</div>

Hi, I have the same problem here, and I cannot open your link. Could you please tell me which link is for `eigs`?

Thanks

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [May 29, 2019, 5:50pm UTC](https://discourse.julialang.org/t/special-method-for-finding-eigenvalues-of-a-large-sparse-symmetric-matrix/9621/18 "2019-05-29T17:50:37Z")

</div>

`eigs` has moved to Arpack.jl. Documentation is here: [Home · Arpack.jl](https://julialinearalgebra.github.io/Arpack.jl/stable/).