Skip to content

Commit 3b95f34

Browse files
mhjensenmhjensen
committed
update
1 parent 3614b54 commit 3b95f34

File tree

4 files changed

+345
-1
lines changed

4 files changed

+345
-1
lines changed

‎doc/src/week9/programs/qft_matrix_simple.png‎

58.2 KB
Loading

‎doc/src/week9/programs/qft_simple.py‎

Lines changed: 344 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,344 @@
1+
"""
2+
Quantum Fourier Transform (QFT) — Simple Functional Implementation
3+
==================================================================
4+
Functions:
5+
- qft_matrix(n) : build the n-qubit QFT matrix
6+
- iqft_matrix(n) : build the n-qubit inverse QFT matrix
7+
- apply_qft(state) : apply QFT to a state vector
8+
- apply_iqft(state) : apply IQFT to a state vector
9+
- basis_state(idx, n) : make a computational basis vector |idx⟩
10+
- random_state(n) : make a random normalised state vector
11+
- run_unitarity_tests(n): run all tests and print + plot results
12+
"""
13+
14+
import numpy as np
15+
import matplotlib
16+
matplotlib.use('Agg')
17+
import matplotlib.pyplot as plt
18+
import matplotlib.gridspec as gridspec
19+
20+
21+
# ─────────────────────────────────────────────────────────────────────────────
22+
# 1. Matrix construction
23+
# ─────────────────────────────────────────────────────────────────────────────
24+
25+
def qft_matrix(n):
26+
"""
27+
Return the 2^n × 2^n QFT matrix.
28+
F[j,k] = exp(+2πi·jk / 2^n) / sqrt(2^n)
29+
"""
30+
dim = 2 ** n
31+
omega = np.exp(2j * np.pi / dim)
32+
idx = np.arange(dim)
33+
return np.power(omega, np.outer(idx, idx)) / np.sqrt(dim)
34+
35+
36+
def iqft_matrix(n):
37+
"""
38+
Return the 2^n × 2^n inverse QFT matrix.
39+
F†[j,k] = exp(−2πi·jk / 2^n) / sqrt(2^n) = conj(F[j,k])
40+
"""
41+
return qft_matrix(n).conj()
42+
43+
44+
# ─────────────────────────────────────────────────────────────────────────────
45+
# 2. Applying the transforms
46+
# ─────────────────────────────────────────────────────────────────────────────
47+
48+
def apply_qft(state):
49+
"""Apply the QFT to a state vector. Returns transformed vector."""
50+
n = int(round(np.log2(len(state))))
51+
return qft_matrix(n) @ np.asarray(state, dtype=complex)
52+
53+
54+
def apply_iqft(state):
55+
"""Apply the inverse QFT to a state vector. Returns transformed vector."""
56+
n = int(round(np.log2(len(state))))
57+
return iqft_matrix(n) @ np.asarray(state, dtype=complex)
58+
59+
60+
# ─────────────────────────────────────────────────────────────────────────────
61+
# 3. Convenience state builders
62+
# ─────────────────────────────────────────────────────────────────────────────
63+
64+
def basis_state(idx, n):
65+
"""Return the computational basis vector |idx⟩ for an n-qubit system."""
66+
v = np.zeros(2 ** n, dtype=complex)
67+
v[idx] = 1.0
68+
return v
69+
70+
71+
def random_state(n, seed=None):
72+
"""Return a Haar-random normalised state vector for an n-qubit system."""
73+
rng = np.random.default_rng(seed)
74+
v = rng.standard_normal(2**n) + 1j * rng.standard_normal(2**n)
75+
return v / np.linalg.norm(v)
76+
77+
78+
# ─────────────────────────────────────────────────────────────────────────────
79+
# 4. Unitarity tests
80+
# ─────────────────────────────────────────────────────────────────────────────
81+
82+
def run_unitarity_tests(n, n_random=500, seed=0, plot=True,
83+
savename_tests="unitarity_tests_simple.png",
84+
savename_matrix="qft_matrix_simple.png"):
85+
"""
86+
Run a full suite of unitarity tests on the n-qubit QFT, print a
87+
pass/fail summary, and optionally save two diagnostic figures.
88+
"""
89+
dim = 2 ** n
90+
F = qft_matrix(n)
91+
Fd = iqft_matrix(n)
92+
I = np.eye(dim)
93+
rng = np.random.default_rng(seed)
94+
95+
# ── 1 & 2: product matrices ────────────────────────────────────────────
96+
FFd_err = np.abs(F @ Fd - I)
97+
FdF_err = np.abs(Fd @ F - I)
98+
99+
# ── 3 & 4: column / row orthonormality ────────────────────────────────
100+
col_err = np.abs(F.conj().T @ F - I)
101+
row_err = np.abs(F @ F.conj().T - I)
102+
103+
# ── 5: norm preservation over random states ────────────────────────────
104+
def _rand_norm():
105+
v = rng.standard_normal(dim) + 1j * rng.standard_normal(dim)
106+
v /= np.linalg.norm(v)
107+
return v
108+
109+
norm_errors = np.array([
110+
abs(np.linalg.norm(apply_qft(_rand_norm())) - 1.0)
111+
for _ in range(n_random)
112+
])
113+
114+
# ── 6: round-trip fidelity ─────────────────────────────────────────────
115+
fidelities = []
116+
for _ in range(n_random):
117+
v = rng.standard_normal(dim) + 1j * rng.standard_normal(dim)
118+
v /= np.linalg.norm(v)
119+
vrt = apply_iqft(apply_qft(v))
120+
fidelities.append(abs(np.vdot(v, vrt)) ** 2)
121+
fidelities = np.array(fidelities)
122+
123+
# ── 7: singular values ─────────────────────────────────────────────────
124+
svs = np.linalg.svd(F, compute_uv=False)
125+
sv_dev = np.max(np.abs(svs - 1.0))
126+
127+
# ── 8: eigenvalues ─────────────────────────────────────────────────────
128+
eigvals = np.linalg.eigvals(F)
129+
ev_dev = np.max(np.abs(np.abs(eigvals) - 1.0))
130+
131+
# ── Print summary ──────────────────────────────────────────────────────
132+
tol = 1e-10
133+
print("=" * 60)
134+
print(f" Unitarity Tests — {n}-qubit QFT (dim = {dim})")
135+
print("=" * 60)
136+
rows = [
137+
("F·F† = I (max error)", FFd_err.max(), tol, True),
138+
("F†·F = I (max error)", FdF_err.max(), tol, True),
139+
("Column orthonormality (max err)", col_err.max(), tol, True),
140+
("Row orthonormality (max err)", row_err.max(), tol, True),
141+
("Norm preservation (max err)", norm_errors.max(), tol, True),
142+
("Round-trip fidelity (min)", fidelities.min(), None, False),
143+
("Singular values ≈ 1 (max dev)", sv_dev, tol, True),
144+
("|eigenvalue| ≈ 1 (max dev)", ev_dev, tol, True),
145+
]
146+
for name, val, threshold, low_is_good in rows:
147+
if threshold is None:
148+
ok = abs(val - 1.0) < 1e-10
149+
else:
150+
ok = val < threshold
151+
print(f" {'✓' if ok else '✗'} {name:<42s} {val:.3e}")
152+
print("=" * 60)
153+
154+
if not plot:
155+
return
156+
157+
# ── Styling ────────────────────────────────────────────────────────────
158+
BG, PANEL, ACCENT = "#0d0f1a", "#131629", "#4fc3f7"
159+
WARM, GREEN, GRID = "#ff7043", "#69f0ae", "#1e2340"
160+
TEXT = "#e8eaf6"
161+
162+
plt.rcParams.update({
163+
'figure.facecolor': BG, 'axes.facecolor': PANEL,
164+
'axes.edgecolor': GRID, 'axes.labelcolor': TEXT,
165+
'xtick.color': TEXT, 'ytick.color': TEXT,
166+
'text.color': TEXT, 'grid.color': GRID,
167+
'grid.linewidth': 0.6,
168+
})
169+
170+
# ════════════════════════════════════════════════════════════════════════
171+
# Figure 1: nine-panel diagnostic grid
172+
# ════════════════════════════════════════════════════════════════════════
173+
fig = plt.figure(figsize=(20, 15), facecolor=BG)
174+
fig.suptitle(f"Unitarity Tests — {n}-qubit QFT (dim = {dim})",
175+
fontsize=17, fontweight='bold', color=TEXT, y=0.98)
176+
gs = gridspec.GridSpec(3, 3, figure=fig,
177+
hspace=0.52, wspace=0.38,
178+
left=0.06, right=0.97, top=0.93, bottom=0.06)
179+
180+
def make_ax(row, col, **kw):
181+
a = fig.add_subplot(gs[row, col], **kw)
182+
a.grid(True, alpha=0.35)
183+
return a
184+
185+
# 1 — F·F† heatmap
186+
a1 = make_ax(0, 0)
187+
im1 = a1.imshow(FFd_err, cmap='inferno', vmin=0,
188+
vmax=max(FFd_err.max(), 1e-16),
189+
interpolation='nearest', aspect='auto')
190+
plt.colorbar(im1, ax=a1, label='|error|')
191+
a1.set_title(f"F·F† − I (max {FFd_err.max():.2e})",
192+
color=TEXT, fontweight='bold')
193+
a1.set_xlabel("column j"); a1.set_ylabel("row i")
194+
195+
# 2 — F†·F heatmap
196+
a2 = make_ax(0, 1)
197+
im2 = a2.imshow(FdF_err, cmap='inferno', vmin=0,
198+
vmax=max(FdF_err.max(), 1e-16),
199+
interpolation='nearest', aspect='auto')
200+
plt.colorbar(im2, ax=a2, label='|error|')
201+
a2.set_title(f"F†·F − I (max {FdF_err.max():.2e})",
202+
color=TEXT, fontweight='bold')
203+
a2.set_xlabel("column j"); a2.set_ylabel("row i")
204+
205+
# 3 — column orthonormality heatmap
206+
a3 = make_ax(0, 2)
207+
im3 = a3.imshow(col_err, cmap='plasma', vmin=0,
208+
vmax=max(col_err.max(), 1e-16),
209+
interpolation='nearest', aspect='auto')
210+
plt.colorbar(im3, ax=a3, label='|error|')
211+
a3.set_title("Column orthonormality error",
212+
color=TEXT, fontweight='bold')
213+
a3.set_xlabel("column j"); a3.set_ylabel("column i")
214+
215+
# 4 — norm preservation histogram
216+
a4 = make_ax(1, 0)
217+
a4.hist(norm_errors, bins=40, color=ACCENT, edgecolor=BG, alpha=0.85)
218+
a4.axvline(norm_errors.mean(), color=WARM, lw=2,
219+
label=f"mean={norm_errors.mean():.2e}")
220+
a4.axvline(norm_errors.max(), color=GREEN, lw=2, ls='--',
221+
label=f"max={norm_errors.max():.2e}")
222+
a4.set_title(f"Norm preservation error\n({n_random} random states)",
223+
color=TEXT, fontweight='bold')
224+
a4.set_xlabel("|‖F|ψ⟩‖ − 1|"); a4.set_ylabel("count")
225+
a4.legend(fontsize=9)
226+
227+
# 5 — round-trip infidelity histogram
228+
a5 = make_ax(1, 1)
229+
infid = 1 - fidelities
230+
a5.hist(infid, bins=40, color=GREEN, edgecolor=BG, alpha=0.85)
231+
a5.axvline(infid.mean(), color=WARM, lw=2,
232+
label=f"mean={infid.mean():.2e}")
233+
a5.set_title("Round-trip infidelity 1−F\n(IQFT∘QFT, random states)",
234+
color=TEXT, fontweight='bold')
235+
a5.set_xlabel("1 − |⟨ψ|IQFT(QFT|ψ⟩)|²"); a5.set_ylabel("count")
236+
a5.legend(fontsize=9)
237+
238+
# 6 — singular values scatter
239+
a6 = make_ax(1, 2)
240+
a6.scatter(range(len(svs)), svs, color=ACCENT, s=14, alpha=0.75, zorder=3)
241+
a6.axhline(1.0, color=GREEN, lw=1.5, ls='--', label='σ = 1 (ideal)')
242+
a6.fill_between(range(len(svs)), 1-1e-10, 1+1e-10,
243+
color=GREEN, alpha=0.12)
244+
a6.set_title(f"Singular values of F\n(max dev: {sv_dev:.2e})",
245+
color=TEXT, fontweight='bold')
246+
a6.set_xlabel("index"); a6.set_ylabel("σ")
247+
a6.legend(fontsize=9)
248+
249+
# 7 — eigenvalues on unit circle
250+
a7 = make_ax(2, 0, aspect='equal')
251+
th = np.linspace(0, 2*np.pi, 400)
252+
a7.plot(np.cos(th), np.sin(th), color=GRID, lw=1.5, zorder=1)
253+
sc = a7.scatter(eigvals.real, eigvals.imag,
254+
c=np.angle(eigvals), cmap='hsv',
255+
s=28, alpha=0.85, zorder=3,
256+
vmin=-np.pi, vmax=np.pi)
257+
plt.colorbar(sc, ax=a7, label='arg(λ) [rad]')
258+
a7.axhline(0, color=GRID, lw=0.8); a7.axvline(0, color=GRID, lw=0.8)
259+
a7.set_title(f"Eigenvalues on unit circle\n(max |λ|−1: {ev_dev:.2e})",
260+
color=TEXT, fontweight='bold')
261+
a7.set_xlabel("Re(λ)"); a7.set_ylabel("Im(λ)")
262+
263+
# 8 — QFT of basis states
264+
a8 = make_ax(2, 1)
265+
for idx in range(min(4, dim)):
266+
v = basis_state(idx, n)
267+
a8.plot(np.abs(apply_qft(v))**2, lw=1.5, alpha=0.8, label=f"|{idx}⟩")
268+
a8.set_title("QFT probability distribution\n(basis states |0⟩…|3⟩)",
269+
color=TEXT, fontweight='bold')
270+
a8.set_xlabel("output basis state"); a8.set_ylabel("probability")
271+
a8.legend(fontsize=9)
272+
273+
# 9 — row orthonormality heatmap
274+
a9 = make_ax(2, 2)
275+
im9 = a9.imshow(row_err, cmap='magma', vmin=0,
276+
vmax=max(row_err.max(), 1e-16),
277+
interpolation='nearest', aspect='auto')
278+
plt.colorbar(im9, ax=a9, label='|error|')
279+
a9.set_title("Row orthonormality error",
280+
color=TEXT, fontweight='bold')
281+
a9.set_xlabel("row j"); a9.set_ylabel("row i")
282+
283+
plt.savefig(savename_tests, dpi=150, bbox_inches='tight', facecolor=BG)
284+
plt.close()
285+
print(f"Saved {savename_tests}")
286+
287+
# ════════════════════════════════════════════════════════════════════════
288+
# Figure 2: QFT matrix visualisation (Re, Im, phase)
289+
# ════════════════════════════════════════════════════════════════════════
290+
fig2, axes = plt.subplots(1, 3, figsize=(18, 6), facecolor=BG)
291+
for ax_i, (data, cmap, title) in zip(axes, [
292+
(F.real, 'RdBu_r', 'Re(F_{jk})'),
293+
(F.imag, 'RdBu_r', 'Im(F_{jk})'),
294+
(np.angle(F), 'hsv', 'arg(F_{jk}) [rad]'),
295+
]):
296+
vabs = max(np.abs(data).max(), 1e-16)
297+
im = ax_i.imshow(data, cmap=cmap, vmin=-vabs, vmax=vabs,
298+
interpolation='nearest', aspect='auto')
299+
plt.colorbar(im, ax=ax_i)
300+
ax_i.set_title(title, color=TEXT, fontweight='bold', fontsize=13)
301+
ax_i.set_xlabel("column k"); ax_i.set_ylabel("row j")
302+
303+
fig2.suptitle(f"QFT Matrix — {n} qubits (dim = {dim})",
304+
fontsize=15, fontweight='bold', color=TEXT, y=1.01)
305+
plt.tight_layout()
306+
plt.savefig(savename_matrix, dpi=150, bbox_inches='tight', facecolor=BG)
307+
plt.close()
308+
print(f"Saved {savename_matrix}")
309+
310+
plt.rcParams.update(plt.rcParamsDefault)
311+
312+
313+
# ─────────────────────────────────────────────────────────────────────────────
314+
# 5. Demo
315+
# ─────────────────────────────────────────────────────────────────────────────
316+
317+
if __name__ == "__main__":
318+
n = 4
319+
dim = 2 ** n
320+
321+
print("=" * 60)
322+
print(f" QFT / IQFT Demo — {n} qubits (dim = {dim})")
323+
print("=" * 60)
324+
325+
# Basis state |5⟩
326+
s5 = basis_state(5, n)
327+
fs5 = apply_qft(s5)
328+
rs5 = apply_iqft(fs5)
329+
fid = abs(np.vdot(s5, rs5)) ** 2
330+
print(f"\n|5⟩ round-trip fidelity : {fid:.15f}")
331+
print(f"QFT(|5⟩) amplitudes : {np.round(fs5, 4)}")
332+
333+
# Random state round-trip
334+
rnd = random_state(n, seed=42)
335+
fid2 = abs(np.vdot(rnd, apply_iqft(apply_qft(rnd)))) ** 2
336+
print(f"\nRandom state round-trip : {fid2:.15f}")
337+
338+
# Confirm IQFT = QFT†
339+
diff = np.max(np.abs(iqft_matrix(n) - qft_matrix(n).conj().T))
340+
print(f"‖IQFT − QFT†‖_max : {diff:.3e} (should be ~0)")
341+
342+
# Full unitarity test suite + plots
343+
print()
344+
run_unitarity_tests(n, n_random=500, seed=0)

‎doc/src/week9/programs/unitarity_tests_simple.png‎

315 KB
Loading

‎doc/src/week9/week9.do.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1254,7 +1254,7 @@ This is just the product representation from earlier, obviously our desired outp
12541254

12551255
!split
12561256
===== QFT example codes =====
1257-
This code sets up the QFT for $n$ qubits.
1257+
This code sets up the QFT for $n$ qubits. Note it is not the way you would encode this on an actual quantum computer
12581258
!bc pycod
12591259
import numpy as np
12601260
def qft(n):

0 commit comments

Comments
 (0)