Tutorial - Introduction to Helmi

Estimating the GHZ fidelity

Preparing the GHZ circuit

# Importing the required modules
import os
import networkx as nx
import numpy as np
import itertools
from iqm.qiskit_iqm import IQMProvider
from iqm.qiskit_iqm.fake_backends import fake_adonis
from iqm.qiskit_iqm.iqm_transpilation import optimize_single_qubit_gates
from qiskit.compiler import transpile
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit.visualization import plot_histogram
from qiskit.result import marginal_counts
from functools import reduce
import matplotlib.pyplot as plt
from scipy.optimize import minimize

# Set up the Helmi backend
HELMI_CORTEX_URL = os.getenv('HELMI_CORTEX_URL')
if not HELMI_CORTEX_URL:
	raise ValueError("Environment variable HELMI_CORTEX_URL is not set")

provider = IQMProvider(HELMI_CORTEX_URL) 
backend = provider.get_backend() # fake_adonis.IQMFakeAdonis()
shots = 1024
n_qubits = 5

print(f"Native operations: {backend.operation_names}")
print(f"Number of qubits: {backend.num_qubits}")
print(f"Coupling map: {backend.coupling_map}")

G = nx.Graph()
G.add_edges_from(backend.coupling_map)
node_labels = {node: f"QB{node + 1}" for node in G.nodes}
nx.draw(G, labels=node_labels, node_color="skyblue", node_size=500, font_size=10)

Let us first prepare the GHZ circuit.

# Simple GHZ circuit
qubits: dict[str, QuantumRegister] = {i: QuantumRegister(1, "QB" + str(i)) for i in range(1, 6)}
qc = QuantumCircuit(*qubits.values())
qc.h(qubits[3])
qc.cx(qubits[3],qubits[2])
qc.cx(qubits[3], qubits[4])
qc.cx(qubits[3], qubits[1])
qc.barrier()
qc.cx(qubits[3],qubits[5])
qc.measure_all() #expected equal majority counts of '00000' and '11111'
qc.draw(output="mpl")

Here is how the circuit will look like if we convert to Helmi’s native gateset.

# Optimising for Helmi -> 
initial_layout: list[int] = [backend.qubit_name_to_index(qubit) for qubit in qubits.values()]
transpiled_circuit = transpile(qc, backend, initial_layout= initial_layout, optimization_level=0)
transpiled_circuit.draw(output="mpl")

As the RZ gates commute with the CZ gates, and the measurement is not affected by the final RZ gate, we can optimise those away to decrease the number of gates we perform.

transpiled_circuit = optimize_single_qubit_gates(transpiled_circuit)
transpiled_circuit = transpile(transpiled_circuit, backend, optimization_level=0)
transpiled_circuit.draw(output="mpl")

job = backend.run(transpiled_circuit, shots=shots)
result = job.result()
counts = result.get_counts()

plot_histogram(counts)

print("Unmitigated GHZ fidelity =", (counts["00000"] + counts["11111"])/shots)

Applying readout mitigation

The code used here was modified from this tutorial from NVIDIA.

To perform readout mitigation with the minimum number of jobs, we will prepare the \(\mid \! 00000 \rangle\) and \(\mid \! 11111 \rangle\) and see which states we actually measure.

First, we observe how the \(\mid \! 00000 \rangle\) state is not strongly affected by state preparation and measurement (SPAM) errors.

qubits: dict[str, QuantumRegister] = {i: QuantumRegister(1, "QB" + str(i)) for i in range(1, 6)}
qc_0s = QuantumCircuit(*qubits.values())
qc_0s.measure_all()

initial_layout: list[int] = [backend.qubit_name_to_index(qubit) for qubit in qubits.values()]
qc_0s = transpile(qc_0s, backend, initial_layout = initial_layout, optimization_level=0)

job_0s = backend.run(qc_0s, shots=shots)
result_0s = job_0s.result()
counts_0s = result_0s.get_counts()
plot_histogram(counts_0s)

