Quantum Programming 2

import numpy as np
import matplotlib.pyplot as plt

from qiskit import QuantumCircuit
from qiskit.quantum_info import Statevector, DensityMatrix, random_unitary
from qiskit.visualization import plot_state_paulivec

from qiskit_aer import AerSimulator

Quantum state tomography

To identify an unknown density matrix $\rho$, we prepare many copies of the state and measure in the tensor product Pauli basis. This allows us to construct $|\rho \rangle \rangle = (1, \text{Tr}(\rho X), \text{Tr}(\rho Y), \text{Tr}(\rho Z), \ldots)^T$.

Define the measurement bases

MEASUREMENT_BASIS = {
    'Z': QuantumCircuit(1),
    'X': QuantumCircuit(1),
    'Y': QuantumCircuit(1)
}
MEASUREMENT_BASIS['X'].h(0)
MEASUREMENT_BASIS['Y'].sdg(0)
MEASUREMENT_BASIS['Y'].h(0)

Write a function to measure in a particular basis

def measure_in_basis(prep_circuit, basis_label, shots=5000):
    """Append measurement basis rotation and measure."""
    circuit = prep_circuit.copy()
    circuit = circuit.compose(MEASUREMENT_BASIS[basis_label], qubits=[0])
    circuit.measure_all()
    sim = AerSimulator()
    # Transpile for the simulator and run
    # t_circ = transpile(circuit, sim)
    job = sim.run(circuit, shots=shots)
    result = job.result()
    counts = result.get_counts()
    # Probability of |0> and |1>
    p0 = counts.get('0', 0)/shots
    p1 = counts.get('1', 0)/shots
    return p0, p1

Reconstruct the density matrix from the measurement results

def reconstruct_density_matrix_1q(prep_circuit, shots=5000):
    """Reconstruct single-qubit density matrix using linear inversion."""
    # Expectation values of Pauli X,Y,Z
    p0z, p1z = measure_in_basis(prep_circuit, 'Z', shots)
    p0x, p1x = measure_in_basis(prep_circuit, 'X', shots)
    p0y, p1y = measure_in_basis(prep_circuit, 'Y', shots)

    expZ = p0z - p1z
    expX = p0x - p1x
    expY = p0y - p1y

    # rho = (I + xX + yY + zZ)/2
    Id = np.eye(2)
    X = np.array([[0,1],[1,0]],dtype=complex)
    Y = np.array([[0,-1j],[1j,0]],dtype=complex)
    Z = np.array([[1,0],[0,-1]],dtype=complex)

    rho = 0.5*(Id + expX*X + expY*Y + expZ*Z)
    return rho

Test on known input states

# Prepare states
prep_circuits = {}

# |0>
prep_circuits['0'] = QuantumCircuit(1)

# |1>
circuit_1 = QuantumCircuit(1)
circuit_1.x(0)
prep_circuits['1'] = circuit_1

# |+>
circuit_plus = QuantumCircuit(1)
circuit_plus.h(0)
prep_circuits['+'] = circuit_plus

# |+i>
circuit_plus_i = QuantumCircuit(1)
circuit_plus_i.h(0)
circuit_plus_i.s(0)
prep_circuits['+i'] = circuit_plus_i

print("Prepared states:", list(prep_circuits.keys()))

prep_circuit = prep_circuits['0']

rho_fit = reconstruct_density_matrix_1q(prep_circuit, shots=5000)
print("Reconstructed density matrix:\n", rho_fit)

ideal_sv = Statevector.from_instruction(prep_circuit)
ideal_dm = DensityMatrix(ideal_sv)

plot_state_paulivec(rho_fit)

plot_state_paulivec(ideal_dm)

prep_circuit = QuantumCircuit(1)
prep_circuit.unitary(random_unitary(2), [0])

rho_fit = reconstruct_density_matrix_1q(prep_circuit, shots=5000)
print("Reconstructed density matrix:\n", rho_fit)

ideal_sv = Statevector.from_instruction(prep_circuit)
ideal_dm = DensityMatrix(ideal_sv)

plot_state_paulivec(rho_fit)

plot_state_paulivec(ideal_dm)

Repeat the above steps for 2-qubit state tomography.

Hint: use the following expression for a 2-qubit density matrix

$\rho = \frac{1}{4}\sum_{i=0}^3 \sum_{j=0}^3 \text{Tr}[\rho (\sigma_i \otimes \sigma_j)] \sigma_i \otimes \sigma_j$,

where $\sigma_0 = I$, $\sigma_1 = X$, $\sigma_2 = Y$, $\sigma_3 = Z$.

# Your code here :)

Quantum process tomography

To identify an unknown 1-qubit quantum channel $\mathcal E$, we use the Pauli transfer matrix. The procedure is as follows:

1. Do state tomography on the states $R_{\mathcal E} |0\rangle\rangle$, $R_{\mathcal E} |1\rangle\rangle$, $R_{\mathcal E} |+\rangle\rangle$, and $R_{\mathcal E} |+i\rangle\rangle$, giving $|\tilde \rho_0\rangle\rangle$, $|\tilde \rho_1\rangle\rangle$, $|\tilde \rho_+\rangle\rangle$, and $|\tilde \rho_{+i}\rangle\rangle$.
2. Form the matrices $V_{\text{in}}$ and $V_{\text{out}}$ by stacking the vectors column-wise

$V_{\text{in}} = (\quad |0\rangle\rangle \quad |1\rangle\rangle \quad |+\rangle\rangle \quad |+i\rangle\rangle \quad)$

$V_{\text{out}} = (\quad |\tilde \rho_0\rangle\rangle \quad |\tilde \rho_1\rangle\rangle \quad |\tilde \rho_+\rangle\rangle \quad |\tilde \rho_{+i}\rangle\rangle \quad)$

3. To get the PTM we invert the equation $V_{\text{out}} = R_{\mathcal E} V_{\text{in}}$, i.e., $R_{\mathcal E} = V_{\text{out}} V_{\text{in}}^{-1}$.

Implement 1-qubit process tomography.

# Your code here :)

Extension: implement 2-qubit process tomography.