I am trying to run a model which simulates a long time series. For each time step, I need to solve the Eigen values and Eigen vectors for one matrix. This matrix is only 4*4, but the computation time is important because the huge number of for-loops of the time steps…
The matrix is sparse, it only has values in diagnal and sub-diagonal elements, like below:
from the LinearAlgebra. However, profiling the whole model codes showing that these two lines:
creat a lot of allocations…because each time the eigvals and eigvecs create new arrays;
takes a lot of time…
So, is there any way to make it a bit faster, and ZERO allocations? I don’t want to use eigvals! to modify it at the original matrix…I could use a preallocated array prealloc_e_val for eigvals. Typically, what I want is: e_val = eigvals(A, prealloc_e_val)
with fastest speed…
Could you please help me with this? Thanks!
Why? What do you do with these at each timestep? (Maybe there is a way to avoid computing them entirely. And since your matrix is non-symmetric, the whole idea of diagonalizing it is often suspect, e.g. what would you do if the matrix were defective?)
The StaticArrays.jl package has optimized methods for working with small fixed-size vectors and matrices. If you have zillions of 4x4 matrices, they will be much more efficient if represented by the SMatrix type from that paper.
However, eigenvectors/values of non-symmetric matrices are not currently supported by StaticArrays.jl.
Note that StaticArrays.jl does not have any methods for eigenvalues of 4x4 matrices and falls back on allocating and re-wrapping, so it’s actually slightly slower than working with a regular array.
julia> using StaticArrays
julia> let A = rand(4, 4)
As = SMatrix{4, 4}(A)
@btime eigvecs($A)
@btime eigvecs($As)
end;
4.100 μs (30 allocations: 5.81 KiB)
4.129 μs (32 allocations: 6.02 KiB)
I don’t know about zero allocations, but it’s possible to use LAPACK in julia. E.g. geev! to find the eigenvalues and the left and right eigenvectors. This might be somewhat faster than calling eigvals/eigvecs. For symmetric matrices, there’s a syev!. There are still some allocations, though.
prealloc .= A # make a copy of A, so it's not destroyed
W, VL, VR = LAPACK.geev!('N', 'V', A) # see ?LAPACK.geev!
There is also a possibility to call the actual unwrapped LAPACK routines via @ccall. Since they don’t allocate, all the workspace must be preallocated. It’s a bit tricky, but absolutely possible. There are quite a few arguments: LAPACK: dgeev
Here’s how to get zero allocations with a direct LAPACK call:
using LinearAlgebra
const dgeev = let
preallocA = fill(0.0, 4, 4)
WR = fill(0.0, 4) # real part of eigval
WI = fill(0.0, 4) # imag part
VL = fill(0.0, 4, 4) # left eigvecs
VR = fill(0.0, 4, 4) # right eigvecs
LWORK = Int32(64)
WORK = Vector{Float64}(undef, LWORK)
INFO = Ref(Int32(0))
function mydgeev(A)
preallocA .= A
@ccall "liblapack".dgeev_(
Ref(Int8('N'))::Ptr{Cchar}, # no left eigvecs
Ref(Int8('V'))::Ptr{Cchar},
Ref{Int32}(size(A, 2))::Ptr{Cint},
preallocA::Ptr{Float64},
Ref{Int32}(size(A, 1))::Ptr{Cint},
WR::Ptr{Float64},
WI::Ptr{Float64},
VL::Ptr{Float64},
Ref{Int32}(size(VL, 1))::Ptr{Cint},
VR::Ptr{Float64},
Ref{Int32}(size(VR, 1))::Ptr{Cint},
WORK::Ptr{Float64},
Ref(LWORK)::Ptr{Cint},
INFO::Ptr{Cint}
)::Cvoid
return WR, WI, VR
end
end
A = rand(4,4)
dgeev(A)
The @ccall to "liblapack".dgeev_ might use a different name on your machine. “liblapack” is the name of the lapack library which is linked into julia, “dgeev_” with the underscore at the end is a convention for some (most?) fortran compilers. The call arguments are standardized. It’s a bit tricky because lapack is written in/for fortran, where all arguments are passed by reference.
The call to dgeev(A) does not allocate. The outputs (WR, WI, VR) are overwritten by the next call to dgeev. So if you store them in some structure, make sure to copy them.
Complex eigenvectors may come in conjugate pairs. See the docs (linked above) for how they are coded in VR.
Note that the dgeev above is not thread safe. If this is a parallel loop, the preallocated buffers (between let and function mydgeev) must be task-local, e.g. stored in Base.task_local_storage() with a get!.
If your matrix is complex, you must call zgeev_ instead, with different arguments.
I have been hesitating to point it out, because I think the bars over your variables mean your matrix is probably complex and that negative sign in the (4,3) position looks like it could be a deal breaker even if a lot of other things are real and positive, but a real tridiagonal T is diagonally similar to a symmetric tridiagonal if t_{j, j+1} t_{j+1, j} > 0 for each sub/super diagonal pair. Something similar holds for a complex matrix if the diagonal is real and the products of the sub/super diagonals are real and positive. On the off chance that this is the case for your matrix, it does open the possibility of using some simpler algorithms that might be faster.
Hi @Oscar_Smith, thanks for your idea! Yes I am also thinking about this, but I didn’t find a solution after I googled on the internet…mostly people are talking about symmetric matrix…
Hi @mstewart thanks for your points! It does make sense, however, the bar means the average between one item and another…nvd, but the final solution does include complex numbers in Eigen values and vectors, so that the final solution of the ODE always takes the real part…
You should just use exp(A * Δt) for matrix exponentials, not eigenvectors, despite what you may have learned in school.
As explained in the famous paper “19 Dubious Ways to Compute a the Exponential of a Matrix”, using eigenvectors is a bad idea to exponentiate non-Hermitian/non-normal matrices, because it might be very inaccurate if the matrix happens to be close to defective.
That being said, if you are solving a system of ODEs, then there are probably much better methods available. See DifferentialEquations.jl.
Meanwhile, I am also curious to see better solution to solve the ODE in above post…I wrote the analytical solution and hard-coded it. However, the most time consuming comes from the Eigen values and vectors calculation (a simple @profileview shows this…)…so I am wondering is there any other way to avoid the eigen vectors calculation? I tried with numerical package suggested by ChatGPT, however, they are slower than the analytical one…
It looks like you are trying to write down the solution for an ODE for just one “timestep”, i.e. for just \Delta t. You should take a step backwards and write down the equations you are actually solving for allt and send those equations to an ODE solver. (A good algorithm will probably not take fixed timestep sizes.)
But it’s hard to give you specific advice without more context. (This is a classic XY problem.)
yes, it’s true. However, this ODE is interacted in a complex model system which consists of empirical/other analytical equations…it’s a huge model, I don’t think it’s a chance to put the equations for the whole time to an ODE solver…because:
the states simulated by this ODE is affected and affect other states from other sub-models…(so-called “interacted”);
currently the big model is write for this only one time step…so that there is a for-loop wrap-up thing out side, and people could run it by specifying the range of the time step…
So you just need to exponentiate the matrix? In that case, StaticArrays.jl is definitely your friend.
julia> let A = rand(4, 4)
As = SMatrix{4, 4}(A)
@btime exp($A)
@btime exp($As)
end;
1.499 μs (18 allocations: 1.58 KiB)
183.229 ns (0 allocations: 0 bytes)
This is way faster than eigensolving.
Also, if what you really want is to calculate exp(t * A) * v, then ExponentialUtilities.jl has tools for doing that even more efficiently than calculating the full matrix exponential:
julia> let A = rand(4, 4)
As = SMatrix{4, 4}(A)
v = rand(4)
vs = SVector{4}(v)
t = 1e-3
# calculate exp(t * A) * v
@btime expv($t, $A, $v)
@btime expv($t, $As, $vs)
end;
1.904 μs (27 allocations: 2.22 KiB)
146.822 ns (0 allocations: 0 bytes)
It’s not clear to me what a lot of this means here, I think the actual problem being solved is still pretty underspecified, but I’ll just mention that DifferentialEquations.jl does support building a differential equation and then remake-ing it with different time-spans.
E.g. you could do something like
using OrdinaryDiffEq, StaticArrays
function f(u, p, t)
(; a, b, c) = p # p stores parameters a, b, and c
# build up some potentially time-dependant matrix
A = SA[a*sin(b*t) b+1
b-1 c*cos(b*t)]
A * u
end
# initial condition
u0 = SA[1.0
0.0]
p = (;a=1.0, b=2.0, c=3.0)
# the timespan we solve over
tspan = (0.0, 1.0)
# construct an ODE problem
prob = ODEProblem(f, u0, tspan, p)
