Initialize a Hamiltonian with a custom eigendecomposition function

When defining a Hamiltonian object, a user-defined eigendecomposition routine can be supplied using a keyword argument EIGS. All the eigendecomposition calls inside the solver will call this function instead of the default one.

For example:

using OpenQuantumTools

# Define a function to construct an eigendecomposition routine for the
# Hamiltonian. The routine should have the signature: (H, t, lvl) -> (w, v).
# The argument of this function is the cache used by the Hamiltonian object.
function build_user_eigen(u_cache)
    EIGS = function(H, t, lvl)
        println("I am the user defined eigendecomposition routine.")
        w, v = eigen(Hermitian(H(t)))
        w[1:lvl], v[:, 1:lvl]
    end
end

H = DenseHamiltonian([(s)->1-s, (s)->s], [σx, σz], EIGS=build_user_eigen)

eigen_decomp(H, 0.5, lvl=2)

I am the user defined eigendecomposition routine.
([-0.7071067811865476, 0.7071067811865476], Complex{Float64}[-0.38268343236
508984 + 0.0im 0.9238795325112866 - 0.0im; 0.9238795325112868 - 0.0im 0.382
6834323650897 - 0.0im])

Constant Hamiltonian

There are two applications for this functionality. First, if the Hamiltonian is constant, one can precalculate that Hamiltonian's eigensystem and build a function that returns those precalculated values. This is particularly helpful for the adiabatic master equation solver.

function build_user_eigen(u_cache)
    # note that to keep the unit consistent, the unit of any value inside the routine should be 1/h
    w, v = eigen(Hermitian(2*π*(σx+σz)))
    EIGS = function(H, t, lvl)
        w[1:lvl], v[:, 1:lvl]
    end
end

H = DenseHamiltonian([(s)->1.0], [σx+σz], EIGS=build_user_eigen)

print(eigen_decomp(H, 0.5, lvl=2))
print(eigen_decomp(H, 0.0, lvl=2))

([-1.4142135623730951, 1.4142135623730951], Complex{Float64}[0.382683432365
0897 + 0.0im -0.9238795325112867 + 0.0im; -0.9238795325112867 + 0.0im -0.38
26834323650897 + 0.0im])([-1.4142135623730951, 1.4142135623730951], Complex
{Float64}[0.3826834323650897 + 0.0im -0.9238795325112867 + 0.0im; -0.923879
5325112867 + 0.0im -0.3826834323650897 + 0.0im])

Sparse Hamiltonian

Another application is to supply special eigendecomposition algorithms for sparse matrices to take advantage of the sparsity.

For example, the default eigendecomposition algorithm for a sparse Hamiltonian is to convert it into dense matrices first and then perform the decomposition.

Hd = standard_driver(4, sp=true);
Hp = two_local_term(rand(3), [[1,2],[2,3],[3,4]], 4, sp=true)
H = SparseHamiltonian([(s)->1-s, (s)->s], [Hd, Hp], unit=:ħ)

# the default eigen_decomposition using the dense matrices algorithm
w, v = eigen_decomp(H, 0.1, lvl=4)

([-0.5734708734977536, -0.30096737258312073, -0.2962565535990648, -0.277406
23690260746], Complex{Float64}[-0.23932287167644808 + 0.0im -0.005131134167
813189 + 0.0im 0.005986080042454323 + 0.0im -0.3173042029285559 + 0.0im; 0.
24733781018649026 + 0.0im 0.029638541287947538 + 0.0im 0.3371592326467173 +
 0.0im -0.027102556462585706 + 0.0im; … ; 0.24733781018649004 + 0.0im -0.02
9638541287947306 + 0.0im -0.337159232646717 + 0.0im 0.027102556462585803 + 
0.0im; -0.23932287167644814 + 0.0im 0.005131134167813223 + 0.0im -0.0059860
80042454208 + 0.0im 0.31730420292855577 + 0.0im])

If the Hamiltonian size becomes large, we can use sparse algorithms provided by Arpack instead. Let's first load Arpack.jl by running:

using Arpack

Next, we can use an Arpack function to replace the default eigendecomposition routine:

function build_user_eigen(u_cache)
    function (H, t, lvl)
        hmat = H(t)
        println("Using sparse matrix algorithm")
        # define all the Arpack routine parameters here
        eigs(hmat, nev = lvl, which=:SR, tol=0.0, maxiter=300)
    end
end

Hd = standard_driver(4, sp=true);
Hp = two_local_term(rand(3), [[1,2],[2,3],[3,4]], 4, sp=true)
H = SparseHamiltonian([(s)->1-s, (s)->s], [Hd, Hp], unit=:ħ, EIGS =build_user_eigen)

eigen_decomp(H, 0.1, lvl=4)

Using sparse matrix algorithm
([-0.573543061746316, -0.3050672399440895, -0.2874370322374558, -0.28610530
626508396], Complex{Float64}[-0.10384629842404558 - 0.2146971342587944im -0
.06404231382566351 + 0.010325370743772974im -0.14186578272135514 + 0.142335
68960019177im 0.05569262969534014 - 0.09933364808438487im; 0.10838705182292
556 + 0.22408491944623096im -0.16652129098615318 + 0.026847781778207985im 0
.01162484664337514 - 0.011663352020027835im -0.2083895156368037 + 0.3716845
644385129im; … ; 0.10838705182292566 + 0.22408491944623113im 0.166521290986
15318 - 0.02684778177820791im -0.01162484664337515 + 0.01166335202002796im 
0.20838951563680333 - 0.3716845644385127im; -0.10384629842404566 - 0.214697
13425879416im 0.06404231382566354 - 0.01032537074377297im 0.141865782721355
22 - 0.14233568960019183im -0.055692629695339835 + 0.09933364808438495im])

---