8. quimb - qiskit MPS Benchmark

Recent studies comparing different quantum backends, such as the one collected in feniqs_lite, have pointed out that \(\texttt{qiskit-aer}\) may be one of the fastest options for executing and sampling a quantum circuit with classical resources, thanks to its underlying Rust/C++ implementation. Nevertheless, in a simple comparison with \(\texttt{quimb}\), we find that this advantage depends strongly on the circuit architecture.

In this example, we explore the contraction and sampling times for two types of circuits: brickwall and random entangling patterns. Carrying out several simulations with the respective MPS classes in \(\texttt{qiskit-aer}\) and \(\texttt{quimb}\), we find that \(\texttt{quimb}\) becomes faster than \(\texttt{qiskit-aer}\) for deep enough circuits built from long-range gates. The data used in this study were computed on a laptop using a single CPU thread and are loaded from the accompanying JSON files.

import json
import os
from collections import Counter

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
from qiskit.visualization import plot_histogram
from qiskit_aer import AerSimulator

from qpe_toolbox.circuit import (
    deserialize_to_qiskit_QuantumCircuit,
    deserialize_to_quimb_Circuit,
    deserialize_to_quimb_CircuitMPS,
    draw_layered_circuit,
    generate_rand_circuit,
    serialize_from_quimb_Circuit,
)

def plot_errorbar(ax, x, y, yerr, mfc, *, fmt="o", label=None):
    return ax.errorbar(
        x,
        y,
        yerr,
        fmt=fmt,
        mfc=mfc,
        label=label,
        markersize=6,
        markeredgecolor="k",
        color="black",
        ecolor="black",
        elinewidth=0.5,
        capsize=6,
        capthick=0.5,
    )

8.1. Circuit Instances

Let us visualize the smallest instances of the circuits explored throughout this example:

circuit_types = ("brick", "rand")
circuit_data = {}
for ctype in circuit_types:
    with open("precomputed_data/" + ctype + "_circ_perfo.json") as fin:
        circuit_data[ctype] = json.load(fin)

for ctype in circuit_types:
    circ = deserialize_to_quimb_Circuit(circuit_data[ctype][ctype + "_16"])
    depth = max([gate.round for gate in circ.gates]) + 1
    fig = draw_layered_circuit(
        circ,
        state_label=r"$0$",
        labels_1qubit=[f"$\\mathrm{{R_y^{{({j})}} }}$" for j in range(1, depth + 1)],
        labels_2qubit=[r"$\mathrm{CNOT}$"] * depth,
        max_depth=depth,
    )

The next cell contains a short code example demonstrating how we compare simulation timings:

(1) we create new circuits or load existing ones with the functions introduced in the Building Circuits example,

(2) we run simulations (circuit execution + sampling) and access the times and samples,

(3) we compare the results for both backends.

The following code also illustrates how the precomputed data were obtained (note that we fixed the number of threads to 1):

# Create small circuit
depth = 3
circ_quimb = generate_rand_circuit(
    n_qubits=5,
    depth=depth,
    one_qubit_gate_label="rx",
    two_qubit_gate_label="cx",
    two_qubit_gate_range=4,
    two_qubit_gate_prob=0.5,
)

# Recast it into the right classes
circ_dict = serialize_from_quimb_Circuit(circ_quimb)
circ_quimb = deserialize_to_quimb_CircuitMPS(
    full_gate_dict=circ_dict, max_bond=2**depth, cutoff=1e-10
)
circ_qiskit = deserialize_to_qiskit_QuantumCircuit(circ_dict, measure=True)

num_samples = 10**4  # pick the number of samples
# Run qiskit
simulator = AerSimulator(
    method="matrix_product_state",
    matrix_product_state_max_bond_dimension=2**depth,
    matrix_product_state_truncation_threshold=1e-10,
    seed_simulator=1,
)
result = simulator.run(circ_qiskit, shots=num_samples).result()
# qiskit inverts the layout order
counts_qiskit = Counter({k[::-1]: v for k, v in result.get_counts().items()})

# Run quimb
counts_quimb = Counter(circ_quimb.sample(C=num_samples, seed=42))

We can compare the bitstring histogram from both simulations as a sanity check:

