5. QPE with Trotterization¶

In this example we perform Quantum Phase Estimation with first and second-order Trotter decompositions of the time-evolution operator \(U(t) = \exp(-i H t)\).

Previously in the Textbook QPE example, we introduced and ran QPE with an exact matrix representation of \(U\); this is only possible for small systems amenable to exact diagonalization. In general, we use a Trotter approximation to exponentiate the Hamiltonian; see the tutorial on Trotter-Suzuki Decomposition for an introduction to Trotter approximants to exponentials of matrices.

We study the energy precision obtained as a function of the number of phase qubits in the QPE circuits and the number of Trotter steps in the time evolution. We also perform some simple resource analysis: we quantify the number of entangling gates and the time required to simulate the circuits with \(\texttt{quimb}\).

import matplotlib.pyplot as plt
import numpy as np
from tqdm import notebook as tqdm

import qpe_toolbox.estimation as qpe
from qpe_toolbox import EXACT
from qpe_toolbox.circuit import count_gates_by_qb, make_circ, make_circMPS
from qpe_toolbox.hamiltonian import do_dmrg, heisenberg_hamiltonian

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

mystyles = [
    {"linestyle": "--", "linewidth": 1, "color": "black"},
    {"linestyle": ":", "marker": "o", "markersize": 10, "color": "tab:blue"},
    {"linestyle": ":", "marker": "s", "markersize": 8, "color": "tab:orange"},
    {"linestyle": ":", "marker": "v", "markersize": 8, "color": "tab:green"},
    {"linestyle": ":", "marker": "*", "markersize": 10, "color": "tab:red"},
    {"linestyle": ":", "marker": "^", "markersize": 8, "color": "tab:brown"},
]

5.1. A First Trotterized Run¶

We consider the 1D nearest-neighbor Heisenberg Hamiltonian with open boundary conditions and take the DMRG ground state as the “exact” reference state.

n_qubits = 4
h_spin = heisenberg_hamiltonian(n_qubits)
exact_energy, psi0_mps = do_dmrg(h_spin)

Then we set the parameters needed to compute the energy with QPE. The QPE circuit output is a phase \(2 \pi \theta\). We need to set an appropriate global phase and total evolution time to make sure we recover the correct energy value from the output \(\theta\) (see the Textbook QPE example).
Note that n_trotter_steps is the number of Trotter steps used to decompose the first controlled time evolution, over the interval \(t\) (evolution_time). As the circuit proceeds, we apply time evolution over an exponentially growing time \(2^k t\) controlled by the \(k\)-th phase qubit; the number of Trotter steps grows accordingly as \(2^k\) so as to keep the Trotter timestep constant.

E_target = exact_energy + 0.2
size_interval = 2
Emax = E_target + size_interval / 2
evolution_time = 2 * np.pi / size_interval
global_phase = Emax * evolution_time

print(f"exact theta = {(Emax - exact_energy) / size_interval}")

n_phase_bits0 = 2
n_trotter_steps0 = 1
print(
    f"Precision in theta = {1 / 2**n_phase_bits0}, precision in energy = {1 / 2**n_phase_bits0 * size_interval}\n"
)

circuit0 = make_circ(n_phase_bits0, psi0_mps)

exact theta = 0.6
Precision in theta = 0.25, precision in energy = 0.5

\(\texttt{quimb}\) represents the QPE circuit as a tensor network:

dt = evolution_time / n_trotter_steps0
_, circ = qpe.qpe_first_stage(
    h_spin, circuit0, evolution_time, dt, global_phase, trotter_order=1
)

phase_reg = list(range(n_phase_bits0))
circ.apply_gates(qpe.iqft_swapped(phase_reg))

circ.psi.draw(
    figsize=(12, 12),
    show_tags=False,
    color={f"I{i}" for i in range(n_qubits + n_phase_bits0)},
    edge_scale=1,
    layout="kamada_kawai",
    edge_color=True,
    legend=False,
)

