top of page
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 suzuki-trotter decomposition of a local Hamiltonian. Our code first demonstrates the use of imaginary time evolution to find the ground state of an infinite (translation-invariant) quantum system in D=1 dimensions.


  • Computational cost: O(d χ^3)


  • Based on a 2-site unit cell (A-B pattern)

  • Makes use of the canonical form for MPS 

Code Language:

TEBD algorithm tutorial (or jump to completed code)

Section 1: The infinite MPS ansatz

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 2-site 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 left-centre-right.

  • 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

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 semi-infinite networks to the left and right of the region of interest. How can such tensors be chosen?


The answer is that the environment tensors σ and μ should be chosen as the dominant fixed points of the contractions shown in Fig.3(i-ii), 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(i-ii) 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 𝛘^2-by-𝛘^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 pre-packaged ARPACK version via the `eigs` function, as demonstrated below.

Section 3: Canonical Form

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 B-A links of the MPS, noting that the A-B 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 B-A 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 eigen-decomposition of the environment tensors σBA and μBA as shown in the code below.  

Section 4: MPS time evolution

We now describe the most important part of the TEBD algorithm: applying the time-evolution gates to the MPS. For simplicity we shall focus on the case of an imaginary-time evolution, although the same approach is used for a real-time 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 time-step, via a Suzuki-Trotter decomposition. Presented below is an example of exponentiating a local Hamiltonian.


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 time-evolution gates across A-B links, and then across B-A links.


  • 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 time-step 𝜏 for the initial iterations, and then reduce to smaller time-steps at later iterations (to improve the final accuracy).  

Anchor 1
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 6e-5

bottom of page