Quantum State Tomography

Quantum tomography is an experimental procedure to reconstruct a description of part of a quantum system from the measurement outcomes of a specific set of experiments. In particular, quantum state tomography reconstructs the density matrix of a quantum state by preparing the state many times and measuring them in a tomographically complete basis of measurement operators.

Note

This tutorial requires the qiskit-aer and qiskit-ibm-runtime packages to run simulations. You can install them with python -m pip install qiskit-aer qiskit-ibm-runtime.

We first initialize a simulator to run the experiments on.

from qiskit_aer import AerSimulator
from qiskit_ibm_runtime.fake_provider import FakePerth

backend = AerSimulator.from_backend(FakePerth())

To run a state tomography experiment, we initialize the experiment with a circuit to prepare the state to be measured. We can also pass in an Operator or a Statevector to describe the preparation circuit.

import qiskit
from qiskit_experiments.framework import ParallelExperiment
from qiskit_experiments.library import StateTomography

# GHZ State preparation circuit
nq = 2
qc_ghz = qiskit.QuantumCircuit(nq)
qc_ghz.h(0)
qc_ghz.s(0)
for i in range(1, nq):
    qc_ghz.cx(0, i)

# QST Experiment
qstexp1 = StateTomography(qc_ghz)
qstdata1 = qstexp1.run(backend, seed_simulation=100).block_for_results()

# Print results
for result in qstdata1.analysis_results():
    print(result)

AnalysisResult
- name: state
- value: DensityMatrix([[ 4.82945830e-01+0.00000000e+00j,
                -2.97321126e-03+9.86498495e-04j,
                 8.37714364e-05+6.19554167e-04j,
                 1.95163931e-03-4.47176844e-01j],
               [-2.97321126e-03-9.86498495e-04j,
                 1.88274415e-02+0.00000000e+00j,
                -2.28713202e-02-1.68160789e-03j,
                 1.12206435e-02+3.05663829e-03j],
               [ 8.37714364e-05-6.19554167e-04j,
                -2.28713202e-02+1.68160789e-03j,
                 3.15487633e-02-1.73472348e-18j,
                -7.36914787e-03-9.18539786e-03j],
               [ 1.95163931e-03+4.47176844e-01j,
                 1.12206435e-02-3.05663829e-03j,
                -7.36914787e-03+9.18539786e-03j,
                 4.66677966e-01+0.00000000e+00j]],
              dims=(2, 2))
- quality: unknown
- extra: <9 items>
- device_components: ['Q0', 'Q1']
- verified: False
AnalysisResult
- name: state_fidelity
- value: 0.9219887412507268
- quality: unknown
- extra: <9 items>
- device_components: ['Q0', 'Q1']
- verified: False
AnalysisResult
- name: positive
- value: True
- quality: unknown
- extra: <9 items>
- device_components: ['Q0', 'Q1']
- verified: False

Tomography Results

The main result for tomography is the fitted state, which is stored as a DensityMatrix object:

state_result = qstdata1.analysis_results("state")
print(state_result.value)

DensityMatrix([[ 4.82945830e-01+0.00000000e+00j,
                -2.97321126e-03+9.86498495e-04j,
                 8.37714364e-05+6.19554167e-04j,
                 1.95163931e-03-4.47176844e-01j],
               [-2.97321126e-03-9.86498495e-04j,
                 1.88274415e-02+0.00000000e+00j,
                -2.28713202e-02-1.68160789e-03j,
                 1.12206435e-02+3.05663829e-03j],
               [ 8.37714364e-05-6.19554167e-04j,
                -2.28713202e-02+1.68160789e-03j,
                 3.15487633e-02-1.73472348e-18j,
                -7.36914787e-03-9.18539786e-03j],
               [ 1.95163931e-03+4.47176844e-01j,
                 1.12206435e-02-3.05663829e-03j,
                -7.36914787e-03+9.18539786e-03j,
                 4.66677966e-01+0.00000000e+00j]],
              dims=(2, 2))

We can also visualize the density matrix:

from qiskit.visualization import plot_state_city
plot_state_city(qstdata1.analysis_results("state").value, title='Density Matrix')

The state fidelity of the fitted state with the ideal state prepared by the input circuit is stored in the "state_fidelity" result field. Note that if the input circuit contained any measurements the ideal state cannot be automatically generated and this field will be set to None.

fid_result = qstdata1.analysis_results("state_fidelity")
print("State Fidelity = {:.5f}".format(fid_result.value))

State Fidelity = 0.92199

Additional state metadata

Additional data is stored in the tomography under the "state_metadata" field. This includes

• eigvals: the eigenvalues of the fitted state

• trace: the trace of the fitted state

• positive: Whether the eigenvalues are all non-negative

