# Seeking faster method for finding top eigenvalues of large dense matrices
**URL:** https://discourse.julialang.org/t/seeking-faster-method-for-finding-top-eigenvalues-of-large-dense-matrices/134231
**Category:** Numerics
**Tags:** question, package
**Created:** [November 30, 2025, 8:23am UTC](https://discourse.julialang.org/t/seeking-faster-method-for-finding-top-eigenvalues-of-large-dense-matrices/134231 "2025-11-30T08:23:37Z")
**Posts on this page:** 19
**Page:** 1
### Author:  [@Udonge](https://discourse.julialang.org/u/Udonge)
#### Post date: [November 30, 2025, 8:23am UTC](https://discourse.julialang.org/t/seeking-faster-method-for-finding-top-eigenvalues-of-large-dense-matrices/134231/1 "2025-11-30T08:23:37Z")
Hello everyone,
I am new to the Julia language, and I’ve encountered a performance issue with an eigenvalue computation. Please feel free to point out any confusing statements or beginner mistakes!
I need to compute the **largest magnitude eigenvalues** (since they are complex) for a set of **large, dense matrices** (size ∼4^8×4^8) with `DoubleFloat` precision.
My current implementation is **extremely slow**!
```
using LinearAlgebra
using DoubleFloats
using GenericLinearAlgebra
using GenericSchur
function biggest_values(S::Matrix{ComplexDF64})::Double64
lambda2 = 1.0
λs = eigvals(S)
mu_pairs = Pair{Float64, ComplexDF64}[]
for λ in λs
diff = λ - 1.0
abs_diff = abs(diff)
if abs_diff < eps()
abs_mu = Inf
else
abs_mu = 1.0 / abs_diff
end
push!(mu_pairs, abs_mu => λ)
end
p = sortperm(mu_pairs, by=p -> abs(p.second), rev=true)
λ2 = mu_pairs[p[2]].second
lambda2 = abs(λ2)
return lambda2
end
```
I am aware that the largest eigenvalue magnitude of the matrices I am dealing with is abs(λ1)=1. I need to find the top two to three eigenvalues by magnitude, e.g., abs(λ2).
Is there an efficient way to calculate these few largest eigenvalues quickly while maintaining `DoubleFloat` precision?
Thank you!
---
### Author:  [@RomeoV](https://discourse.julialang.org/u/RomeoV)
#### Post date: [November 30, 2025, 9:08am UTC](https://discourse.julialang.org/t/seeking-faster-method-for-finding-top-eigenvalues-of-large-dense-matrices/134231/2 "2025-11-30T09:08:46Z")
Probably you should replace the `eigvals(S)` computation, which computes all 4^8 eigenvalues, by an iterative method that only computes the top two, for instance [Power method · IterativeSolvers.jl](https://iterativesolvers.julialinearalgebra.org/stable/eigenproblems/power_method/#Implementation-details).
Besides the algorithmic discussion, I’d also suggest you first try with regular Float64 eltypes and without `GenericLinearAlgebra` to check if that runs “fast enough”. Your extra requirements will almost certainly be slower than this “normal” baseline.
---
### Author:  [@Udonge](https://discourse.julialang.org/u/Udonge)
#### Post date: [November 30, 2025, 10:50am UTC](https://discourse.julialang.org/t/seeking-faster-method-for-finding-top-eigenvalues-of-large-dense-matrices/134231/3 "2025-11-30T10:50:10Z")
Thank you! This method was very effective and fast! But if I want to find the third and fourth largest eigenvalue magnitudes, this method no longer seems to work. Is there an effective method that can solve them?
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [November 30, 2025, 12:08pm UTC](https://discourse.julialang.org/t/seeking-faster-method-for-finding-top-eigenvalues-of-large-dense-matrices/134231/4 "2025-11-30T12:08:29Z")
> [@Udonge](#):
>
> But if I want to find the third and fourth largest eigenvalue magnitudes, this method no longer seems to work. Is there an effective method that can solve them?
Use the Arnoldi iterative method, a Krylov method which can find several eigenvalues.
(The basic power method is the first iterative method that everyone learns in linear algebra, but it is mostly superseded by the superior Arnoldi method. Essentially, Arnoldi accelerates the power method by keeping a “memory” of the space of all the vectors in the iteration and searches this whole “Krylov subspace” for the best eigenvector estimates.)
Another thing you could do is to compute the eigenvalues in `Float64` precision first, which is much faster, and then use these as initial guesses for shift-invert steps, or even Newton steps, though I don’t know if there is a good implementation of a Newton eigensolver for general nonsymmetric matrices lying around in Julia right now. (QuadGK.jl has a Newton eigensolver that it uses for arbitrary-precision eigensolves, starting with `Float64` initial guesses, but its implementation is only for symmetric matrices.)
---
### Author:  [@abraemer](https://discourse.julialang.org/u/abraemer)
#### Post date: [November 30, 2025, 2:14pm UTC](https://discourse.julialang.org/t/seeking-faster-method-for-finding-top-eigenvalues-of-large-dense-matrices/134231/5 "2025-11-30T14:14:02Z")
> [@stevengj](#):
>
> Use the Arnoldi iterative method, a Krylov method which can find several eigenvalues.
As for how to this: I would just try to stick your matrix into KrylovKit.eigsolve:
> **[Eigenvalue problems · KrylovKit.jl](https://jutho.github.io/KrylovKit.jl/stable/man/eig/#KrylovKit.eigsolve)**
>
> Documentation for KrylovKit.jl.
That should also work with double floats I think.
---
### Author:  [@lkdvos](https://discourse.julialang.org/u/lkdvos)
#### Post date: [December 1, 2025, 6:20am UTC](https://discourse.julialang.org/t/seeking-faster-method-for-finding-top-eigenvalues-of-large-dense-matrices/134231/6 "2025-12-01T06:20:38Z")
Unfortunately KrylovKit currently only handles arbitrary vector and operator types, but still more or less requires BLAS floats as the underlying scalars. Under the hood, we still are depending on LAPACK for some of the factorizations of the (smaller) matrices, but there are some ongoing efforts to extend that. No due dates yet though, unfortunately…
---
### Author:  [@Udonge](https://discourse.julialang.org/u/Udonge)
#### Post date: [December 1, 2025, 6:52am UTC](https://discourse.julialang.org/t/seeking-faster-method-for-finding-top-eigenvalues-of-large-dense-matrices/134231/7 "2025-12-01T06:52:52Z")
> [@stevengj](#):
>
> Use the Arnoldi iterative method, a Krylov method which can find several eigenvalues.
Thank you for your help! Now I’m trying to use `ArnoldiMethod.jl` to handle my problem with `Float64`, now the speed is much faster than the old one. But the final result from numerical calculation (~10^-4) is ~10^3 larger than theoretical prediction (~10^-7), so I guess I have to try better methods to get high precision eigenvalues.
> [@stevengj](#):
>
> Another thing you could do is to compute the eigenvalues in `Float64` precision first, which is much faster, and then use these as initial guesses for shift-invert steps, or even Newton steps, though I don’t know if there is a good implementation of a Newton eigensolver for general nonsymmetric matrices lying around in Julia right now. ([QuadGK.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/QuadGK) has a Newton eigensolver that it uses for arbitrary-precision eigensolves, starting with `Float64` initial guesses, but its implementation is only for symmetric matrices.)
That’s a good idea! I’ll try it, thank you for all your suggestions!
---
### Author:  [@Udonge](https://discourse.julialang.org/u/Udonge)
#### Post date: [December 1, 2025, 6:58am UTC](https://discourse.julialang.org/t/seeking-faster-method-for-finding-top-eigenvalues-of-large-dense-matrices/134231/8 "2025-12-01T06:58:27Z")
> [@abraemer](#):
>
> I would just try to stick your matrix into KrylovKit.eigsolve
Thank you for your recommend!
But this package is not support with `DoubleFloat` precision ([issue](https://github.com/Jutho/KrylovKit.jl/issues/137)). Now I’m trying to find other package supported `DoubleFloat`. Your suggestion is very valuable for me.
---
### Author:  [@Udonge](https://discourse.julialang.org/u/Udonge)
#### Post date: [December 1, 2025, 7:03am UTC](https://discourse.julialang.org/t/seeking-faster-method-for-finding-top-eigenvalues-of-large-dense-matrices/134231/9 "2025-12-01T07:03:48Z")
> [@lkdvos](#):
>
> Unfortunately KrylovKit currently only handles arbitrary vector and operator types, but still more or less requires BLAS floats as the underlying scalars. Under the hood, we still are depending on LAPACK for some of the factorizations of the (smaller) matrices, but there are some ongoing efforts to extend that. No due dates yet though, unfortunately…
Thank you for the reminder! Now I’m trying to find other package supported `DoubleFloat`. I hope to see this support implemented in `KrylovKit.jl` soon.
---
### Author:  [@hz-xiaxz](https://discourse.julialang.org/u/hz-xiaxz)
#### Post date: [December 1, 2025, 9:01am UTC](https://discourse.julialang.org/t/seeking-faster-method-for-finding-top-eigenvalues-of-large-dense-matrices/134231/10 "2025-12-01T09:01:39Z")
also check out [GitHub - dzhang314/MultiFloats.jl: Fast, SIMD-accelerated extended-precision arithmetic for Julia](https://github.com/dzhang314/MultiFloats.jl)
which claims a 2x speedup.
Also, I would think the ArnoldiMethods.jl problem can be solved by adding up `atol` precision?
---
### Author:  [@Udonge](https://discourse.julialang.org/u/Udonge)
#### Post date: [December 1, 2025, 11:17am UTC](https://discourse.julialang.org/t/seeking-faster-method-for-finding-top-eigenvalues-of-large-dense-matrices/134231/11 "2025-12-01T11:17:06Z")
> [@hz-xiaxz](#):
>
> also check out [GitHub - dzhang314/MultiFloats.jl: Fast, SIMD-accelerated extended-precision arithmetic for Julia](https://github.com/dzhang314/MultiFloats.jl)
Thank you for your recommend! The speedup from `MultiFloats.jl` is pretty important for my calculations. Now the main issue is finding a Arnold-like package supported higher precision…
> [@hz-xiaxz](#):
>
> Also, I would think the [ArnoldiMethods.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/ArnoldiMethods) problem can be solved by adding up `atol` precision?
Yes, I have try to use `ArnoldiMethods.jl` with `eps()` precision, but the numerical results (10^-8) are still ~10^2 greater than theoretical prediction values (~10^-10):
```
using ArnoldiMethod
using LinearAlgebra
function biggest_values(S::Matrix{ComplexF64})::Float64
lambda2 = 1.0
n = size(S, 1)
decomp, _ = partialschur(construct_linear_map(S-I(n)); nev = 2, which = :LM, tol = eps(), mindim = 4, maxdim = 8)
λs_inv, _ = partialeigen(decomp)
λs = 1 ./ λs_inv
lambda2 = abs(λs[2] + 1.0) - abs(λs[1] + 1.0) + 1.0
return lambda2
end
```
Here I have used the fact: the largest eigenvalue is abs(lambda1)=1 from theory, so I add a correction term `1.0 - abs(λs[1] + 1.0)` to lambda2. But this is still not enough. So I’m trying to use `doublefloat` now. Thank you for all your suggestions!
---
### Author:  [@hz-xiaxz](https://discourse.julialang.org/u/hz-xiaxz)
#### Post date: [December 2, 2025, 1:58pm UTC](https://discourse.julialang.org/t/seeking-faster-method-for-finding-top-eigenvalues-of-large-dense-matrices/134231/12 "2025-12-02T13:58:43Z")
> [@Udonge](#):
>
> `function biggest_values(S::Matrix{ComplexF64})::Float64`
Please remove the type annotation, which limits the matrix precision to ComplexF64, and the return type annotation. In Julia, these types could be correctly inferred, restricting these may lose accuracy with respect to DoubleFloats.
It might be beneficial to see the DoubleFloat.jl usage example in [GitHub - JuliaLinearAlgebra/ArnoldiMethod.jl: The Arnoldi Method with Krylov-Schur restart, natively in Julia.](https://github.com/JuliaLinearAlgebra/ArnoldiMethod.jl)
---
### Author:  [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)
#### Post date: [December 3, 2025, 11:52pm UTC](https://discourse.julialang.org/t/seeking-faster-method-for-finding-top-eigenvalues-of-large-dense-matrices/134231/13 "2025-12-03T23:52:07Z")
It appears that the eigenvalues of interest to you are close together. This is challenging for the Arnoldi schemes, so one usually does shift-invert to make the cluster easier to resolve accurately, but that also requires an expensive factorization in high precision.
An alternative is to get rough estimates from the Arnoldi scheme and use iterative refinement for the cluster (like Arnoldi this just uses matrix-vector products, not large-scale factorization). This is now straightforward with the IterativeRefinement.jl package.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [December 4, 2025, 1:40am UTC](https://discourse.julialang.org/t/seeking-faster-method-for-finding-top-eigenvalues-of-large-dense-matrices/134231/14 "2025-12-04T01:40:09Z")
> [@stevengj](#):
>
> Another thing you could do is to compute the eigenvalues in `Float64` precision first, which is much faster, and then use these as initial guesses for shift-invert steps, or even Newton steps, though I don’t know if there is a good implementation of a Newton eigensolver for general nonsymmetric matrices lying around in Julia right now. ([QuadGK.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/QuadGK) has a Newton eigensolver that it uses for arbitrary-precision eigensolves, starting with `Float64` initial guesses, but its implementation is only for symmetric matrices.)
Just for reference, I wrote an implementation of this scheme for general nonsymmetric matrices here: [fast bignum eigensolver exploiting Newtons' method from Float64? · Issue #167 · JuliaLinearAlgebra/GenericLinearAlgebra.jl · GitHub](https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl/issues/167)
However, because it requires an O(m^3) Hessenberg factorization, it is still probably too slow for your case.
---
### Author:  [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)
#### Post date: [December 4, 2025, 11:39am UTC](https://discourse.julialang.org/t/seeking-faster-method-for-finding-top-eigenvalues-of-large-dense-matrices/134231/15 "2025-12-04T11:39:47Z")
It might be helpful if you told us a little more about your problem, and why e.g. Float32 is not enough for you.
Between the lines, you suggested that you know from theoretical considerations that your matrix has largest eigenvalue 1. Is your matrix maybe stochastic? Other potentially useful properties?
Then, you care about the next couple of leading eigenvalues. Do you have a spectral gap between the eigenvalues you care about and the ones you don’t care about? Both empirically and theoretically.
Then, is your issue that eigenvalues you care about almost collide?
If so, why do you believe that eigenvalues / eigenspaces are the thing to care about?
If not, then why do you care about knowing the eigenvalues with high precision?
Why do you believe that DoubleFloats will give you better results than Float64, if you only know your input data up to Float64 precision?
Or is your setting rather: You know that your leading eigenvalue is 1; and all other eigenvalues have smallish absolute value. You don’t particularly care about all the other eigenvalues, but you do care about the spectral radius that remains after removing the leading eigenvalue (this happens to be the second eigenvalue).
If you had a stochastic matrix from a markov chain, this would tell you the convergence rate to the ergodic measure thingy, and would nicely explain why you care whether the second eigenvalue has absolute value ~1e-5 or ~1e-8. (I don’t know how to solve this specific question, but I’m pretty sure people have thought about it a lot; and using “generic” methods would suck, because you have a separation of scales between the leading eigenvalue and the remaining stuff you care about, which is very bad for floating point numbers!)
---
### Author:  [@Udonge](https://discourse.julialang.org/u/Udonge)
#### Post date: [December 4, 2025, 12:29pm UTC](https://discourse.julialang.org/t/seeking-faster-method-for-finding-top-eigenvalues-of-large-dense-matrices/134231/16 "2025-12-04T12:29:22Z")
> [@Ralph\_Smith](#):
>
> An alternative is to get rough estimates from the Arnoldi scheme and use iterative refinement for the cluster (like Arnoldi this just uses matrix-vector products, not large-scale factorization). This is now straightforward with the [IterativeRefinement.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/IterativeRefinement) package.
Thank you for your suggestions! I’ll try this package latter. By using `ArnoldiMethod.jl` and shift-invert I can get good estimates now.
---
### Author:  [@Udonge](https://discourse.julialang.org/u/Udonge)
#### Post date: [December 4, 2025, 12:36pm UTC](https://discourse.julialang.org/t/seeking-faster-method-for-finding-top-eigenvalues-of-large-dense-matrices/134231/17 "2025-12-04T12:36:45Z")
> [@stevengj](#):
>
> Just for reference, I wrote an implementation of this scheme for general nonsymmetric matrices here: [fast bignum eigensolver exploiting Newtons’ method from Float64? · Issue #167 · JuliaLinearAlgebra/GenericLinearAlgebra.jl · GitHub](https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl/issues/167)
Thank you for your helpful reference! I’ll try your method later.
The speed is the secondary important question now, I hope this method gives more precise eigenvalues.
---
### Author:  [@Udonge](https://discourse.julialang.org/u/Udonge)
#### Post date: [December 4, 2025, 12:57pm UTC](https://discourse.julialang.org/t/seeking-faster-method-for-finding-top-eigenvalues-of-large-dense-matrices/134231/18 "2025-12-04T12:57:00Z")
> [@foobar\_lv2](#):
>
> Between the lines, you suggested that you know from theoretical considerations that your matrix has largest eigenvalue 1. Is your matrix maybe stochastic? Other potentially useful properties?
YES, the matrix which I handled is a random matrix with some physical settings, they are not unitary, but I can solve its largest eigenvalue with abs(λ)=1 corresponding to steady state. And all other eigenvalues with abs(λ)\<1 because this is a dissipation process.
> [@foobar\_lv2](#):
>
> Then, you care about the next couple of leading eigenvalues. Do you have a spectral gap between the eigenvalues you care about and the ones you don’t care about? Both empirically and theoretically.
YES, they have spectral gaps between each other from theoretical considerations, magnitudes ~2^-d where d is the dimension of matrix.
> [@foobar\_lv2](#):
>
> Then, is your issue that eigenvalues you care about almost collide?
YES, the spectral gap I care about is ~10^-6-10^-12.
> [@foobar\_lv2](#):
>
> If so, why do you believe that eigenvalues / eigenspaces are the thing to care about?
Well, because the only way to get these gaps I know is calculating the differences between eigenvalues. If you have any suggestions about that, they’re valuable for me!
> [@foobar\_lv2](#):
>
> Why do you believe that DoubleFloats will give you better results than Float64, if you only know your input data up to Float64 precision?
The random matrix I use is created in `DoubleFloat` precision, but I didn’t display them here. I apologize for any confusions caused by this.
> [@foobar\_lv2](#):
>
> Or is your setting rather: You know that your leading eigenvalue is 1; and all other eigenvalues have smallish absolute value. You don’t particularly care about all the other eigenvalues, but you do care about the spectral radius that remains after removing the leading eigenvalue (this happens to be the second eigenvalue).
YES, this is what I wanna. I know the leading eigenvalue is exactly 1, and spectral gaps are pretty small. And the most important thing I care about is the gap between abs(λ1)=1 and abs(λ2).
---
### Author:  [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)
#### Post date: [December 4, 2025, 1:49pm UTC](https://discourse.julialang.org/t/seeking-faster-method-for-finding-top-eigenvalues-of-large-dense-matrices/134231/19 "2025-12-04T13:49:11Z")
Some more clarification questions:
1. Is your matrix stochastic? (you answered “Yes, it is either stochastic or has some other useful physical property”, without giving more hints on the useful physical property)
2. Do you care about \lambda\_3, at all? You were not quite clear on whether you only care about the spectral radius after removing \lambda\_1, or also about lower eigenvalues.
3. If you care about \lambda\_3, is there any theoretical reason why |\lambda\_2| \neq |\lambda\_3|?
> [@Udonge](#):
>
> magnitudes ~2^-d where d is the dimension of matrix
Since d~2^16, double float won’t cut any cakes on resolving that; you would need an eps of `2^-65000`.
I think your setting is something like `A = v*v' + B`, where `v` is your stationary distribution, and B is your leftover that you care about; however, `v` is on scale one, while `B` is close to zero. So even if you knew `v` and `B`, you would lose most info if you computed `BB = (B + v*v') - v*v'`, due to numerical errors.
The ideal answer would be to never have `A` exist, and directly generate `v`, `v'` and `B`.
Generally, I would think that using generic eigenvalue techniques is bad here. Quick googling found a bunch papers/algorithms that directly tackle “estimate/compute |\lambda\_2| for stochastic matrix”, aka “compute diffusion/convergence rate for markov process”; you care especially about the case where diffusion/convergence is very fast, i.e. |\lambda\_2| \ll 1, as opposed to the case of slow convergence 1-|\lambda\_2|\ll 1 (I think that was your case, right? Just tell us the rough size of |\lambda\_2|, explicitly, in Float64 or even Float32. Both settings will probably want entirely different approaches / algorithms!).
I think you should write out a more detailed problem description for us/yourself, and do a little research on specific algorithms for this issue.