plot_histogram(
    [counts_qiskit, counts_quimb],
    legend=["qiskit", "quimb"],
    bar_labels=False,
    figsize=(12, 5),
)

8.2. Contraction and Sampling Times in quimb

In the following we present a non-exhaustive study of the dependence of the MPS (layer-by-layer) contraction time \(t_{\mathrm{contr}}\) and sampling time \(t_{\mathrm{sampl}}\) in \(\texttt{quimb}\). We track how these times grow with the number of qubits n_qubits (length of the MPS) and the depth; note that the depth of the circuit upper bounds the bond dimension \(\chi\) of the fully contracted MPS as \(\chi \leq 2^{\rm depth}\).

We will compare the results for both types of quantum circuit instances introduced above. Since these classes of circuits are both random (in the sense that the parametrization of the one-qubit gates is random), a thorough study would require averaging over circuit architectures and parametrizations. Such a study is beyond the scope of this example; we will see that solid comparisons can be drawn even without averaging over this randomness.

The conditions of the simulations include:

• A maximum bond dimension max_bond \(\chi=2^{\rm depth}\), such that the simulation should not introduce any norm truncation due to bond dimension capping.

• A cutoff parameter, set to \(\varepsilon=10^{-10}\); this cutoff controls the tail of singular values that is cut out after a local gate application.

A more natural prescription for the study than the one chosen here would be to set a maximum accumulated error in the norm \(\varepsilon_{\mathrm{total}}\), and study the growth of \(t_{\mathrm{contr}}\), \(t_{\mathrm{sampl}}\) and \(\chi\) for that fixed threshold in n_qubits and depth. We avoid this prescription due to the difficulty of fixing the maximum error across different backends.

8.2.1. Scaling with the Number of Qubits

In the following plot we present our results in two columns: the left for brickwall circuits, and the right for random circuits. The chosen depths are depth=1, 3, 5, and the \(\texttt{quimb}\) class CircuitPermMPS (labeled PermMPS in the plots) is also used for random circuits with depth=1:

mycolors = ("tab:red", "tab:orange", "tab:olive", "tab:purple", "tab:pink", "tab:cyan")
cases = ("quimb MPS contr", "quimb MPS sampl 1", "quimb MPS sampl 10")
cases_perm = ("quimb PermMPS contr", "quimb PermMPS sampl 1", "quimb PermMPS sampl 10")
dum = np.array([np.nan])
circuit_sizes = 2 ** np.arange(4, 9)
x_lines = np.linspace(circuit_sizes[0] / 2, circuit_sizes[-1] * 2)
ylims = ((5e-4, 1.0), (1e-3, 1.0), (1e-3, 10.0))
fig, axes = plt.subplots(3, 2, sharex=True, figsize=(5 * 2, 5 * 3 / 1.61))
for idx, ax in np.ndenumerate(axes):
    for instance in circuit_data[circuit_types[idx[1]]].values():
        for i, c in enumerate(cases):
            values = instance[f"depth{2 * idx[0] + 1}"][c]
            plot_errorbar(
                ax, instance["n_qubits"], np.mean(values), np.std(values), mycolors[i]
            )
    ax.set_xlim(circuit_sizes[0] / 2, circuit_sizes[-1] * 2)
    ax.set_xscale("log", base=2)
    ax.set_yscale("log")
    ax.set_xticks(circuit_sizes)
    ax.get_xaxis().set_major_formatter(plt.ScalarFormatter())
    ax.grid(which="both", linestyle="-", linewidth=0.25, alpha=0.5)
    ax.set_ylim(ylims[idx[0]])
    for m in range(16):
        ax.plot(x_lines, 10 ** (m / 2) / 8e5 * x_lines, "-", color="grey", alpha=0.1)
for instance in circuit_data[circuit_types[1]].values():
    for i, c in enumerate(cases_perm):
        values = instance["depth1"][c]
        plot_errorbar(
            axes[0, 1],
            instance["n_qubits"] * 1.15,
            np.mean(values),
            np.std(values),
            mycolors[i + 3],
            fmt="s",
        )
for i in range(3):
    axes[i, 0].set_ylabel(rf"$\text{{depth}}={2 * i + 1}$")
    axes[i, 1].set_yticklabels([])
    plot_errorbar(axes[0, 1], dum, dum, dum, mycolors[i], label=cases[i][6:])
    plot_errorbar(
        axes[0, 1], dum, dum, dum, mycolors[i + 3], fmt="s", label=cases_perm[i][6:]
    )
