|
| 1 | +""" |
| 2 | +Fourier series structure of data re-uploading PQC. |
| 3 | +
|
| 4 | +Reference: Schuld et al., "The effect of data encoding on the expressive |
| 5 | +power of variational quantum-classical algorithms", PRA 103, 032430 (2021). |
| 6 | +
|
| 7 | +Circuit: L-layer single-qubit re-uploading. |
| 8 | + Per layer: Ry(θ_l) · Rz(φ_l) · Rz(x) |
| 9 | + Observable: <X> |
| 10 | +
|
| 11 | +Rz(x) is the data-encoding gate. Rz(φ_l+x) = Rz(φ_l)·Rz(x) introduces a |
| 12 | +phase between the |0⟩ and |1⟩ amplitudes that depends on x; the off-diagonal |
| 13 | +observable <X> is sensitive to this phase, making the output a Fourier series |
| 14 | +in x. With L layers, the accessible frequencies are integers {-L,...,0,...,L}, |
| 15 | +giving a truncated Fourier series of degree L. |
| 16 | +
|
| 17 | +All SymbolCircuit construction must happen BEFORE tc.set_backend("jax"). |
| 18 | +""" |
| 19 | + |
| 20 | +import sympy |
| 21 | +import jax |
| 22 | +import jax.numpy as jnp |
| 23 | +import numpy as np |
| 24 | +import tensorcircuit as tc |
| 25 | + |
| 26 | +# ── helpers ──────────────────────────────────────────────────────────────────── |
| 27 | + |
| 28 | + |
| 29 | +def make_jax_fn(expr, syms): |
| 30 | + return sympy.lambdify(syms, expr, modules=[jnp, "math"]) |
| 31 | + |
| 32 | + |
| 33 | +def build_reuploading(L, x_sym): |
| 34 | + """ |
| 35 | + Single-qubit L-layer re-uploading: [Ry(θ_l) · Rz(φ_l) · Rz(x)]^L. |
| 36 | + Rz(φ_l)·Rz(x) = Rz(φ_l+x): trainable phase + data encoding combined. |
| 37 | + Returns (SymbolCircuit, param_syms) where param_syms = [θ_0,...,φ_0,...]. |
| 38 | + """ |
| 39 | + thetas = [sympy.Symbol(f"t{l}", real=True) for l in range(L)] |
| 40 | + phis = [sympy.Symbol(f"p{l}", real=True) for l in range(L)] |
| 41 | + sc = tc.SymbolCircuit(1) |
| 42 | + for l in range(L): |
| 43 | + sc.ry(0, theta=thetas[l]) |
| 44 | + sc.rz(0, theta=phis[l]) # trainable phase |
| 45 | + sc.rz(0, theta=x_sym) # data encoding |
| 46 | + return sc, thetas + phis |
| 47 | + |
| 48 | + |
| 49 | +def fourier_coeffs_sym(expr, x_sym, L): |
| 50 | + """ |
| 51 | + Symbolically extract Fourier coefficients via trig orthogonality: |
| 52 | + a_k = (1/π) ∫₀^{2π} f(x) cos(kx) dx (a_0 uses 1/2π) |
| 53 | + b_k = (1/π) ∫₀^{2π} f(x) sin(kx) dx |
| 54 | + Returns (a_dict, b_dict). Results are rewritten in trig form. |
| 55 | + """ |
| 56 | + pi = sympy.pi |
| 57 | + |
| 58 | + def _trig(e): |
| 59 | + return sympy.trigsimp(sympy.expand(e.rewrite(sympy.cos))) |
| 60 | + |
| 61 | + a, b = {}, {} |
| 62 | + a[0] = _trig(sympy.integrate(expr, (x_sym, 0, 2 * pi)) / (2 * pi)) |
| 63 | + for k in range(1, L + 1): |
| 64 | + a[k] = _trig( |
| 65 | + sympy.integrate(expr * sympy.cos(k * x_sym), (x_sym, 0, 2 * pi)) / pi |
| 66 | + ) |
| 67 | + b[k] = _trig( |
| 68 | + sympy.integrate(expr * sympy.sin(k * x_sym), (x_sym, 0, 2 * pi)) / pi |
| 69 | + ) |
| 70 | + return a, b |
| 71 | + |
| 72 | + |
| 73 | +def dft_power(expr, x_sym, params, param_vals, N=512): |
| 74 | + """ |
| 75 | + Evaluate expr numerically at N evenly-spaced x in [0, 2π], |
| 76 | + then compute the DFT power spectrum. |
| 77 | + """ |
| 78 | + param_dict = dict(zip(params, param_vals)) |
| 79 | + expr_x = expr.subs(param_dict) |
| 80 | + f_x = sympy.lambdify(x_sym, expr_x, modules="numpy") |
| 81 | + x_grid = np.linspace(0, 2 * np.pi, N, endpoint=False) |
| 82 | + y = np.array([float(np.real(f_x(xi))) for xi in x_grid]) |
| 83 | + ck = np.fft.rfft(y) / N |
| 84 | + return np.arange(len(ck)), np.abs(ck) ** 2 |
| 85 | + |
| 86 | + |
| 87 | +# ══ build all circuits BEFORE tc.set_backend ══════════════════════════════════ |
| 88 | + |
| 89 | +x = sympy.Symbol("x", real=True) |
| 90 | +L_MAX = 3 |
| 91 | + |
| 92 | +circuits = {} # L → (sc, params, expr_X) |
| 93 | +for L in range(1, L_MAX + 1): |
| 94 | + sc, params = build_reuploading(L, x) |
| 95 | + # <X> = expectation of Pauli X on qubit 0 |
| 96 | + expr_X = sc.expectation_ps(x=[0]) |
| 97 | + circuits[L] = (sc, params, expr_X) |
| 98 | + |
| 99 | + |
| 100 | +# ── 1. Symbolic expressions ──────────────────────────────────────────────────── |
| 101 | + |
| 102 | +print("=== Symbolic <X>(x; θ, φ) ===") |
| 103 | +print() |
| 104 | +for L in range(1, L_MAX + 1): |
| 105 | + _, params, expr = circuits[L] |
| 106 | + fully = sympy.trigsimp(sympy.expand(expr)) |
| 107 | + n_add = len(fully.as_ordered_terms()) |
| 108 | + print(f"L={L}: {n_add} additive terms") |
| 109 | + print(f" {fully}") |
| 110 | +print() |
| 111 | + |
| 112 | + |
| 113 | +# ── 2. Symbolic Fourier coefficients for L=1 ────────────────────────────────── |
| 114 | +# L=1: <X> = sin(t0)·cos(x + p0) = sin(t0)cos(p0)·cos(x) − sin(t0)sin(p0)·sin(x) |
| 115 | +# Both a_1 and b_1 are controlled independently via t0 (amplitude) and p0 (phase). |
| 116 | + |
| 117 | +print("=== Symbolic Fourier coefficients, L=1 ===") |
| 118 | +_, params1, expr1 = circuits[1] |
| 119 | +a, b = fourier_coeffs_sym(expr1, x, L=1) |
| 120 | +print(f" a_0 = {a[0]}") |
| 121 | +print(f" a_1 = {a[1]} (cos x coefficient)") |
| 122 | +print(f" b_1 = {b[1]} (sin x coefficient)") |
| 123 | + |
| 124 | +# Verify: frequency-2 coefficient is zero for L=1 |
| 125 | +a2_check = sympy.simplify( |
| 126 | + sympy.integrate(expr1 * sympy.cos(2 * x), (x, 0, 2 * sympy.pi)) / sympy.pi |
| 127 | +) |
| 128 | +print(f" a_2 = {a2_check} (should be 0 — freq-2 inaccessible at L=1)") |
| 129 | +print() |
| 130 | + |
| 131 | + |
| 132 | +# ── 3. Numerical Fourier power spectrum (fixed random params) ───────────────── |
| 133 | +# Peaks above frequency L must be zero regardless of parameter values. |
| 134 | + |
| 135 | +print("=== Fourier power spectrum (DFT, fixed random params) ===") |
| 136 | +rng = np.random.default_rng(0) |
| 137 | +for L in range(1, L_MAX + 1): |
| 138 | + _, params, expr = circuits[L] |
| 139 | + param_vals = rng.uniform(0.3, 2.0, size=len(params)) |
| 140 | + freqs, power = dft_power(expr, x, params, param_vals) |
| 141 | + shown = " ".join(f"k={k}: {power[k]:.4f}" for k in range(L + 3)) |
| 142 | + print(f" L={L}: {shown}") |
| 143 | +print() |
| 144 | + |
| 145 | + |
| 146 | +# ── 4. JAX gradient fitting ──────────────────────────────────────────────────── |
| 147 | +# Target: sin(2x) — a pure second harmonic. |
| 148 | +# L=1: b_2 = 0 identically → cannot represent sin(2x) → MSE stays at ~0.5. |
| 149 | +# L=2: b_2 free → can represent sin(2x) → MSE → 0. |
| 150 | +# L=3: same, with more params. |
| 151 | +# |
| 152 | +# Analytically (from the derivation): |
| 153 | +# L=2 achieves sin(2x) exactly with: |
| 154 | +# t1=0, t0=π/2, p0+p1=−π/2 → <X> = sin(2x) |
| 155 | + |
| 156 | +tc.set_backend("jax") |
| 157 | + |
| 158 | +target_fn = lambda xv: jnp.sin(2.0 * xv) |
| 159 | + |
| 160 | + |
| 161 | +def fit(L, n_steps=2000, lr=0.05, seed=0): |
| 162 | + _, params, expr = circuits[L] |
| 163 | + all_syms = [x] + params |
| 164 | + f_jax = make_jax_fn(expr, all_syms) |
| 165 | + |
| 166 | + # batch evaluation: vmap over x, broadcast params |
| 167 | + def predict(pv, xb): |
| 168 | + return jax.vmap(lambda xi: jnp.real(f_jax(xi, *pv)))(xb) |
| 169 | + |
| 170 | + @jax.jit |
| 171 | + def loss(pv, xb, yb): |
| 172 | + return jnp.mean((predict(pv, xb) - yb) ** 2) |
| 173 | + |
| 174 | + grad_loss = jax.jit(jax.grad(loss)) |
| 175 | + |
| 176 | + # try multiple seeds, keep best |
| 177 | + best_mse, best_pv = jnp.inf, None |
| 178 | + x_train = jnp.linspace(0.0, 2.0 * jnp.pi, 80) |
| 179 | + y_train = target_fn(x_train) |
| 180 | + for s in range(3): |
| 181 | + key = jax.random.PRNGKey(seed + s) |
| 182 | + pv = jax.random.uniform(key, (len(params),), minval=0.1, maxval=2.0) |
| 183 | + for _ in range(n_steps): |
| 184 | + pv = pv - lr * grad_loss(pv, x_train, y_train) |
| 185 | + mse = float(loss(pv, x_train, y_train)) |
| 186 | + if mse < best_mse: |
| 187 | + best_mse, best_pv = mse, pv |
| 188 | + |
| 189 | + return best_mse, np.array(best_pv) |
| 190 | + |
| 191 | + |
| 192 | +print("=== Gradient fitting of target sin(2x) ===") |
| 193 | +print(" (MSE ≈ 0.5 = Var[sin(2x)] means the circuit cannot represent the target)") |
| 194 | +for L in range(1, L_MAX + 1): |
| 195 | + mse, pv = fit(L, n_steps=2000) |
| 196 | + print(f" L={L} ({len(pv)} params): final MSE = {mse:.6f}") |
| 197 | +print() |
| 198 | + |
| 199 | + |
| 200 | +# ── 5. Analytical verification for L=2 ──────────────────────────────────────── |
| 201 | +# Show the exact parameter values that make <X> = sin(2x), then cross-check |
| 202 | +# with the numerical tc circuit. |
| 203 | + |
| 204 | +print("=== Analytical solution: L=2, <X> = sin(2x) ===") |
| 205 | +# t1=0, t0=π/2, p0+p1=−π/2 → pick p0=0, p1=−π/2 |
| 206 | +t0_sol, t1_sol = np.pi / 2, 0.0 |
| 207 | +p0_sol, p1_sol = 0.0, -np.pi / 2 |
| 208 | + |
| 209 | +_, params2, expr2 = circuits[2] |
| 210 | +all_syms2 = [x] + params2 |
| 211 | +f2 = make_jax_fn(expr2, all_syms2) |
| 212 | + |
| 213 | +x_check = jnp.linspace(0.0, 2.0 * jnp.pi, 200) |
| 214 | +pv_sol = jnp.array([t0_sol, t1_sol, p0_sol, p1_sol]) |
| 215 | +y_model = jax.vmap(lambda xi: jnp.real(f2(xi, *pv_sol)))(x_check) |
| 216 | +y_target = target_fn(x_check) |
| 217 | +mse_sol = float(jnp.mean((y_model - y_target) ** 2)) |
| 218 | +print(f" params: t0={t0_sol:.4f}, t1={t1_sol:.4f}, p0={p0_sol:.4f}, p1={p1_sol:.4f}") |
| 219 | +print(f" MSE against sin(2x): {mse_sol:.2e} (numerical zero)") |
| 220 | + |
| 221 | +# Cross-check with tc numerical circuit (bind x too, since it's also a gate param) |
| 222 | +sc2, params2, _ = circuits[2] |
| 223 | +param_dict_sol = dict(zip(params2, [t0_sol, t1_sol, p0_sol, p1_sol])) |
| 224 | +x_sample = 0.7 |
| 225 | +c_num = sc2.to_circuit({**param_dict_sol, x: x_sample}) |
| 226 | +ref = float(jnp.real(c_num.expectation_ps(x=[0]))) |
| 227 | +model_val = float(jnp.real(f2(jnp.array(x_sample), *pv_sol))) |
| 228 | +print(f" at x=0.7: lambdify = {model_val:.8f}, tc.Circuit = {ref:.8f}") |
0 commit comments