import random import numpy as np import math import matplotlib.pyplot as plt import time def generate_random_measurements(num_rounds, system_size): """ Generate a list of measurement rounds for a quantum system of 'system_size' qubits. Each measurement round is a list of (pauli, outcome) for each qubit. - 'pauli' is a string in {'X', 'Y', 'Z'}. - 'outcome' is an integer in {+1, -1}. """ full_measurement = [] for _ in range(num_rounds): single_round = [] for _ in range(system_size): pauli = random.choice(["X", "Y", "Z"]) outcome = random.choice([+1, -1]) single_round.append((pauli, outcome)) full_measurement.append(single_round) return full_measurement

def single_qubit_estimator(pauli, outcome): """ Build a 2x2 single-qubit estimator matrix. For a given measurement in basis pauli with outcome (+1 or -1), the estimator is defined as (3*projector - I)/2. """ # Define Pauli matrices and identity 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) I = np.eye(2, dtype=complex) if pauli == "Z": # outcome +1 -> |0><0|; outcome -1 -> |1><1| projector = np.array([[1, 0], [0, 0]]) if outcome == +1 else np.array([[0, 0], [0, 1]]) elif pauli == "X": # |+> = (|0>+|1>)/sqrt2, |-> = (|0>-|1>)/sqrt2 plus = 1/np.sqrt(2) * np.array([1, 1], dtype=complex) minus = 1/np.sqrt(2) * np.array([1, -1], dtype=complex) projector = np.outer(plus, plus.conj()) if outcome == +1 else np.outer(minus, minus.conj()) elif pauli == "Y": # |i> = (|0>+i|1>)/sqrt2, |-i> = (|0>-i|1>)/sqrt2 i_state = 1/np.sqrt(2) * np.array([1, 1j], dtype=complex) minus_i_state = 1/np.sqrt(2) * np.array([1, -1j], dtype=complex) projector = np.outer(i_state, i_state.conj()) if outcome == +1 else np.outer(minus_i_state, minus_i_state.conj()) else: raise ValueError("Unknown Pauli basis") return (3 * projector - I) / 2 def build_snapshot(single_measurement): """ Given a single measurement round (a list of (pauli, outcome) for each qubit), construct the full snapshot operator as a tensor product of single-qubit estimators. Note: This returns a 2^n x 2^n matrix (feasible only for small n). """ snapshot = np.array([[1.0 + 0.0j]]) for (pauli, outcome) in single_measurement: snapshot = np.kron(snapshot, single_qubit_estimator(pauli, outcome)) return snapshot def build_all_snapshots(full_measurement): """ Build a list of snapshots from all measurement rounds. Each snapshot is obtained from one measurement round. """ snapshots = [] for round_data in full_measurement: snapshots.append(build_snapshot(round_data)) return snapshots

# Set parameters for demonstration (small system for feasibility) num_rounds = 10 # number of measurement rounds system_size = 5 # number of qubits; increase cautiously as 2^n scaling applies # Generate random measurement data measurement_data = generate_random_measurements(num_rounds, system_size) print("Generated Measurement Data:") for idx, round_data in enumerate(measurement_data, start=1): print(f"Round {idx}: {round_data}") # Build snapshots from measurement data snapshots = build_all_snapshots(measurement_data) # Estimate purity from the snapshots purity_estimate = estimate_purity_from_snapshots(snapshots) print("\nEstimated Purity (tr(rho^2)) =", purity_estimate)