Next, we can observe how the \(\mid \! 11111 \rangle\) is much more strongly affected by SPAM error. Part of this comes from the \(\mid \! 1 \rangle\) state relaxing to \(\mid \! 0 \rangle\), however, the exact ways in which Helmi is calibrated also plays an important role here.

qubits: dict[str, QuantumRegister] = {i: QuantumRegister(1, "QB" + str(i)) for i in range(1, 6)}
qc_1s = QuantumCircuit(*qubits.values())
for qubit in range(1,6):
    qc_1s.x(qubits[qubit])
qc_1s.measure_all()

initial_layout: list[int] = [backend.qubit_name_to_index(qubit) for qubit in qubits.values()]
qc_1s = transpile(qc_1s, backend, initial_layout = initial_layout, optimization_level=0)

job_1s = backend.run(qc_1s, shots=shots)
result_1s = job_1s.result()
counts_1s = result_1s.get_counts()
plot_histogram(counts_1s)

local_states = ["0" * n_qubits, "1" * n_qubits]
results = {"00000": result_0s, "11111": result_1s}
for state in local_states:
    res = dict(list(results[state].get_counts().items()))
    print(f"{state} becomes {res}")

possible_counts = [
    dict(list(results[state].get_counts().items())) for state in local_states
]
matrices = []

for k in range(n_qubits):
    matrix = np.zeros([2, 2], dtype=float)
    marginalized_counts = []
    total_shots = []
    for i in local_states:
        marginal_cts = marginal_counts(results[i], indices = [k]).get_counts()
        marginalized_counts.append(marginal_cts)
        total_shots.append(sum(marginal_cts.values()))

    # matrix[i][j] is the probability of counting i for expected j
    for i in range(2):
        for j in range(2):
            matrix[i][j] = marginalized_counts[j].get(str(i),
                                                      0) / total_shots[j]
    matrices.append(matrix)

We have just prepared the confusion matrices for each of the \(5\) qubits using just \(2\) jobs, which can be used to calculate the confusion matrix for the whole system. This does assume that measurement or doing gates on one qubit does not affect other gates. We could create these confusion matrices using \(2n_{\text{qubits}}\), or even more properly using \(2^{n_{\text{qubits}}}\) jobs, however, just these \(2\) jobs should suffice.

for i, matrix in enumerate(matrices):
    print(f"Confusion matrix for QB{i+1}:")
    print(matrix, "\n")

labels = list(map(list, itertools.product([0, 1], repeat=n_qubits)))
states = list(map(lambda label: "".join(map(str, label)), labels))

Let us have another look at the counts from the GHZ circuit. We will perform readout error mitigation on it, and hopefully, improve the GHZ state fidelity.

new_counts = dict(counts.items())
noisy_counts = np.array(
    [new_counts.get(state, 0) for i, state in enumerate(states)])
noisy_counts

def find_closest_distribution(empirical_dist):
    """
    Find the closest distribution to an empirical distribution by minimizing the L1 norm.

    Args:
        empirical_dist: Empirical distribution that you want to find the closest distribution to.

    Returns:
        Closest distribution to `empirical_dist`
    """

    def objective(x):
        return np.linalg.norm(empirical_dist - x, ord=1)

    # Constraint: all elements of p must be positive, and the distribution must sum to 1
    cons = (
        {
            "type": "ineq",
            "fun": lambda p: p
        },
        {
            "type": "eq",
            "fun": lambda p: np.sum(p) - 1
        },
    )
    bnds = [(0, 1) for _ in range(len(empirical_dist))]
    initial_value = np.random.uniform(size=len(empirical_dist))

    res = minimize(
        objective,
        initial_value,
        method="SLSQP",
        options={"maxiter": 1000},
        bounds=bnds,
        constraints=cons,
    )

    return res.x

def get_counts_from_distribution(n_qubits, size, dist):
    """
    Generates samples based on a given distribution and returns the counts of each sample value.

    Args:
        n_qubits: The number of qubits in the quantum circuit.
        dist: The probability distribution from which samples are drawn.

    Returns:
        An array of counts for each possible value in the distribution. The array has a length of 2^n_qubits.
    """
    samples = np.random.choice(np.arange(2**n_qubits), size=size, p=dist)
    values, counts = np.unique(samples, return_counts=True)
    res = np.zeros(2**n_qubits, dtype=int)
    res[values] = counts
    return res

