Note

This page was generated from tutorials/noise/8_tomography.ipynb.

Quantum Tomography

Introduction

Quantum tomography is an experimental procedure to reconstruct a description of part of quantum system from the measurement outcomes of a specific set of experiments. In Qiskit we implement the following types of tomography:

1. Quantum state tomography: Given a state-preparation circuit that prepares a system in a state, reconstruct a description of the density matrix \(\rho\) of the actual state obtained in the system.

2. Quantum process tomography: Given a circuit, reconstruct a description of the quantum channel \(\mathcal{E}\) that describes the circuit’s operator when running on the system.

3. Quantum gate set tomography: Performs process tomography on a set of gates in a self-consistent manner, meaning quantum noises on gates used by the tomography process itself is also taken into account.

This notebook gives examples for how to use the ignis.verification.tomography modules.

# Needed for functions
import numpy as np
import time
from copy import deepcopy

# Import Qiskit classes
import qiskit
import qiskit.quantum_info as qi
from qiskit import QuantumRegister, QuantumCircuit, ClassicalRegister, Aer
from qiskit.providers.aer import noise
from qiskit.compiler import assemble

# Tomography functions
from qiskit.ignis.verification.tomography import state_tomography_circuits, StateTomographyFitter
from qiskit.ignis.verification.tomography import process_tomography_circuits, ProcessTomographyFitter
from qiskit.ignis.verification.tomography import gateset_tomography_circuits, GatesetTomographyFitter
import qiskit.ignis.mitigation.measurement as mc

# Auxiliary methods
from qiskit.quantum_info import Choi, Kraus
from qiskit.extensions import HGate, XGate

Initial examples

2-Qubit state tomography Example

In the below example we want to perform state tomography on a 2Q Bell state between qubits 3 and 5. To make the reference circuit we generate the expected statevector using statevector_simulator between qubits 0 and 1.

# Create the expected statevector
q2 = QuantumRegister(2)
bell = QuantumCircuit(q2)
bell.h(q2[0])
bell.cx(q2[0], q2[1])
print(bell)

target_state_bell = qi.Statevector.from_instruction(bell)
print(target_state_bell)

      ┌───┐
q0_0: ┤ H ├──■──
      └───┘┌─┴─┐
q0_1: ─────┤ X ├
           └───┘
Statevector([0.70710678+0.j, 0.        +0.j, 0.        +0.j,
             0.70710678+0.j],
            dims=(2, 2))

# Create the actual circuit
q2 = QuantumRegister(6)
bell = QuantumCircuit(q2)
bell.h(q2[3])
bell.cx(q2[3], q2[5])
print(bell)

q1_0: ──────────

q1_1: ──────────

q1_2: ──────────
      ┌───┐
q1_3: ┤ H ├──■──
      └───┘  │
q1_4: ───────┼──
           ┌─┴─┐
q1_5: ─────┤ X ├
           └───┘

Here we are going to generate and run the state tomography circuits. By only passing in the 2 registers we want to measure the state tomography will only run on that reduced \(2^2\) Hilbert space. However, if we pass the whole register in the state tomography module will try and fit the full \(2^6\) space.

# Generate circuits and run on simulator
t = time.time()

# Generate the state tomography circuits.
qst_bell = state_tomography_circuits(bell, [q2[3], q2[5]])

# Execute
job = qiskit.execute(qst_bell, Aer.get_backend('qasm_simulator'), shots=5000)
print('Time taken:', time.time() - t)

# Fit result
tomo_fitter_bell = StateTomographyFitter(job.result(), qst_bell)

/home/computertreker/git/qiskit/qiskit/.tox/docs/lib/python3.7/site-packages/qiskit/ignis/verification/tomography/basis/circuits.py:468: DeprecationWarning: The QuantumCircuit.__iadd__() method is being deprecated. Use the compose() (potentially with the inplace=True argument) and tensor() methods which are more flexible w.r.t circuit register compatibility.
  prep += circuit