• positive_delta: the deviation from positivity given by 1-norm of negative eigenvalues.

If trace rescaling was performed this dictionary will also contain a raw_trace field containing the trace before rescaling. Futhermore, if the state was rescaled to be positive or trace 1 an additional field raw_eigvals will contain the state eigenvalues before rescaling was performed.

state_result.extra

{'trace': 1.000000000000001,
 'eigvals': array([0.92223018, 0.05310093, 0.02466889, 0.        ]),
 'raw_eigvals': array([ 0.92309311,  0.05396386,  0.02553181, -0.00258878]),
 'rescaled_psd': True,
 'fitter_metadata': {'fitter': 'linear_inversion',
  'fitter_time': 0.004037380218505859},
 'conditional_probability': 1.0,
 'positive': True,
 'experiment': 'StateTomography',
 'run_time': None}

To see the effect of rescaling, we can perform a “bad” fit with very low counts:

# QST Experiment
bad_data = qstexp1.run(backend, shots=10, seed_simulation=100).block_for_results()
bad_state_result = bad_data.analysis_results("state")

# Print result
print(bad_state_result)

# Show extra data
bad_state_result.extra

AnalysisResult
- name: state
- value: DensityMatrix([[ 0.53783545+0.00000000e+00j, -0.04854462+6.17528333e-02j,
                 0.11056488+6.61864453e-03j, -0.05424073-3.64589567e-01j],
               [-0.04854462-6.17528333e-02j,  0.04892982+0.00000000e+00j,
                -0.04210178+9.31473198e-03j,  0.01549988+1.05529334e-02j],
               [ 0.11056488-6.61864453e-03j, -0.04210178-9.31473198e-03j,
                 0.06532003+3.46944695e-18j, -0.07894382-8.08559588e-02j],
               [-0.05424073+3.64589567e-01j,  0.01549988-1.05529334e-02j,
                -0.07894382+8.08559588e-02j,  0.3479147 -1.38777878e-17j]],
              dims=(2, 2))
- quality: unknown
- extra: <9 items>
- device_components: ['Q0', 'Q1']
- verified: False

{'trace': 1.0000000000000007,
 'eigvals': array([0.85745346, 0.14254654, 0.        , 0.        ]),
 'raw_eigvals': array([ 0.95939749,  0.24449057,  0.02949267, -0.23338073]),
 'rescaled_psd': True,
 'fitter_metadata': {'fitter': 'linear_inversion',
  'fitter_time': 0.0029458999633789062},
 'conditional_probability': 1.0,
 'positive': True,
 'experiment': 'StateTomography',
 'run_time': None}

Tomography Fitters

The default fitters is linear_inversion, which reconstructs the state using dual basis of the tomography basis. This will typically result in a non-positive reconstructed state. This state is rescaled to be positive-semidefinite (PSD) by computing its eigen-decomposition and rescaling its eigenvalues using the approach from Ref. [1].

There are several other fitters are included (See API documentation for details). For example, if cvxpy is installed we can use the cvxpy_gaussian_lstsq() fitter, which allows constraining the fit to be PSD without requiring rescaling.

try:
    import cvxpy

    # Set analysis option for cvxpy fitter
    qstexp1.analysis.set_options(fitter='cvxpy_gaussian_lstsq')

    # Re-run experiment
    qstdata2 = qstexp1.run(backend, seed_simulation=100).block_for_results()

    state_result2 = qstdata2.analysis_results("state")
    print(state_result2)
    print("\nextra:")
    for key, val in state_result2.extra.items():
        print(f"- {key}: {val}")

except ModuleNotFoundError:
    print("CVXPY is not installed")

AnalysisResult
- name: state
- value: DensityMatrix([[ 0.4669578 +0.00000000e+00j,  0.01682927-5.66719099e-03j,
                 0.0174351 -3.37508009e-04j,  0.00874897-4.37373576e-01j],
               [ 0.01682927+5.66719099e-03j,  0.0278597 +0.00000000e+00j,
                 0.00880431+2.42351490e-04j, -0.00813464-8.99591120e-03j],
               [ 0.0174351 +3.37508009e-04j,  0.00880431-2.42351490e-04j,
                 0.02394919+0.00000000e+00j,  0.0029724 -6.80490582e-03j],
               [ 0.00874897+4.37373576e-01j, -0.00813464+8.99591120e-03j,
                 0.0029724 +6.80490582e-03j,  0.4812333 +0.00000000e+00j]],
              dims=(2, 2))
- quality: unknown
- extra: <9 items>
- device_components: ['Q0', 'Q1']
- verified: False