axes[0, 0].set_title("brickwall circuits")
axes[0, 1].set_title("random circuits")
axes[0, 1].legend(loc="lower right", fontsize=8)
axes[2, 0].set_xlabel("number of qubits")
axes[2, 1].set_xlabel("number of qubits")
fig.text(
    0.04, 0.5, "wall-clock time (s)", va="center", rotation="vertical", fontsize=16
)
fig.tight_layout(rect=[0.05, 0.05, 1, 0.95])

Both the contraction and sampling times vary linearly with the number of qubits. Indeed, the circuit classes treated here involve a number of gates proportional to the number of qubits in the register, therefore the circuit contraction process repeats the gate application procedure a number of times proportional to \(n_{\rm qubits}\). This behavior is reflected by the scaling of the red markers. To highlight the scaling, a set of parallel power-law lines \(\propto n_{\rm qubits}\) is added to the grid as a guide to the eye.

Similarly, sampling requires computing probability marginals qubit by qubit to fix each of the \(n_{\rm qubits}\) bits of the outcome. Therefore, the scaling of sampling times with the number of qubits is also \(\propto n_{\rm qubits}\).

Note, though, that sampling a single bitstring is not equivalent to sampling 10 bitstrings in one go: the sample method from \(\texttt{quimb}\) caches the marginal distributions between calls; this means that sampling gets progressively cheaper per bitstring. This effect is reflected by the gap between the orange markers (1 sample) and the yellow markers (10 samples): the latter indicate the average sampling time for drawing 10 bitstrings (thus we are timing sample(C=10)/10 rather than sample(C=1)).

The discussion so far focused on the CircuitMPS class from \(\texttt{quimb}\). Nevertheless, there exists another MPS class that tries to leverage the non-locality of long-range circuits: CircuitPermMPS. The goal of this class can be understood by diving a bit into the procedure for applying a long-range gate: it requires permuting the two involved qubits through the 1-dimensional layout until they are nearest neighbors, followed by the local application of the two-body gate. While CircuitMPS would restore the initial layout after the gate application, CircuitPermMPS leaves this new layout and tracks the intermediate trajectories of the qubits. This class may reduce wall-clock times when entangling gates are sparse. However, our circuit instances are not sparse enough to show a clear advantage. This can be seen for depth=1 in the right column (random entangling pattern). Despite these results, we expect that a study with varying two_qubit_gate_prob and two_qubit_gate_range in the function generate_rand_circuit would reveal a regime where CircuitPermMPS outperforms CircuitMPS for sparse circuits (low two_qubit_gate_prob, medium-low two_qubit_gate_range).

8.2.2. Scaling with the Bond Dimension

In the following plot we present our results in two columns: the left for brickwall circuits, and the right for random circuits. The chosen circuit sizes are \(n_{\rm qubits} \in \{64, 128, 256\}\):

x_range = np.linspace(1, 256)
n_qubits_values = (64, 128, 256)
depths = range(1, 8)
ylims = ((1e-2, 5.0), (1e-2, 1e2), (1e-2, 1.0), (1e-2, 5.0))
mycolors = ("tab:blue", "tab:orange", "tab:green")
fig, axes = plt.subplots(2, 2, sharex=True, figsize=(10, 5 * 2 / 1.61))
for iN, n_qubits in enumerate(n_qubits_values):
    plot_errorbar(
        axes[0, 0], dum, dum, dum, mycolors[iN], label=rf"n_qubits={n_qubits}"
    )
    for j, ctype in enumerate(circuit_types):
        for depth in depths:
            data = circuit_data[ctype][f"{ctype}_{n_qubits}"]
            bond_dim = data[f"depth{depth}"]["quimb MPS bondim"]
            values = (
                data[f"depth{depth}"]["quimb MPS contr"],
                data[f"depth{depth}"]["quimb MPS sampl 1"],
            )
            for k, ax in enumerate(axes[:, j]):
                plot_errorbar(
                    ax,
                    bond_dim,
                    np.mean(values[k]),
                    np.std(values[k]),
                    mycolors[iN],
                )
                if depth == max(depths):
                    y = x_range ** (3 - k) * np.mean(values[k]) / bond_dim ** (3 - k)
                    ax.plot(x_range, y, "-", color="grey", alpha=0.1)