/home/computertreker/git/qiskit/qiskit/.tox/docs/lib/python3.7/site-packages/qiskit/circuit/quantumcircuit.py:942: DeprecationWarning: The QuantumCircuit.extend() method is being deprecated. Use the compose() (potentially with the inplace=True argument) and tensor() methods which are more flexible w.r.t circuit register compatibility.
  return self.extend(rhs)
/home/computertreker/git/qiskit/qiskit/.tox/docs/lib/python3.7/site-packages/qiskit/ignis/verification/tomography/basis/circuits.py:478: DeprecationWarning: The QuantumCircuit.__add__() method is being deprecated.Use the compose() method which is more flexible w.r.t circuit register compatibility.
  circ = prep + meas
/home/computertreker/git/qiskit/qiskit/.tox/docs/lib/python3.7/site-packages/qiskit/circuit/quantumcircuit.py:933: DeprecationWarning: The QuantumCircuit.combine() method is being deprecated. Use the compose() method which is more flexible w.r.t circuit register compatibility.
  return self.combine(rhs)

Time taken: 0.4530761241912842

The fitter will output a density matrix ordered according to how we passed in the registers to state_tomography_circuits.

# Perform the tomography fit
# which outputs a density matrix
rho_fit_bell = tomo_fitter_bell.fit(method='lstsq')
F_bell = qi.state_fidelity(rho_fit_bell, target_state_bell)
print('State Fidelity: F = {:.5f}'.format(F_bell))

State Fidelity: F = 0.99382

Repeat the Example with Measurement Noise

#Add measurement noise
noise_model = noise.NoiseModel()
for qubit in range(6):
    read_err = noise.errors.readout_error.ReadoutError([[0.75, 0.25],[0.1,0.9]])
    noise_model.add_readout_error(read_err,[qubit])

#generate the calibration circuits
meas_calibs, state_labels = mc.complete_meas_cal(qubit_list=[3,5])

backend = Aer.get_backend('qasm_simulator')
job_cal = qiskit.execute(meas_calibs, backend=backend, shots=15000, noise_model=noise_model)
job_tomo = qiskit.execute(qst_bell, backend=backend, shots=15000, noise_model=noise_model)

meas_fitter = mc.CompleteMeasFitter(job_cal.result(),state_labels)

tomo_bell = StateTomographyFitter(job_tomo.result(), qst_bell)

#no correction
rho_bell = tomo_bell.fit(method='lstsq')
F_bell = qi.state_fidelity(rho_bell, target_state_bell)
print('State fidelity (no correction): F = {:.5f}'.format(F_bell))

#correct data
correct_tomo_results = meas_fitter.filter.apply(job_tomo.result(), method='least_squares')
tomo_bell_mit = StateTomographyFitter(correct_tomo_results, qst_bell)
rho_fit_bell_mit = tomo_bell_mit.fit(method='lstsq')
F_bell_mit = qi.state_fidelity(rho_fit_bell_mit, target_state_bell)
print('State fidelity (w/ correction): F = {:.5f}'.format(F_bell_mit))

State fidelity (no correction): F = 0.57282
State fidelity (w/ correction): F = 0.98599

1-qubit process tomography example

# Process tomography of a Hadamard gate
q = QuantumRegister(1)
circ = QuantumCircuit(q)
circ.h(q[0])

# Get the ideal unitary operator
target_unitary = qi.Operator(circ)

# Generate process tomography circuits and run on qasm simulator
qpt_circs = process_tomography_circuits(circ, q)
job = qiskit.execute(qpt_circs, Aer.get_backend('qasm_simulator'), shots=4000)

# Extract tomography data so that counts are indexed by measurement configuration
qpt_tomo = ProcessTomographyFitter(job.result(), qpt_circs)
qpt_tomo.data