def purity_estimation_error(num_runs, num_rounds, system_size): """ Runs the purity estimation pipeline num_runs times for a given number of measurement rounds (num_rounds) and returns the mean and standard deviation of the purity estimates. """ estimates = [] for _ in range(num_runs): # Generate synthetic measurement data measurement_data = generate_random_measurements(num_rounds, system_size) # Build snapshots from the measurements snapshots = build_all_snapshots(measurement_data) # Estimate purity from the snapshots purity = estimate_purity_from_snapshots(snapshots) estimates.append(purity) return np.mean(estimates), np.std(estimates) # Parameters for the experiment system_size = 3 # For demonstration; full matrices scale as 2^n so keep n small. num_runs = 20 # Number of independent runs to estimate variance # List of different numbers of measurement rounds to test measurement_rounds_list = [10, 20, 50, 100, 200, 500, 1000, 5000] mean_purities = [] std_purities = [] print("Assessing purity estimation error:") for M in measurement_rounds_list: mean_val, std_val = purity_estimation_error(num_runs, M, system_size) mean_purities.append(mean_val) std_purities.append(std_val) print(f"M = {M:3d}: Mean Purity = {mean_val:.4f}, Std Dev = {std_val:.4f}") # Plotting the results plt.errorbar(measurement_rounds_list, mean_purities, yerr=std_purities, fmt='o-', capsize=5) plt.xlabel("Number of Measurement Rounds (M)") plt.ylabel("Estimated Purity (tr(ρ²))") plt.title("Purity Estimation Error vs. Number of Measurements") plt.grid(True, which="both", ls="--") plt.xscale("log") plt.show()

def recommended_num_snapshots(subsystem_size, p2_estimate, epsilon, delta): """ subsystem_size: |AB| p2_estimate: an estimate or upper bound for p2 = tr((rho^TA)^2) epsilon, delta: desired accuracy & confidence (Pr[ |hat{p2} - p2| > epsilon ] <= delta) Returns the recommended M from Lemma 2. """ # Term 1 term1 = (2**subsystem_size) * p2_estimate / (epsilon**2 * delta) # Term 2 term2 = (2**(1.5 * subsystem_size)) / (epsilon * math.sqrt(delta)) return int(math.ceil(8 * max(term1, term2))) M = recommended_num_snapshots(3, 0.125, 0.01, 0.01) print(M)

from qiskit import QuantumCircuit, transpile from qiskit_aer import AerSimulator def apply_identity(qc, qubit): # For Z basis, no gate is needed. pass def apply_hadamard(qc, qubit): qc.h(qubit) def apply_y_measurement(qc, qubit): qc.sdg(qubit) # S† gate qc.h(qubit)

def prepare_random_measurement_circuit(state_preparation, system_size): """ Given a Qiskit QuantumCircuit 'state_preparation' that prepares the unknown state, append a layer of random local Clifford unitaries chosen from {"Z", "X", "Y"}. Here: - "Z": do nothing, - "X": apply H, - "Y": apply S† then H. Then measure all qubits. Returns the full circuit (with measurements) and a list of chosen bases (one per qubit). """ qc = state_preparation.copy() chosen_bases = [] for q in range(system_size): basis = random.choice(["Z", "X", "Y"]) chosen_bases.append(basis) if basis == "Z": apply_identity(qc, q) elif basis == "X": apply_hadamard(qc, q) elif basis == "Y": apply_y_measurement(qc, q) qc.measure_all() return qc, chosen_bases

