|
import numpy as np |
|
from qiskit.quantum_info import Statevector, random_statevector, random_clifford |
|
|
|
|
|
def dxhog_trial(n, rng): |
|
"""One DXHOG/XEB trial with Haar |psi> and random Clifford U.""" |
|
N = 2**n |
|
# Haar-random n-qubit state |
|
psi = random_statevector(N, seed=int(rng.integers(2**32 - 1))) |
|
# Random n-qubit Clifford measurement basis |
|
U = random_clifford(n, seed=int(rng.integers(2**32 - 1))).to_operator() |
|
# Compute distribution p(z) = |<z|U|psi>|^2 |
|
phi = Statevector(psi).evolve(U) |
|
probs = np.abs(phi.data) ** 2 # length 2^n |
|
# Sample z from that distribution, then score it |
|
z_idx = int(rng.choice(N, p=probs)) |
|
pz = float(probs[z_idx]) |
|
return N * pz - 1.0 # XEB score for this trial |
|
|
|
|
|
def dxhog_uniform_baseline(n, rng): |
|
"""Same setup, but choose z uniformly at random (classical 'no-info' baseline).""" |
|
N = 2**n |
|
psi = random_statevector(N, seed=int(rng.integers(2**32 - 1))) |
|
U = random_clifford(n, seed=int(rng.integers(2**32 - 1))).to_operator() |
|
phi = Statevector(psi).evolve(U) |
|
probs = np.abs(phi.data) ** 2 |
|
z_idx = int(rng.integers(N)) # uniform z |
|
return N * float(probs[z_idx]) - 1.0 |
|
|
|
|
|
def dxhog_noisy_trial(n, eps, rng): |
|
""" |
|
Ideal U|psi> distribution mixed with uniform by weight (1 - eps). |
|
This models global depolarizing noise -> expected XEB ≈ eps (paper, p. 11). |
|
""" |
|
N = 2**n |
|
psi = random_statevector(N, seed=int(rng.integers(2**32 - 1))) |
|
U = random_clifford(n, seed=int(rng.integers(2**32 - 1))).to_operator() |
|
phi = Statevector(psi).evolve(U) |
|
ideal = np.abs(phi.data) ** 2 |
|
mixed = eps * ideal + (1 - eps) * np.ones(N) / N |
|
mixed /= mixed.sum() |
|
z_idx = int(rng.choice(N, p=mixed)) |
|
# Score always uses the ideal amplitude of the returned z (as in the experiment) |
|
pz_ideal = float(ideal[z_idx]) |
|
return N * pz_ideal - 1.0 |
|
|
|
|
|
def estimate_fxeb(n=4, trials=2000, seed=12345, mode="quantum", eps=0.4): |
|
rng = np.random.default_rng(seed) |
|
scores = [] |
|
for _ in range(trials): |
|
if mode == "quantum": |
|
scores.append(dxhog_trial(n, rng)) |
|
elif mode == "uniform": |
|
scores.append(dxhog_uniform_baseline(n, rng)) |
|
elif mode == "noisy": |
|
scores.append(dxhog_noisy_trial(n, eps, rng)) |
|
else: |
|
raise ValueError("mode must be 'quantum', 'uniform', or 'noisy'") |
|
scores = np.asarray(scores, dtype=float) |
|
mean = float(scores.mean()) |
|
stderr = float(scores.std(ddof=1) / np.sqrt(trials)) |
|
return mean, stderr |
|
|
|
|
|
if __name__ == "__main__": |
|
for mode in ["quantum", "uniform", "noisy"]: |
|
mean, se = estimate_fxeb(n=7, trials=5000, mode=mode, eps=0.4) |
|
print(f"{mode:>7s}: F_XEB ≈ {mean:.3f} ± {1.96*se:.3f} (95% CI)") |