Quantum Programming 1¶

Visualizing noise processes¶

We will code some examples to visualize the effect of different single-qubit noise processes on the Bloch sphere. We begin from depolarizing noise. Read through the following code, try to understand it and run it on your machine.

We begin from some imports.

In [ ]:
from qiskit import QuantumCircuit
from qiskit.circuit.library import IGate
from qiskit.quantum_info import random_statevector
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit_aer import AerSimulator
from qiskit_aer.noise import NoiseModel, errors, amplitude_damping_error, pauli_error
from qiskit_aer.primitives import EstimatorV2 as Estimator
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

We define a custom identity gate which will introduce the noise in our simulation.

In [ ]:
i_gate_dep = IGate(label="depolarizing")
i_error_dep = errors.depolarizing_error(0.2, 1)

depolarizing_noise_model = NoiseModel()
depolarizing_noise_model.add_all_qubit_quantum_error(i_error_dep, i_gate_dep.label)

qc_dep = QuantumCircuit(1) 

qc_dep.append(i_gate_dep,[0])
qc_dep.draw(output = "mpl")

Any single qubit density matrix can be written $\rho = \frac{1}{2}\mathbb{I} +\frac{1}{2}\left( r_x X + r_y Y + r_z Z\right)$.
Show that $r_j$ is the expectation value of the measurement of the corresponding Pauli matrix given the state $\rho$.

In [ ]:
from qiskit.quantum_info import Pauli
observables = [Pauli("X"),Pauli("Y"),Pauli("Z")]

We want to estimate the mean values of Pauli observables. This can be achieved in qiskit with the Estimator primitives.

In [ ]:
# Instantiate a new estimator
exact_estimator = Estimator()
# The circuit needs to be transpiled to the AerSimulator target 
# (optimization level is set to zero to prevent optimization which might remove our noisy gate)
pass_manager = generate_preset_pass_manager(0, AerSimulator())

According to qiskit's documentation "primitives provide you with the ability to efficiently define vectorized workloads by using a data structure known as a Primitive Unified Bloc (PUB). These PUBs are the fundamental unit of work a QPU needs to execute these workloads. They are used as inputs to the run() method for the Sampler and Estimator primitives, which execute the defined workload as a job".
In particular, we want to specify which circuit we want to run and which observables to estimate.

In [ ]:
# define PUB and run estimation
pub = (qc_dep, observables)
job = exact_estimator.run([pub])
In [ ]:
# Extract the result for the 0th pub (this example only has one pub).
result = job.result()
 
# print the expectation values
print(result[0].data.evs)

We now define another estimator which will simulate the noisy version of the circuit.

In [ ]:
noisy_estimator_dep = Estimator(options=dict(backend_options=dict(noise_model=depolarizing_noise_model)))
job = noisy_estimator_dep.run([pub])
result = job.result()
pub_result = result[0]
print(pub_result.data.evs)

Let's make the above functions that can take arbitrary input states.

In [ ]:
def estimate_exact_evs(input_state):
    
    qc_dep = QuantumCircuit(1) 
    qc_dep.initialize(input_state)
    qc_dep.append(i_gate_dep,[0])
    # display(qc_dep.draw(output = "mpl"))

    
    # Instantiate a new estimator
    exact_estimator = Estimator()
    # The circuit needs to be transpiled to the AerSimulator target (zero to prevent optimization)
    pass_manager = generate_preset_pass_manager(0, AerSimulator())
    # define PUB and run estimation
    pub = (qc_dep, observables)
    job = exact_estimator.run([pub])
    # Extract the result for the 0th pub (this example only has one pub).
    result = job.result()
    
    return result[0].data.evs
