A First Glimpse At IBM's Quantum Open Science Price Challenge

It's a surprise to see what the actual problem is

IBM just announced its second Quantum Open Science Prize. They ask for a solution to a quantum simulation problem. They want us to simulate a Heisenberg model Hamiltonian for a three-particle system on IBM Quantum’s 7-qubit Jakarta system using Trotterization.

Even though IBM explains what a Heisenberg model Hamiltonian and Trotterization are, it all appears mysterious, unless you are a physicist or a quantum computing senior.

Fortunately, they also provide a working code example of what they expect. So, let’s see what we computer scientists can learn about the challenge.

First, you need a working Jupyter environment (see my previous post). Second, download the Jupyter notebook (source). We are working on a copy of it because we will make some changes. You can find mine here.

So, let’s go through the code very briefly. First, I stripped away anything unnecessary, such as the classical explanation in section 1.

So, we start with some imports. In line 6, I added the Aer package.

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 16})  # enlarge matplotlib fonts

# Import qubit states Zero (|0>) and One (|1>), and Pauli operators (X, Y, Z)
from qiskit.opflow import Zero, One, I, X, Y, Z
from qiskit import QuantumCircuit, QuantumRegister, IBMQ, execute, transpile, Aer
from qiskit.providers.aer import QasmSimulator
from qiskit.tools.monitor import job_monitor
from qiskit.circuit import Parameter

# Import state tomography modules
from qiskit.ignis.verification.tomography import state_tomography_circuits, StateTomographyFitter
from qiskit.quantum_info import state_fidelity

# Suppress warnings
import warnings
warnings.filterwarnings('ignore')

The original code tries to connect to your IBM account and the Jakarta backend in the next step. We skip this and continue directly with the following few cells defining the custom trotterization gate. It is worth noticing that these cells are surrounded with the comments “YOUR TROTTERIZATION GOES HERE — START (beginning of example)” and “FINISH (end of example).”

Apparently, IBM wants us to replace this part with our solution. But, for now, let’s stick to their example.

# YOUR TROTTERIZATION GOES HERE -- START (beginning of example)

# Parameterize variable t to be evaluated at t=pi later
t = Parameter('t')

# Build a subcircuit for XX(t) two-qubit gate
XX_qr = QuantumRegister(2)
XX_qc = QuantumCircuit(XX_qr, name='XX')

XX_qc.ry(np.pi/2,[0,1])
XX_qc.cnot(0,1)
XX_qc.rz(2 * t, 1)
XX_qc.cnot(0,1)
XX_qc.ry(-np.pi/2,[0,1])

# Convert custom quantum circuit into a gate
XX = XX_qc.to_instruction()

# Build a subcircuit for YY(t) two-qubit gate
YY_qr = QuantumRegister(2)
YY_qc = QuantumCircuit(YY_qr, name='YY')

YY_qc.rx(np.pi/2,[0,1])
YY_qc.cnot(0,1)
YY_qc.rz(2 * t, 1)
YY_qc.cnot(0,1)
YY_qc.rx(-np.pi/2,[0,1])

# Convert custom quantum circuit into a gate
YY = YY_qc.to_instruction()

# Build a subcircuit for ZZ(t) two-qubit gate
ZZ_qr = QuantumRegister(2)
ZZ_qc = QuantumCircuit(ZZ_qr, name='ZZ')

ZZ_qc.cnot(0,1)
ZZ_qc.rz(2 * t, 1)
ZZ_qc.cnot(0,1)

# Convert custom quantum circuit into a gate
ZZ = ZZ_qc.to_instruction()

# Combine subcircuits into a single multiqubit gate representing a single trotter step
num_qubits = 3

Trot_qr = QuantumRegister(num_qubits)
Trot_qc = QuantumCircuit(Trot_qr, name='Trot')

for i in range(0, num_qubits - 1):
    Trot_qc.append(ZZ, [Trot_qr[i], Trot_qr[i+1]])
    Trot_qc.append(YY, [Trot_qr[i], Trot_qr[i+1]])
    Trot_qc.append(XX, [Trot_qr[i], Trot_qr[i+1]])

