How to efficiently exponentiate a complex (symmetric) tridiagonal matrix?

General Usage Performance
, , ,
1

Hi, I’m trying to efficiently obtain the exponential of a complex symmetric tridiagonal matrix (SymTridiagonal).

Unfortunately, exp is neither defined for SymTridiagonal nor Tridiagonal. I believe this is because eigen is not defined for those matrices (although eigen is defined for real-valued symmetric tridiagonal matrices, although not complex ones). An MWE:

julia> using LinearAlgebra

julia> n = 3
3

julia> St = SymTridiagonal(-2ones(n), ones(n-1))
3×3 SymTridiagonal{Float64,Array{Float64,1}}:
 -2.0   1.0    ⋅ 
  1.0  -2.0   1.0
   ⋅    1.0  -2.0

julia> T = Tridiagonal(ones(n-1), -2ones(n), ones(n-1))
3×3 Tridiagonal{Float64,Array{Float64,1}}:
 -2.0   1.0    ⋅ 
  1.0  -2.0   1.0
   ⋅    1.0  -2.0

julia> exp(St)  # similar error for exp(complex(St))
ERROR: MethodError: no method matching exp(::SymTridiagonal{Float64,Array{Float64,1}})
Closest candidates are:
  exp(::BigFloat) at mpfr.jl:619
  exp(::Missing) at math.jl:1167
  exp(::Complex{Float16}) at math.jl:1115
  ...
Stacktrace:
 [1] top-level scope at REPL[5]:1

julia> exp(T)  # similar error for exp(complex(T))
ERROR: MethodError: no method matching exp(::Tridiagonal{Float64,Array{Float64,1}})
Closest candidates are:
  exp(::BigFloat) at mpfr.jl:619
  exp(::Missing) at math.jl:1167
  exp(::Complex{Float16}) at math.jl:1115
  ...
Stacktrace:
 [1] top-level scope at REPL[6]:1

julia> eigen(St)
Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}
values:
3-element Array{Float64,1}:
 -3.414213562373093
 -1.9999999999999987
 -0.5857864376269051
vectors:
3×3 Array{Float64,2}:
 -0.5       -0.707107     0.5
  0.707107  -9.42055e-16  0.707107
 -0.5        0.707107     0.5