axes = axes.flatten()
for i, ax in enumerate(axes):
    ax.set_xlim(x_range[0], x_range[-1])
    ax.set_xticks(2 ** np.arange(2, 8))
    ax.set_xscale("log", base=2)
    if i > 1:
        ax.set_xlabel("bond dimension")
    ax.get_xaxis().set_major_formatter(plt.FuncFormatter(lambda x, _: f"{int(x)}"))
    ax.set_ylim(ylims[i])
    ax.set_yscale("log")
    ax.grid(which="both", linestyle="-", linewidth=0.25, alpha=0.5)
axes[0].set_title("brickwall circuits")
axes[0].set_ylabel("contraction")
axes[1].set_title("random circuits")
axes[2].set_ylabel("sampling")
axes[0].legend(loc="lower right", fontsize=8)
fig.text(
    0.04, 0.5, "wall-clock time (s)", va="center", rotation="vertical", fontsize=16
)
fig.tight_layout(rect=[0.05, 0.05, 1, 0.95])

To understand the scaling of the contraction and sampling times with the bond dimension, we first recall the expected algorithmic cost in each case.

The contraction procedure involves (as mentioned earlier) the application of a two-body gate onto a pair of local tensors on an MPS, which at a given depth may have a bond dimension \(\chi\leq 2^{\rm depth}\). The algorithmic cost for this contraction is \(\propto \chi^3\), as can easily be checked by dimension arguments. This scaling only becomes visible for the largest bond dimensions in the right column for the random entangling pattern circuits, which at the same time saturate the bound between the circuit depth and the bond dimension: \(\chi=2^{\rm depth}\). As a guide to the eye, we have added the cubic power law that passes through the points with highest \(\chi\).

For the case of sampling a single bitstring, the scaling can be trickier: sampling a single bitstring requires constructing the diagonal of a reduced density matrix, which can be done with cost \(\propto \chi^2\) if we consider only contracting the bra and ket parts of the MPS around the sites where the bitstring outcome is to be fixed. This is actually the scaling observed in the lower row of results for all the instances (see quadratic power law in grey passing through the points with highest \(\chi\)). Note, however, that before contracting the MPS bra and ket, the MPS must be canonicalized. If the cost of such a process were to be included within the sampling time, then the cost would scale as \(\propto \chi^3\), due to the cost of performing rank-revealing matrix decompositions.

8.3. Comparing quimb and qiskit-aer

For the \(\texttt{qiskit-aer}\) backend, the contraction and sampling times cannot be easily split apart; this means that the simulator only returns \(t_{\mathrm{total}}(n_{\mathrm{sampl}})\), which does not necessarily decompose as \(t_{\mathrm{contr}}+n_{\mathrm{sampl}}\cdot t_{\mathrm{sampl}}\). The simplest way to estimate \(t_{\mathrm{contr}}\) would be to repeat the simulation for a fixed circuit architecture and a list of \(n_{\mathrm{sampl}}\) values, and to deduce \(t_{\mathrm{contr}}\) and \(t_{\mathrm{sampl}}\) from a \(n_{\mathrm{sampl}}\to 0\) extrapolation. Nevertheless, when we explored this behavior for a few instances, we found that there exist different regimes where the extrapolation would follow a particular scaling power law in \(n_{\mathrm{sampl}}\):

n_qubits = 32
num_samples = np.rint(np.logspace(np.log10(1), np.log10(5000), num=40)).astype(int)
num_repetitions = 5
depth = 3
qc = deserialize_to_qiskit_QuantumCircuit(
    circuit_data["brick"][f"brick_{n_qubits}"], max_depth=depth + 1, measure=True
)
dict_total_qiskit = {}
for shots in num_samples:
    batch_total_qiskit = np.empty(num_repetitions)
    for r in range(num_repetitions):
        simulator = AerSimulator(
            method="matrix_product_state",
            matrix_product_state_max_bond_dimension=2**depth,
            matrix_product_state_truncation_threshold=1e-10,
            mps_log_data=True,
        )
        result = simulator.run(qc, shots=shots).result()
        batch_total_qiskit[r] = result.time_taken

    dict_total_qiskit[shots] = {
        "mean": np.mean(batch_total_qiskit),
        "std": np.std(batch_total_qiskit),
    }