# Convert custom quantum circuit into a gate
Trot_gate = Trot_qc.to_instruction()

# YOUR TROTTERIZATION GOES HERE -- FINISH (end of example)

We don’t yet look into the details of the trotterization, but we continue with the quantum circuit. In the following code, we define the overall quantum circuit (qc) and generate the state tomography circuits to evaluate the fidelity of the simulation. Fidelity is the performance score that we aim to optimize. Fidelity ranges from 0.0 to 1.0, with 1.0 being the best result we could achieve.

# The final time of the state evolution
target_time = np.pi

# Number of trotter steps
trotter_steps = 4  ### CAN BE >= 4

# Initialize quantum circuit for 3 qubits
qr = QuantumRegister(7)
qc = QuantumCircuit(qr)

# Prepare initial state (remember we are only evolving 3 of the 7 qubits on jakarta qubits (q_5, q_3, q_1) corresponding to the state |110>)
qc.x([3,5])  # DO NOT MODIFY (|q_5,q_3,q_1> = |110>)

# Simulate time evolution under H_heis3 Hamiltonian
for _ in range(trotter_steps):
    qc.append(Trot_gate, [qr[1], qr[3], qr[5]])

# Evaluate simulation at target_time (t=pi) meaning each trotter step evolves pi/trotter_steps in time
qc = qc.bind_parameters({t: target_time/trotter_steps})

# Generate state tomography circuits to evaluate fidelity of simulation
st_qcs = state_tomography_circuits(qc, [qr[1], qr[3], qr[5]])

# Display circuit for confirmation
# st_qcs[-1].decompose().draw()  # view decomposition of trotter gates
st_qcs[-1].draw()  # only view trotter gates

Finally, the above code displays a drawing of our circuit.

Our trotterization gate appears four times in the figure. This is what we define in line 5 — the trotter_steps. The comment says that it can be greater than four. Let's keep this in mind because it is one parameter we can work with.

We’re almost ready to execute the circuit. We only need to make up for not setting up the backend earlier. So, we use the Aer package that lets us choose a local simulation backend. Specifically, we use the qasm_simulator—a noiseless simulation backend.

shots = 8192
reps = 8

# WE USE A NOISELESS SIMULATION HERE
backend = Aer.get_backend('qasm_simulator')

jobs = []
for _ in range(reps):
    # execute
    job = execute(st_qcs, backend, shots=shots)
    print('Job ID', job.job_id())
    jobs.append(job)

When you run this cell, you’ll get some output like this.

Job ID 58e33e7a-b539-4275-bd7c-f349f74635e1
Job ID 6fd2d9bc-a945-452b-8f14-233396574582
Job ID 3a09070b-b4e6-4441-b433-32ad7d1362ee
Job ID 6c548328-42e8-4821-8133-256ef596adb1
Job ID 74bad2f3-d821-4d55-bf12-efd884e3ec43
Job ID 577b6965-1dbd-4728-ae91-c4c88868ec53
Job ID 8f5bcc1a-df17-4ea0-952b-a936dce8b921
Job ID 260e575d-6be6-4b3a-81ca-195ab78bccd7

We see that we executed the job eight times — the number of replications.

Finally, let’s evaluate the circuit.

# Compute the state tomography based on the st_qcs quantum circuits and the results from those ciricuits
def state_tomo(result, st_qcs):
    # The expected final state; necessary to determine state tomography fidelity
    target_state = (One^One^Zero).to_matrix()  # DO NOT MODIFY (|q_5,q_3,q_1> = |110>)
    # Fit state tomography results
    tomo_fitter = StateTomographyFitter(result, st_qcs)
    rho_fit = tomo_fitter.fit(method='lstsq')
    # Compute fidelity
    fid = state_fidelity(rho_fit, target_state)
    return fid