def single_qubit_estimator(basis, outcome): """ Given a measurement outcome from a single qubit with a chosen basis (from {"Z","X","Y"}) and outcome (+1 or -1), return the corresponding single-qubit estimator. The estimator is defined as (3 * projector - I)/2, where the projector is obtained by "rotating back" the computational-basis projector using the same Clifford unitary. """ I = np.eye(2, dtype=complex) # The Clifford that was applied to map the chosen basis to Z is as follows: if basis == "Z": U = np.array([[1, 0], [0, 1]], dtype=complex) elif basis == "X": U = 1/np.sqrt(2) * np.array([[1, 1], [1, -1]], dtype=complex) elif basis == "Y": # The circuit applies H * Sdg, so the inverse is S * H S = np.array([[1, 0], [0, 1j]], dtype=complex) # S = diag(1, i) Hmat = 1/np.sqrt(2) * np.array([[1, 1], [1, -1]], dtype=complex) U = S.dot(Hmat) else: raise ValueError("Unknown basis") # In the computational basis, assign outcome +1 to |0><0| and -1 to |1><1|. proj = np.array([[1, 0], [0, 0]], dtype=complex) if outcome == +1 else np.array([[0, 0], [0, 1]], dtype=complex) # Rotate back: effective projector = U† * proj * U. projector = U.conj().T.dot(proj).dot(U) return (3 * projector - I) def build_snapshot(single_measurement): """ Given a single measurement round (a list of (basis, outcome) for each qubit), return the full snapshot operator as a tensor product of the single-qubit estimators. """ snapshot = np.array([[1.0+0.0j]]) for basis, outcome in single_measurement: snapshot = np.kron(snapshot, single_qubit_estimator(basis, outcome)) return snapshot def build_all_snapshots(full_measurement): """ Given full measurement data (a list of rounds), build the list of snapshots. """ snapshots = [] for round_data in full_measurement: snapshots.append(build_snapshot(round_data)) return snapshots def estimate_purity_from_snapshots(snapshots): """ Estimate purity using the formula: P_hat = (1 / [M*(M-1)]) * sum_{i != j} trace( snapshots[i] @ snapshots[j] ) where snapshots is a list of M snapshot operators, each a 2^n x 2^n matrix. """ M = len(snapshots) if M < 2: raise ValueError("Need at least 2 snapshots to estimate purity.") total = 0.0 for i in range(M): for j in range(M): if i != j: total += np.trace(snapshots[i].dot(snapshots[j])).real return total / (M*(M-1)) def estimate_purity_biased(snapshots): """ A *biased* purity estimator that sums over ALL pairs (i,j), including i=j, then divides by M^2. """ M = len(snapshots) total = 0.0 for i in range(M): for j in range(M): total += np.trace(snapshots[i].dot(snapshots[j])).real return total / (M*M)

# For demonstration, we prepare a 3-qubit GHZ state. system_size = 3 state_circuit = QuantumCircuit(system_size) state_circuit.h(0) for q in range(system_size-1): state_circuit.cx(q, q+1) # Create an AerSimulator instance. simulator = AerSimulator() def get_random_measurement_round(state_circuit, system_size, simulator): """ For one measurement round, create a new circuit from the state-preparation circuit, append a new random measurement layer (using the three Clifford unitaries), and run it with 1 shot. Returns a measurement round in classical shadows format (list of (basis, outcome) per qubit). """ # Create a new circuit copy from state_circuit. qc, chosen_bases = prepare_random_measurement_circuit(state_circuit, system_size) compiled_qc = transpile(qc, simulator) job = simulator.run(compiled_qc, shots=1) result = job.result() counts = result.get_counts() # Expect one outcome per shot. bitstring = list(counts.keys())[0] # Qiskit returns bitstrings in little-endian; reverse if needed. measurement_round = [] for bit, basis in zip(bitstring[::-1], chosen_bases): outcome = +1 if bit == '0' else -1 measurement_round.append((basis, outcome)) return measurement_round def get_measurement_rounds(state_circuit, system_size, simulator, M): """ Generate M measurement rounds, each with its own random measurement layer. """ rounds = [] for _ in range(M): rounds.append(get_random_measurement_round(state_circuit, system_size, simulator)) return rounds