../../_images/b04e36aadd1f9c01acb5bcf9981a07eb80ad271d641cbbb0a3d2fd3bc7bc1d07.png

We compute the energy and check that the error is below the theoretical bound:

traces, energy = qpe.qpe_energy(
    h_spin,
    circuit0,
    n_trotter_steps0,
    E_target,
    size_interval,
    trotter_order=1,
    verbosity=1,
)

print(f"\nBest guess = {energy:.4f} with probability {traces['prob']:.4f}")
print(f"exact energy = {exact_energy:.4f}")

# check error is below QPE upper bound
print(
    f"error = {abs(exact_energy - energy):.4f} < error bound = {size_interval / 2**n_phase_bits0:.4f}"
)

10 |2> 0.5  0.6581

01 |1> 0.25 0.1777

00 |0> 0.0  0.0936

11 |3> 0.75 0.0706

Best guess = -1.4160 with probability 0.6581
exact energy = -1.6160
error = 0.2000 < error bound = 0.5000

5.2. Influence of Phase Qubits and Trotter Steps¶

We run several QPE circuits, varying the number of phase qubits and Trotter steps, and measure the energy. We also record the gate count and simulation time, and compare two modes of circuit simulation available with \(\texttt{quimb}\): the Circuit mode and the CircuitMPS mode. For the Circuit mode we choose a greedy contraction-path optimizer.

Let us first run the circuits (this may take a couple of minutes):

res = {"durations_tn": [], "durations_mps": [], "energies": [], "entangling_gates": []}
trotter_order = 2
nphase_list = np.array([1, 2, 3, 4, 5])
ns_list = [1, 2, 3, 4, EXACT]
max_tn_nphase = 5
max_tn_nsteps = 3

for n_trotter_steps in tqdm.tqdm(ns_list):
    duration_tn = []
    duration_mps = []
    energies = []
    gates_count = []
    for n_phase_bits in tqdm.tqdm(nphase_list, leave=False):
        # contracting the circuit in MPS mode is much more efficient than generic TN
        # deeper, wider circuits can be classically simulated
        initial_circMPS = make_circMPS(n_phase_bits, psi0_mps)
        traces, energy = qpe.qpe_energy(
            h_spin,
            initial_circMPS,
            n_trotter_steps,
            E_target,
            size_interval,
            trotter_order=trotter_order,
        )
        energies.append(energy)
        count = count_gates_by_qb(traces["gates_count"])
        gates_count.append(count["2qb"] + count["3+qb"])
        duration_mps.append(traces["ctimes"][-1])

        if (n_trotter_steps is not EXACT) and (
            n_phase_bits < max_tn_nphase or n_trotter_steps < max_tn_nsteps
        ):
            # using generic tensor network contraction to simulate a quantum circuit
            # is usually very expensive. Only try for small numbers of qubits.
            initial_circ = make_circ(n_phase_bits, psi0_mps)
            traces_tn, energy_tn = qpe.qpe_energy(
                h_spin,
                initial_circ,
                n_trotter_steps,
                E_target,
                size_interval,
                trotter_order=trotter_order,
                optimize="greedy",
            )
            duration_tn.append(traces_tn["ctimes"][-1])
            assert abs(energy_tn - energy) < 1e-6

    res["durations_tn"].append(duration_tn)
    res["energies"].append(np.array(energies))
    res["durations_mps"].append(duration_mps)
    res["entangling_gates"].append(gates_count)

5.2.1. Error Versus Number of Phase Qubits¶

Let us see how the precision improves when adding more qubits in the phase register, for different numbers of Trotter steps:

