Example codes: TEBD for MPS
Please look at the readme page if you have not done so already. Here we present an implementation of the time evolving block decimation (TEBD) algorithm, which implements real or imaginary time evolution of a matrix product state (MPS). This method is based on a suzukitrotter decomposition of a local Hamiltonian. Our code first demonstrates the use of imaginary time evolution to find the ground state of an infinite (translationinvariant) quantum system in D=1 dimensions.

Computational cost: O(d χ^3)

Based on a 2site unit cell (AB pattern)

Makes use of the canonical form for MPS
Code Language:
MATLAB
Julia
Python
TEBD algorithm tutorial (or jump to completed code)
Section 1: The infinite MPS ansatz
Fig.1:
The TEBD algorithm makes use of an infinite MPS, which describes a quantum state on an infinite 1D lattice.

we use an MPS with a 2site unit cell, composed of `A` and `B` tensors, as this is the minimum unit cell compatible with TEBD. An MPS with a larger unit cell could be used if the periodicity of the quantum state that we are trying to describe is greater than 2 sites.

for convenience we allow weights `sBA` and `sAB`, which are restricted to be diagonal matrices, to be positioned on links between the MPS tensors.

by convention, the indices on MPS tensors are ordered leftcentreright.

we call the local dimension as `d` (i.e. the dimension of the Hilbert space associated to each lattice site) and the MPS bond dimension as `𝛘` (although this is also called `m` in some of the literature.
Section 2: Contraction of an infinite MPS
Fig.2:
Key in any tensor network algorithm is the evaluation of local expectation values from the network, like that of <ΨoABΨ> shown in Fig.2(i). This is complicated in the infinite MPS, as we cannot directly contract an infinite network. In practice, as shown Fig. 2(ii), this difficulty can circumvented through the introduction of environment tensors σ and μ, which should be chosen to effectively account for the semiinfinite networks to the left and right of the region of interest. How can such tensors be chosen?
Fig.3:
The answer is that the environment tensors σ and μ should be chosen as the dominant fixed points of the contractions shown in Fig.3(iii), which are known as the (left/right) MPS transfer operators. In practice, a simple method to find the environment tensors is to first initialize random tensors, and to then apply the transfer operators iteratively until they are converged to within a desired tolerance, as demonstrated below.
Notice that Fig.3(iii) is equivalent to an eigenvalue problem: the σ and μ tensors that we are trying to compute are the dominant eigenvectors of the left/right transfer operators (if interpreted as 𝛘^2by𝛘^2 matrices). The iterative method that we employed above is what is known as a power method. Although a valid approach, more efficient iterative methods are known (i.e. which require fewer iterations to converge), such as the Lanczos method. Fortunately we do not have to write our own implementation of the Lanczos method; instead we can employ a prepackaged ARPACK version via the `eigs` function, as demonstrated below.
Fig.4:
Next we describe how to bring the MPS into canonical form via an appropriate gauge transformation. As discussed in tutorials 3 and 4, use of the canonical form simplifies network contractions and also allow for optimal truncations of tensors within a network. Use of the canonical for is thus a key component of the TEBD algorithm.

we shall begin by orthogonalizing only the BA links of the MPS, noting that the AB links can afterwards be handled in an identical manner.

we shall use the `direct orthogonalization` approach discussed in tutorial 3, where the change of gauge is made via X and Y matrices (and their inverses) as shown in Fig.4(i).

Fig.4(ii) represents the contraction of MPS with itself across a BA link. The X and Y gauge change matrices are chosen such that the environment tensors σBA and μBA are be proportional to the identity after the change of gauge as shown.

we find the gauge change matrices via the eigendecomposition of the environment tensors σBA and μBA as shown in the code below.
Section 4: MPS time evolution
Fig.5:
We now describe the most important part of the TEBD algorithm: applying the timeevolution gates to the MPS. For simplicity we shall focus on the case of an imaginarytime evolution, although the same approach is used for a realtime evolution. Starting from an initial MPS Ψ> we wish to construct the MPS for the time evolved state exp(tH)Ψ>. As shown in Fig.5, the exponential of a local Hamiltonian H can be approximated as a product of local gates exp(𝜏h), where 𝜏 is a small timestep, via a SuzukiTrotter decomposition. Presented below is an example of exponentiating a local Hamiltonian.
Fig.6:
As shown in Fig.6, a key step in TEBD is the contraction of a gate into a pair of MPS tensors, followed by splitting back into a new pair of MPS tensors.

the adjacent link weights sBA should be absorbed into the composite tensor before truncation. This ensures that the truncation will minimize the global error as discussed in tutorial 3.

the splitting of the composite tensor is accomplished via a truncated SVD as discussed extensively in tutorial 2.

after the splitting, the adjacent link weights sBA should be again factored out, which can be accomplished by applying their inverses.

some care should be taken to ensure that link weights sBA are not too small, which can otherwise lead to numerical instability due to finite precision arithmetic.
Section 5: Putting everything together
The TEBD algorithm begins with a randomly initialized MPS and then iterates many of the steps described above . Specifically, each iteration involves:

solve for left/right environment tensors (using a Lanzcos method).

bring MPS into canonical form.

apply timeevolution gates across AB links, and then across BA links.
Notes:

solving for the left/right environment tensors and bringing the MPS into canonical form is computationally expensive to do every iteration. In practice it is more efficient to perform multiple steps of applying the time evolution gates between solving for the left/right environment tensors.

it is often most efficient to start with a larger timestep 𝜏 for the initial iterations, and then reduce to smaller timesteps at later iterations (to improve the final accuracy).
Time evolution of MPS (function):
Initialization script:
'mainTEBD' benchmark:
Method: TEBD for an MPS of bond dimension χ = 32
Test problem: 1D quantum XX model (on an infinite lattice)
Running time: approx 70 seconds
Quantities computed: ground energy density
Typical results:
Error in ground energy density (TEBD): approx 6e5