def purity_estimation_error_qiskit(num_runs, num_rounds, state_circuit, system_size, simulator): unbiased_estimates = [] biased_estimates = [] for _ in range(num_runs): rounds = get_measurement_rounds(state_circuit, system_size, simulator, num_rounds) snapshots = build_all_snapshots(rounds) if len(snapshots) < 2: continue p_unbiased = estimate_purity_from_snapshots(snapshots) p_biased = estimate_purity_biased(snapshots) unbiased_estimates.append(p_unbiased) biased_estimates.append(p_biased) return (np.mean(unbiased_estimates), np.std(unbiased_estimates), np.mean(biased_estimates), np.std(biased_estimates)) # For a pure GHZ state, the true purity is 1. true_purity = 1.0 num_runs = 5 measurement_rounds_list = [10, 20, 50, 100, 200, 500] unbiased_means = [] unbiased_stds = [] biased_means = [] biased_stds = [] print("Running Qiskit-based measurement simulation (GHZ state):") print("True purity for GHZ state: {:.3f}".format(true_purity)) for M in measurement_rounds_list: ub_mean, ub_std, b_mean, b_std = purity_estimation_error_qiskit(num_runs, M, state_circuit, system_size, simulator) unbiased_means.append(ub_mean) unbiased_stds.append(ub_std) biased_means.append(b_mean) biased_stds.append(b_std) print(f"M = {M:3d}: Unbiased = {ub_mean:.4f} ± {ub_std:.4f}, Biased = {b_mean:.4f} ± {b_std:.4f}") # --- Part 5: Plot the Results --- plt.figure(figsize=(8,6)) plt.errorbar(measurement_rounds_list, unbiased_means, yerr=unbiased_stds, fmt='o-', capsize=5, label='Unbiased Estimator') plt.errorbar(measurement_rounds_list, biased_means, yerr=biased_stds, fmt='s--', capsize=5, label='Biased Estimator') plt.axhline(true_purity, color='red', linestyle='--', label=f"True Purity ({true_purity:.3f})") plt.xscale("log") plt.xlabel("Number of Measurement Rounds (M)", fontsize=12) plt.ylabel("Estimated Purity (tr(ρ²))", fontsize=12) plt.title("Purity Estimation (GHZ state): Biased vs. Unbiased", fontsize=14) plt.grid(True, which="both", linestyle="--", alpha=0.6) plt.legend(fontsize=12) plt.tight_layout() plt.show()