extra:
- trace: 1.0000000014421637
- eigvals: [0.9123321  0.04807894 0.02622124 0.01336772]
- raw_eigvals: [0.9123321  0.04807894 0.02622124 0.01336772]
- rescaled_psd: False
- fitter_metadata: {'fitter': 'cvxpy_gaussian_lstsq', 'cvxpy_solver': 'SCS', 'cvxpy_status': ['optimal'], 'psd_constraint': True, 'trace_preserving': True, 'fitter_time': 0.019433259963989258}
- conditional_probability: 1.0
- positive: True
- experiment: StateTomography
- run_time: None

Parallel Tomography Experiment

We can also use the ParallelExperiment class to run subsystem tomography on multiple qubits in parallel.

For example if we want to perform 1-qubit QST on several qubits at once:

from math import pi
num_qubits = 5
gates = [qiskit.circuit.library.RXGate(i * pi / (num_qubits - 1))
         for i in range(num_qubits)]

subexps = [
    StateTomography(gate, physical_qubits=(i,))
    for i, gate in enumerate(gates)
]
parexp = ParallelExperiment(subexps)
pardata = parexp.run(backend, seed_simulation=100).block_for_results()

for result in pardata.analysis_results():
    print(result)

AnalysisResult
- name: state
- value: DensityMatrix([[ 0.96679687+0.j, -0.02832031+0.j],
               [-0.02832031+0.j,  0.03320313+0.j]],
              dims=(2,))
- quality: unknown
- extra: <9 items>
- device_components: ['Q0']
- verified: False
AnalysisResult
- name: state_fidelity
- value: 0.9667968750000001
- quality: unknown
- extra: <9 items>
- device_components: ['Q0']
- verified: False
AnalysisResult
- name: positive
- value: True
- quality: unknown
- extra: <9 items>
- device_components: ['Q0']
- verified: False
AnalysisResult
- name: state
- value: DensityMatrix([[0.84960938+0.j     , 0.00390625+0.34375j],
               [0.00390625-0.34375j, 0.15039063+0.j     ]],
              dims=(2,))
- quality: unknown
- extra: <9 items>
- device_components: ['Q1']
- verified: False
AnalysisResult
- name: state_fidelity
- value: 0.9902791158617665
- quality: unknown
- extra: <9 items>
- device_components: ['Q1']
- verified: False
AnalysisResult
- name: positive
- value: True
- quality: unknown
- extra: <9 items>
- device_components: ['Q1']
- verified: False
AnalysisResult
- name: state
- value: DensityMatrix([[0.515625  +0.j     , 0.01855469+0.46875j],
               [0.01855469-0.46875j, 0.484375  +0.j     ]],
              dims=(2,))
- quality: unknown
- extra: <9 items>
- device_components: ['Q2']
- verified: False
AnalysisResult
- name: state_fidelity
- value: 0.9687500000000002
- quality: unknown
- extra: <9 items>
- device_components: ['Q2']
- verified: False
AnalysisResult
- name: positive
- value: True
- quality: unknown
- extra: <9 items>
- device_components: ['Q2']
- verified: False
AnalysisResult
- name: state
- value: DensityMatrix([[ 0.17773438+0.j        , -0.01464844+0.32617188j],
               [-0.01464844-0.32617188j,  0.82226562+0.j        ]],
              dims=(2,))
- quality: unknown
- extra: <9 items>
- device_components: ['Q3']
- verified: False
AnalysisResult
- name: state_fidelity
- value: 0.9585145534256518
- quality: unknown
- extra: <9 items>
- device_components: ['Q3']
- verified: False
AnalysisResult
- name: positive
- value: True
- quality: unknown
- extra: <9 items>
- device_components: ['Q3']
- verified: False
AnalysisResult
- name: state
- value: DensityMatrix([[ 0.03417969+0.j        , -0.00976562-0.03808594j],
               [-0.00976562+0.03808594j,  0.96582031+0.j        ]],
              dims=(2,))
- quality: unknown
- extra: <9 items>
- device_components: ['Q4']
- verified: False
AnalysisResult
- name: state_fidelity
- value: 0.9658203124999989
- quality: unknown
- extra: <9 items>
- device_components: ['Q4']
- verified: False
AnalysisResult
- name: positive
- value: True
- quality: unknown
- extra: <9 items>
- device_components: ['Q4']
- verified: False

View component experiment analysis results:

for i, expdata in enumerate(pardata.child_data()):
    state_result_i = expdata.analysis_results("state")
    fid_result_i = expdata.analysis_results("state_fidelity")

    print(f'\nPARALLEL EXP {i}')
    print("State Fidelity: {:.5f}".format(fid_result_i.value))
    print("State: {}".format(state_result_i.value))

References

[1]

J Smolin, JM Gambetta, G Smith, Phys. Rev. Lett. 108, 070502 (2012), open access.

See also

• API documentation: StateTomography