fig, ax = plt.subplots()
ax.plot(
    nphase_list,
    size_interval / 2**nphase_list,
    label="$\\Delta/2^{n_{\\rm phase~bits}}$",
    **mystyles[0],
)
ax.plot(
    nphase_list,
    np.abs(res["energies"][-1] - exact_energy),
    label="exact time evolution",
    **mystyles[1],
)
for i, n_trotter_steps in enumerate(ns_list[:-1]):
    ax.plot(
        nphase_list,
        np.abs(res["energies"][i] - exact_energy),
        label=f"{n_trotter_steps} Trotter steps",
        **mystyles[i + 2],
    )
ax.legend(loc="lower right", bbox_to_anchor=(1.5, 0.4))
ax.set_ylim(0, 0.25)
ax.set_xticks(range(1, 6))
ax.set_xlabel("number of phase qubits")
ax.set_ylabel("energy error");

../../_images/037c9045c665d81f362f6716a61c306d0e28a6c98e83fc51031ae6f1a44de25a.png

\(\Delta=\) size_interval is the size of the search window \([E_{\rm target} - \Delta/2, E_{\rm target} + \Delta/2]\) where we expect to find \(E_0\), see the Textbook QPE example. The black dashed curve \(\Delta/2^{n_{\rm phase~bits}}\) is the upper bound on the QPE error.

For n_trotter_steps >= 3, the Trotter error becomes negligible and the curves overlap with the exact-time-evolution curve.

5.2.2. Error Versus Number of Trotter Steps¶

We now plot the error as a function of the number of Trotter steps \(n_{\rm steps}\) for a given number of phase qubits:

energies = np.array([res["energies"][i][-1] for i in range(len(ns_list) - 1)])

fig, ax = plt.subplots()
xfit = np.linspace(1, 5, 41)
ax.plot(
    xfit,
    abs(energies[0] - exact_energy) * (ns_list[1] / xfit) ** 2,
    linestyle="--",
    color="r",
    label="$\\mathcal{O}(t^3/n_{\\rm steps}^2)$",
)
ax.axhline(
    y=size_interval / 2 ** nphase_list[-1],
    linestyle=":",
    color="k",
    label="$\\Delta/2^{n_{\\rm phase~bits}}$",
)
ax.plot(
    ns_list[:-1],
    np.abs(energies - exact_energy),
    label=f"{nphase_list[-1]} phase qubits",
    **mystyles[1],
)
ax.legend(loc="upper right", framealpha=1)
ax.set_xticks(range(1, 5))
ax.set_xlim(0.8, 4.2)
ax.set_ylim(0, 0.14)
ax.set_xlabel("number of Trotter steps")
ax.set_ylabel("energy error");

../../_images/99eb94fe0183b9b6da2cf15ebf6f8c365eb1a08fc2426713126461c1cd26d192.png

The red dashed line shows the second-order Trotter error bound \(t^3 / n_{\rm steps}^2\). The black dashed line shows the theoretical QPE error bound \(\Delta / 2^{n_{\rm phase~bits}}\). Here \(2\) Trotter steps are sufficient to reach the QPE precision for this number of phase qubits.

5.2.3. Gate Count and Computation Time¶

Let us now visualize the gate count. We only count entangling gates, i.e. multi-qubit gates. The Trotterized time-evolution operator is implemented with CNOT gates and one-qubit gates (rotations). In QPE, the time evolution is controlled on the phase register. Hence, each CNOT gate acquires an additional control and becomes a CCNOT (Toffoli) gate, and the rotation gates become controlled rotations. The number of Trotter steps in the controlled time-evolution sequence grows exponentially with the number of phase qubits, and so does the number of entangling gates.

fig, (ax_n, ax_t) = plt.subplots(1, 2, figsize=(12, 4))

for i, n_trotter_steps in enumerate(ns_list[:-1]):
    ax_n.plot(
        nphase_list,
        res["entangling_gates"][i],
        label=f"{n_trotter_steps} steps",
        **mystyles[i + 2],
    )
    ax_t.plot(
        nphase_list,
        res["durations_mps"][i],
        label=f"{n_trotter_steps} steps MPS",
        **mystyles[i + 2],
    )

