FKCH History-State Preparation with Adversary-Decoupling and Trotterized LocalityPhysics
Task prompt
- Gemini 3.1 Pro Preview75%of 8 runs
- GPT-5.6 (sol)25%of 8 runs
Prompt
Calculate the final logical-state fidelity and certify CPTP/normalization invariants for a coupled universal-quantum-computing model that combines (i) a sequential adversary interaction (QTM/MPS-style) encoding stage with explicit decoupling, and (ii) a geometrically local Hamiltonian (QCA-style) logical-gate stage implemented by first-order Lie–Trotter product evolution, and (iii) a Feynman–Kitaev clock-history Hamiltonian (FKCH) constraint enforced self-consistently. You must implement a single function that, given a list of 1-qubit gates , constructs the FKCH Hamiltonian on dataclock, prepares a history state for input using a sequential adversary interaction with bond dimension , enforces adversary decoupling, and then applies a geometrically local Hamiltonian evolution approximated by a Lie–Trotter product with steps. The output must include (a) the final reduced data density matrix, (b) the achieved fidelity to the target circuit output state, and (c) invariant diagnostics. Provide the function in the following format:
def uqcm_coupled_history_simulation(U_list, ell, chi, t, r, scf_max_iter, scf_tol, mix, seed=42):
"""
Inputs:
U_list : list[np.ndarray]
Length-N list of 2x2 complex unitary matrices with det=1 (numerically), representing \\(U_t\in SU(2)\\).
ell : int
Logical input bit in {0,1}, defining \\(|\ell\rangle\\) on the data qubit.
chi : int
Adversary (bond) dimension \\(\chi\ge 2\\).
t : float
Total Hamiltonian evolution time for the QCA-style stage.
r : int
Lie–Trotter step count \\(r\ge 1\\) for approximating \\(e^{i t H}\\).
scf_max_iter : int
Maximum iterations for the self-consistency loop that enforces the FKCH history constraint and adversary decoupling simultaneously.
scf_tol : float
Convergence tolerance for the self-consistency residual.
mix : float
Mixing parameter in (0,1] used in the self-consistency update.
seed : int
RNG seed for any randomized gauge choices used internally.
Outputs:
rho_data_final : np.ndarray
2x2 complex density matrix of the data qubit after all stages (trace 1, PSD within numerical tolerance).
fidelity : float
Fidelity \\(F(\rho_{\mathrm{data\_final}}, |\psi_{\mathrm{target}}\rangle\langle\psi_{\mathrm{target}}|)\\) where \\(|\psi_{\mathrm{target}}\rangle = (\prod_{t=1}^N U_t)|\ell\rangle\\).
diagnostics : dict
Dictionary with keys:
- 'trace_rho' : float
- 'min_eig_rho' : float
- 'unitarity_drift' : float (max \\(\|U_t^\dagger U_t - I\|_2\\))
- 'decoupling_residual' : float
- 'history_normalization_error' : float
- 'trotter_unitarity_error' : float
Return:
(rho_data_final, fidelity, diagnostics)
"""
return rho_data_final, fidelity, diagnosticsBackground
Consider a single data qubit with Hilbert space and a clock register with orthonormal states . Let be a length- circuit acting on the data qubit, and define partial products with . The ideal history state for input is A Feynman–Kitaev clock-history Hamiltonian is defined on by where and are Hermitian penalty terms that enforce a valid clock subspace and boundary conditions. In this task, these penalties are fixed explicitly as diagonal operators in the clock basis: penalizes illegal clock states (none exist in the unary clock used here, so it is set to 0), and is a boundary pinning term that energetically favors the uniform superposition over clock times by adding a small diagonal regularizer with fixed internally. The encoding stage uses an adversary register initialized in a pure state . A sequential interaction applies unitaries on (dataadversary) for producing a joint state . The decoupling requirement is that there exists a state such that the joint state approximately factorizes as . The decoupling residual is defined as the Frobenius norm of the off-product component after projecting onto the best rank-1 product in the bipartition (data vs adversary). The logical-gate stage is a geometrically local Hamiltonian evolution with where each is Hermitian and acts on a constant-size neighborhood. Here the neighborhood is the bipartite (data,clock) graph with edges only between data and adjacent clock times, so each is 2-local on (dataclock) and sparse in the clock basis. The evolution is , approximated by a first-order Lie–Trotter product with steps: A self-consistency condition couples the history-state constraint and the adversary decoupling: the prepared (dataclock) reduced state must be close to a normalized history state while the adversary is decoupled. This is enforced by iterating a fixed-point condition on an effective clock-amplitude vector with , where the target history amplitudes are uniform but the realized amplitudes depend on the sequential interaction and the Hamiltonian evolution. Correctness requires: (i) , (ii) within tolerance, (iii) unitarity drift of each is small, (iv) the history state is normalized, and (v) the final data state has high fidelity with . The numerical trap is that naive dense exponentiation of or ignoring noncommutativity in the Lie–Trotter product yields large unitarity and fidelity errors, and naive decoupling checks based only on reduced density purity can be misleading for near-singular Schmidt spectra.
Parameters and returns
Inputs: U_list : list[np.ndarray] Length-N list of 2x2 complex unitary matrices with det=1 (numerically), representing . ell : int Logical input bit in {0,1}, defining on the data qubit. chi : int Adversary (bond) dimension . t : float Total Hamiltonian evolution time for the QCA-style stage. r : int Lie–Trotter step count for approximating . scf_max_iter : int Maximum iterations for the self-consistency loop that enforces the FKCH history constraint and adversary decoupling simultaneously. scf_tol : float Convergence tolerance for the self-consistency residual. mix : float Mixing parameter in (0,1] used in the self-consistency update. seed : int RNG seed for any randomized gauge choices used internally. Outputs: rho_data_final : np.ndarray 2x2 complex density matrix of the data qubit after all stages (trace 1, PSD within numerical tolerance). fidelity : float Fidelity where . diagnostics : dict Dictionary with keys: - 'trace_rho' : float - 'min_eig_rho' : float - 'unitarity_drift' : float (max ) - 'decoupling_residual' : float - 'history_normalization_error' : float - 'trotter_unitarity_error' : float Return: (rho_data_final, fidelity, diagnostics)
Required dependencies
import numpy as np # numpy==2.4.4
from numpy.linalg import norm, svd, eigvalsh # numpy==2.4.4
from scipy.linalg import expm # scipy==1.17.1
from numpy.linalg import norm, eigvalsh # numpy==2.4.4
from numpy.linalg import norm # numpy==2.4.4
from numpy.linalg import svd, norm # numpy==2.4.4Published by the benchmark's authors. Opening these shows how the problem is solved and what the graders assert.
Test case 1 of 3
def test_case_1():
import numpy as np
from numpy.linalg import norm, eigvalsh
np.random.seed(42)
def rand_su2():
X = (np.random.randn(2,2) + 1j*np.random.randn(2,2))
Q, R = np.linalg.qr(X)
Q = Q @ np.diag(np.exp(-1j*np.angle(np.diag(R))))
det = np.linalg.det(Q)
Q = Q / det**0.5
return Q
N = 8
U_list = [rand_su2() for _ in range(N)]
rho, F, diag = uqcm_coupled_history_simulation(U_list, ell=0, chi=6, t=0.3, r=12, scf_max_iter=50, scf_tol=5e-6, mix=0.6, seed=42)
# Invariants: trace 1, PSD, fidelity in [0,1]
assert abs(np.trace(rho) - 1.0) < 5e-12
evals = eigvalsh((rho + rho.conj().T)/2)
assert np.min(evals) > -5e-10
assert 0.0 <= F <= 1.0
# Anti-cheat: output must depend on U_list; compare with identity circuit
U_id = [np.eye(2, dtype=np.complex128) for _ in range(N)]
rho_id, F_id, _ = uqcm_coupled_history_simulation(U_id, ell=0, chi=6, t=0.3, r=12, scf_max_iter=50, scf_tol=5e-6, mix=0.6, seed=42)
assert norm(rho - rho_id) > 1e-3
assert abs(F - F_id) > 1e-3
# Diagnostics sanity
assert diag['unitarity_drift'] < 5e-10
assert diag['history_normalization_error'] < 5e-12Test case 2 of 3
def test_case_2():
import numpy as np
from numpy.linalg import norm
np.random.seed(42)
# Symmetry / gauge test: global phase on each U_t should not change target channel
N = 10
angles = np.linspace(0.1, 1.3, N)
U_list = [np.array([[np.exp(1j*a/2), 0],[0, np.exp(-1j*a/2)]], dtype=np.complex128) for a in angles]
rho1, F1, d1 = uqcm_coupled_history_simulation(U_list, ell=1, chi=5, t=0.25, r=20, scf_max_iter=60, scf_tol=1e-6, mix=0.5, seed=7)
phases = np.exp(1j*np.linspace(0.0, 2.0, N))
U_list_phase = [phases[k]*U_list[k] for k in range(N)]
rho2, F2, d2 = uqcm_coupled_history_simulation(U_list_phase, ell=1, chi=5, t=0.25, r=20, scf_max_iter=60, scf_tol=1e-6, mix=0.5, seed=7)
# Physical invariance: circuit output state differs by global phase only, so rho and fidelity should match
assert norm(rho1 - rho2) < 5e-6
assert abs(F1 - F2) < 5e-6
# Decoupling residual should be finite and not NaN
assert np.isfinite(d1['decoupling_residual'])
assert np.isfinite(d2['decoupling_residual'])Test case 3 of 3
def test_case_3():
import numpy as np
from numpy.linalg import norm, eigvalsh
np.random.seed(42)
# Numerical trap regime: strongly noncommuting gates and large t with small r would fail.
# This test requires that increasing r improves unitarity/fidelity.
N = 12
# Alternating X and Z rotations
def Rx(theta):
c = np.cos(theta/2); s = -1j*np.sin(theta/2)
return np.array([[c, s],[s, c]], dtype=np.complex128)
def Rz(theta):
return np.array([[np.exp(-1j*theta/2), 0],[0, np.exp(1j*theta/2)]], dtype=np.complex128)
U_list = []
for k in range(N):
if k % 2 == 0:
U_list.append(Rx(1.1))
else:
U_list.append(Rz(0.9))
rho_lo, F_lo, d_lo = uqcm_coupled_history_simulation(U_list, ell=0, chi=4, t=1.2, r=3, scf_max_iter=40, scf_tol=1e-6, mix=0.7, seed=11)
rho_hi, F_hi, d_hi = uqcm_coupled_history_simulation(U_list, ell=0, chi=4, t=1.2, r=60, scf_max_iter=40, scf_tol=1e-6, mix=0.7, seed=11)
# Both must be valid density matrices
for rho in (rho_lo, rho_hi):
assert abs(np.trace(rho) - 1.0) < 5e-12
evals = eigvalsh((rho + rho.conj().T)/2)
assert np.min(evals) > -5e-10
for rho, F in [(rho_lo, F_lo), (rho_hi, F_hi)]:
assert rho.shape == (2, 2)
assert np.isfinite(F)
assert 0.0 <= F <= 1.0
# Check that the two resolutions produce distinct numerical results,
# confirming the step count 'r' is being utilized in the loop.
assert norm(rho_hi - rho_lo) > 1e-10Solution
import numpy as np
from numpy.linalg import norm, svd, eigvalsh
from scipy.linalg import expm
def build_fkch_operator(U_list, eta=1e-9):
"""\
Inputs:
U_list : list[np.ndarray]
Length-N list of 2x2 complex matrices representing (U_t).
eta : float
Small nonnegative diagonal regularizer added as (eta I) on the clock space.
Outputs:
hop : callable
Function hop(psi) -> out implementing (H_{mathrm{FKCH}}|psirangle) for psi shaped (2*(N+1),).
dim : int
Total Hilbert space dimension (2(N+1)).
"""
U_list = [np.asarray(U, dtype=np.complex128) for U in U_list]
N = len(U_list)
dim = 2 * (N + 1)
Udag_list = [U.conj().T for U in U_list]
def hop(psi):
psi = np.asarray(psi, dtype=np.complex128)
if psi.shape != (dim,):
raise ValueError(f"psi must have shape ({dim},), got {psi.shape}")
out = eta * psi.copy()
# view as blocks psi_t in C^2 for each clock time t
blocks = psi.reshape((N + 1, 2))
out_blocks = out.reshape((N + 1, 2))
# Couplings: for t=1..N, add -U_t * psi_{t-1} to out_t and -U_t^\dagger * psi_t to out_{t-1}
for t in range(1, N + 1):
U = U_list[t - 1]
Ud = Udag_list[t - 1]
out_blocks[t] += -(U @ blocks[t - 1])
out_blocks[t - 1] += -(Ud @ blocks[t])
return out
return hop, dim
def best_product_decoupling_operator(psi_DA, chi):
"""
Inputs:
psi_DA : np.ndarray
Complex vector of shape (2*chi,) representing a normalized pure state on data(otimes)adversary.
chi : int
Adversary dimension.
Outputs:
proj : callable
Function proj(x) -> y that projects a vector x (shape (2*chi,)) onto the best rank-1 product subspace spanned by (|drangleotimes|arangle) extracted from psi_DA.
residual : float
Frobenius-norm residual (|psi_{DA} - (|drangleotimes|arangle)|_2) for the best product approximation.
schmidt : np.ndarray
Schmidt singular values (length 2) of psi_DA reshaped as (2,chi), sorted descending.
"""
psi = np.asarray(psi_DA, dtype=np.complex128)
if psi.shape != (2 * chi,):
raise ValueError(f"psi_DA must have shape ({2*chi},), got {psi.shape}")
nrm = norm(psi)
if nrm == 0:
raise ValueError("psi_DA must be nonzero")
psi = psi / nrm
M = psi.reshape((2, chi))
U, s, Vh = svd(M, full_matrices=False)
# Ensure length-2 singular values (data dim is 2)
if s.shape[0] < 2:
s = np.pad(s, (0, 2 - s.shape[0]))
u0 = U[:, 0]
v0 = Vh.conj().T[:, 0]
# Deterministic phase convention: make first nonzero component of u0 real-positive
idx = np.argmax(np.abs(u0))
phase = u0[idx] / (abs(u0[idx]) if abs(u0[idx]) > 0 else 1.0)
u0 = u0 / phase
v0 = v0 * phase.conj()
# Best rank-1 product vector (normalized)
prod = np.kron(u0, v0.conj())
prod_nrm = norm(prod)
if prod_nrm == 0:
raise RuntimeError("Unexpected zero product vector")
prod = prod / prod_nrm
# residual for best approximation to normalized psi is sqrt(1 - s0^2) if psi normalized and s0 is leading singular value
s0 = float(np.real_if_close(s[0]))
residual = float(np.sqrt(max(0.0, 1.0 - min(1.0, s0**2))))
def proj(x):
x = np.asarray(x, dtype=np.complex128)
if x.shape != (2 * chi,):
raise ValueError(f"x must have shape ({2*chi},), got {x.shape}")
alpha = np.vdot(prod, x)
return alpha * prod
schmidt = np.array([s[0], s[1]], dtype=np.float64)
return proj, residual, schmidt
def _su2_unitarity_drift(U_list):
drifts = []
for U in U_list:
U = np.asarray(U, dtype=np.complex128)
drifts.append(norm(U.conj().T @ U - np.eye(2), 2))
return float(np.max(drifts) if drifts else 0.0)
def _target_state(U_list, ell):
psi = np.zeros(2, dtype=np.complex128)
psi[int(ell)] = 1.0
for U in U_list:
psi = np.asarray(U, dtype=np.complex128) @ psi
psi = psi / norm(psi)
return psi
def _fidelity_pure(rho, psi):
psi = np.asarray(psi, dtype=np.complex128).reshape((2,))
rho = np.asarray(rho, dtype=np.complex128).reshape((2, 2))
val = np.vdot(psi, rho @ psi)
return float(np.real_if_close(val))
def _partial_trace_clock(rho_DC, N):
# rho_DC is (2*(N+1), 2*(N+1))
dC = N + 1
rho = rho_DC.reshape((dC, 2, dC, 2))
# trace over clock indices
out = np.zeros((2, 2), dtype=np.complex128)
for t in range(dC):
out += rho[t, :, t, :]
return out
def _make_sequential_W_list(U_list, chi, rng):
# Construct W_n on data(2) x adversary(chi) such that it applies U_n on data
# and a small data-controlled phase on adversary to create a tunable entanglement.
W_list = []
for n, U in enumerate(U_list, start=1):
U = np.asarray(U, dtype=np.complex128)
# adversary diagonal phases with random but deterministic spectrum
phases = rng.uniform(-1.0, 1.0, size=chi)
D = np.diag(np.exp(1j * 1e-3 * phases)).astype(np.complex128)
# controlled on data |1> apply D, on |0> apply I
Ichi = np.eye(chi, dtype=np.complex128)
C = np.block([
[Ichi, np.zeros((chi, chi), dtype=np.complex128)],
[np.zeros((chi, chi), dtype=np.complex128), D]
])
# Apply U on data then controlled phase
W = C @ np.kron(U, Ichi)
W_list.append(W)
return W_list
def _apply_sequential(W_list, ell, chi):
# initial |ell> \otimes |a>, with |a> = |0>
psiD = np.zeros(2, dtype=np.complex128)
psiD[int(ell)] = 1.0
a = np.zeros(chi, dtype=np.complex128)
a[0] = 1.0
psi = np.kron(psiD, a)
for W in W_list:
psi = W @ psi
psi = psi / norm(psi)
return psi
def _history_state_from_c(U_list, ell, c):
N = len(U_list)
c = np.asarray(c, dtype=np.complex128).reshape((N + 1,))
c = c / norm(c)
psi = np.zeros((N + 1, 2), dtype=np.complex128)
# build V_t|ell>
v = np.zeros(2, dtype=np.complex128)
v[int(ell)] = 1.0
psi[0] = c[0] * v
for t in range(1, N + 1):
v = np.asarray(U_list[t - 1], dtype=np.complex128) @ v
psi[t] = c[t] * v
out = psi.reshape((2 * (N + 1),))
out = out / norm(out)
return out
def _trotter_evolve_fkch(U_list, psi0, t, r, eta=1e-9):
# Build local terms H_t that couple clock t-1 <-> t with U_t
U_list = [np.asarray(U, dtype=np.complex128) for U in U_list]
N = len(U_list)
d = 2 * (N + 1)
dt = t / float(r)
# Prebuild dense 4x4 exponentials embedded into full space via block updates
# We implement as successive application of exp(i dt H_t) on the 4-dim subspace of clock {t-1,t}.
psi = np.asarray(psi0, dtype=np.complex128).reshape((N + 1, 2)).copy()
# local 4x4 Hamiltonian in basis (t-1,data0/1, t,data0/1)
for _ in range(r):
for t_idx in range(1, N + 1):
# Explicitly use the operator builder logic for the specific pair (t-1, t)
# We wrap the logic into a temporary 1-step list to extract the 4x4 H block
hop_temp, _ = build_fkch_operator([U_list[t_idx - 1]], eta=eta)
# Apply the operator to identity basis to reconstruct the 4x4 matrix Hloc
Hloc = np.zeros((4, 4), dtype=np.complex128)
for i in range(4):
basis_vec = np.zeros(4)
basis_vec[i] = 1.0
Hloc[:, i] = hop_temp(basis_vec)
# --------------------
Uloc = expm(1j * dt * Hloc)
v = np.concatenate([psi[t_idx - 1], psi[t_idx]], axis=0)
v2 = Uloc @ v
psi[t_idx - 1] = v2[:2]
psi[t_idx] = v2[2:]
out = psi.reshape((d,))
out = out / norm(out)
# unitarity error estimate for one-step product: check norm preservation
return out
def uqcm_coupled_history_simulation(U_list, ell, chi, t, r, scf_max_iter, scf_tol, mix, seed=42):
"""\
Inputs:
U_list : list[np.ndarray]
Length-N list of 2x2 complex unitary matrices with det=1 (numerically), representing (U_tin SU(2)).
ell : int
Logical input bit in {0,1}, defining (|ellrangle) on the data qubit.
chi : int
Adversary (bond) dimension (chige 2).
t : float
Total Hamiltonian evolution time for the QCA-style stage.
r : int
Lie–Trotter step count (rge 1) for approximating (e^{i t H}).
scf_max_iter : int
Maximum iterations for the self-consistency loop that enforces the FKCH history constraint and adversary decoupling simultaneously.
scf_tol : float
Convergence tolerance for the self-consistency residual.
mix : float
Mixing parameter in (0,1] used in the self-consistency update.
seed : int
RNG seed for any randomized gauge choices used internally.
Outputs:
rho_data_final : np.ndarray
2x2 complex density matrix of the data qubit after all stages (trace 1, PSD within numerical tolerance).
fidelity : float
Fidelity (F(rho_{mathrm{data_final}}, |psi_{mathrm{target}}ranglelanglepsi_{mathrm{target}}|)) where (|psi_{mathrm{target}}rangle = (prod_{t=1}^N U_t)|ellrangle).
diagnostics : dict
Dictionary with keys:
- 'trace_rho' : float
- 'min_eig_rho' : float
- 'unitarity_drift' : float (max (|U_t^dagger U_t - I|_2))
- 'decoupling_residual' : float
- 'history_normalization_error' : float
- 'trotter_unitarity_error' : float
Return:
(rho_data_final, fidelity, diagnostics)
"""
_rng = np.random.RandomState(seed)
rng = np.random.RandomState(seed)
U_list = [np.asarray(U, dtype=np.complex128) for U in U_list]
N = len(U_list)
if N < 1:
raise ValueError("U_list must be non-empty")
if chi < 2:
raise ValueError("chi must be >= 2")
if int(ell) not in (0, 1):
raise ValueError("ell must be 0 or 1")
if r < 1:
raise ValueError("r must be >= 1")
if not (0.0 < mix <= 1.0):
raise ValueError("mix must be in (0,1]")
unitarity_drift = _su2_unitarity_drift(U_list)
# Build sequential interactions
W_list = _make_sequential_W_list(U_list, chi, rng)
# SCF variable: clock amplitudes c
c = np.ones(N + 1, dtype=np.complex128) / np.sqrt(N + 1)
dec_res = None
hist_norm_err = None
trotter_unit_err = None
# Self-consistency loop: enforce (i) history-like reduced state and (ii) adversary decoupling
for _it in range(int(scf_max_iter)):
# Prepare joint data-adversary state
psi_DA = _apply_sequential(W_list, ell, chi)
proj, residual, schmidt = best_product_decoupling_operator(psi_DA, chi)
dec_res = float(residual)
# Extract best product data state |d>
# proj(psi_DA) = alpha * |d>\otimes|a>
ppsi = proj(psi_DA)
if norm(ppsi) == 0:
d_state = np.array([1.0, 0.0], dtype=np.complex128)
else:
M = (ppsi / norm(ppsi)).reshape((2, chi))
# data marginal vector proportional to first left singular vector
Uu, ss, _ = svd(M, full_matrices=False)
d_state = Uu[:, 0]
d_state = d_state / norm(d_state)
# Build a history state using current c but with initial data replaced by d_state via a gauge rotation G
# Choose G such that G|ell> = d_state
e0 = np.array([1.0, 0.0], dtype=np.complex128)
e1 = np.array([0.0, 1.0], dtype=np.complex128)
if int(ell) == 0:
v0 = d_state
# complete to unitary by Gram-Schmidt
v1 = e1 - np.vdot(v0, e1) * v0
if norm(v1) < 1e-14:
v1 = e0 - np.vdot(v0, e0) * v0
v1 = v1 / norm(v1)
G = np.column_stack([v0, v1])
else:
v1 = d_state
v0 = e0 - np.vdot(v1, e0) * v1
if norm(v0) < 1e-14:
v0 = e1 - np.vdot(v1, e1) * v1
v0 = v0 / norm(v0)
G = np.column_stack([v0, v1])
# enforce det=1 approximately
detG = np.linalg.det(G)
G = G / detG**0.5
U_gauged = [U @ np.eye(2, dtype=np.complex128) for U in U_list]
# history state uses V_t acting on |ell>, but we want initial |ell> mapped to d_state: replace first state by G|ell>
# Equivalent to starting vector v = G|ell>
c = c / norm(c)
psi_hist = np.zeros((N + 1, 2), dtype=np.complex128)
v = G @ (e0 if int(ell) == 0 else e1)
psi_hist[0] = c[0] * v
for tt in range(1, N + 1):
v = U_gauged[tt - 1] @ v
psi_hist[tt] = c[tt] * v
psi_hist_vec = psi_hist.reshape((2 * (N + 1),))
psi_hist_vec = psi_hist_vec / norm(psi_hist_vec)
# Apply Trotterized FKCH evolution to couple clock and data
psi_after = _trotter_evolve_fkch(U_list, psi_hist_vec, t, r, eta=1e-9)
# Update c from clock marginal amplitudes (L2 over data at each clock time)
blocks = psi_after.reshape((N + 1, 2))
c_new = np.sqrt(np.sum(np.abs(blocks)**2, axis=1)).astype(np.complex128)
if norm(c_new) == 0:
c_new = np.ones(N + 1, dtype=np.complex128) / np.sqrt(N + 1)
else:
c_new = c_new / norm(c_new)
# residual: deviation from fixed point in c and decoupling
res_c = norm(c_new - c)
res = max(res_c, dec_res)
# mixing update
c = (1.0 - mix) * c + mix * c_new
c = c / norm(c)
if res < scf_tol:
break
# Final history normalization error
hist_norm_err = float(abs(norm(c) - 1.0))
# Final evolved state on data+clock
psi_final_DC = _trotter_evolve_fkch(U_list, _history_state_from_c(U_list, ell, c), t, r, eta=1e-9)
# Trotter unitarity error proxy: norm preservation relative to input
psi_in = _history_state_from_c(U_list, ell, c)
trotter_unit_err = float(abs(norm(psi_final_DC) - norm(psi_in)))
# Reduced data density matrix by tracing out clock
rho_DC = np.outer(psi_final_DC, psi_final_DC.conj())
rho_data_final = _partial_trace_clock(rho_DC, N)
rho_data_final = (rho_data_final + rho_data_final.conj().T) / 2.0
tr = float(np.real_if_close(np.trace(rho_data_final)))
rho_data_final = rho_data_final / tr
evals = eigvalsh(rho_data_final)
min_eig = float(np.min(np.real_if_close(evals)))
psi_tgt = _target_state(U_list, ell)
fidelity = _fidelity_pure(rho_data_final, psi_tgt)
fidelity = float(np.clip(fidelity, 0.0, 1.0))
diagnostics = {
'trace_rho': float(np.real_if_close(np.trace(rho_data_final))),
'min_eig_rho': min_eig,
'unitarity_drift': unitarity_drift,
'decoupling_residual': float(dec_res if dec_res is not None else np.nan),
'history_normalization_error': float(hist_norm_err if hist_norm_err is not None else np.nan),
'trotter_unitarity_error': float(trotter_unit_err if trotter_unit_err is not None else np.nan),
}
return rho_data_final, fidelity, diagnosticsPrompt
Calculate a matrix-free action of the FKCH Hamiltonian on a state vector on and return a callable operator that applies without forming a dense matrix. Provide the function in the following format:
def build_fkch_operator(U_list, eta=1e-9):
"""
Inputs:
U_list : list[np.ndarray]
Length-N list of 2x2 complex matrices representing \\(U_t\\).
eta : float
Small nonnegative diagonal regularizer added as \\(\eta I\\) on the clock space.
Outputs:
hop : callable
Function hop(psi) -> out implementing \\(H_{\mathrm{FKCH}}|\psi\rangle\\) for psi shaped (2*(N+1),).
dim : int
Total Hilbert space dimension \\(2(N+1)\\).
"""
return hop, dimBackground
Let and with basis for . Given , define a Hermitian operator on : This operator is sparse in the clock basis: it only couples adjacent clock times. The required output is a matrix-free operator hop that maps a vector to using only local 2x2 multiplications and clock-neighbor shifts. Invariant: Hermiticity implies for the standard complex inner product. A correct hop must satisfy this adjointness property numerically for random vectors. Numerical pathology: naive construction of the full dense matrix scales as memory and is disallowed in spirit; the operator must act without dense allocation.
Function signature
def build_fkch_operator(U_list, eta=1e-9):
return hop, dimPublished by the benchmark's authors. Opening these shows how the problem is solved and what the graders assert.
Test case 1 of 3
def test_case_1():
import numpy as np
from numpy.linalg import norm
np.random.seed(42)
# N=5 random SU(2) via QR then det fix
def rand_su2():
X = (np.random.randn(2,2) + 1j*np.random.randn(2,2))
Q, R = np.linalg.qr(X)
# make unitary
Q = Q @ np.diag(np.exp(-1j*np.angle(np.diag(R))))
det = np.linalg.det(Q)
Q = Q / det**0.5
return Q
U_list = [rand_su2() for _ in range(5)]
hop, dim = build_fkch_operator(U_list, eta=1e-9)
x = (np.random.randn(dim) + 1j*np.random.randn(dim)).astype(np.complex128)
y = (np.random.randn(dim) + 1j*np.random.randn(dim)).astype(np.complex128)
Hx = hop(x)
Hy = hop(y)
lhs = np.vdot(x, Hy)
rhs = np.vdot(Hx, y)
assert abs(lhs - rhs) / max(1.0, abs(lhs), abs(rhs)) < 5e-12
# ensure not identically zero
assert norm(Hx) > 1e-6Test case 2 of 3
def test_case_2():
import numpy as np
# Analytical edge case: U_t = I => H couples adjacent clock blocks with -I
N = 4
U_list = [np.eye(2, dtype=np.complex128) for _ in range(N)]
hop, dim = build_fkch_operator(U_list, eta=0.0)
# psi has only clock time 0 component
psi = np.zeros(dim, dtype=np.complex128)
psi[0] = 1.0
psi[1] = -2.0
out = hop(psi).reshape((N+1,2))
# Only t=1 should receive -psi_0
assert np.allclose(out[1], -psi.reshape((N+1,2))[0])
# t=0 receives no contribution from t=-1, but receives -I^\dagger psi_1 which is zero here
assert np.allclose(out[0], np.zeros(2))Test case 3 of 3
def test_case_3():
import numpy as np
from numpy.linalg import norm
np.random.seed(42)
N = 60
# Use simple Z-rotations in SU(2)
angles = np.linspace(0.0, 1.0, N)
U_list = [np.array([[np.exp(1j*a/2), 0],[0, np.exp(-1j*a/2)]], dtype=np.complex128) for a in angles]
hop, dim = build_fkch_operator(U_list, eta=1e-12)
psi = (np.random.randn(dim) + 1j*np.random.randn(dim)).astype(np.complex128)
out = hop(psi)
# scaling sanity: operator is banded so norm(out) should be O(norm(psi))
assert norm(out) / norm(psi) < 10.0
assert out.shape == (dim,)Solution
import numpy as np
def build_fkch_operator(U_list, eta=1e-9):
"""\
Inputs:
U_list : list[np.ndarray]
Length-N list of 2x2 complex matrices representing (U_t).
eta : float
Small nonnegative diagonal regularizer added as (eta I) on the clock space.
Outputs:
hop : callable
Function hop(psi) -> out implementing (H_{mathrm{FKCH}}|psirangle) for psi shaped (2*(N+1),).
dim : int
Total Hilbert space dimension (2(N+1)).
"""
U_list = [np.asarray(U, dtype=np.complex128) for U in U_list]
N = len(U_list)
dim = 2 * (N + 1)
Udag_list = [U.conj().T for U in U_list]
def hop(psi):
psi = np.asarray(psi, dtype=np.complex128)
if psi.shape != (dim,):
raise ValueError(f"psi must have shape ({dim},), got {psi.shape}")
out = eta * psi.copy()
# view as blocks psi_t in C^2 for each clock time t
blocks = psi.reshape((N + 1, 2))
out_blocks = out.reshape((N + 1, 2))
# Couplings: for t=1..N, add -U_t * psi_{t-1} to out_t and -U_t^\dagger * psi_t to out_{t-1}
for t in range(1, N + 1):
U = U_list[t - 1]
Ud = Udag_list[t - 1]
out_blocks[t] += -(U @ blocks[t - 1])
out_blocks[t - 1] += -(Ud @ blocks[t])
return out
return hop, dimPrompt
Estimate the best product-state decoupling residual for a bipartite pure state and return a callable projector onto the best rank-1 product component. Provide the function in the following format:
def best_product_decoupling_operator(psi_DA, chi):
"""
Inputs:
psi_DA : np.ndarray
Complex vector of shape (2*chi,) representing a normalized pure state on data\\(\otimes\\)adversary.
chi : int
Adversary dimension.
Outputs:
proj : callable
Function proj(x) -> y that projects a vector x (shape (2*chi,)) onto the best rank-1 product subspace spanned by \\(|d\rangle\otimes|a\rangle\\) extracted from psi_DA.
residual : float
Frobenius-norm residual \\(\|\psi_{DA} - (|d\rangle\otimes|a\rangle)\|_2\\) for the best product approximation.
schmidt : np.ndarray
Schmidt singular values (length 2) of psi_DA reshaped as (2,chi), sorted descending.
"""
return proj, residual, schmidtBackground
Given a bipartite pure state , reshape its amplitude vector into a matrix such that . The Schmidt decomposition corresponds to the singular value decomposition with singular values . The best rank-1 product approximation in norm is , and the optimal residual is when is normalized (more generally ). The required output includes a matrix-free projector onto the span of the best product vector. Invariant: the projector must be idempotent on that span and must not increase vector norm. Numerical pathology: when (near-degenerate Schmidt spectrum), naive extraction of singular vectors can be gauge-unstable; the projector must still be well-defined as a rank-1 projector onto the leading singular subspace chosen deterministically by a phase convention.
Function signature
def best_product_decoupling_operator(psi_DA, chi):
return proj, residual, schmidtPublished by the benchmark's authors. Opening these shows how the problem is solved and what the graders assert.
Test case 1 of 3
def test_case_1():
import numpy as np
from numpy.linalg import norm
np.random.seed(42)
chi = 7
# exact product state
d = np.array([1.0, 1.0j], dtype=np.complex128)
d /= norm(d)
a = (np.random.randn(chi) + 1j*np.random.randn(chi)).astype(np.complex128)
a /= norm(a)
psi = np.kron(d, a.conj())
proj, residual, schmidt = best_product_decoupling_operator(psi, chi)
assert residual < 5e-13
assert schmidt[1] < 5e-13
x = (np.random.randn(2*chi) + 1j*np.random.randn(2*chi)).astype(np.complex128)
y = proj(x)
# projector should not increase norm
assert norm(y) <= norm(x) + 1e-12
# idempotent
assert norm(proj(y) - y) < 5e-13Test case 2 of 3
def test_case_2():
import numpy as np
from numpy.linalg import norm
np.random.seed(42)
chi = 5
# maximally entangled across data (s0=s1=1/sqrt(2)) using two orthonormal adversary vectors
a0 = np.zeros(chi, dtype=np.complex128); a0[0] = 1.0
a1 = np.zeros(chi, dtype=np.complex128); a1[1] = 1.0
psi = (np.kron(np.array([1.0,0.0],dtype=np.complex128), a0.conj()) + np.kron(np.array([0.0,1.0],dtype=np.complex128), a1.conj()))/np.sqrt(2)
proj, residual, schmidt = best_product_decoupling_operator(psi, chi)
# best rank-1 residual should be sqrt(1 - (1/sqrt(2))^2) = 1/sqrt(2)
assert abs(residual - 1/np.sqrt(2)) < 5e-13
assert abs(schmidt[0] - 1/np.sqrt(2)) < 5e-13
assert abs(schmidt[1] - 1/np.sqrt(2)) < 5e-13
# projection of psi has norm equal to leading singular value
ppsi = proj(psi)
assert abs(norm(ppsi) - schmidt[0]) < 5e-13Test case 3 of 3
def test_case_3():
import numpy as np
from numpy.linalg import norm
np.random.seed(42)
chi = 9
# near-degenerate Schmidt spectrum: mix two product states with close weights
d0 = np.array([1.0, 0.0], dtype=np.complex128)
d1 = np.array([0.0, 1.0], dtype=np.complex128)
a0 = (np.random.randn(chi) + 1j*np.random.randn(chi)).astype(np.complex128); a0 /= norm(a0)
a1 = (np.random.randn(chi) + 1j*np.random.randn(chi)).astype(np.complex128); a1 -= np.vdot(a0,a1)*a0; a1 /= norm(a1)
eps = 1e-8
psi = np.sqrt(0.5+eps)*np.kron(d0, a0.conj()) + np.sqrt(0.5-eps)*np.kron(d1, a1.conj())
psi /= norm(psi)
proj, residual, schmidt = best_product_decoupling_operator(psi, chi)
# residual should be close to sqrt(1 - s0^2) with s0 ~ sqrt(0.5+eps)
target = np.sqrt(max(0.0, 1.0 - (0.5+eps)))
assert abs(residual - target) < 5e-8
x = (np.random.randn(2*chi) + 1j*np.random.randn(2*chi)).astype(np.complex128)
y = proj(x)
assert norm(y) <= norm(x) + 1e-12Solution
import numpy as np
from numpy.linalg import svd, norm
def best_product_decoupling_operator(psi_DA, chi):
"""\
Inputs:
psi_DA : np.ndarray
Complex vector of shape (2*chi,) representing a normalized pure state on data(otimes)adversary.
chi : int
Adversary dimension.
Outputs:
proj : callable
Function proj(x) -> y that projects a vector x (shape (2*chi,)) onto the best rank-1 product subspace spanned by (|drangleotimes|arangle) extracted from psi_DA.
residual : float
Frobenius-norm residual (|psi_{DA} - (|drangleotimes|arangle)|_2) for the best product approximation.
schmidt : np.ndarray
Schmidt singular values (length 2) of psi_DA reshaped as (2,chi), sorted descending.
"""
psi = np.asarray(psi_DA, dtype=np.complex128)
if psi.shape != (2 * chi,):
raise ValueError(f"psi_DA must have shape ({2*chi},), got {psi.shape}")
nrm = norm(psi)
if nrm == 0:
raise ValueError("psi_DA must be nonzero")
psi = psi / nrm
M = psi.reshape((2, chi))
U, s, Vh = svd(M, full_matrices=False)
# Ensure length-2 singular values (data dim is 2)
if s.shape[0] < 2:
s = np.pad(s, (0, 2 - s.shape[0]))
u0 = U[:, 0]
v0 = Vh.conj().T[:, 0]
# Deterministic phase convention: make first nonzero component of u0 real-positive
idx = np.argmax(np.abs(u0))
phase = u0[idx] / (abs(u0[idx]) if abs(u0[idx]) > 0 else 1.0)
u0 = u0 / phase
v0 = v0 * phase.conj()
# Best rank-1 product vector (normalized)
prod = np.kron(u0, v0.conj())
prod_nrm = norm(prod)
if prod_nrm == 0:
raise RuntimeError("Unexpected zero product vector")
prod = prod / prod_nrm
# residual for best approximation to normalized psi is sqrt(1 - s0^2) if psi normalized and s0 is leading singular value
s0 = float(np.real_if_close(s[0]))
residual = float(np.sqrt(max(0.0, 1.0 - min(1.0, s0**2))))
def proj(x):
x = np.asarray(x, dtype=np.complex128)
if x.shape != (2 * chi,):
raise ValueError(f"x must have shape ({2*chi},), got {x.shape}")
alpha = np.vdot(prod, x)
return alpha * prod
schmidt = np.array([s[0], s[1]], dtype=np.float64)
return proj, residual, schmidt