# Compute tomography fidelities for each repetition
fids = []
for job in jobs:
    fid = state_tomo(job.result(), st_qcs)
    fids.append(fid)
    
print('state tomography fidelity = {:.4f} \u00B1 {:.4f}'.format(np.mean(fids), np.std(fids)))

state tomography fidelity = 0.0003 ± 0.0002

We see a devastating state tomography fidelity of almost 0.

So, let’s see what we can do easily. Since we will run and execute the code a few times, let’s write a wrapper function. The function run_and_evaluate_with_steps takes the number of trotterization steps and the backend to run on as arguments.

def run_and_evaluate_with_steps(steps, backend):
    
    # The final time of the state evolution
    target_time = np.pi

    # Number of trotter steps
    trotter_steps = steps  ### CAN BE >= 4

    # Initialize quantum circuit for 3 qubits
    qr = QuantumRegister(7)
    qc = QuantumCircuit(qr)

    # Prepare initial state (remember we are only evolving 3 of the 7 qubits on jakarta qubits (q_5, q_3, q_1) corresponding to the state |110>)
    qc.x([3,5])  # DO NOT MODIFY (|q_5,q_3,q_1> = |110>)

    # Simulate time evolution under H_heis3 Hamiltonian
    for _ in range(trotter_steps):
        qc.append(Trot_gate, [qr[1], qr[3], qr[5]])

    # Evaluate simulation at target_time (t=pi) meaning each trotter step evolves pi/trotter_steps in time
    qc = qc.bind_parameters({t: target_time/trotter_steps})

    # Generate state tomography circuits to evaluate fidelity of simulation
    st_qcs = state_tomography_circuits(qc, [qr[1], qr[3], qr[5]])

    shots = 8192
    reps = 8

    jobs = []
    for _ in range(reps):
        # execute
        job = execute(st_qcs, backend, shots=shots)
        print('Job ID', job.job_id())
        jobs.append(job)

    
    # Compute tomography fidelities for each repetition
    fids = []
    for job in jobs:
        fid = state_tomo(job.result(), st_qcs)
        fids.append(fid)

    print('state tomography fidelity = {:.4f} \u00B1 {:.4f}'.format(np.mean(fids), np.std(fids)))

We can now run the whole code with different numbers of trotter steps. Let’s try it with eight steps.

run_and_evaluate_with_steps(8, Aer.get_backend('qasm_simulator'))

state tomography fidelity = 0.8543 ± 0.0013

The fidelity is not too bad, is it? So, if increasing the number of trotter steps has such a significant effect, why don’t we try it with even more steps?

run_and_evaluate_with_steps(12, Aer.get_backend('qasm_simulator'))

state tomography fidelity = 0.9687 ± 0.0006

With twelve trotter steps, the fidelity is almost 0.97. That’s practically perfect. So, the exemplary code seems to work pretty well. So, where’s the problem?

IBM asks us to run the circuit on their Jakarta device. Thus far, we simulated a perfect quantum computer. But we know that actual quantum computers are noisy and error-prone. Therefore, let’s simulate an actual device. Qiskit provides the test.mock package that offers simulators corresponding to the behavior of the real devices, such as FakeJakarta.

from qiskit.test.mock import FakeJakarta
device = FakeJakarta()
run_and_evaluate_with_steps(12, device)

state tomography fidelity = 0.1442 ± 0.0021

When we run the code with twelve trotter steps on a noisy simulation of the Jakarta device, the fidelity drops to 0.14.

You may want to play with the number of trotter steps. But, you won’t get any fidelity better than 0.2.

So, the actual problem is not to simulate a Heisenberg model Hamiltonian for a three-particle system using Trotterization. Instead, the problem is to do this on a real 7-qubit device.

When we look at the exemplary trotterization gate, we can see that it uses only three qubits. The other four qubits remain unused. The actual challenge is to use these additional four qubits to immunize the circuit against noise and errors.

You're interested in quantum computing and machine learning...

...But you don't know how to get started?