{(('Zp',), ('X',)): {'0': 4000},
 (('Zp',), ('Y',)): {'1': 2028, '0': 1972},
 (('Zp',), ('Z',)): {'0': 1994, '1': 2006},
 (('Zm',), ('X',)): {'1': 4000},
 (('Zm',), ('Y',)): {'1': 2082, '0': 1918},
 (('Zm',), ('Z',)): {'0': 1981, '1': 2019},
 (('Xp',), ('X',)): {'1': 2041, '0': 1959},
 (('Xp',), ('Y',)): {'1': 2041, '0': 1959},
 (('Xp',), ('Z',)): {'0': 4000},
 (('Yp',), ('X',)): {'1': 1976, '0': 2024},
 (('Yp',), ('Y',)): {'1': 4000},
 (('Yp',), ('Z',)): {'1': 1989, '0': 2011}}

# Tomographic reconstruction
t = time.time()
choi_fit_lstsq = qpt_tomo.fit(method='lstsq')

print('Fit time:', time.time() - t)
print('Average gate fidelity: F = {:.5f}'.format(qi.average_gate_fidelity(choi_fit_lstsq, target=target_unitary)))

Fit time: 0.004435539245605469
Average gate fidelity: F = 0.99310

1-qubit process tomography of two-qubit swap gate

We will prepare qubit-0 and measure qubit-1 so the reconstructed channel should be an identity.

# Process tomography of a Hadamard gate
q = QuantumRegister(2)
circ = QuantumCircuit(q)
circ.swap(q[0], q[1])

# Generate process tomography circuits and run on qasm simulator
# We use the optional prepared_qubits kwarg to specify that the prepared qubit was different to measured qubit
qpt_circs = process_tomography_circuits(circ, q[1], prepared_qubits=q[0])
job = qiskit.execute(qpt_circs, Aer.get_backend('qasm_simulator'), shots=2000)

# Extract tomography data so that counts are indexed by measurement configuration
qpt_tomo = ProcessTomographyFitter(job.result(), qpt_circs)
qpt_tomo.data

{(('Zp',), ('X',)): {'1': 958, '0': 1042},
 (('Zp',), ('Y',)): {'0': 993, '1': 1007},
 (('Zp',), ('Z',)): {'0': 2000},
 (('Zm',), ('X',)): {'0': 1034, '1': 966},
 (('Zm',), ('Y',)): {'0': 1020, '1': 980},
 (('Zm',), ('Z',)): {'1': 2000},
 (('Xp',), ('X',)): {'0': 2000},
 (('Xp',), ('Y',)): {'0': 998, '1': 1002},
 (('Xp',), ('Z',)): {'0': 987, '1': 1013},
 (('Yp',), ('X',)): {'1': 1030, '0': 970},
 (('Yp',), ('Y',)): {'0': 2000},
 (('Yp',), ('Z',)): {'1': 1016, '0': 984}}

# Tomographic reconstruction

t = time.time()
choi_fit = qpt_tomo.fit(method='lstsq')
print('Fit time:', time.time() - t)
print('Average gate fidelity: F = {:.5f}'.format(qi.average_gate_fidelity(choi_fit)))

Fit time: 0.004413127899169922
Average gate fidelity: F = 0.98710

Advances examples

Generating and fitting random states

We now test the functions on the state generated by a circuit consisting of a layer of random single qubit unitaries u3.

def random_u_tomo(nq, shots):

    def rand_angles():
        return tuple(2 * np.pi * np.random.random(3) - np.pi)
    q = QuantumRegister(nq)
    circ = QuantumCircuit(q)
    for j in range(nq):
        circ.u(*rand_angles(), q[j])
    target_state = qi.Statevector.from_instruction(circ)

    qst_circs = state_tomography_circuits(circ, q)
    job = qiskit.execute(qst_circs, Aer.get_backend('qasm_simulator'),
                         shots=shots)
    tomo_data = StateTomographyFitter(job.result(), qst_circs)
    rho_fit = tomo_data.fit(method='lstsq')

    print('F = {:.5f}'.format(qi.state_fidelity(rho_fit, target_state)))

for j in range(5):
    print('Random single-qubit unitaries: set {}'.format(j))
    random_u_tomo(3, 5000)

Random single-qubit unitaries: set 0
F = 0.99656
Random single-qubit unitaries: set 1
F = 0.99467
Random single-qubit unitaries: set 2
F = 0.99799
Random single-qubit unitaries: set 3
F = 0.99561
Random single-qubit unitaries: set 4
F = 0.99730

