A quick introduction to HSMPO

The Hybrid Stabilizer MPO representation of a quantum circuit combines the efficiency of stabilizers with the flexibility of a tensor network representation. The effectivness of such a representation is more pronounced when the amount of non-Clifford gates in the circuit is small, or when we don’t have too much entanglement stored in the final MPO representation. Let’s code some examples of quantum circuits to practice with this formalism.

A first simple circuit

To start, we need to recall the quantum circuit simulation interface required by Qibo. To prepare a quantum circuit we will need a few Qibo’s modules

import numpy as np
from qibo import Circuit, gates

def construct_circuit(nqubits):
    """Simple function to construct a dummy quantum circuit."""
    # We first define an empty circuit of nqubits
    circ = Circuit(nqubits)
    # We then add some gates
    for _ in range(3):
        for q in range(nqubits):
            circ.add(gates.H(q))
            # We add some random rotations
            if np.random.uniform(0,1) <= 0.8:
                circ.add(gates.RY(q, theta=np.random.uniform(0., 2 * np.pi)))
        for q in range(nqubits):
            # We add some entanglement
            circ.add(gates.CNOT(q0=q%nqubits, q1=(q+1)%nqubits))
    return circ

Now we know how to prepare a quantum circuit using Qibo. Let’s construct its HSMPO.

from mpstab import HSMPO

# Set the number of qubits
nqubits = 9
# Construct a circuit
circ = construct_circuit(nqubits)

# Print the circuit structure
print(circ)

# Construct its hybrid stabilizer MPO
surrogate = HSMPO(circ)

0: ─H─RY─o───────────────X─H─RY─o───────────────X─H─RY─o───────────────X─
1: ─H─RY─X─o─────────────|─H────X─o─────────────|─H─RY─X─o─────────────|─
2: ─H─RY───X─o───────────|─H─RY───X─o───────────|─H─RY───X─o───────────|─
3: ─H─RY─────X─o─────────|─H────────X─o─────────|─H─RY─────X─o─────────|─
4: ─H─RY───────X─o───────|─H─RY───────X─o───────|─H─RY───────X─o───────|─
5: ─H─RY─────────X─o─────|─H─RY─────────X─o─────|─H─RY─────────X─o─────|─
6: ─H─RY───────────X─o───|─H─RY───────────X─o───|─H─RY───────────X─o───|─
7: ─H─RY─────────────X─o─|─H─RY─────────────X─o─|─H─RY─────────────X─o─|─
8: ─H─RY───────────────X─o─H──────────────────X─o─H─RY───────────────X─o─

Now, we may be interested in computing the expectation value of some observable over the state prepared by executing the given quantum circuit.

# Defining an observable
obs_str = "I" * int(nqubits / 2) + "X" + "I" * int(nqubits / 2)
print(f"Chosen observable: {obs_str}")

mpstab_expval = surrogate.expectation(observable=obs_str)
print(f"Computed expectation value: {mpstab_expval}")

Chosen observable: IIIIXIIII
Computed expectation value: 0.06246974547469883

Let’s check if this value matched a statevector simulator, like the ones implemented in Qibo!

# We have implemented a function to convert strings into Qibo observables
from mpstab.utils import obs_string_to_qibo_hamiltonian

qibo_hamiltonian = obs_string_to_qibo_hamiltonian(obs_str)
# Compute expectation value
qibo_expval = qibo_hamiltonian.expectation_from_state(
    circ().state()
)
print(qibo_expval)

0.06246974547470002

This match is expected, since the HSMPO representation is exact in this simulation regime. We can always decide to lighten the MPO representation of the circuit by setting a maximum bond dimension in our tensor network representation! If the approximation is too extreme, we expect the final expectation value carrying some error.

approx_surrogate = HSMPO(circ, max_bond_dimension=2)
approx_mpstab_expval = approx_surrogate.expectation(obs_str)

print(approx_mpstab_expval)

0.21156020226400218

And the error decreases if we set an higher bond dimension.

approx_surrogate = HSMPO(circ, max_bond_dimension=8)
approx_mpstab_expval = approx_surrogate.expectation(obs_str)

print(approx_mpstab_expval)

0.06256210325489374

Simulating a Clifford circuit

When we simulate a circuit which is composed solely by Clifford gates, then the simulation capabilities become the same of a pure stabilizers simulator, allowing us to scale up to hundreds of qubits. As a simple example, let’s prepare the GHZ state.

def ghz_circuit(nqubits):
    """Prepare the GHZ circuit."""
    circ = Circuit(nqubits)
    circ.add(gates.H(0))
    for q in range(nqubits - 1):
        circ.add(gates.CNOT(q, q + 1))
    return circ

# Let's compute the expval of "ZZZ....ZZZ" for 500 qubits
%time
ghz_qubits = 500

ghz_surr = HSMPO(ghz_circuit(ghz_qubits))
ghz_expval = ghz_surr.expectation("Z" * ghz_qubits)
print(ghz_expval)

CPU times: user 2 µs, sys: 1 µs, total: 3 µs
Wall time: 4.29 µs
1.0

In summary, mpstab offers its maximum benefit in regimes where we have a good amount of Clifford gates or a reasonable amount of entanglement stored in our final MPO.