# Function to draw the confusion matrix
def plot_cmat(mat):
    fig, ax = plt.subplots()
    n = len(mat)
    im2 = ax.matshow(mat, cmap=plt.cm.Reds, vmin=0, vmax=1.0)
    ax.set_yticks(np.arange(n))
    ax.set_xticks(np.arange(n))
    ax.set_yticklabels(n * [""])
    ax.set_xticklabels(n * [""])
    ax.set_title(r"Confusion Matrix", fontsize=16)
    ax.set_xlabel("Prepared State")
    ax.xaxis.set_label_position("top")
    ax.set_ylabel("Measured State")
    fig.colorbar(im2, ax=ax)
    plt.show()

First, we invert the confusion matrix for each individual qubit. Next, we approximate the inverse confusion matrix for Helmi using the tensor product. Now, we can simply apply this to our “noisy” results and we should have mitigated the effects of the readout errors.

As this is a mathematical operation, we can end up with negative counts, which lack any physical meaning. We can find another distribution with all positive values that is similar to our mitigated distribution. We sample this distribution probabilistically to obtain our corrected counts.

pinv_confusion_matrices = [np.linalg.pinv(m) for m in matrices]
A_pinv = reduce(np.kron, pinv_confusion_matrices)
mitigated = np.array(np.dot(A_pinv, noisy_counts), dtype=int)
print(f"Mitigated counts:\n{mitigated}")

if not np.all(mitigated >= 0):
    positive_dist = find_closest_distribution(mitigated / shots)
    mitigated = get_counts_from_distribution(n_qubits, shots, positive_dist)
    print(f"\nCorrected for negative counts:\n{mitigated}")

A_joint = reduce(np.kron, matrices)
plot_cmat(A_joint)

As we can see, the mitigated fidelity is higher than the unmitigated fidelity, at the cost of two more jobs run.

plot_histogram([counts, dict({bin(i)[2:].zfill(5): x for i, x in enumerate(mitigated)})], legend=["Unmitigated", "Mitigated"])

print("Unmitigated GHZ fidelity =", (counts["00000"] + counts["11111"])/shots)
print("Corrected GHZ fidelity =", (mitigated[0] + mitigated[-1])/shots)

Multiple Quantum Coherences

Theory

Multiple Quantum Coherences (MQC) offers an alternative way to estimate the GHZ fidelity. It allows us to calculate a lower and upper bound on the GHZ fidelity, as well as calculating the exact GHZ fidelity if run alongside a GHZ circuit.

MQC works by utilising phase kickback, a highly versatile tool used in many algorithms to “kick” a phase from a target qubit “back” to a control qubit. First, we prepare a \(N\)-qubit GHZ state, and then we apply a phase \(\phi\) to all \(N\) qubits. Then, we “undo” the GHZ state, and this “kicksback” the phase to the control qubit, applying a phase shift of \(N\phi\) to the control qubit.

To help explain this more easily, I will quote the paper directly here.

Starting from the \(N\)-qubit ground state: \(\mid \! GS \rangle = \mid \! 000..00 \rangle\), apply a Hadamard gate on qubit \(0\) followed by a sequence of CX gates. Ideally this brings the system into the GHZ state: \(\mid \! GHZ \rangle = \tfrac 1{\sqrt 2} \left( \mid \! 000..00 \rangle + \mid \! 111..11 \rangle\right)\)

Apply a collective rotation given by the unitary \(U_{\phi}\) on all qubits. This amounts to adding a phase \(N\phi\) to the GHZ state: \(\tfrac 1{\sqrt 2} \left( \mid \! 000..00 \rangle + e^{-iN\phi}\mid \! 111..11 \rangle\right)\)

Disentangle the GHZ state by performing the CX gate sequence in reverse order. The amplified phase is mapped onto qubit \(0\): \(\tfrac 1{\sqrt 2} \left( \mid \! 000..00 \rangle + e^{-iN\phi}\mid \! 111..11 \rangle\right) \otimes \mid \! 00..00 \rangle\)

