Tutorial - Introduction to Helmi

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)

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

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

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?

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)))