9. MPO to circuit transpilation

This notebook reproduces the algorithms of Causer et al., which itself builds on the method of Vidal (note that we reference the first version because it is substantially different from the later versions and contains the appropriate information to follow the procedure).

Other teams built on top of it, adding variants and state preparation, in Kanno et al., 2024 and Kanno et al., 2026.

---

The goal of this notebook is to illustrate the transpilation of an MPO unitary operator (\(\mathrm{M}_{\mathrm{ref}}\) representing the unitary time evolution induced by a Hamiltonian during a time \(\Delta t\)) into a nearest-neighbor brickwall circuit that can be run on a QPU (\(U_{\mathrm{bw}}\)):

We impose an Ansatz circuit made out of three rows of entangling gates on even-odd-even links. Note that in some references the prescription for increasing the depth/layer counting is just a row of even or of odd entangling gates, while other references consider that two consecutive rows even-odd constitute a layer

In order to translate the circuit that best reproduces the action of the unitary MPO, we maximize the overlap between the two unitaries. If this cost function is maximized, then \(U_{\mathrm{bw}}\) will act on another state or operator in the same way as \(\mathrm{M}_{\mathrm{ref}}\) does.

In this case, the cost function can be computed as a fully contracted tensor network with a cylindrical topology. An illustration of such a cost function for the previously introduced ansatz circuit is:

import os

os.environ["NUMBA_NUM_THREADS"] = "1"
os.environ["OMP_NUM_THREADS"] = "1"
os.environ["MKL_NUM_THREADS"] = "1"

import matplotlib.pyplot as plt
import numpy as np
import quimb.tensor as qtn

from qpe_toolbox.circuit import init_cost_tn, transpile_mpo_to_circuit
from qpe_toolbox.hamiltonian import Hamiltonian, trotter_approx_as_MPO
from qpe_toolbox.tensor import state_preparation_mpo

plt.rcParams.update({"font.size": 12})

9.1. Dynamics Induced by the Next-Nearest-Neighbor Ising Model

L = 11
gx, gzz, gz1z = 0.2, 0.5, 0.1
terms_NNIM = []
for x in range(L):
    terms_NNIM.append((gx, "x", [x]))
for x in range(L - 1):
    terms_NNIM.append((gzz, "zz", [x, x + 1]))
for x in range(L - 2):
    terms_NNIM.append((gz1z, "zz", [x, x + 2]))

ham_NNIM = Hamiltonian(terms_NNIM, L)

trotter_mpo_ham_NNIM = trotter_approx_as_MPO(
    ham_NNIM,
    0.5,
    trotter_order=4,
    cutoff=1e-12,
    max_bond=64,
    verbosity=1,
)

Building 4th order Trotter
  Building 2nd order Trotter (1st and 3rd layers)
  Building 2nd order Trotter
    Building 1st order Trotter (1st half)
    Building 1st order Trotter (2nd half)

  Building 2nd order Trotter (2nd layer)
  Building 2nd order Trotter
    Building 1st order Trotter (1st half)
    Building 1st order Trotter (2nd half)
  Multiplying the 3 MPO layers
    Final bond dimension: 12

We now initialize the cost function.

depth = 3

for boundary_bool in [False, True]:
    cost_tn = init_cost_tn(
        ref_mpo=trotter_mpo_ham_NNIM,
        depth=depth,
        param_scaling=1e-1,
        closed=boundary_bool,
        rng=np.random.default_rng(42),
    )
    cost_tn.draw([f"ROUND_{i}" for i in range(depth)], show_inds=True, show_tags=False)

Now we optimize the overlap as suggested in the context of algorithms for entanglement renormalization: to solve the Constrained Linear problem where the isometry to be found decomposes in a given circuit structure, we solve a series of Unconstrained Linear problems for each of the gates forming the circuit. The latter is known to have an analytical solution.

The procedure goes as follows: first, we select a given qubit (say, site 0) and target the optimization of all the gates exposed to that qubit. In the following picture we highlight in yellow the set of gates we are referring to:

Second, we contract the subnetwork that is not being updated. That contraction is the right environment of qubit 0, \(R_0\); such an environment can be built recursively, starting from the right environment of qubit 1, \(R_1\), and so on until the last qubit:

By knowing the first right environment (that of the second-to-last qubit, \(R_{\mathrm{n_qubits-1}}\), which coincides with the last tensor of \(\mathrm{M}_{\mathrm{ref}}\) ) and the tensors required to transfer from one environment to another, we can build all the right environments for the first sweep. Similarly, this can be done for left environments:

Once the environment structure is found, we need to find the optimal gate to update a particular two-qubit unitary. To do so, we contract all the tensors around it. For example, say we want to update the first unitary at depth 0 between qubits 0 and 1, \(U^{(0)}_{0,1}\); we contract all the tensors around it into its environment \(E^{(0)}_{0,1}\) (this \(E\) environment contains the \(L\) and \(R\) of each site, together with the corresponding tensor from the reference unitary at that site).

The environment \(E\) is not necessarily a unitary object. To re-introduce unitarity, we can compute the SVD decomposition of it and remove the singular values. This substitution \(E^{(0)}_{0,1} \longrightarrow U^{\dagger (0) \mathrm{opt}}_{0,1}=u v^\dagger\) is indeed the exact analytical solution of the Unconstrained Linear problem:

Note that the single-gate update above is the same update used by the generic tn_fit routine in the circuit_preparation_opt tutorial’s Local Optimization section. tn_fit handles an arbitrary tensor-network topology at the cost of fully re-contracting each tensor’s environment from scratch on every sweep. On the other hand, this notebook is specialized for a 1D chain and caches the left/right partial contractions (contracted_envs), updating only the one environment adjacent to each optimized gate.

Causer et al. find that the model is prone to get stuck on local minima, even when starting from different initial circuits. We will check this by running the same optimization with different seeds of the Ansatz:

rtol = 1e-5
n_sweeps_max = 200
n_seeds = 4
depth = 5
overlaps5 = np.empty(n_seeds)
ss = np.random.SeedSequence(12)

for i, seed in enumerate(ss.spawn(n_seeds)):
    cost_tn, contracted_envs, overlaps5[i] = transpile_mpo_to_circuit(
        trotter_mpo_ham_NNIM,
        depth,
        rtol,
        n_sweeps_max,
        param_scaling=1e-1,
        closed=True,
        rng=np.random.default_rng(seed),
    )
    print(f"seed key {seed.spawn_key}: overlap = {overlaps5[i]:.6f}")

seed key (0,): overlap = 0.999351

seed key (1,): overlap = 0.999314

seed key (2,): overlap = 0.999363

seed key (3,): overlap = 0.999411

We observe some dependence on the initial state; however, for such a simple problem the effect is small and the algorithm converges to very high fidelity. Causer et al. overcome the local minimum issue by designing a circuit Ansatz that looks like the second-order Trotter expansion of the circuit, where some SWAPs are held fixed and only the remaining gates need to be optimized.

Let us now consider different depths:

depths = np.arange(1, 6)
n_seeds = 4
overlaps = np.zeros((depths.size, n_seeds))
ss = np.random.SeedSequence(42)
rngs = [np.random.default_rng(seed) for seed in ss.spawn(n_seeds)]

for i_depth in range(depths.size - 1):
    for i_seed in range(n_seeds):
        cost_tn, contracted_envs, overlaps[i_depth, i_seed] = transpile_mpo_to_circuit(
            trotter_mpo_ham_NNIM,
            depths[i_depth],
            rtol,
            n_sweeps_max,
            param_scaling=1e-1,
            closed=True,
            rng=rngs[i_seed],
        )
    print(f"depth = {depths[i_depth]}: overlaps = {overlaps[i_depth]}")

overlaps[-1] = overlaps5
print(f"depth = {depths[-1]}: overlaps = {overlaps5}")

depth = 1: overlaps = [0.98876218 0.98876218 0.98876218 0.98876218]