In the following plot we detect several regimes with different potential extrapolation slopes:

fig, ax = plt.subplots(1, 1, figsize=(6, 6 / 1.61))
for key, instance in dict_total_qiskit.items():
    plot_errorbar(ax, key, instance["mean"], instance["std"], "tab:cyan")
ax.set_ylabel("contr. + sampl. time (s)", fontsize=12)
ax.set_xlabel("number of shots", fontsize=12)
ax.set_xscale("log")
ax.set_yscale("log")
ax.grid(which="both", linestyle="-", linewidth=0.25)

We therefore restrict ourselves to comparing the total wall-clock time to contract each circuit and sample 10 bitstrings in both backends:

n_qubits_values = 2 ** np.arange(4, 9)
fig, axes = plt.subplots(3, 2, sharex=True, figsize=(5 * 2, 5 * 3 / 1.61))
for idx, ax in np.ndenumerate(axes):
    depth = idx[0] + 4
    for instance in circuit_data[circuit_types[idx[1]]].values():
        contr_quimb = np.asarray(instance[f"depth{depth}"]["quimb MPS contr"])
        sampl_quimb = np.asarray(instance[f"depth{depth}"]["quimb MPS sampl 10"])
        total_quimb = contr_quimb + sampl_quimb
        total_qiskit = np.asarray(instance[f"depth{depth}"]["qiskit total"])
        x = instance["n_qubits"]
        plot_errorbar(ax, x, np.mean(total_quimb), np.std(total_quimb), "tab:orange")
        plot_errorbar(ax, x, np.mean(total_qiskit), np.std(total_qiskit), "tab:cyan")

    ax.set_xlim(n_qubits_values[0] / 2, n_qubits_values[-1] * 2)
    ax.set_xscale("log", base=2)
    ax.set_yscale("log")
    ax.set_xticks(n_qubits_values)
    ax.get_xaxis().set_major_formatter(plt.ScalarFormatter())
    ax.grid(which="both", linestyle="-", linewidth=0.25)
    if idx[1] == 0:
        ax.set_ylabel(rf"$\text{{depth}}={depth}$")
    if idx[0] == 2:
        ax.set_xlabel("number of qubits")
plot_errorbar(axes[0, 0], dum, dum, dum, "tab:orange", label="quimb")
plot_errorbar(axes[0, 0], dum, dum, dum, "tab:cyan", label="qiskit-aer")
axes[0, 0].legend(loc="lower right", fontsize=12)
axes[0, 0].set_title("brickwall circuits", fontsize=12)
axes[0, 1].set_title("random circuits", fontsize=12)
fig.text(
    0.04, 0.5, "wall-clock time (s)", va="center", rotation="vertical", fontsize=12
)
fig.tight_layout(rect=[0.05, 0.05, 1, 0.95])

Interestingly, while for all sizes and depths of the chosen brickwall circuits \(\texttt{qiskit-aer}\) dominates with speedups of almost an order of magnitude, the same cannot be said for random circuits with long-range entangling gates. We observe that for deep enough instances (\({\rm depth}\geq 4\)), \(\texttt{qiskit-aer}\) and \(\texttt{quimb}\) already perform comparably; for the deepest random circuits explored here, \(\texttt{quimb}\) surpasses \(\texttt{qiskit-aer}\) by an order of magnitude.

8.4. Other Sampling Strategies in quimb

Although we have found counterexamples to results such as those gathered in feniqs_lite, we are well aware that \(\texttt{quimb}\) can offer even more by introducing more flexible quantum circuit classes (like Circuit) and methods suited for different architectures:

• sample: uses light-cone simplification and caching of the marginal distributions for chosen groupings of qubits.

• sample_gate_by_gate: achieves a cost for sampling similar to that of computing single amplitudes, but increases the number of tensor networks to be contracted. This method was introduced in 2021.

• sample_chaotic: assumes that only subgroups of qubits are correlated, while the rest of the sampled bitstring results are essentially random. This strategy was successful at classically replicating one of the most prominent quantum advantage claims.