import numpy as np import random import matplotlib.pyplot as plt from qiskit import QuantumCircuit, transpile from qiskit_aer import AerSimulator # --- Helper functions for the factorized approach --- def unitary_for_basis(basis): """ Returns the 2x2 unitary that maps the computational (Z) basis to the eigenbasis of the given Pauli. For "Z": Identity. For "X": Hadamard. For "Y": For our protocol, we assume the circuit applies S† then H, so the effective transformation is U = H * S†. """ if basis == "Z": return np.array([[1, 0], [0, 1]], dtype=complex) elif basis == "X": return 1/np.sqrt(2) * np.array([[1, 1], [1, -1]], dtype=complex) elif basis == "Y": Sdg = np.array([[1, 0], [0, -1j]], dtype=complex) H = 1/np.sqrt(2) * np.array([[1, 1], [1, -1]], dtype=complex) return H.dot(Sdg) else: raise ValueError("Unknown basis") def state_from_measurement(basis, outcome): """ Returns the state vector (in the computational basis) corresponding to a measurement outcome in the given basis. For "Z": outcome +1 -> |0>, outcome -1 -> |1> For "X": outcome +1 -> |+> = (|0>+|1>)/√2, outcome -1 -> |-> = (|0>-|1>)/√2 For "Y": outcome +1 -> |i> = (|0>+i|1>)/√2, outcome -1 -> |-i> = (|0>-i|1>)/√2 """ if basis == "Z": return np.array([1,0], dtype=complex) if outcome == +1 else np.array([0,1], dtype=complex) elif basis == "X": return 1/np.sqrt(2)*np.array([1,1], dtype=complex) if outcome==+1 else 1/np.sqrt(2)*np.array([1,-1], dtype=complex) elif basis == "Y": return 1/np.sqrt(2)*np.array([1,1j], dtype=complex) if outcome==+1 else 1/np.sqrt(2)*np.array([1,-1j], dtype=complex) else: raise ValueError("Unknown basis") def T_value_for_qubit(basis_i, outcome_i, basis_j, outcome_j): """ Computes the per-qubit contribution: T_k = 9 * |<ψ_i|ψ_j>|^2 - 4, where |ψ> = U(basis) |o>, with U(basis) from unitary_for_basis. If the same basis is chosen: - outcomes matching: returns 5, - outcomes differing: returns -4. When the basis are different: - outcomes returns: 0.5. """ if basis_i == basis_j: return 5.0 if outcome_i == outcome_j else -4.0 else: psi_i = unitary_for_basis(basis_i).dot(state_from_measurement(basis_i, outcome_i)) psi_j = unitary_for_basis(basis_j).dot(state_from_measurement(basis_j, outcome_j)) A = np.vdot(psi_i, psi_j) return 9*(abs(A)**2) - 4 def trace_product_from_rounds(round_i, round_j): """ Given two measurement rounds (each a list of (basis, outcome) for each qubit), compute the product over qubits: ∏_{k=1}^{N} T_k, which is the same as Tr[hat(ρ)^(i) hat(ρ)^(j)] by the factorization. """ prod = 1.0 for (basis_i, outcome_i), (basis_j, outcome_j) in zip(round_i, round_j): prod *= T_value_for_qubit(basis_i, outcome_i, basis_j, outcome_j) return prod def estimate_purity_unbiased_factorized(measurement_rounds): """ Unbiased estimator for purity (tr(ρ²)) using the factorized approach: P = (1/(M*(M-1))) * sum_{i≠j} ∏_{k=1}^{N} T_k(i,j) """ M = len(measurement_rounds) if M < 2: raise ValueError("Need at least 2 rounds") total = 0.0 for i in range(M): for j in range(M): if i != j: total += trace_product_from_rounds(measurement_rounds[i], measurement_rounds[j]) return total / (M*(M-1)) def estimate_purity_biased_factorized(measurement_rounds): """ Biased estimator that sums over all pairs (including i=j): P = (1/M²) * sum_{i,j} ∏_{k=1}^{N} T_k(i,j) """ M = len(measurement_rounds) if M < 1: raise ValueError("Need at least 1 round") total = 0.0 for i in range(M): for j in range(M): total += trace_product_from_rounds(measurement_rounds[i], measurement_rounds[j]) return total / (M*M) # --- Qiskit-Based Measurement Generation --- def prepare_random_measurement_circuit(state_preparation, system_size): """ Given a state-preparation circuit, append a random measurement layer. For each qubit, randomly choose a basis from {"Z", "X", "Y"}. - For "Z": do nothing, - For "X": apply H, - For "Y": apply Sdg then H. Then measure all qubits. Returns the new circuit and the list of chosen bases. """ qc = state_preparation.copy() chosen_bases = [] for q in range(system_size): basis = random.choice(["Z", "X", "Y"]) chosen_bases.append(basis) if basis == "Z": pass elif basis == "X": qc.h(q) elif basis == "Y": qc.sdg(q) qc.h(q) qc.measure_all() return qc, chosen_bases simulator = AerSimulator() def get_random_measurement_round(state_circuit, system_size, simulator): """ For one shot, create a new circuit (copying the state-prep circuit), append a random measurement layer (new for each shot), run it, and convert the result into the classical shadows round format. """ qc, chosen_bases = prepare_random_measurement_circuit(state_circuit, system_size) compiled_qc = transpile(qc, simulator) job = simulator.run(compiled_qc, shots=1) result = job.result() counts = result.get_counts() # There should be one outcome. bitstring = list(counts.keys())[0] round_data = [] for bit, basis in zip(bitstring[::-1], chosen_bases): outcome = +1 if bit == '0' else -1 round_data.append((basis, outcome)) return round_data def get_measurement_rounds(state_circuit, system_size, simulator, M): rounds = [] for _ in range(M): rounds.append(get_random_measurement_round(state_circuit, system_size, simulator)) return rounds # --- Prepare a 3-Qubit GHZ State --- system_size = 3 state_circuit = QuantumCircuit(system_size) state_circuit.h(0) for q in range(system_size-1): state_circuit.cx(q, q+1) # --- Experiment: Compute Purity Estimates for Different Numbers of Rounds --- def purity_estimation_error_qiskit(num_runs, num_rounds, state_circuit, system_size, simulator): unbiased_estimates = [] biased_estimates = [] for _ in range(num_runs): rounds = get_measurement_rounds(state_circuit, system_size, simulator, num_rounds) p_unbiased = estimate_purity_unbiased_factorized(rounds) p_biased = estimate_purity_biased_factorized(rounds) unbiased_estimates.append(p_unbiased) biased_estimates.append(p_biased) return (np.mean(unbiased_estimates), np.std(unbiased_estimates), np.mean(biased_estimates), np.std(biased_estimates)) true_purity = 1.0 # GHZ is pure. num_runs = 5 measurement_rounds_list = [10, 20, 50, 100, 200, 500] unbiased_means = [] unbiased_stds = [] biased_means = [] biased_stds = [] print("Factorized purity estimation for GHZ state (true purity = 1):") for M in measurement_rounds_list: ub_mean, ub_std, b_mean, b_std = purity_estimation_error_qiskit(num_runs, M, state_circuit, system_size, simulator) unbiased_means.append(ub_mean) unbiased_stds.append(ub_std) biased_means.append(b_mean) biased_stds.append(b_std) print(f"M = {M:3d}: Unbiased = {ub_mean:.4f} ± {ub_std:.4f}, Biased = {b_mean:.4f} ± {b_std:.4f}") plt.figure(figsize=(8,6)) plt.errorbar(measurement_rounds_list, unbiased_means, yerr=unbiased_stds, fmt='o-', capsize=5, label='Unbiased Factorized') plt.errorbar(measurement_rounds_list, biased_means, yerr=biased_stds, fmt='s--', capsize=5, label='Biased Factorized') plt.axhline(true_purity, color='red', linestyle='--', label=f"True Purity ({true_purity})") plt.xscale("log") plt.xlabel("Number of Measurement Rounds (M)", fontsize=12) plt.ylabel("Estimated Purity (tr(ρ²))", fontsize=12) plt.title("Factorized Purity Estimation (GHZ state): Biased vs. Unbiased", fontsize=14) plt.grid(True, which="both", linestyle="--", alpha=0.6) plt.legend(fontsize=12) plt.tight_layout() plt.show()