depth = 2: overlaps = [0.98882621 0.98882578 0.98885796 0.98884728]

depth = 3: overlaps = [0.98888903 0.99400312 0.99404461 0.98887494]

depth = 4: overlaps = [0.99703867 0.99705265 0.99919038 0.99920812]
depth = 5: overlaps = [0.99935092 0.99931353 0.99936344 0.99941055]

fig, ax = plt.subplots()
ax.plot(depths, 1 - overlaps, ls="", color="k", marker="o", ms=4)
ax.plot(depths, 1 - overlaps.max(axis=1), "rv-", ms=6)
ax.set_ylim(0, 0.012)
ax.grid(visible=True, alpha=0.3)
ax.set_xlabel("depth")
ax.set_ylabel("1 - overlap")
ax.set_xticks(depths)
ax.set_title(f"MPO transpilation L = {L}");

9.2. State Preparation

With the same cost function we can target also the problem of finding a good preparation circuit for some initial state in the form of an MPS. We just need to build \(\text{M}_{\mathrm{ref}}\) to represent the transition from an empty quantum register into the target MPS \(|0 \rangle^{\otimes L} \langle \Psi_{\mathrm{ref}} |\).

Note that in this case, the cost function will unfold into a square tensor network with open boundary conditions, since the outer product used to build the reference MPO is just an encoding of the tensors.

# Let us pick a target state: the ground state of the next-nearest-neighbor Ising Hamiltonian defined above, in MPS form.
ham_NNIM_mpo = ham_NNIM.to_mpo()
p0 = qtn.MPS_rand_state(L, bond_dim=2, seed=42)
dmrg = qtn.DMRG2(ham_NNIM_mpo, p0=p0)
dmrg.solve(max_sweeps=16, bond_dims=64, verbosity=1, cutoffs=1e-12)
GS = dmrg.state

# build an MPO by taking the outer product with an empty register
GS_mpo = state_preparation_mpo(state_mps=GS)

1, R, max_bond=(2/64), cutoff:1e-12

Energy: -4.4333778147802345 ... not converged.

2, R, max_bond=(4/64), cutoff:1e-12

Energy: -4.433465142686813 ... converged!

# optimize the reference MPO using the same routine as above
n_sweeps_max = 20
rtol = 1e-6
state_depths = np.arange(1, 6)
n_seeds = 4
state_overlaps = np.empty((state_depths.size, n_seeds))
for i, depth in enumerate(state_depths):
    for seed in range(n_seeds):
        rng = np.random.default_rng(seed)
        cost_tn, contracted_envs, state_overlaps[i, seed] = transpile_mpo_to_circuit(
            GS_mpo, depth, rtol, n_sweeps_max, param_scaling=1e-1, closed=True, rng=rng
        )
    print(f"depth {depth}: overlaps = {state_overlaps[i]}")

depth 1: overlaps = [0.99943381 0.99943381 0.99943381 0.99943381]

depth 2: overlaps = [0.99984926 0.99985794 0.99986034 0.9998426 ]

depth 3: overlaps = [0.99988782 0.99988774 0.99988699 0.99987628]

depth 4: overlaps = [0.99988088 0.99989366 0.99988383 0.99989819]

depth 5: overlaps = [0.99989283 0.99989003 0.99988806 0.9998929 ]

The problem is simpler than MPO transpilation above, so it is unsurprising that we reach very high fidelity with lower runtimes. We could raise it slightly more with a lower rtol, at the cost of more iterations / longer runtimes.

fig, ax = plt.subplots()
ax.plot(state_depths, 1 - state_overlaps, ls="", color="k", marker="o", ms=4)
ax.plot(state_depths, 1 - state_overlaps.max(axis=1), "rv-", ms=6)
ax.set_xlabel("depth")
ax.set_ylabel("1 - overlap")
ax.grid(visible=True, alpha=0.3)
ax.set_xticks(state_depths)
ax.set_ylim(0, 6e-4)
ax.set_title(f"State Preparation L = {L}");