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

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!

Try `eigs`.

1 Like

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

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

@ 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.

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?

You could try simeaultanious inverse iteration with Rayleigh shifts.

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.

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.

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, changing to `N=500` to increase the matrix size.)

3 Likes

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.

1 Like

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.

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.

4 Likes

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)

sorry correction more like .03% fill

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

1 Like

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

Thanks

`eigs` has moved to Arpack.jl. Documentation is here: https://julialinearalgebra.github.io/Arpack.jl/stable/.

1 Like