Read out the amplified phase by measuring the probability of the system returning to its initial state: \(\mid \! GS \rangle \)> The measured signal of this protocol is given by \(S_\phi = \left| \langle 000..00 \! \mid \! U^{\dagger}_{GHZ} U_{\phi} U_{GHZ} \! \mid \! 000..00 \rangle \right|^2 = Tr(\rho_\phi\rho)\) where \(\rho = U_{GHZ} \! \mid \! GS \rangle \langle GS \! \mid \! U^{\dagger}_{GHZ}\), \(U_{GHZ} = U_{CZ}H_0\), and \(\rho_\phi = U_{\phi}\rho U^{\dagger}_{\phi}\).

If we run this on a noiseless quantum computer, we expect \(S_{\phi} = \tfrac 12(1 + \cos(N\phi))\). We can then measure \(S_{\phi}\) for multiple angles \(\phi\), and then perform a Fourier transform to find which frequency our \(S_{\phi}\) actually corresponds to.

\[I_q=\mathcal{N}^{-1}\left|\sum_\phi e^{i q \phi} S_\phi\right|\]

Ideally, it should be completely dependent on \(N\), and have no reliance on any of the other frequencies - however, due to noise, we shall get a spread of frequencies with varying amplitudes. We can then use these to estimate the lower and upper bounds of the GHZ fidelity.

\[2 \sqrt{I_N} \leq F \leq \sqrt{I_0 / 2}+\sqrt{I_N}\]

Preparing the MQC circuits

First, let us define a parametrised circuit for the MQC experiment. Developed based on the original paper.

def mqc_circuit(angle: float):
    qubits: dict[str, QuantumRegister] = {i: QuantumRegister(1, "QB" + str(i)) for i in range(1, 6)}
    mqc = QuantumCircuit(*qubits.values())
    mqc.h(qubits[3])
    mqc.cx(qubits[3],qubits[2])
    mqc.cx(qubits[3], qubits[4])
    mqc.cx(qubits[3], qubits[1])
    mqc.cx(qubits[3],qubits[5])
    for qubit in qubits.values():
        mqc.x(qubit) # We will comment out this line later
        mqc.rz(angle, qubit)
    mqc.cx(qubits[3],qubits[5])
    mqc.cx(qubits[3],qubits[2])
    mqc.cx(qubits[3], qubits[4])
    mqc.cx(qubits[3], qubits[1])
    mqc.h(qubits[3])
    mqc.measure_all()
    return mqc

Next, we define which \(\phi\) to sweep over, and how many shots to execute for each circuit. We need a minimum of \(12\) experiments here, so that our Fourier transform can detect frequencies upto \(6\).

n_exp = 2 * n_qubits + 2
angles = [2 * np.pi * j / n_exp for j in range(n_exp)]
shots = 16_384

The parameterized circuits are collected in a list mqc_circuits.

mqc_circuits = []

for angle in angles:
    mqc = mqc_circuit(angle)    
    mqc_circuits.append(mqc)

Let us observe how the untranspiled circuit looks.

mqc_circuits[1].draw(output='mpl')

Before running, we transpile each circuit to the native gate set and map the logical qubits to the physical qubits.

initial_layout: list[int] = [backend.qubit_name_to_index(qubit) for qubit in ["QB1", "QB2", "QB3", "QB4", "QB5"]]

fidelity_circuits = transpile(mqc_circuits, backend, initial_layout= initial_layout, optimization_level=2)
fidelity_circuits = optimize_single_qubit_gates(fidelity_circuits)
fidelity_circuits = transpile(fidelity_circuits, backend, initial_layout= initial_layout, optimization_level=0)

Let’s have a look at the transpiled circuit.

fidelity_circuits[1].draw(output='mpl')

With Ideal Backend

fake_backend = fake_adonis.IQMFakeAdonis()

fake_jobs = fake_backend.run(fidelity_circuits, shots=shots)

fake_results = fake_jobs.result()
fake_counts = fake_results.get_counts()

plot_histogram(fake_counts[1])