In [ ]:
def estimate_noisy_evs(input_state):
    
    qc_dep = QuantumCircuit(1) 
    qc_dep.initialize(input_state)
    qc_dep.append(i_gate_dep,[0])
    # display(qc_dep.draw(output = "mpl"))

    # define PUB and run estimation
    pub = (qc_dep, observables)
    # job = exact_estimator.run([pub])
    noisy_estimator_dep = Estimator(options=dict(backend_options=dict(noise_model=depolarizing_noise_model)))
    job = noisy_estimator_dep.run([pub])
    result = job.result()
    pub_result = result[0]

    return pub_result.data.evs

We can now run them on many random input states.

In [ ]:
# Prepare lists to store Bloch vectors
initial_blochs = []
output_blochs = []
states = [[1,0],[0,1]]
for j in range(100):
# for j in range(2):
    input_state = random_statevector(2)
    # input_state = states[j]
    # print(f"state: {input_state}")
    initial = estimate_exact_evs(input_state)
    output = estimate_noisy_evs(input_state)
    
    # print(initial)
    # print(output)

    initial_blochs.append(initial)
    output_blochs.append(output)

We can visualize the effect of noise on the states by plotting input and output states on the Bloch sphere with the following cell.

In [ ]:
# Convert lists to numpy arrays for easier plotting
initial_blochs = np.array(initial_blochs)
output_blochs = np.array(output_blochs)

def plot_all_bloch_vectors(initial_blochs, output_blochs):
    # Plot all Bloch vectors on the same sphere
    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(111, projection='3d')
    
    # Plot initial states
    ax.scatter(initial_blochs[:, 0], initial_blochs[:, 1], initial_blochs[:, 2], c='b', label='Initial States')
    
    # Plot output states
    ax.scatter(output_blochs[:, 0], output_blochs[:, 1], output_blochs[:, 2], c='r', label='Output States')
    
    # Draw the Bloch sphere
    u = np.linspace(0, 2 * np.pi, 100)
    v = np.linspace(0, np.pi, 100)
    x = np.outer(np.cos(u), np.sin(v))
    y = np.outer(np.sin(u), np.sin(v))
    z = np.outer(np.ones(np.size(u)), np.cos(v))
    ax.plot_surface(x, y, z, color='g', alpha=0.1)
    
    # Set labels and title
    ax.set_xlabel("X")
    ax.set_ylabel("Y")
    ax.set_zlabel("Z")
    ax.set_title("Bloch Vectors Before and After Noise Channel")
    ax.legend()
    
    # Show the plot
    plt.tight_layout()
    plt.show()
    
    # Print average shrinkage factor
    initial_norms = np.linalg.norm(initial_blochs, axis=1)
    output_norms = np.linalg.norm(output_blochs, axis=1)
    avg_shrinkage = np.mean(output_norms / initial_norms)
    print(f"Average shrinkage factor: {avg_shrinkage:.4f}")

plot_all_bloch_vectors(initial_blochs, output_blochs)

To gain more intuition, repeat the above for pure input states chosen from the XY, XZ and YZ planes (display all of the new input-output states in the same plot).

In [ ]:
# YOUR ANSWER HERE...TAKE AS MANY CELLS AS NEEDED

Repeat the steps above to visualize amplitude damping, phase damping, bit flip and phase flip errors.

In [ ]:
In [ ]:
# Prepare lists to store Bloch vectors
initial_blochs = []
output_blochs = []
# states = [[1,0],[0,1]]
for j in range(100):
# for j in range(2):
    input_state = random_statevector(2)
    # input_state = states[j]
    # print(f"state: {input_state}")
    initial = estimate_exact_evs(input_state)
    output = estimate_noisy_evs_amp(input_state) # write this function
    
    # print(initial)
    # print(output)

    initial_blochs.append(initial)
    output_blochs.append(output)

$T_1$ and $T_2$ experiments¶

Suppose a qubit is subject to amplitude damping noise with strength $p$ and dephasing noise with strength $q$.

In [ ]:
p = 0.02
q = 0.01

