Skip to main content

Operator function & Krylov time evolution

Evolving quantum states is a common routine in numerical methods. There are many ways to do this, the most common and accessible being full exact diagonalization and evolving an initial state using knowledge of the full spectrum of a Hamiltonian HH,

∣ψ(t)⟩=∑mcme−iEmt∣Em⟩|\psi(t)\rangle = \sum_m c_m e^{-iE_m t} |E_m\rangle

where EkE_k and ∣Ek⟩|E_k\rangle are the eigenvalues and eigenvectors of H^{\hat{H}}. For a spin 12\frac{1}{2} Hamiltonian this limits most computers to 14−1814-18 sites due to large memory costs. The advantage however is that this is numerically exact. One can exploit symmetries and commonly get this number up to 24−2624-26 sites. In this document we will assume this numerical exactness is the goal of the user, so we will leave approximation schemes like time evolution with matrix product states for another tutorial. How can we go beyond these system sizes while staying numerically accurate? The answer to this is combining aleph's symbolic operator framework with Krylov time evolution. Krylov evolution is for example discussed here but we will produce the main arguments below.

Krylov evolution​

Let H^\hat{H} be the Hamiltonian that generates our dynamics. Our goal is to approximate the process,

∣ψ(t+δ)⟩=e−iδH^∣ψ(t)⟩. |\psi(t+\delta) \rangle = e^{-i \delta \hat{H}} |\psi(t) \rangle.

To do this we build a Krylov subspace of size NN, {∣ψ(t)⟩,H^∣ψ(t)⟩,H^2∣ψ(t)⟩,…,H^N−1∣ψ(t)⟩}\{\ket{\psi(t)}, \hat{H} \ket{\psi(t)},\hat{H}^2 \ket{\psi(t)}, \ldots, \hat{H}^{N-1} \ket{\psi(t)} \}. We want to build an orthonormal basis out of this subspace which is simple using a modified Gram-Schmidt procedure or performing a QR factorization on a matrix with these vectors as the columns of a matrix. The resulting basis we will call ∣vi⟩\ket{v_i}, vectors are normalized and orthogonal to each other. This allows us to define the projector,

P=∑j=0N−1∣vj⟩⟨vj∣. P = \sum_{j=0}^{N-1} \ket{v_j} \bra{v_j}.

or conveniently we may fill a matrix with our orthonormal basis, V=[∣v0⟩,∣v1⟩,…vN−1]V = [ |v_0\rangle, |v_1\rangle, \ldots v_{N-1} ], which allows us to write P=VV†P = V V^\dagger.

From here we approximate the time evolution as,

∣ψ(t+δ)⟩≈P†e−iδH^P∣ψ(t)⟩, |\psi(t+\delta) \rangle \approx P^\dagger e^{-i \delta \hat{H}} P |\psi(t) \rangle, ∣ψ(t+δ)⟩≈Ve−δhV†∣ψ(t)⟩, |\psi(t+\delta) \rangle \approx V e^{-\delta h} V^\dagger|\psi(t)\rangle,

where h=V†H^Vh = V^\dagger \hat{H} V is a tri-diagonal matrix. Since NN is typically chosen to be small hh can be efficiently diagonalized and this expression can be simplified. Or alternatively this expression can already be efficiently evaluated with aleph's linear algebra kit.

The error on this approximation scheme is remarkably small, one can look at the bounds for exact arithmetic to see it is exponentially small in the Krylov subspace dimension. Let W=Emax−EminW = E_{\text{max}} - E_{\text{min}} be the spectral width. Then the error in exact arithmetic is O(δN)O(\delta^N) so long as Wδ≤N\sqrt{W\delta} \leq N. To learn more see theorem 2 of this reference.

operator_function​

Operator function allows you to approximate any function of an operator expression with a Krylov method. Let A^\hat{A} be any operator expression generated with aleph, then operator_function allows you to efficiently evaluate operators like,

∣ψ⟩=f(A^)∣ϕ⟩ |\psi\rangle = f(\hat{A}) |\phi\rangle

The syntax for this is simple, you need to give the operator_function an operator, a lambda function to evaluate and a size to the Krylov subspace.

var Nkrylov = 10;
dt = 0.1
// the function captures the time step dt
var expiHt = fun[dt](x) { return exp(-1i * dt * x) }
var U = operator_function(any_operator, expiHt, ["krylov_dimension" : Nkrylov]);