# Plot TN timings on top
for i, n_trotter_steps in enumerate(ns_list[:-1]):
    ax_t.plot(
        nphase_list[nphase_list < max_tn_nphase]
        if n_trotter_steps >= max_tn_nsteps
        else nphase_list,
        res["durations_tn"][i],
        label=f"{n_trotter_steps} steps TN",
        marker="P",
        markersize=8,
        linestyle="--",
        markeredgecolor="k",
        color=mystyles[i + 2]["color"],
    )

ax_n.legend(fontsize=12)
ax_t.legend()
ax_n.set_xlabel("number of phase qubits")
ax_n.set_ylabel("number of entangling gates")
ax_t.set_xlabel("number of phase qubits")
ax_t.set_ylabel("computation time (seconds)");

../../_images/9b7eb869051d2234f064ddb91a010f07f89fee256e2e25d0d99121f25a1c15c7.png

The CircuitMPS mode is much more efficient than the Circuit mode (actual timing depends on the contraction order found by the optimizer). The computation time is directly correlated with the number of entangling gates, which grows exponentially with the number of phase qubits.
We quickly reach computation times of tens of seconds due to the exponentially growing circuit depth. Recall that for such small systems, the Hamiltonian can be exactly diagonalized in a fraction of a second on any laptop. For QPE to gain an advantage over exact diagonalization or DMRG, larger systems, with more than 30 qubits and strong correlations, must be considered.
Note that a comparison between \(\texttt{quimb}\)’s CircuitMPS and \(\texttt{qiskit}\) MPS can be found in the Performance MPS example.

5.3. Resource Analysis¶

For larger numbers of phase qubits, computation becomes costly: in the following we do not run the simulation but only generate the list of gate instructions.

_, evolution_time_resource, global_phase_resource = qpe.set_search_window(
    h_spin, E_target, size_interval
)

nphase_list_resource = [1, 3, 5, 7, 9, 11]
n_trotter_steps_resource = 4

entangling_gates = []
for n_phase_bits in nphase_list_resource:
    initial_circ = make_circ(n_phase_bits, psi0_mps)
    dt = evolution_time_resource / n_trotter_steps_resource

    traces_resource, res_resource = qpe.qpe_sample(
        h_spin,
        initial_circ,
        evolution_time_resource,
        dt,
        global_phase_resource,
        trotter_order=2,
        run_simulation=False,
    )
    count_qb = count_gates_by_qb(traces_resource["gates_count"])
    entangling_gates.append(count_qb["2qb"] + count_qb["3+qb"])

fig, ax = plt.subplots()
ax.semilogy(nphase_list_resource, entangling_gates, "-o")
ax.set_xlabel("number of phase qubits")
ax.set_ylabel("number of entangling gates")
fig.suptitle(
    f"QPE Heisenberg {n_qubits} spins, Trotter order 2 with $dt=t/{n_trotter_steps_resource}$"
);

../../_images/33c787e7ec1d2f908f6c236601b374a4ea4ed01c5f4c2d87b3072b479c3dd4e1.png

In the resource analysis mode (when run_simulation=False) the output is a list of quimb.tensor.Gate objects storing the details of the quantum circuit gates.

print("(label, params, qubits, controls)")
print(*res_resource[:5], sep="\n")

(label, params, qubits, controls)
<Gate(label=H, params=[], qubits=(0,), round=0)>
<Gate(label=H, params=[], qubits=(1,), round=0)>
<Gate(label=H, params=[], qubits=(2,), round=0)>
<Gate(label=H, params=[], qubits=(3,), round=0)>
<Gate(label=H, params=[], qubits=(4,), round=0)>

As a final remark, note that in the process of releasing qpe-toolbox we became aware of recent related work on numerical simulations of textbook QPE on a \(3\)-qubit Heisenberg Hamiltonian using qiskit: arxiv:2602.22349.