After running the circuit, we should only have the states \(\mid \! 00000 \rangle\) and \(\mid \! 00100 \rangle\).

fake_outcomes = [fake_count['00000'] / (fake_count['00000'] + fake_count['00100']) for fake_count in fake_counts] # (fake_count['00000'] + fake_count['00100'])

fig, ax = plt.subplots(1, 1, figsize=(10, 6))
all_angles = np.linspace(0, 2 * np.pi, 10_000)
expected_fidelity = [(1 + np.cos(n_qubits * j)) / 2 for j in all_angles]
ax.plot(angles, fake_outcomes, 'o', label='Fake', color='red')
ax.plot(all_angles, expected_fidelity, '-', label='Expected', color='black')
ax.set_xlabel('Angle $(\phi)$')
ax.set_ylabel('$S_{\phi}$')
ax.legend(loc='lower right')
plt.show()

I_0 = 0
I_n = 0
for i, angle in enumerate(angles):
    I_0 += fake_outcomes[i]
    I_n += np.exp(1j * n_qubits * angle) * fake_outcomes[i]
I_0 = np.abs(I_0 / n_exp)
I_n = np.abs(I_n / n_exp)
print("Lower bound for MQC fidelity = " + str(2 * np.sqrt(I_n)))
print("Upper bound for MQC fidelity = " + str(np.sqrt(I_0 / 2) + np.sqrt(I_n)))

Simulations vs Actual Hardware

Let us now rerun this, but comment out the mqc.x(qubit) line at the very start of the circuit. Similarly, try changing the n_exp to a higher value. What changes do you observe?

Spoiler warning

We should not observe any significant changes in the fidelity estimates for the simulator. However, we will now see how this changes when we run on actual hardware.

With Helmi Backend

jobs = backend.run(fidelity_circuits, shots=shots)

jobs.status()

mqc_results = jobs.result()
mqc_counts = mqc_results.get_counts()

mqc_results.timestamps

plot_histogram(mqc_counts[1])

outcomes = [mqc_count['00000'] / (mqc_count['00000'] + mqc_count['00100'])  for mqc_count in mqc_counts] # (mqc_count['00000'] + mqc_count['00100'])

fig, ax = plt.subplots(1, 1, figsize=(10, 6))
all_angles = np.linspace(0, 2 * np.pi, 10_000)
expected_fidelity = [(1 + np.cos(n_qubits * i)) / 2 for i in all_angles]
ax.plot(angles, fake_outcomes, 'o', label='Simulated', color='red')
ax.plot(angles, outcomes, 'o', label='Experimental', color='blue')
ax.plot(all_angles, expected_fidelity, '-', label='Expected', color='black')
ax.set_xlabel('Angle $(\phi)$')
ax.set_ylabel('$S_{\phi}$')
ax.legend(loc='lower right')
plt.show()

I_0 = 0
I_n = 0
for i, angle in enumerate(angles):
    I_0 += outcomes[i]
    I_n += np.exp(1j * n_qubits * angle) * outcomes[i]
I_0 = np.abs(I_0 / n_exp)
I_n = np.abs(I_n / n_exp)
print("Lower bound for MQC fidelity = " + str(2 * np.sqrt(I_n)))
print("Upper bound for MQC fidelity = " + str(np.sqrt(I_0 / 2) + np.sqrt(I_n)))

What we observe is that removing the \(X\) gates results in a noticeable phase shift in the results, even though it should technically not matter. The presenter’s hypothesis is that this phenomenon is related to qubit relaxation, which causes the state \(\mid \! 00000\rangle\) to dominate over \(\mid \! 11111\rangle\), even though they should be even. Consequently, applying the \(X\) gates inverts this distribution, and further relaxation should lead to a more balanced ratio of these states. An analogy can be drawn to the Hahn-Echo experiment, where a \(Z\) gate can help cancel out some of the dephasing effects.

Furthermore, running more experiments provides a higher fidelity. The presenter hypothesises that this can be attributed to the noise after the Fourier transform being spread out over more frequencies, leading to improved results. We are not changing the experiment - just collecting additional data to help average out errors.