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?
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.
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.
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?
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.)
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 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.
That’s a good idea! I’ll try it, thank you for all your suggestions!
Thank you for your recommend!
But this package is not support with DoubleFloat precision (issue). Now I’m trying to find other package supported DoubleFloat. Your suggestion is very valuable for me.
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…
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!
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 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.
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!)
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.
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.
YES, they have spectral gaps between each other from theoretical considerations, magnitudes ~2^-d where d is the dimension of matrix.
YES, the spectral gap I care about is ~10^-6-10^-12.
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!
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.
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).
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)
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.
If you care about \lambda_3, is there any theoretical reason why |\lambda_2| \neq |\lambda_3|?
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.