import numpy as np import matplotlib.pyplot as plt import random from qiskit import QuantumCircuit, transpile from qiskit_aer import AerSimulator # Encode the bases as integers basis_map = {"Z": 0, "X": 1, "Y": 2} def get_measurement_round_array(state_circuit, system_size, simulator): """ Build a circuit from the GHZ state circuit plus a random measurement layer, run it with 1 shot, and convert the result into a single round encoded as [ (basis_code, outcome), ... ]. We'll store it in a shape (N,2) NumPy array for convenience. """ qc = state_circuit.copy() chosen_bases = [] for q in range(system_size): basis = random.choice(["Z", "X", "Y"]) chosen_bases.append(basis) if basis == "Z": pass elif basis == "X": qc.h(q) elif basis == "Y": qc.sdg(q) qc.h(q) qc.measure_all() compiled_qc = transpile(qc, simulator) job = simulator.run(compiled_qc, shots=1) result = job.result() counts = result.get_counts() bitstring = list(counts.keys())[0] # Convert to shape (N,2) array: each row is [basis_code, outcome]. arr = np.zeros((system_size, 2), dtype=np.int32) for i, (bit, basis) in enumerate(zip(bitstring[::-1], chosen_bases)): outcome = +1 if bit == '0' else -1 arr[i, 0] = basis_map[basis] arr[i, 1] = outcome return arr def get_measurement_rounds_array(state_circuit, system_size, simulator, M): """ Generate M measurement rounds, each is shape (N,2). We'll stack them into (M,N,2). """ rounds_array = [] for _ in range(M): arr = get_measurement_round_array(state_circuit, system_size, simulator) rounds_array.append(arr) return np.stack(rounds_array, axis=0) # shape (M, N, 2) def factorized_trace_matrix(rounds_array): """ Vectorized computation of T_k for each pair (i,j) of measurement rounds. Input: rounds_array of shape (M, N, 2). - rounds_array[i, k, 0] = basis code for qubit k in round i - rounds_array[i, k, 1] = outcome for qubit k in round i Output: an (M, M) matrix "prod_matrix" where prod_matrix[i,j] = ∏_{k=1}^N T_k(i,j). We do this by accumulating the product over qubits with NumPy broadcasting. """ M, N, _ = rounds_array.shape prod_matrix = np.ones((M, M), dtype=np.float64) # For each qubit for k in range(N): basis_k = rounds_array[:, k, 0] # shape (M,) outcome_k = rounds_array[:, k, 1] # shape (M,) # same_basis is shape (M, M) same_basis = (basis_k[:, None] == basis_k[None, :]) same_outcome = (outcome_k[:, None] == outcome_k[None, :]) # If same basis & same outcome => T=5, same basis & diff outcome => T=-4 T_same = np.where(same_outcome, 5.0, -4.0) # If basis differ => T=0.5 T_k = np.where(same_basis, T_same, 0.5) prod_matrix *= T_k return prod_matrix def estimate_purity_unbiased_factorized_vectorized(rounds_array): """ Unbiased estimator: P = (1 / [M*(M-1)]) sum_{i != j} ∏_{k=1}^N T_k(i,j). """ prod_matrix = factorized_trace_matrix(rounds_array) M = rounds_array.shape[0] total = prod_matrix.sum() - np.trace(prod_matrix) return total / (M*(M-1)) def estimate_purity_biased_factorized_vectorized(rounds_array): """ Biased estimator: P = (1 / M^2) sum_{i,j} ∏_{k=1}^N T_k(i,j). """ prod_matrix = factorized_trace_matrix(rounds_array) M = rounds_array.shape[0] return prod_matrix.sum() / (M*M) # Now define an experiment with a GHZ state system_size = 3 ghz_circuit = QuantumCircuit(system_size) ghz_circuit.h(0) for q in range(system_size-1): ghz_circuit.cx(q, q+1) simulator = AerSimulator() def run_experiment(num_runs, M, ghz_circuit, system_size, simulator): """ For each run, generate M measurement rounds, convert to an array, compute unbiased & biased purity. Return the average & std dev across runs. """ unbiased_vals = [] biased_vals = [] for _ in range(num_runs): rounds_arr = get_measurement_rounds_array(ghz_circuit, system_size, simulator, M) ub = estimate_purity_unbiased_factorized_vectorized(rounds_arr) b = estimate_purity_biased_factorized_vectorized(rounds_arr) unbiased_vals.append(ub) biased_vals.append(b) return (np.mean(unbiased_vals), np.std(unbiased_vals), np.mean(biased_vals), np.std(biased_vals)) # Let's run it for various M import math true_purity = 1.0 num_runs = 5 measurement_rounds_list = [10, 20, 50, 100, 200, 500] unbiased_means = [] unbiased_stds = [] biased_means = [] biased_stds = [] for M in measurement_rounds_list: ub_mean, ub_std, b_mean, b_std = run_experiment(num_runs, M, ghz_circuit, system_size, simulator) unbiased_means.append(ub_mean) unbiased_stds.append(ub_std) biased_means.append(b_mean) biased_stds.append(b_std) print(f"M={M:4d}: Unbiased = {ub_mean:.4f} ± {ub_std:.4f}, Biased = {b_mean:.4f} ± {b_std:.4f}") plt.figure(figsize=(8,6)) plt.errorbar(measurement_rounds_list, unbiased_means, yerr=unbiased_stds, fmt='o-', capsize=5, label='Unbiased Factorized') plt.errorbar(measurement_rounds_list, biased_means, yerr=biased_stds, fmt='s--', capsize=5, label='Biased Factorized') plt.axhline(true_purity, color='red', linestyle='--', label='True Purity = 1') plt.xscale("log") plt.xlabel("Number of Measurement Rounds (M)", fontsize=12) plt.ylabel("Estimated Purity (tr(ρ²))", fontsize=12) plt.title("Vectorized Factorized Purity Estimation for GHZ State", fontsize=14) plt.grid(True, which="both", linestyle="--", alpha=0.6) plt.legend(fontsize=12) plt.tight_layout() plt.show()