BQN's Quantum Noise

Preamble

We will implement and test a Schrödinger-style1 quantum simulator in the BQN2 programming language. Initially, we import the necessary system functions and define a 1-modifier for handling complex-valued matrix products. Next, we define a namespace containing various quantum gates: 

Sin‿Cos‿GCD ← •math
U ← •rand.Range
_cp ← {(-´𝔽¨)⋈(+´𝔽¨)⟜⌽}

g ← {
  IM ← (⊢⋈≢⥊⟜0⊢)¨
  x‿id‿h ⇐ (⌽‿⊢{𝕎𝕩}¨<=⌜˜↕2) ∾○IM <-⌾(1‿1⊸⊑)2‿2⥊÷√2
  swap‿cnot ⇐ IM ⟨1‿2, 2‿3⟩ {⌽⌾(𝕨⊸⊏)𝕩}¨ <=⌜˜↕4
  P ⇐ {⟨=⌜˜↕4,  4‿4⥊0⟩ {𝕨⌾(3‿3⊸⊑)𝕩}¨˜ Sin⊸⋈⟜Cos𝕩}
}

Interpreter

The (282 chars3) quantum interpreter is based on references arXiv:1711.02086 and arXiv:1608.03355. For simplicity, we always measure at the end of the execution: 

Q ← {𝕨𝕊nq‿qc:
  wf ← 1⌾⊑⊸⋈0⥊˜2⋆nq
  M ← +˝∘×⎉1‿∞ _cp
  K ← {1𝕊𝕩:𝕩; 𝕨𝕊1:𝕨; 𝕨∾∘×⟜<_cp𝕩}
  L ← {𝕩K⍟𝕨1}{K´⟨g.id𝔽˜nq-𝕨+⌊2⋆⁼≠⊑𝕩, 𝕩, 𝕨𝔽g.id⟩}
  T ← ∾↕∘≠{(⊢∾1↓⌽)𝕨↓↕𝕨⊑⟜𝕩˜•_while_>𝕨⊑𝕩}¨<
  A ← {1=≠⊑𝕨 ? 𝕩M˜⊑⊸L´𝕨;
    𝕩M´⟨0⊸L´𝕨⟩(⊢∾∾⟜⌽)L⟜g.swap¨T𝕨(⌽∘⊢∾¬∘∊/⊣)⟜⊑˜↕nq
  }
  »⊸<∨`0>𝕨-`+○(ט)´wf A´ ⌽qc
}

Physical context

The interpreter’s logic is straightforward: it evolves the wavefunction \(\Psi_0^n\) in the full Hilbert space, whose dimension is \(2^n\) for \(n\) qubits. Each quantum gate is extended to this space through Kronecker lifting, performed by the function L:

\begin{equation*} L_U = I \otimes \cdots U \cdots \otimes I \end{equation*}

Note that the function L performs more than a plain Kronecker embedding. Before constructing \(L_U\), the affected qubits are permuted to become index-adjacent in the tensor order. This permutation is decomposed into adjacent swaps, applied as unitary transpositions (T). The system then evolves under the sequence of lifted gates, defined by A:

\begin{equation*} \Psi_m^n = \Psi_0^n \prod_i^m L_{U_i} \end{equation*}

Measurement is implemented by sampling from the CDF of the squared amplitudes of the state vector. Statistically, this corresponds to drawing a discrete random variable \(q_j\) whose cumulative distribution is:

\begin{equation*} F(q_j) = \sum_{i \le j}p(q_i) \end{equation*}

Given a uniform random variable \(U \sim \text{Uniform}(0, 1)\), we sample by inversion:

\begin{equation*} \text{Pr}(U \in (F(q_{j-1}), F(q_{j})]) = \text{Pr}(F(q_{j-1}) < U \le F(q_{j})) = F(q_{j}) - F(q_{j-1}) = p(q_{j}) \end{equation*}

Shor's algorithm

As a test case, we employ the quantum circuit of Shor's algorithm for the number fifteen and base eleven, following references arXiv:1804.03719 and arXiv:2306.09122. The resulting compiled circuit uses five qubits, three of which serve as control. To enhance statistical accuracy, the experiment is repeated multiple times. 

n‿a‿nq‿r ← ⟨15, 11, 5, 0 U˜ 2⋆10⟩

sc ← ⟨
  ⟨0⟩‿g.h ⋄ ⟨1⟩‿g.h ⋄ ⟨2⟩‿g.h
  ⟨2, 3⟩‿g.cnot ⋄ ⟨2, 4⟩‿g.cnot ⋄ ⟨1⟩‿g.h
  ⟨⟨1, 0⟩, g.P π÷2⟩ ⋄ ⟨0⟩‿g.h
  ⟨⟨1, 2⟩, g.P π÷4⟩ ⋄ ⟨⟨0, 2⟩, g.P π÷2⟩ ⋄ ⟨2⟩‿g.h
⟩

C ← {n (⊣≡×´∘GCD) +‿-{𝕩𝕎1}¨ <a⋆(≠÷2×⊑∘⍒) 0⌾⊑+˝∘‿(2⋆nq-2)⥊𝕩}

The function C performs the classical post-processing part of the algorithm. Given the measurement registers, it reconstructs the candidate period \(r\) and checks the Shor certificate: it is true iff both \(\gcd(a^{\,r/2}\pm 1,N)\) are non-trivial and their product equals \(N\).

Let’s run the quantum simulation and complete the reasoning in the classical domain; let us kill the ignorabimus. Wir müssen wissen, wir werden wissen!4 

C +´Q⟜nq‿sc¨r

1

Compare the result with that obtained using an IBM's Eagle QPU!

Epilogue

This was, of course, a conceptually simple simulation, the minimal instance of Shor’s algorithm that has been verified experimentally. The interpreter itself remains in its basic, non-optimized form. Yet for me, the goal was different: to stress-test BQN on a genuinely technical task. I enjoyed distilling the essence of a long and intricate implementation into a few concise expressions, like trimming a bonsai.

In that process, the language revealed its character. Why the hieroglyphs, you may ask? The tacit and functional style woven from combinators, turns programming into a quiet algebraic game. The resulting program feels less written than discovered, its complexity compressed into a handful of primitives. BQN is the epitome of minimalism's power, and the primitive statistics below tell the same story in numbers:

Primitive's stats

The prog string contains the full source code. We used: 

prog (+´⊸≍⟜≠∊)˜ ⊑¨•primitives

⟨ 43 64 ⟩

With this distribution: 

prog (⍉·⌽˘∘∨⊣≍˘˜≠∘⊣↑/⁼∘⊐)˜"•↩←⇐·𝕊𝕏𝕎𝔽𝔾𝕤𝕩𝕨𝕣𝕗𝕘"∾⊑¨•primitives

┌─                                                                                                                                                                                                                                                                                                                                 
╵ '𝕩' '𝕨' '←' '⊑' '˜' '¨' '´' '⟜' '⊸' '∘' '⌽' '∾' '<' '+' '⊢' '↕' '-' '⥊' '⌾' '⋆' '÷' '×' '𝔽' '⋈' '≠' '=' '𝕊' '⌜' '⇐' '•' '˝' '>' '𝕎' '○' '⊣' '↓' '`' '⎉' '⍟' '⍒' '⌊' '⊏' '≢' '≡' '∨' '√' '∊' '⁼' '»' '¬' '/' '𝕤' '𝕣' '𝕘' '𝕗' '𝕏' '𝔾' '⚇' '◶' '⎊' '⍷' '⍋' '⍉' '⌈' '⊘' '⊔' '⊒' '⊐' '≥' '≤' '≍' '∧' '↩' '↑' '˙' '˘' '·' '«' '|' '!'  
  17  15  15  11  11  10  9   8   8   8   7   7   7   7   6   6   6   5   5   5   5   5   4   4   4   4   3   3   3   3   3   3   2   2   2   2   2   1   1   1   1   1   1   1   1   1   1   1   1   1   1   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0    
                                                                                                                                                                                                                                                                                                                                  ┘

BQN is also fast:

Benchmarks

Here is a comparison with the reference, not particularly optimized Common Lisp code running also 1024 samples: 

Benchmark 1: cbqn -f ../bqn/q.bqn
  Time (mean ± σ):     806.9 ms ± 118.0 ms    [User: 800.6 ms, System: 2.9 ms]
  Range (min … max):   752.2 ms … 1017.9 ms    5 runs
 
Benchmark 2: sbcl --script ../../supp/perf_qi/q.lisp
  Time (mean ± σ):     37.015 s ±  0.184 s    [User: 37.271 s, System: 0.223 s]
  Range (min … max):   36.826 s … 37.286 s    5 runs
 
Summary
  cbqn -f ../bqn/q.bqn ran
   45.87 ± 6.71 times faster than sbcl --script ../../supp/perf_qi/q.lisp

And here is a full program's profile. All time is spent in the Kronecker and matrix products: 

)profile C +´Q⟜nq‿sc¨r

Got 5598 samples
GC: 12 samples
(REPL): 5567 samples:
     1│Q ← {𝕨𝕊nq‿qc:
     1│  wf ← 1⌾⊑⊸⋈0⥊˜2⋆qb
      │  M ← +˝∘×⎉1‿∞ _cp
  2977│  K ← {1𝕊𝕩:𝕩; 𝕨𝕊1:𝕨; 𝕨∾∘×⟜<_cp𝕩}
    86│  L ← {𝕩K⍟𝕨1}{K´⟨g.id𝔽˜nq-𝕨+⌊2⋆⁼≠⊑𝕩, 𝕩, 𝕨𝔽g.id⟩}
    26│  T ← ∾↕∘≠{(⊢∾1↓⌽)𝕨↓↕𝕨⊑⟜𝕩˜•_while_>𝕨⊑𝕩}¨<
   284│  A ← {1=≠⊑𝕨 ? 𝕩M˜⊑⊸L´𝕨;
  2186│    𝕩M´⟨0⊸L´𝕨⟩(⊢∾∾⟜⌽)L⟜g.swap¨T𝕨(⌽∘⊢∾¬∘∊/⊣)⟜⊑˜↕nq
      │  }
     6│  »⊸<∨`0>𝕨-`+○(ט)´wf A´ ⌽qc
      │}
(REPL): 19 samples:
      │Sin‿Cos‿GCD ← •math
      │U ← •rand.Range
    19│_cp ← {(-´𝔽¨)⋈(+´𝔽¨)⟜⌽}
      │  
      │g ← {
      │  IM ← (⊢⋈≢⥊⟜0⊢)¨
      │  x‿id‿h ⇐ (⌽‿⊢{𝕎𝕩}¨<=⌜˜↕2) ∾○IM <-⌾(1‿1⊸⊑)2‿2⥊÷√2
      │  swap‿cnot ⇐ IM ⟨1‿2, 2‿3⟩ {⌽⌾(𝕨⊸⊏)𝕩}¨ <=⌜˜↕4
      │  P ⇐ {⟨=⌜˜↕4,  4‿4⥊0⟩ {𝕨⌾(3‿3⊸⊑)𝕩}¨˜ Sin⊸⋈⟜Cos𝕩}
      │}

Try running the simulation in the BQN repl and explore it! If you are an Emacs user, the org-mode computational notebook in the blog's repository provides the best experience.

One last thing: if you are particularly pedantic like me, you have by now realized that the Shor's algorithm test above is a surprisingly weak unit test. Because it relies on extracting global periodic structures, it can easily mask severe, localized topological routing bugs. To formally verify our tensor geometry, we must test unitary reversibility:

Uncomputation test 
UC ← {𝕊nq‿depth:
  fs ← ⟨g.h⋈˜U ⋄ g.cnot⋈˜2⊸•rand.Deal ⋄ 2⊸•rand.Deal⋈·g.P 2×π×U∘0⟩⊏˜U⟜3 depth
  (⍉¨⋈⟜-´)⌾(1⊸⊑)¨∘⌽⊸∾{𝕏nq}¨fs
}

By generating a chaotic sequence of parameterized gates and appending its exact conjugate-transpose inverse in reverse chronological order, the total circuit evaluates to the identity matrix. If our simulator's geometry is mathematically flawless, applying this circuit will perfectly uncompute the system and restore the initial state vector unmodified: 

{(1⌾⊑0⥊˜2⋆𝕩) ≡ 0 U⊸Q𝕩⋈UC𝕩‿20} 10

1

Footnotes:

1

Although conceptually straightforward, a Hilbert space of size 2⋆n makes this type of simulation a true computational challenge. For an efficient implementation, see arXiv:1710.05867.

2

This post's title is a playful recursive acronym that employs quantum computing terminology, without any specific significance beyond that.

3

I golfed this down a fair bit, though true Kolmogorov complexity isn't the goal of course. I'm not fond of that •_while_ loop, but my only alternative is the under {𝕩𝕊⍟≢⊏⟜𝕩⌾((/𝕩<↕≠𝕩)⊸⊏)𝕩}. Edit: duh {𝕩𝕊⍟≢⊒˜⊸⌊⊸⊏𝕩}.

4

Hilbert's radio address in 1930.