Below we will do a fully worked example.

Example J1-J2 model​

In this example we will use an operator_function to time evolve J1-J2 system starting in a N'eel state. The Hamiltonian is given by (taking periodic boundaries),

H^=∑m=0L−1J1(XmXm+1+YmYm+1+γ1ZmZm+1)+J2(XmXm+2+YmYm+2+γ2ZmZm+2). \hat{H} = \sum_{m=0}^{L-1} J_1(X_m X_{m+1} + Y_m Y_{m+1} + \gamma_1 Z_m Z_{m+1}) + J_2(X_m X_{m+2} + Y_m Y_{m+2} + \gamma_2 Z_m Z_{m+2}).

And the initial state is taken to be,

∣ψ⟩=12(∣↑↓…↑↓⟩+∣↓↑…↓↑⟩). |\psi\rangle = \frac{1}{\sqrt{2}} \left(|\uparrow \downarrow \ldots \uparrow \downarrow\rangle + | \downarrow \uparrow \ldots \downarrow \uparrow \rangle \right).

We will track the evolutions of two observables,

A^0=Z0Z1,A^1=Z0Z2. \hat{A}_0 = Z_0 Z_1, \enspace \hat{A}_1 = Z_0Z_2.

Our goal will be to approximate the average value of these observables in time.

// system size
var L = 18;
// krylov subspace size
var Nkrylov = 6;
// time step
var dt = 0.1;
// number of time steps
// make this bigger to see equilibration
var num_steps = 10;
// model params
var J1 = 1.0;
var J2 = -0.5;
var gamma1 = 0.5;
var gamma2 = 0.5;

// operator summation for our Hamiltonian
var ham = operator_sum(as_complex)

var s1; var s2;
var s3;
// we are doing periodic boundaries
for(m : [0..L])
{
    s1 = m;
    s2 = (m+1)%L
    s3 = (m+2)%L
    // XXPYY is a faster subroutine to do XX+YY calculations
    ham += J1*XXPYY(s1,s2) + (J1*gamma1)*ZZ(s1,s2) + J2*XXPYY(s1,s3) + (J2*gamma2)*ZZ(s1,s3);
}

//print(ham)

// Use operator function to approximate time evolution
// error is exponentially suppressed in Nkrylov 
var expiHt = fun[dt](x) { return exp(-1i * dt * x) }
var U = operator_function(ham, expiHt, ["krylov_dimension" : Nkrylov]);

// we can track these observables for example
var Obs0 = ZZ(0,1);
var Obs1 = ZZ(0,2);

// building a dense vector for our state (Neel state)
var state = zeros([2**L]) // -> 1D complex matrix
var psi_prod = Qbit(L) // -> product state with all bits 0
psi_prod.flip([0..L][::2]) // -> flip only the even ones
state[psi_prod] = 1 // -> adding|101010..> 
psi_prod.translate() // -> shift bits
state[psi_prod] = 1; // -> adding |010101...>
state.normalize(); // normalizing

// container to collect results
var data0 = zeros([num_steps], as_array | as_real);
var data1 = zeros([num_steps], as_array | as_real);
// initial conditions
data0[0] = expval(Obs0,state).real;
data1[0] = expval(Obs1,state).real;

print("${0}, ${data0[0] }, ${data1[0]}")
// time evolving using our operator function
for(m:[1..num_steps])
{
    state = U*state;
    data0[m] = expval(Obs0,state).real;
    data1[m] = expval(Obs1,state).real;
    print("${m*dt}, ${data0[m]}, ${data1[m]}")
}


print("-----------------------------")
// average values
print("${num_steps*dt}, ${data0.mean}, ${data1.mean}")

What's next?​

  • Write a routine that checks to see if the system has equilibrated. This can be done by tracking how much our time average changes if we take different cuts of continuous regions in time and average over the observables expectation value.
  • How big of a system can you simulate? L=18L=18 is close to instant, can you do L=24?L=24? What about L=28L=28? How big of a computer would you need to do, to do L=32L=32?
  • For longer runs or larger sizes, experiment with the job scheduler instead of doing direct execution.
  • How would I check if the equilibrium expectation values we see are thermal expectation values?
  • Can you use the ZZN and ZZNN operators to speedup this calculation?

Documentation Contributors

Jonathon Riddell