5-Qubit Bell State

# Create a state preparation circuit
q5 = QuantumRegister(5)
bell5 = QuantumCircuit(q5)
bell5.h(q5[0])
for j in range(4):
    bell5.cx(q5[j], q5[j + 1])

# Get ideal output state
target_state_bell5 = qi.Statevector.from_instruction(bell5)

# Generate circuits and run on simulator
t = time.time()
qst_bell5 = state_tomography_circuits(bell5, q5)
job = qiskit.execute(qst_bell5, Aer.get_backend('qasm_simulator'), shots=5000)

# Extract tomography data so that counts are indexed by measurement configuration
tomo_bell5 = StateTomographyFitter(job.result(), qst_bell5)
print('Time taken:', time.time() - t)

Time taken: 3.46025013923645

t = time.time()
rho_fit_bell5 = tomo_bell5.fit(method='lstsq')
print('Time taken:', time.time() - t)
print('State fidelity: F = {:.5f}'.format(qi.state_fidelity(rho_fit_bell5, target_state_bell5)))

Time taken: 2.87139630317688
State fidelity: F = 0.99398

2-Qubit Conditional State Tomography

In this example, we have a three-qubit system where one of the qubits will be an ancilla for performing state tomography, i.e. only perform tomography when the third qubit is in the state “1”. The circuit is setup in such a way that after conditional tomography we will get a Bell state on the first two qubits.

First make a 3Q GHZ state with no classical measurements.

# Create the actual circuit
q2 = QuantumRegister(3)
ghz = QuantumCircuit(q2)
ghz.h(q2[0])
ghz.cx(q2[0], q2[1])
ghz.cx(q2[1], q2[2])
ghz.h(q2[2])
print(ghz)

       ┌───┐
q11_0: ┤ H ├──■────────────
       └───┘┌─┴─┐
q11_1: ─────┤ X ├──■───────
            └───┘┌─┴─┐┌───┐
q11_2: ──────────┤ X ├┤ H ├
                 └───┘└───┘

Here we are going to generate and run the state tomography circuits. Only pass the registers we want to perform state tomography on. The code will generate a new classical register for only those measurements.

qst_ghz = state_tomography_circuits(ghz, [q2[0],q2[1]])
print(qst_ghz[0])

       ┌───┐                ░ ┌───┐┌─┐
q11_0: ┤ H ├──■─────────────░─┤ H ├┤M├───
       └───┘┌─┴─┐           ░ ├───┤└╥┘┌─┐
q11_1: ─────┤ X ├──■────────░─┤ H ├─╫─┤M├
            └───┘┌─┴─┐┌───┐ ░ └───┘ ║ └╥┘
q11_2: ──────────┤ X ├┤ H ├─░───────╫──╫─
                 └───┘└───┘ ░       ║  ║
c10: 2/═════════════════════════════╩══╩═
                                    0  1

Now make a copy of this circuit (we will need it for the fitter) and make a new circuit with an ancilla measurement attached (this is what will be run):

#Make a copy without the ancilla register
qst_ghz_no_anc = deepcopy(qst_ghz)
ca = ClassicalRegister(1)
for qst_ghz_circ in qst_ghz:
    qst_ghz_circ.add_register(ca)
    qst_ghz_circ.measure(q2[2],ca[0])

#Run in Aer
job = qiskit.execute(qst_ghz, Aer.get_backend('qasm_simulator'), shots=10000)
raw_results = job.result()

Before sending the results to the state tomography fitter we must strip the register for the Q2 measurement and only keep the results when that register is 1.

new_result = deepcopy(raw_results)

for resultidx, _ in enumerate(raw_results.results):
    old_counts = raw_results.get_counts(resultidx)
    new_counts = {}

    #change the size of the classical register
    new_result.results[resultidx].header.creg_sizes = [new_result.results[resultidx].header.creg_sizes[0]]
    new_result.results[resultidx].header.clbit_labels = new_result.results[resultidx].header.clbit_labels[0:-1]
    new_result.results[resultidx].header.memory_slots = 2

    for reg_key in old_counts:
        reg_bits = reg_key.split(' ')
        if reg_bits[0]=='1':
            new_counts[reg_bits[1]]=old_counts[reg_key]

        new_result.results[resultidx].data.counts = new_counts