# Create a noise model with amplitude damping and phase flips
noise_model = NoiseModel()
amp_damp = amplitude_damping_error(p)
phase_damp = pauli_error([('Z', q), ('I', 1-q)])

# Compose the errors and add to the noise model
combined_error = amp_damp.compose(phase_damp)
noise_model.add_all_qubit_quantum_error(combined_error, ['id'])

To estimate $T_1$, we need to run an inversion experiment

In [ ]:
gate_time = 1e-6  # 1 microsecond per delay step

# Function to create inversion experiment
def create_t1_circuit(delay_time):
    qc = QuantumCircuit(1, 1)
    qc.x(0)  # Prepare qubit in |1⟩ state (inversion)
    
    # Implement delay as a sequence of identity gates
    num_delays = int(delay_time / gate_time)
    for _ in range(num_delays):
        qc.id(0)
    
    qc.measure(0, 0)
    return qc

We can run the simulation using qiskit aer

In [ ]:
# Define delay times for the experiment (in microseconds)
delay_times = np.linspace(0, 150e-6, 30)  # 0 to 150 μs

# Run experiments
simulator = AerSimulator(noise_model=noise_model)
excited_state_probs = []

print("Running inversion recovery experiments...")
for delay in delay_times:
    qc = create_t1_circuit(delay)
    job = simulator.run(qc, shots=2048)
    result = job.result()
    counts = result.get_counts()
    
    # Calculate probability of measuring |1⟩
    prob_1 = counts.get('1', 0) / 2048
    excited_state_probs.append(prob_1)
print("Done")

Next we fit to $p(1) = A + B \exp (-t/T_1)$

In [ ]:
# Fit to exponential decay: p(1) = A + B * exp(-t/T1)
def exponential_decay(t, A, T1, B):
    return A + B * np.exp(-t / T1)

# Perform curve fitting to extract T1
popt, pcov = curve_fit(
    exponential_decay, 
    delay_times, 
    excited_state_probs,
    p0=[0.0, 50e-6, 1.0],  # Initial guess
    bounds=([-0.1, 1e-9, 0], [0.1, 1e-3, 1.1])  # Reasonable bounds
)

A_fit, T1_extracted, B_fit = popt
T1_error = np.sqrt(pcov[1, 1])

# Calculate theoretical T1 for comparison
T1_theoretical = gate_time / p

Finally we plot the results

In [ ]:
# Plot results
plt.figure(figsize=(10, 6))
plt.scatter(delay_times * 1e6, excited_state_probs, label='Experimental data', 
            color='blue', alpha=0.6, s=50)
plt.plot(delay_times * 1e6, exponential_decay(delay_times, *popt), 
         'r-', linewidth=2, label=f'Exponential fit: T₁ = {T1_extracted*1e6:.2f} ± {T1_error*1e6:.2f} μs')
plt.axhline(y=0.0, color='gray', linestyle='--', alpha=0.5, label='prob = 0')
plt.xlabel('Delay Time (μs)', fontsize=12)
plt.ylabel('Probability of |1⟩', fontsize=12)
plt.title('T₁ Inversion Experiment', fontsize=14)
plt.legend(fontsize=10)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

print("\n" + "="*60)
print("RESULTS")
print("="*60)
print(f"Extracted T₁ from experiment: {T1_extracted * 1e6:.2f} ± {T1_error * 1e6:.2f} μs")
print(f"Theoretical T₁ (from p={p}): {T1_theoretical * 1e6:.2f} μs")
print(f"Relative error: {abs(T1_extracted - T1_theoretical)/T1_theoretical * 100:.1f}%")
print("\nFit parameters:")
print(f"  A (offset) = {A_fit:.3f}")
print(f"  B (initial population) = {B_fit:.3f}")
print("="*60)

Now write qiskit code to simulate a Ramsey experiment and use the results to estimate the coherence time $T_2$.

In [ ]:
# YOUR ANSWER HERE...TAKE AS MANY CELLS AS NEEDED
In [ ]: