Samarth0710/qksd-tfim-qiskit-demo

Quantum Krylov Subspace Diagonalization (QKSD) on the TFIM

Real-time-evolution Quantum Krylov diagonalization for the ground-state energy of a 4-qubit transverse-field Ising model at criticality (J = h = 1, open boundaries), implemented from scratch in Qiskit 2.x.

References:

• Kirby, Analysis of quantum Krylov algorithms with errors, arXiv:2401.01246.
• Klymko, Williams-Young, Sud, de Jong, Tubman, Real-Time Krylov Theory, arXiv:2208.01063.
• Cortes & Gray, QKSD for ground and excited states, arXiv:2109.06868.
• Lee, Choi, Huh, Izmaylov, Efficient strategies for reducing QKSD sampling error, arXiv:2409.02504.

Algorithm

Pick reference state |φ_0⟩ = |+⟩^⊗4 (overlap with the true ground state ~ 0.81). Build the non-orthogonal Krylov basis

V = [|φ_0⟩, U(Δt)|φ_0⟩, U(2Δt)|φ_0⟩, U(3Δt)|φ_0⟩],   U(t) = e^(−iHt).

Form the (D × D) Toeplitz matrices

S_jk = ⟨φ_0|U((k−j)Δt)|φ_0⟩
H_jk = ⟨φ_0|U((k−j)Δt) H|φ_0⟩

and solve the generalized eigenvalue problem H c = E S c with thresholded SVD on S (truncate eigenvalues < ε λ_max(S)) to handle ill-conditioning. The smallest E approximates E_0.

Each matrix element is measured by a Hadamard test: an ancilla in |+⟩, a controlled-W where W is either U(mΔt) (for S) or U(mΔt) P_α (for each Pauli term P_α in H), and read out ⟨Z_anc⟩ → real part, ⟨Y_anc⟩ → imaginary part. The script measures all 7 Pauli terms × 7 distinct lags × 2 (real / imag) per shot count.

Time evolution is implemented with PauliEvolutionGate (LieTrotter synthesis); for this model the Trotter representation is exact to machine precision.

Three modes compared

mode	matrix elements	purpose
A	exact statevector	algorithm-only test (Krylov truncation only)
B	Hadamard test, exact (no shots)	sanity-check the Hadamard-test wiring
C	Hadamard test, finite shots	full quantum simulation

Results

Exact ground-state energy of the TFIM(N=4, J=h=1): E_0 = −4.75877.

method	E_min	error
classical exact diag	−4.75877	0
Krylov A (classical, D=4)	−4.74338	+0.0154
Krylov B (exact Hadamard, D=4)	−4.74338	+0.0154
Krylov C (sampled, 8000 sh./elem.)	−5.0958	−0.337 ← non-variational
Krylov C (sampled, 8000 sh., eps=0.1, kept=2)	−4.62	+0.14 (variational)

Convergence with Krylov dimension D (mode A)

D	E_min	error
1	−4.0000	0.759
2	−4.5907	0.168
3	−4.7005	0.058
4	−4.7434	0.015

Roughly exponential convergence — exactly what real-time Krylov is known to give for gapped Hamiltonians.

What goes wrong under shot noise

S is genuinely ill-conditioned (eigenvalues 0.029, 0.24, 0.45, 3.28 → cond ≈ 110). With finite shots, S_C is a noisy version of S_A; even at 8000 shots per element the noise on the smallest eigenvalue (~0.03) is comparable to its magnitude. The regularized GEVP then explores a slightly wrong subspace and can give a non-variational eigenvalue below the true ground state — this is exactly the pathology discussed in arXiv:2401.01246 and 2409.02504. Aggressive truncation (eps = 0.1, keep only 2 directions) restores variationality but sacrifices resolution.

The shots scan also shows that simply throwing more shots is not enough: even at 32 000 shots per element, the best |error| stays ~0.12. The condition number of S is the bottleneck, not raw statistical precision. Mitigations the literature proposes: thresholded GEVP (Epperly et al.), measurement-efficient construction of Krylov matrices (Zhang et al. 2301.13353), measurement-allocation strategies (Lee et al. 2409.02504), or symmetry-aware schemes (Cortes & Gray).

Files

qksd_tfim.py
artifacts/
├── results.json           # all matrix elements, energies, S spectra
├── run.log                # full execution log (everything printed above + per-element values)
├── bar.png                # 4-bar comparison
├── convergence_D.png      # |E−E_exact| vs D (log-y)
├── shots_scan.png         # |E−E_exact| vs shots (log-log, mode C)
└── S_matrices.png         # |S| heatmaps for modes A, B, C
README.md

Run

pip install "qiskit~=2.4" qiskit-aer scipy matplotlib
python qksd_tfim.py

Tested with qiskit==2.4.1, qiskit-aer==0.17.2. Runtime ~50 s on a CPU sandbox.

Bugs I hit while building this (so they're documented)

1. Hand-rolled Trotter circuit had wrong sign convention — gave a spectral-norm error of ~2 (independent of reps!) which made modes B and C wildly wrong on the first run. Fix: use Qiskit's PauliEvolutionGate + LieTrotter.
2. Setting AerSimulator(seed_simulator=…) at construction reuses the same RNG state on every .run() call, biasing aggregated estimators. Fix: pass a fresh seed_simulator= per .run() call (already carried over from the previous shadows demo).