julia> eigen(T)  # similar error for eigen(complex(St)) and eigen(complex(T))
ERROR: MethodError: no method matching eigen!(::Tridiagonal{Float64,Array{Float64,1}}; permute=true, scale=true, sortby=LinearAlgebra.eigsortby)
Closest candidates are:
  eigen!(::SymTridiagonal{#s664,V} where V<:AbstractArray{#s664,1} where #s664<:Union{Float32, Float64}) at /build/julia-98cBbp/julia-1.4.1+dfsg/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/tridiag.jl:232 got unsupported keyword arguments "permute", "scale", "sortby"
  eigen!(::SymTridiagonal{#s664,V} where V<:AbstractArray{#s664,1} where #s664<:Union{Float32, Float64}, ::UnitRange) at /build/julia-98cBbp/julia-1.4.1+dfsg/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/tridiag.jl:235 got unsupported keyword arguments "permute", "scale", "sortby"
  eigen!(::SymTridiagonal{#s664,V} where V<:AbstractArray{#s664,1} where #s664<:Union{Float32, Float64}, ::Real, ::Real) at /build/julia-98cBbp/julia-1.4.1+dfsg/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/tridiag.jl:240 got unsupported keyword arguments "permute", "scale", "sortby"
  ...
Stacktrace:
 [1] eigen(::Tridiagonal{Float64,Array{Float64,1}}; permute::Bool, scale::Bool, sortby::typeof(LinearAlgebra.eigsortby)) at /build/julia-98cBbp/julia-1.4.1+dfsg/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/eigen.jl:237
 [2] eigen(::Tridiagonal{Float64,Array{Float64,1}}) at /build/julia-98cBbp/julia-1.4.1+dfsg/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/eigen.jl:235
 [3] top-level scope at REPL[8]:1

This might seem like a corner case, but SymTridiagonal matrices can be quite useful for solving 1D problems in Quantum Mechanics. An example application follows:

using LinearAlgebra

D²(n, Δx) = SymTridiagonal(-2ones(n), ones(n-1)) / Δx^2

function hamiltonian(x, V; m=1, ħ=1)
    # Assuming x is equally spaced
    Δx = x[2] - x[1]

    T = -(ħ^2 / 2m) * D²(length(x), Δx)
    return T + V
end

ħ = 1
ω₀ = 1
m = 1

t = 1
L = 5
Δx =0.1

x = -(L-Δx):Δx:(L-Δx)
V = Diagonal(m * ω₀^2 * x.^2 / 2);
H = hamiltonian(x, V, m=m);

typeof(H)  # => SymTridiagonal{Float64,Array{Float64,1}}
exp(-im * H * t/ ħ)
# =>
# ERROR: MethodError: no method matching exp(::SymTridiagonal{Complex{Float64},Array{Complex{Float64},1}})
# Closest candidates are:
#   exp(::BigFloat) at mpfr.jl:619
#   exp(::Missing) at math.jl:1167
#   exp(::Complex{Float16}) at math.jl:1115
#   ...
# Stacktrace:
#  [1] top-level scope at REPL[74]:1

I can do exp(Array(-im * H * t/ ħ)), but then I would lose the nice structure of H. Is there a more efficient way of doing this?

2

In general, it is not safe (not numerically stable) to use eigen-decomposition methods for exponentiating non-Hermitian matrices. The issue is that if you have an arbitrary non-Hermitian matrix, even an arbitrary complex-symmetric tridiagonal matrix, it may be very close to defective, which means the eigenvectors may form an ill-conditioned basis.

In consequence, for non-Hermitian matrices one falls back on a variety of other methods. If I remember correctly, a typical method involves writing A^p = (A^{p/n})^n for a sufficiently large n and then using a Padé approximant for A^{p/n}.

In principle, it might be possible to specialize some of these latter methods for complex-symtridiagonal matrices, but since such matrices aren’t that common I don’t know whether anyone has gone to the trouble. The good news is that Julia is a good language in which to implement such specialized algorithms if you have the knowledge and the time.

5 Likes
3

Do you need the matrix exponential, or the action of it on a vector? Since it is a Hamiltonian from quantum mechanics that you’re dealing with, you can just use one of the many approaches to propagate wavefunctions that are available in the literature. The Hamiltonian you have is Hermitian and purely real, the complexity comes about from the -im*t in the propagator, but this does not change any properties of the underlying Hamiltonian, except rotating the spectrum by -pi/2, i.e. all the real eigenvalues are now purely imaginary.

The simplest, but still stable (A-stable in fact) method that I know of is Crank–Nicolson:

\exp(-\mathrm{i}Ht) \approx [I + \mathrm{i} H t/2]^{-1}[I - \mathrm{i} H t/2]

which has a local error of t^3 and a global error of t^2. It amounts to multiplication by tridiagonal matrix, and solving by another [whose factorization you can easily precompute, factorize(I + im*H*t/2)]. This is a Padé approximant to the exponential function, as @stevengj mentioned.

There are also Krylov methods that can help you compute the action, see e.g. ExponentialUtilities.jl, but beware that they struggle when your grid spacing becomes too fine.

4 Likes
4

If that is the only reason for the non-Hermitian-ness, then you are have an easy task — the matrix A=iH is anti-Hermitian, not Hermitian, but it is still normal with orthogonal eigenvectors.

More simply, if you have a matrix A where H=-iA is Hermitian, then simply diagonalize H=Q\Lambda Q^* and use e^A = e^{iH}=Qe^{i\Lambda}Q^*.

Furthermore, if A is tridiagonal, then H is tridiagonal and Hermitian. If H is not already real, there is a trivial change of basis to make the problem real symmetric tridiagonal. So you can always use the efficient real-SymTridiagonal eigensolver.

2 Likes
5

Precisely, this is used all the time in QM; you usually build up the Krylov subspace using the real Hamiltonian, and then compute the subspace exponential by scaling the Ritz values (~eigenvalues) by -im*t, which works precisely because of the transformation you mentioned.

I might add that Crank–Nicolson can handle mild non-Hermitian Hamiltonians as well, such as the implementation of absorbing boundary conditions using a complex absorbing potential at the end of the box.

1 Like
6

For the same reason, we’ve discussed including an optimized cis(::Hermitian) function in the LinearAlgebra stdlib, in order to facilitate e^{iH}= \operatorname{cis}(H) computation for Hermitian matrices.

Should be a straightforward PR if someone wants to take a crack at it. (Update: done in Julia 1.7.)

3 Likes
7

Thank you all for the excellent suggestions!

@jagot, I realized I only need the time propagation, not the exponential explicitly. As such, your suggestion on the Crank–Nicolson method fits my particular use case perfectly, thanks!

@stevengj Somehow, I missed the decomposition you mentioned. This might be the best approach in the general case. In particular, the cis(::Hermitian) idea got me interested; I’ll make a PR for it soon.

2 Likes