Towards Solving IBM's Quantum Open Science Prize
Today, we meet the entry criteria to participate
How do you find out whether to pursue a career in quantum computing? Almost three months ago, I claimed the best way is to work on IBM's Quantum Open Science Prize.
And this is exactly what we did. In my previous posts, we worked through the nitty-gritty details of that challenge. IBM wants us to simulate a Heisenberg model Hamiltonian for a three-particle system on IBM Quantum's 7-qubit Jakarta system using Trotterization. First, we learned that the real challenge is not simulating a Heisenberg model Hamiltonian for a three-particle system. Instead, it is about how to do it on a noisy device. So, we set up our development environment, looked at how to use a real quantum computer, and introduced ourselves to quantum error mitigation methods.
We chose the Clifford Data Regression (CDR) method. In CDR, we use quantum circuits that we can simulate classically to train a noise model. This noise model lets us predict the noise-free value from a noisy one.
We implemented this method with Qiskit and Mitiq, and the results were promising, both in a local simulation and on a real quantum computer. Unfortunately, we recognized that we couldn't use the CDR method unless we rewrite the way IBM assesses the performance of our solution. We refrained from doing that because it could disqualify our submission. Instead, we decided to change our approach.
Therefore, we looked at how to mitigate the noise to meet IBM's expectations and [mitigated the noise on the measurement level].
So, we're now ready to put it all together.
Conceptually, we follow the approach of the CDR method.
Generate training data
Train a noise model
Predict the noise-free measurements from noisy ones
Generate Training Data
Our training data consists of pairs of measurements. Each pair consists of a noisy and noise-free measurement of the same quantum circuit. To create these training data, we need respective runtime environments.
A noise-free (classical) simulator
A noisy quantum device or a noisy simulator
Today, we use Qiskit's QasmSimulator as both environments. First, we initialize the QasmSimulator Python class without any parameters to get the noise-free simulator. Second, we use our IBM Quantum account to create a simulated backend based on the Jakarta device's noise profile.
| # IBM account | |
| from qiskit.providers.aer import QasmSimulator | |
| # load IBMQ Account data | |
| from qiskit import IBMQ | |
| # Noiseless simulated backend | |
| sim = QasmSimulator() | |
| # replace TOKEN with your API token string (https://quantum-computing.ibm.com/lab/docs/iql/manage/account/ibmq) | |
| IBMQ.save_account("TOKEN", overwrite=True) | |
| account = IBMQ.load_account() | |
| provider = IBMQ.get_provider(hub='ibm-q-community', group='ibmquantumawards', project='open-science-22') | |
| jakarta = provider.get_backend('ibmq_jakarta') | |
| # Simulated backend based on ibmq_jakarta's device noise profile | |
| sim_noisy_jakarta = QasmSimulator.from_backend(provider.get_backend('ibmq_jakarta')) |
In the next step, we need to prepare the quantum training circuits. The challenge here is to use circuits that are representative of the problem but classically simulatable. The original CDR method achieves representability by building upon a linear relationship between the noisy and the noise-free expectation value of a circuit. Further, it uses circuits that consist of Clifford gates. These are gates that are classically simulatable.
Unfortunately, we can't use the expectation values, but we need to mitigate the noise on the measurement level. We can't build upon a linear relationship between the noisy and the noise-free values. Furthermore, the problem with the Heisenberg Hamiltonian simulation is its potential size. So, even though it consists entirely of Clifford gates, the overall circuit becomes intractable for any classical computer given a specific size. So, while a three-particle Heisenberg Hamiltonian is a piece of cake, a fifty-particle Heisenberg Hamiltonian is impossible to simulate classically.
Consequently, we need to rethink our training data. We start with the trotterization algorithm provided by IBM. This algorithm contains at least four trotterization gates. In a previous post, we learned that the more trotterization gates, the better the performance without noise. However, the more trotterization gates we use, the more sensitive to noise the algorithm becomes. Without mitigation, this negative effect of noise exceeds the performance improvement. But with mitigation, it may make sense to use more than the minimum number of four trotterization gates.
For training, however, even the four trotterization gates are too big. We need shorter circuits. Fortunately, each trotter gate consists of three subcircuits, the XX, YY, and ZZ gate. So, the plan is to construct training circuits that contain only two out of the three subcircuits in a single trotterization gate.
The following function provides these circuits for us.
| # Trotterized circuit from IBM material | |
| from qiskit.circuit import Parameter | |
| from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister | |
| def get_circuit(measure, trotter_steps, X=True, Y=True, Z=True): | |
| # 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): | |
| if Z: | |
| Trot_qc.append(ZZ, [Trot_qr[i], Trot_qr[i+1]]) | |
| if Y: | |
| Trot_qc.append(YY, [Trot_qr[i], Trot_qr[i+1]]) | |
| if X: | |
| 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) | |
| # The final time of the state evolution | |
| target_time = np.pi | |
| # Number of trotter steps | |
| #trotter_steps = 8 ### CAN BE >= 4 | |
| # Initialize quantum circuit for 3 qubits | |
| qr = QuantumRegister(7) | |
| cr = ClassicalRegister(3) | |
| qc = QuantumCircuit(qr, cr) if measure is True else 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]]) | |
| if not X or not Y or not Z: | |
| break | |
| # 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}) | |
| if measure: | |
| qc.measure([1,3,5], cr) | |
| return qc |
In general, this is the unchanged trotterization circuit of IBM. The get_circuit function provides a parameterizable interface to create the training circuits from it. The idea is that we can customize which subcircuits we want to include through the parameters X, Y, and Z. Unless we include all three subcircuits, we add only a single trotterization gate. For instance, if we have the XX and the YY gates, the function creates a circuit with these gates for a single trotter step. It is essential to mention that we need to use the intended number of trotterization steps because we use this number as a parameterization for qubit rotations. So, the circuits with a different number of trotter steps are entirely different.
We are now ready to create the training circuits. Let's use 12 trotter steps. You can play with different numbers of steps if you like.
We create three sets of circuits. These are circuits with
XXandYYgatesYYandZZgatesZZandXXgates
Of course, you can also play with the combinations of gates.
Furthermore, we do not create a single circuit each, but 27! The reason is the quantum state tomography IBM uses to assess the performance of our simulation. Quantum state tomography recreates a quantum state by looking at the quantum system from all exclusive angles. For a three-qubit system, these are 3^3=27 angles. Each angle results in a slightly different quantum circuit. But these slight differences significantly affect the measurements. Essentially, we treat each of these circuits separately.
So, in the following step, we create 27 circuits for each specified combination.
| from qiskit.ignis.verification.tomography import state_tomography_circuits | |
| trotters = 12 | |
| train_st_qcs_xy = state_tomography_circuits(get_circuit(False, trotters, True, True, False), [1,3,5]) | |
| train_st_qcs_yz = state_tomography_circuits(get_circuit(False, trotters, False, True, True), [1,3,5]) | |
| train_st_qcs_zx = state_tomography_circuits(get_circuit(False, trotters, True, False, True), [1,3,5]) |
We create two helper functions before we can execute all these circuits (see this previous post). In a nutshell, the sorted_counts function sorts the measurement counts from an executed circuit and ensures that all possible states exist, even if their count is zero.
The get_modifiers function takes a circuit, executes it on our noise-free and the noisy simulator, and calculates a list of modifiers. Our circuit can result in one out of eight states. These range from 000 to 111. There is one modifier for each of these values. So, we get eight modifiers per circuit. A modifier is a number that turns the noisy count into the correct noise-free count when we multiply by it.
| from collections import Counter | |
| from functools import reduce | |
| from qiskit import assemble, execute, transpile | |
| def sorted_counts(counts): | |
| complete = dict(reduce(lambda a, b: a.update(b) or a, [{'000': 0, '001': 0, '010': 0, '011': 0, '100': 0, '101': 0, '110': 0, '111': 0}, counts], Counter())) | |
| return {k: v for k, v in sorted(complete.items(), key=lambda item: item[0])} | |
| def get_modifiers(qc, shots = 4096, display=True): | |
| t_qc_sim = transpile(qc, sim) | |
| noiseless_result = sim.run( assemble(t_qc_sim), shots=shots).result() | |
| noiseless_counts = sorted_counts(noiseless_result.get_counts()) | |
| t_qc = transpile(qc, sim_noisy_jakarta) | |
| qobj = assemble(t_qc) | |
| counts = sorted_counts(sim_noisy_jakarta.run(qobj, shots=shots).result().get_counts()) | |
| zipped = list(zip(noiseless_counts.values(), counts.values())) | |
| modifier = list(map(lambda pair: pair[0]/pair[1] if pair[1] > 0 else 1, zipped)) | |
| if display is True: | |
| print("noisy: ", counts) | |
| print("nose-free: ", noiseless_counts) | |
| print("modifier: ", modifier) | |
| print("\n") | |
| return modifier |
We turn all our circuits into the respective modifiers by applying the get_modifiers function. The variables modifiers_xy, modifiers_yz, and modifiers_zx are lists of the 27 circuits each. Each element is a list of eight modifiers.
| modifiers_xy = list(map(lambda qc: get_modifiers(qc, display=False), train_st_qcs_xy)) | |
| modifiers_yz = list(map(lambda qc: get_modifiers(qc, display=False), train_st_qcs_yz)) | |
| modifiers_zx = list(map(lambda qc: get_modifiers(qc, display=False), train_st_qcs_zx)) | |
| modifiers_zipped = list(zip(modifiers_xy, modifiers_yz, modifiers_zx)) | |
| def mult(tup): | |
| zipped = zip(*tup) | |
| return list(map(lambda x: (x[0]*x[1]*x[2])**(trotters/2), zipped)) | |
| mods = list(map(mult, modifiers_zipped)) |
We zip these into a single list (modifiers_zipped). The result is a list of 27 circuits containing three-item tuples of 8 modifiers.
In the next step, we calculate the final modifiers for each circuit. We do this in the mult function. This function takes the tuple mentioned above of modifiers for a circuit. By zipping these tuples, we get a list of eight items with the corresponding three modifiers each. Effectively, we sorted the modifiers to have all of them in a list that corresponds to a specific count.
We calculate the final modifier as the product of all modifiers in such a list to the half of the trotter steps.
To summarize, we created three small, classically simulatable circuits representing the overall circuit whose measurements we want to mitigate. For each of these three circuits, we calculate a list of modifiers that tell us the factor of how much the noisy count deviates from the noiseless. We combine all modifiers corresponding to a count by multiplying them.
Altogether, each of the small circuits has two-thirds of one trotter step. So, these three small circuits combined denote two trotter steps. To calculate a modifier that corresponds to the overall circuit, we need to do this multiplication multiple times. That is the number of trotter steps divided by two.
Finally, we do not have a single overall circuit but 27 circuits because we aim to use state tomography.
Let's see how well our final modifiers work. We reuse two more helpers we created in previous posts. These are the OwnResult Python class and the state_tomo function that calculates the final state tomography based on that class (see this post).
| from qiskit.result import Result | |
| from qiskit.circuit.quantumcircuit import QuantumCircuit | |
| from qiskit.exceptions import QiskitError | |
| from qiskit.result.counts import Counts | |
| class OwnResult(Result): | |
| def __init__(self, result): | |
| self._result = result | |
| self._counts = {} | |
| def get_counts(self, experiment=None): | |
| if experiment is None: | |
| exp_keys = range(len(self._result.results)) | |
| else: | |
| exp_keys = [experiment] | |
| dict_list = [] | |
| for key in exp_keys: | |
| exp = self._result._get_experiment(key) | |
| try: | |
| header = exp.header.to_dict() | |
| except (AttributeError, QiskitError): # header is not available | |
| header = None | |
| if "counts" in self._result.data(key).keys(): | |
| if header: | |
| counts_header = { | |
| k: v | |
| for k, v in header.items() | |
| if k in {"time_taken", "creg_sizes", "memory_slots"} | |
| } | |
| else: | |
| counts_header = {} | |
| # CUSTOM CODE STARTS HERE ####################### | |
| dict_list.append(Counts( | |
| self._counts[str(key)] if str(key) in map(lambda k: str(k), self._counts.keys()) else self._result.data(key)["counts"] | |
| , **counts_header)) | |
| # CUSTOM CODE ENDS HERE ####################### | |
| elif "statevector" in self._result.data(key).keys(): | |
| vec = postprocess.format_statevector(self._result.data(key)["statevector"]) | |
| dict_list.append(statevector.Statevector(vec).probabilities_dict(decimals=15)) | |
| else: | |
| raise QiskitError(f'No counts for experiment "{repr(key)}"') | |
| # Return first item of dict_list if size is 1 | |
| if len(dict_list) == 1: | |
| return dict_list[0] | |
| else: | |
| return dict_list | |
| def set_counts(self, counts, experiment=None): | |
| self._counts[str(experiment) if experiment is not None else "0"] = counts |
| from qiskit.opflow import Zero, One | |
| from qiskit.ignis.verification.tomography import StateTomographyFitter | |
| from qiskit.quantum_info import state_fidelity | |
| # Compute the state tomography based on the st_qcs quantum circuits and the results from those ciricuits | |
| def state_tomo(result, st_qcs, modifiers, mitigate=False): | |
| # 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>) | |
| own_res = OwnResult(result) | |
| idx = 0 | |
| if mitigate: | |
| for experiment in st_qcs: | |
| exp_keys = [experiment] | |
| for key in exp_keys: | |
| exp = result._get_experiment(key) | |
| counts = sorted_counts(result.get_counts(key)) | |
| mitigated = {item[0]: item[1]*modifiers[idx][i] for i, item in enumerate(counts.items())} | |
| #print("original: ", sorted_counts(result.get_counts(key))) | |
| #print("mitigated: ", sorted_counts(mitigated)) | |
| #print("\n") | |
| own_res.set_counts(mitigated, key) | |
| idx = idx + 1 | |
| # Fit state tomography results | |
| tomo_fitter = StateTomographyFitter(own_res if mitigate else result, st_qcs) | |
| rho_fit = tomo_fitter.fit(method='lstsq') | |
| # Compute fidelity | |
| fid = state_fidelity(rho_fit, target_state) | |
| return fid |
| shots=4096 | |
| st_qcs = state_tomography_circuits(get_circuit(False, trotters, True, True, True), [1,3,5]) | |
| noisy_job = execute(st_qcs, sim_noisy_jakarta, shots=shots) | |
| noisefree_job = execute(st_qcs, sim, shots=shots) | |
| noisy_fid = state_tomo(noisy_job.result(), st_qcs, mods, mitigate=False) | |
| noisefree_fid = state_tomo(noisefree_job.result(), st_qcs, mods, mitigate=False) | |
| mitigated_fid = state_tomo(noisy_job.result(), st_qcs, mods, mitigate=True) | |
| print('noisy state tomography fidelity = {:.4f} \u00B1 {:.4f}'.format(np.mean([noisy_fid]), np.std([noisy_fid]))) | |
| print('noise-free state tomography fidelity = {:.4f} \u00B1 {:.4f}'.format(np.mean([noisefree_fid]), np.std([noisefree_fid]))) | |
| print('mitigated state tomography fidelity = {:.4f} \u00B1 {:.4f}'.format(np.mean([mitigated_fid]), np.std([mitigated_fid]))) |
noisy state tomography fidelity = 0.2861 ± 0.0000
noise-free state tomography fidelity = 0.9699 ± 0.0000
mitigated state tomography fidelity = 0.7518 ± 0.0000The overall results show a significant improvement. With 12 trotter steps, the best possible state tomography fidelity is about 0.97. This is what a noise-free quantum computer would achieve.
However, due to the unmitigated (simulated) noise, the state tomography drops to 0.29. But once we mitigate the noise with our modifiers, we get a state tomography fidelity of 0.75. The entry criteria to participate in IBM's Quantum Open Science Prize is the fidelity of 0.3. So, I'd say we're pretty good.
Looking at the statistics discloses the effectiveness of our approach. We reduce around 68% of the noise.
| unmitigated = np.mean([noisy_fid]) | |
| ideal = np.mean([noisefree_fid]) | |
| mitigated = np.mean([mitigated_fid]) | |
| error_unmitigated = abs(unmitigated-ideal) | |
| error_mitigated = abs(mitigated-ideal) | |
| print("Error (unmitigated):", error_unmitigated) | |
| print("Error (mitigated):", error_mitigated) | |
| print("Relative error (unmitigated):", (error_unmitigated/ideal)) | |
| print("Relative error (mitigatedR):", error_mitigated/ideal) | |
| print(f"Error reduction: {(error_unmitigated-error_mitigated)/error_unmitigated :.1%}.") |
Error (unmitigated): 0.6838963295732885
Error (mitigated): 0.21813650699579556
Relative error (unmitigated): 0.7050850779029165
Relative error (mitigatedR): 0.22489489909174665
Error reduction: 68.1%.You're interested in quantum computing and machine learning...
...But you don't know how to get started?
