Skip to content

Latest commit

 

History

History

25-dd-sequence-alignment

Folders and files

NameName
Last commit message
Last commit date

parent directory

..
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

README.md

Quantum shift-search for short-read barcode demultiplexing: amplitude amplification over alignment offsets with classical-data and quantum-data oracles, simulated with decision diagrams up to 79 qubits

Abstract. Quantum sequence-alignment proposals cast "where does this barcode sit in this read" as a Grover search over alignment offsets and demonstrate it on 5–10 qubits. We build both oracle styles the literature blurs together and push them past the statevector limit with decision-diagram (DD) simulation. Oracle A compiles the read and barcode into the circuit as constants and marks the matching offsets directly: it needs only the ⌈log₂ #offsets⌉-qubit offset register and is a classical lookup table in disguise. Oracle B loads the read into a 2·n_r-qubit register, aligns it with a barrel shifter of controlled SWAPs driven by the offset register, compares against the barcode into per-position flags, phases exact matches or ≤ 1-mismatch patterns, and uncomputes. On 108 seeded instances (barcodes of 4–16 bases in reads of 7–23 bases, 0 or 1 planted mismatch) both oracles reproduce sin²((2k+1)θ) to 10⁻¹⁵ on the qcore statevector backend where it fits (≤ 22 qubits) and to within shot noise (max deviation 0.014 at 4 096 shots) on MQT DDSIM up to 79 qubits — where every DD run took 1–32 ms, because the data register is a computational-basis state and the DD stays linear in the qubit count. Oracle B costs 420–5 418 T gates and 48–264 Fredkin gates per Grover iteration for 1–3 iterations, i.e. 0.4–4.6 ms per read at a 10⁶ logical-T/s clock, against 37–67 ns per read for Python's str.find and 75–260 ns for a numpy mismatch-tolerant scan: a 10⁸-read lane takes 4–7 s classically and 10⁴–10⁵ s quantum-mechanically before any error-correction overhead. The √N Grover factor is irrelevant here (N = 4–16 offsets, k_opt = 1–3); what the experiment does establish is that DD simulators, not statevectors, are the right tool for validating oracle circuits of this kind at realistic register sizes.

1. Background & gap

Grover-style sequence alignment appears in several proposals — MAFFT-inspired quantum shift alignment (arXiv:2602.23848), arXiv:2506.22775, arXiv:2603.22245, arXiv:2507.04111, QuASeR's QUBO alignment (Sarkar et al. arXiv:2004.06663), quantum pattern matching (Niroula & Nam, "A quantum algorithm for string matching", npj QI 2021, arXiv:2012.11276) — and all are validated on ≤ 10-qubit statevector simulations, where the read register cannot even be represented. Decision-diagram simulation (Zulehner & Wille, arXiv:1707.00865; MQT DDSIM, Burgholzer et al.) represents states with few distinct amplitudes compactly, which is exactly the situation in an oracle circuit whose data register is classical-valued: the width of the register costs nothing. Barcode/UMI demultiplexing (Illumina index reads; cutadapt) is a real instance with a tiny alphabet, short patterns and a mismatch tolerance of one base.

Gap: nobody has separated the "classical-data" and "quantum-data" oracle constructions, verified them beyond 24 qubits, or costed the per-read quantum query against the classical scan it would replace.

2. Problem formulation

Alphabet {A, C, G, T} → 2 bits per base (A = 00, C = 01, G = 10, T = 11). Read R of n_r bases, barcode Q of n_q bases, offset s ∈ {0 … 2^{n_o} − 1} with n_r = n_q + 2^{n_o} − 1 so that every offset-register value is a legal alignment start. Match(s) ⇔ Hamming(R[s : s + n_q], Q) ≤ m, m ∈ {0, 1}. Grover over the offset register: N = 2^{n_o}, M = #matching offsets, k_opt = round(π/(4 arcsin √(M/N)) − ½).

Registers (oracle.py): O (n_o qubits, bit i = qubit i), D (2 n_r qubits, base j on qubits 2j, 2j+1), F (n_q match flags, oracle B), A (ladder ancillas for multi-controlled Z: max(n_q, n_o) − 2).

  • Oracle A (classical-data): X-conjugated multi-controlled Z on the bit pattern of each matching offset. Qubits n_o + max(n_o − 2, 0).
  • Oracle B (quantum-data): load |R⟩ into D (X gates); barrel shifter — for offset bit b, a controlled cyclic left shift of D by 2^b bases, decomposed into gcd-cycles of adjacent transpositions, each base transposition = 2 Fredkin (CSWAP) gates; comparator — for each barcode position j, X-conjugate D_j to Q_j's code and CCX into F_j; threshold — multi-controlled Z on F for the all-ones pattern (exact) or the n_q + 1 disjoint patterns with ≤ 1 zero (≤ 1 mismatch); uncompute comparator and shifter (all gates self-inverse). Diffusion acts on O only. Qubits n_o + 2n_r + n_q + max(n_q, n_o) − 2.

3. Method

  • Instances (dna.py): seeded random reads with the barcode planted at a random offset with 0 or 1 substitution; M counted classically (extra chance matches are included in M).
  • Statevector: qcore aer for circuits ≤ 24 qubits; P(success) = probability mass of the offset register on matching offsets, for k = 0 … k_opt + 1, compared with theory.
  • Decision diagrams: the qcore circuit is exported with Circuit.to_qasm2(), loaded by Qiskit and run on MQT DDSIM's qasm_simulator (4 096 shots) at k_opt; success fraction and wall-clock recorded; the qcore→QASM bridge is the only place a non-qcore simulator is used.
  • Resources: gates, CSWAP, CCX, depth and T-count (7 per Toffoli/Fredkin) of one Grover iteration from the IR.
  • Classical baselines: Python str.find per read and a vectorised numpy sliding-window Hamming scan (≤ 1 mismatch) over 20 000 random reads; time per read and the 10⁸-read extrapolation.
  • Quantum per-read cost model: k_opt × T/iteration at 10⁶ logical T gates per second (no error-correction or T-factory overhead).

4. Experimental setup

item value
(n_q, n_o) (4,2) (4,3) (6,2) (6,3) (8,2) (8,3) (8,4) (12,3) (16,3) → n_r = 7 … 23 bases
seeds × mismatch 3 × {0, 1} → 108 (instance, oracle) rows
statevector aer, ≤ 24 qubits (oracle A always; oracle B only at (4,2) = 22 qubits)
DDSIM all rows ≤ 80 qubits, 4 096 shots, 120 s total budget (119.4 s unused)
classical timing 20 000 reads per configuration
hardware Apple Silicon laptop, Python 3.13, mqt.ddsim 2.6.0
wall-clock, full run 7 s
wall-clock, --quick ≈ 3 s

5. Results

5.1 Correctness and simulator agreement

  • aer vs theory: max |P(k) − sin²((2k+1)θ)| = 1.3 × 10⁻¹⁵ over all statevector rows and all k.
  • DDSIM vs theory at k_opt: max deviation 0.0137 at 4 096 shots (binomial standard error ≈ 0.004–0.007), over 108 rows up to 79 qubits.
  • Seeds where a chance second match appears (n_q = 4, n_o = 3, m = 1: M = 2, 3, 1) are handled by the M-dependent k_opt and land on theory (0.844–1.000 as predicted).

success vs k

5.2 Resources per Grover iteration (seed 0; the ≤ 1-mismatch threshold adds n_q multi-controlled phases)

n_q n_o n_r oracle qubits CSWAP CCX T / iteration (m = 0) T / iteration (m = 1) depth (m = 0 / 1) k_opt P(k_opt) theory
4 2 7 A 2 0 0 0 0 6 / 6 1 1.000
4 2 7 B 22 48 12 420 532 48 / 71 1 1.000
4 3 11 B 31 120 14 938 1 050 94 / 118 2 0.945
6 3 13 B 39 144 22 1 162 1 498 114 / 171 2 0.945
8 3 15 B 47 168 30 1 386 2 058 130 / 237 2 0.945
8 4 23 B 64 352 32 2 688 3 360 238 / 346 3 0.961
12 3 19 B 63 216 46 1 834 3 514 166 / 421 2 0.945
16 3 23 B 79 264 62 2 282 5 418 202 / 669 2 0.945
any 3 – A 4 0 4 28 28–42 10–17 2 0.945

The barrel shifter dominates oracle B (CSWAP count = 2 × Σ_b (n_r − gcd(n_r, 2^b))); the ≤ 1-mismatch threshold multiplies the phase stage by n_q + 1.

resources

5.3 Decision diagrams vs statevectors

oracle B qubits 22 31 39 47 63–64 79
aer statevector, per circuit 0.25 s – – – – –
DDSIM, k_opt circuit, 4 096 shots 2–4 ms 4–7 ms 6 ms 7–8 ms 10–32 ms 17–21 ms

DD time grows with gate count, not with qubit count: the state is |offset superposition⟩ ⊗ |R⟩ ⊗ |flags⟩ with at most 2^{n_o} distinct branches, so the DD has O(qubits + 2^{n_o}) nodes at every point of the circuit.

DDSIM time vs qubits

5.4 Per-read cost: classical vs quantum

n_r / n_q str.find per read numpy ≤ 1-mismatch scan per read 10⁸ reads classical (find / scan) oracle B: k_opt × T per read at 10⁶ logical T/s 10⁸ reads quantum
7 / 4 37 ns 120 ns 3.7 s / 12 s 420 0.42 ms 4.2 × 10⁴ s
9 / 6 49 ns 104 ns 4.9 s / 10 s 588 0.59 ms 5.9 × 10⁴ s
11 / 8 45 ns 75 ns 4.5 s / 7.5 s 756 0.76 ms 7.6 × 10⁴ s
23 / 16 67 ns 258 ns 6.7 s / 26 s 2 × 2 282 = 4 564 4.6 ms 4.6 × 10⁵ s

6. Discussion & limitations

  1. Oracle A is not a quantum algorithm. Its circuit is the answer: the matching offsets are computed classically and compiled into multi-controlled phases on a 2–6-qubit register. Any demonstration that reports only this register's size has measured the cost of writing down the solution.
  2. Oracle B is the honest construction, and it is cheap to validate — with the right simulator. At 79 qubits the statevector would need 2^79 amplitudes; DDSIM ran it in 21 ms because the read register never leaves the computational basis. This is the useful methodological result: oracle circuits with classical-valued data registers should be validated with decision diagrams at realistic register widths, not truncated to fit a statevector.
  3. There is no speedup to be had. The search space is the offset register: N = n_r − n_q + 1 ≤ 16 here and ≤ 300 for a 300-base read, so √N buys a factor of ≤ 17 in oracle calls, each oracle call costs 10²–10⁴ T gates, and demultiplexing is embarrassingly parallel over 10⁸ reads. Classical scanning of a lane takes seconds; the quantum per-read cost is milliseconds per read before error correction (which adds 10³–10⁴×). The literature's qubit-count comparisons hide this because they never price the oracle.
  4. What would change the picture: a single search over a huge offset space (a full-genome alignment, N ~ 3 × 10⁹) — but then the data register is the genome (6 × 10⁹ qubits with this encoding) and loading it is the QRAM problem.
  5. Limitations. Offset register sized to a power of two (reads padded to n_q + 2^{n_o} − 1); Hamming distance only (no indels); CSWAP/CCX T-count of 7 without routing overhead; DDSIM shot statistics only (no DD node counts exported by the backend — timings are the proxy); no noise.

7. Reproduce

cd experiments/25-dd-sequence-alignment
uv sync
uv run pytest -q                 # swap network, barrel shifter, oracle phases, theory agreement, DDSIM == aer, mismatch scan
uv run python run.py --quick     # (4,2) and (4,3), ~3 s
uv run python run.py             # 108 rows up to 79 qubits, ~7 s

Outputs: results/results.json, results/success_vs_k.png, results/resources_vs_read_length.png, results/ddsim_time_vs_qubits.png.

References

  • Grover. A fast quantum mechanical algorithm for database search. STOC 1996. arXiv:quant-ph/9605043.
  • Niroula, Nam. A quantum algorithm for string matching. npj Quantum Inf. 7, 37 (2021). arXiv:2012.11276.
  • Sarkar et al. QuASeR: quantum accelerated de novo DNA sequence reconstruction. PLoS ONE 2021. arXiv:2004.06663.
  • Quantum sequence alignment 2025–2026: arXiv:2602.23848, arXiv:2506.22775, arXiv:2603.22245, arXiv:2507.04111.
  • Zulehner, Wille. Advanced simulation of quantum computations. IEEE TCAD 38(5), 2019. arXiv:1707.00865.
  • Burgholzer, Bauer, Wille. Hybrid Schrödinger–Feynman simulation with decision diagrams. QCE 2021; MQT DDSIM, mqt.readthedocs.io.
  • Martin. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal 17(1), 2011.
  • Aaronson. Read the fine print. Nat. Phys. 11, 291 (2015).