tomo_bell = StateTomographyFitter(new_result, qst_ghz_no_anc)
# Perform the tomography fit
# which outputs a density matrix
rho_fit_bell = tomo_bell.fit(method='lstsq')

np.around(rho_fit_bell, 3)

array([[ 0.504+0.j   ,  0.004+0.002j,  0.001-0.002j, -0.498+0.004j],
       [ 0.004-0.002j,  0.001+0.j   , -0.   +0.j   , -0.003+0.001j],
       [ 0.001+0.002j, -0.   -0.j   ,  0.001+0.j   , -0.002-0.003j],
       [-0.498-0.004j, -0.003-0.001j, -0.002+0.003j,  0.495+0.j   ]])

Gate set tomography

1-Qubit gate set tomography Examples

The main difference between gate set tomography and process tomography is that in gate set tomography, the input consists of a gate set basis: A set of gates that are both used in the initialization/measurement phase of the tomography, and are being reconstructed.

Qiskit supplies a default gateset basis; in order to use this gateset basis in order to reconstruct another gate, this gate should be added to the basis. We use the following method to simplify the process:

from qiskit.ignis.verification.tomography.basis import default_gateset_basis

def collect_tomography_data(shots=10000,
                            noise_model=None,
                            gateset_basis='Standard GST'):
    backend_qasm = Aer.get_backend('qasm_simulator')
    circuits = gateset_tomography_circuits(gateset_basis=gateset_basis)
    qobj = assemble(circuits, shots=shots)
    result = backend_qasm.run(qobj, noise_model=noise_model).result()
    fitter = GatesetTomographyFitter(result, circuits, gateset_basis)
    return fitter

def gate_set_tomography(gate, noise_model=None):
    basis = default_gateset_basis()
    basis.add_gate(gate)
    fitter = collect_tomography_data(shots=10000, noise_model=noise_model, gateset_basis=basis)
    result_gates = fitter.fit()
    result_gate = result_gates[gate.name]
    return Choi(result_gate)

Noiseless 1-qubit gate set tomography

target_unitary = qi.Operator(HGate())
t = time.time()
channel_fit = gate_set_tomography(HGate())
print('fit time:', time.time() - t)
print('Average gate fidelity: F = {:.5f}'.format(qi.average_gate_fidelity(channel_fit, target_unitary)))

fit time: 1.4129164218902588
Average gate fidelity: F = 0.99378

import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright

/home/computertreker/git/qiskit/qiskit/.tox/docs/lib/python3.7/site-packages/qiskit/aqua/__init__.py:86: DeprecationWarning: The package qiskit.aqua is deprecated. It was moved/refactored to qiskit-terra For more information see <https://github.com/Qiskit/qiskit-aqua/blob/main/README.md#migration-guide>
  warn_package('aqua', 'qiskit-terra')

Version Information

Qiskit Software	Version
qiskit-terra	0.18.2
qiskit-aer	0.8.2
qiskit-ignis	0.6.0
qiskit-ibmq-provider	0.16.0
qiskit-aqua	0.9.5
qiskit	0.29.1
qiskit-nature	0.2.2
qiskit-finance	0.3.0
qiskit-optimization	0.2.3
qiskit-machine-learning	0.2.1
System information
Python	3.7.12 (default, Nov 22 2021, 14:57:10) [GCC 11.1.0]
OS	Linux
CPUs	32
Memory (Gb)	125.71650314331055
Tue Jan 04 11:14:35 2022 EST

This code is a part of Qiskit

© Copyright IBM 2017, 2022.

This code is licensed under the Apache License, Version 2.0. You may
obtain a copy of this license in the LICENSE.txt file in the root directory
of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.

Any modifications or derivative works of this code must retain this
copyright notice, and modified files need to carry a notice indicating
that they have been altered from the originals.