QSO User Guide

QSO is a quantum simulator written in pure Rust, with a Python package and a command-line tool. It simulates quantum circuits exactly, as pure states or as density matrices with noise, and runs variational and textbook algorithms on top of them. It also simulates quantum walks on arbitrary graphs and samples photons in linear interferometers. Its signature geometry is Metatron’s Cube, the complete graph K₁₃ on the 13 centres of the Fruit of Life.

The rule of the project is no claim without a test. Every algorithm is compared with an independent reference: a closed form, a published value, a dense-matrix computation, a second method or another simulator. The chapter How QSO is validated describes the checks, and the technical note reports them in full.

What this guide covers

The tutorials are Jupyter notebooks from docs/tutorials. They are executed every time this site is built, so their output always comes from the current code.

Examples use the Python package, with Rust equivalents where the interfaces differ. The two share one implementation. Every code block in this guide is run when the site is built.

Versions

This guide describes the code on the main branch. Released versions are on PyPI and crates.io; the changelog lists what changed in each release.

Installation

Python

pip install metatron-qso

One abi3 wheel covers CPython 3.9 and newer. Prebuilt wheels exist for these platforms:

• Linux x86_64 and aarch64, glibc and musl;
• macOS on Intel and Apple silicon;
• Windows x86_64 and arm64.

The only runtime dependency is NumPy. On other platforms pip builds the package from source, which needs Rust 1.89 or newer.

import metatron_qso as q

print(q.__version__, q.MAX_QUBITS)

Rust

cargo add metatron-qso-rs

The crate is imported as metatron_qso. It has two features:

Without parallel, everything runs on one thread. The results are bit-for-bit identical either way: parallel sums use a fixed chunking that does not depend on the number of threads.

Command-line tool

cargo install metatron-qso-rs --features cli
qso graph metatron

From source

No C or C++ toolchain is needed, only Rust, plus maturin for the Python package:

git clone https://github.com/LashSesh/qso && cd qso
cargo test -p metatron-qso-rs
pip install maturin
maturin develop --release -m metatron_qso_py/Cargo.toml

Conventions

Qubit order

Qubit 0 is the least significant bit of a basis index, as in Qiskit. The amplitude at index k belongs to the basis state with bit q of k on qubit q. Bit strings are printed with qubit n−1 first, so "001" means qubit 0 is 1.

import metatron_qso as q

state = q.simulate(q.Circuit(3).x(0))
assert abs(state.amplitudes[1]) == 1.0      # index 1 = 0b001
assert q.bitstring(1, 3) == "001"

Cirq and PennyLane order qubits the other way. The cross-framework comparison reverses the qubits when it compares their states with QSO’s.

Gates and time evolution

Randomness

All randomness goes through one portable, seeded generator (xoshiro256++). Functions that sample or restart take a seed. The same seed gives the same result on every platform and for every thread count.

Errors

Invalid input raises ValueError in Python and returns QsoError in Rust, never a panic. This covers, among others, out-of-range qubits, wrong parameter counts, non-finite values, unsupported OpenQASM and oversized requests. The fuzz tests check this for the parsers.

import metatron_qso as q

try:
    q.Circuit(2).cx(0, 2)
except ValueError as e:
    print(e)

Limits

Circuits and states

Building circuits

A Circuit holds a fixed number of qubits and a list of gates. Builder methods return the circuit, so calls chain:

import metatron_qso as q

ghz = q.Circuit(3).h(0).cx(0, 1).cx(1, 2)
print(len(ghz), ghz.depth, ghz.gate_counts())

The gate set:

• Single-qubit gates: H, X, Y, Z, S, S†, T, T†, √X, √X†.
• Rotations: RX, RY, RZ, P and U.
• Two-qubit gates: SWAP, RXX, RYY and RZZ.
• Controlled gates: CX, CY, CZ, CH, CRX, CRY, CRZ, CP, CCX, CSWAP, and multi-controlled X and Z.
• Custom unitaries on up to six qubits, with any number of controls.

import numpy as np
import metatron_qso as q

c = q.Circuit(3)
c.ccx(0, 1, 2).mcz([0, 1, 2]).rzz(0, 2, 0.4)
iswap = np.array([[1, 0, 0, 0], [0, 0, 1j, 0], [0, 1j, 0, 0], [0, 0, 0, 1]])
c.unitary(iswap, targets=[1, 2], controls=[0])

Simulating

simulate starts from |0…0⟩ and returns a StateVector:

import numpy as np
import metatron_qso as q

state = q.simulate(q.Circuit(3).h(0).cx(0, 1).cx(1, 2))
assert np.allclose(state.probabilities(), [0.5, 0, 0, 0, 0, 0, 0, 0.5])
print(state.sample(1000, seed=1))              # {'000': ..., '111': ...}
print(state.entanglement_entropy([0]))         # 1.0 bit

Other starting states are available as StateVector.zero, plus, basis, random and from_amplitudes. state.apply(circuit) continues from any state. States support inner products, fidelities and reduced density matrices.

Observables

Observables are sums of Pauli strings with real coefficients. The parser accepts expressions like the ones in papers:

import metatron_qso as q

h = q.PauliSum("0.5*Z0 Z1 - 1.2 X0 + 0.3*Y2*Y3 - 1.05")
print(h.num_terms, h.constant)
state = q.simulate(q.Circuit(4).h(0))
print(h.expectation(state))
energy, ground = h.ground_state()            # Lanczos
print(energy)

Sums can be added, subtracted and scaled. The models transverse_field_ising, heisenberg_xxz, maxcut_cost and h2_sto3g build common Hamiltonians. Expectation values are computed without building a matrix.

OpenQASM 2.0

Circuits export to OpenQASM 2.0, using only gates of the original qelib1.inc. Gates outside that set are decomposed exactly. Import reads the same gates plus Qiskit’s extensions:

import metatron_qso as q

c = q.Circuit(2).h(0).cx(0, 1).rz(1, 0.25)
text = c.to_qasm(measure=True)
again = q.Circuit.from_qasm(text)
assert abs(q.simulate(c).fidelity(q.simulate(again)) - 1) < 1e-12

Import rejects, with an error, the following:

• custom gate definitions, reset and classical conditions;
• parameters that are not finite;
• programs larger than 2¹⁴ qubits or 2²² gates.

In Rust

The Rust interface mirrors the Python one:

use metatron_qso::prelude::*;

fn main() -> Result<()> {
    let mut ghz = Circuit::new(3);
    ghz.h(0).cx(0, 1).cx(1, 2);
    let state = StateVector::from_circuit(&ghz, &[])?;
    let counts = state.sample(1000, &mut Rng::seed_from_u64(1));
    println!("{counts:?}");
    let zz: PauliSum = "Z0 Z2".parse()?;
    assert!((zz.expectation(&state)? - 1.0).abs() < 1e-12);
    println!("{}", to_qasm2(&ghz, &[], true)?);
    Ok(())
}

Parameters and gradients

Symbolic parameters

new_parameter returns a parameter p[k] of the circuit. A gate angle can be any affine expression a·p[k] + b. The parameter values theta are passed when the circuit is simulated, so one circuit serves every evaluation:

import metatron_qso as q

c = q.Circuit(2)
a, b = c.new_parameters(2)
c.ry(0, a).cx(0, 1).rz(1, 2.0 * b + 0.5)
print(c.num_parameters)                       # 2
state = q.simulate(c, [0.3, 0.1])
fixed = c.bind([0.3, 0.1])                    # the same circuit with numbers

Gradients

expectation_and_gradient returns ⟨ψ(θ)|H|ψ(θ)⟩ and its gradient with respect to all parameters:

import numpy as np
import metatron_qso as q

c = q.Circuit(2)
a, b = c.new_parameters(2)
c.ry(0, a).cx(0, 1).rz(1, 2.0 * b + 0.5).rx(0, b)
h = q.PauliSum("0.5*Z0 Z1 - X0 + 0.2 Y1")
theta = [0.3, 0.1]

value, adjoint = q.expectation_and_gradient(c, h, theta)                    # default
_, shift = q.expectation_and_gradient(c, h, theta, method="parameter_shift")
_, fd = q.expectation_and_gradient(c, h, theta, method="finite_difference")
assert np.allclose(adjoint, shift, atol=1e-12)
assert np.allclose(adjoint, fd, atol=1e-6)

With 100 parameters on 10 qubits, the adjoint method is about 50 times faster than parameter shift (BENCHMARKS.md). Shared parameters and affine angles are handled by the chain rule: in the example, b enters two gates with factors 2 and 1.

Optimisers

run_vqe and MaxCutQaoa.optimize accept method= with one of:

• lbfgs (default; L-BFGS with a strong-Wolfe line search);
• adam;
• gradient_descent;
• nelder_mead;
• spsa.

In Rust, optimize::minimize works with any objective that implements the Objective trait:

use metatron_qso::optimize::{Method, Objective, Options, minimize};
use metatron_qso::prelude::*;

/// Rosenbrock's function with its gradient.
struct Rosenbrock;

impl Objective for Rosenbrock {
    fn dim(&self) -> usize {
        2
    }
    fn value(&self, x: &[f64]) -> f64 {
        (1.0 - x[0]).powi(2) + 100.0 * (x[1] - x[0] * x[0]).powi(2)
    }
    fn value_and_gradient(&self, x: &[f64]) -> Option<(f64, Vec<f64>)> {
        let g = vec![
            -2.0 * (1.0 - x[0]) - 400.0 * x[0] * (x[1] - x[0] * x[0]),
            200.0 * (x[1] - x[0] * x[0]),
        ];
        Some((self.value(x), g))
    }
}

fn main() -> Result<()> {
    let r = minimize(&Rosenbrock, &[-1.2, 1.0], &Method::lbfgs(), &Options::default())?;
    assert!((r.x[0] - 1.0).abs() < 1e-6 && (r.x[1] - 1.0).abs() < 1e-6);
    println!("{} iterations", r.iterations);
    Ok(())
}

Variational algorithms

VQE

run_vqe minimises ⟨ψ(θ)|H|ψ(θ)⟩ over the parameters of an ansatz. Each run starts from several seeded random points, which are evaluated in parallel, and uses exact gradients:

import metatron_qso as q

h2 = q.h2_sto3g()                                   # 2 qubits, parity mapping
result = q.run_vqe(h2["hamiltonian"], q.real_amplitudes(2, 1), restarts=4)
total = result["energy"] + h2["nuclear_repulsion"]
assert abs(total - (-1.137306)) < 1e-6              # full CI of the basis

Three ansatz families are available:

• real_amplitudes: RY layers with CX entanglers.
• efficient_su2: RY and RZ layers with CX entanglers.
• hardware_efficient: choose the rotations and the entanglement pattern (linear, circular or full).

Any parametrised Circuit works as an ansatz. For larger models, compare with the exact ground energy from Lanczos:

import metatron_qso as q

h = q.transverse_field_ising(6, 1.0, 0.8)
exact = h.ground_state_energy()
vqe = q.run_vqe(h, q.efficient_su2(6, 4), restarts=4, seed=1)
print(f"VQE {vqe['energy']:.6f}, exact {exact:.6f}")

QAOA for MaxCut

MaxCutQaoa simulates QAOA for weighted MaxCut in a diagonal basis, which is much faster than a generic circuit. It also finds the exact optimum by enumeration, so every result comes with its approximation ratio:

import metatron_qso as q

ring = q.MaxCutQaoa(q.Graph.cycle(10), p=1).optimize(restarts=8, seed=0)
assert abs(ring["approximation_ratio"] - 0.75) < 1e-9   # 3/4 on rings at p = 1
print(ring["most_likely_bitstring"])

for r in q.optimize_qaoa_layerwise(q.Graph.metatron_cube(), 3):
    print(r["p"], round(r["approximation_ratio"], 5))

optimize_qaoa_layerwise seeds each depth with the INTERP extrapolation of the previous optimum (Zhou et al. 2020), so the ratio improves monotonically. The validated reference values are listed in the technical note:

Classifiers

VariationalClassifier is a data re-uploading classifier for points in the plane. GraphClassifier classifies graphs. Each graph defines its own cut Hamiltonian, the model applies the same trainable QAOA layers to every graph, and the normalised expected cut is the feature:

import metatron_qso as q

def cubic(n, count, seed, want_triangle):
    graphs, s = [], seed
    while len(graphs) < count:
        g = q.Graph.random_regular(n, 3, seed=s)
        s += 1
        if (g.triangle_count() > 0) == want_triangle:
            graphs.append(g)
    return graphs

train = cubic(12, 10, 1, True) + cubic(12, 10, 500, False)
labels = [1] * 10 + [0] * 10
model = q.GraphClassifier(1)
report = model.fit(train, labels, restarts=3, seed=0)
print(report["accuracy"])          # one layer sees triangles: 1.0

After p layers, the model sees each edge’s neighbourhood up to distance p, which is the QAOA light cone. Depth therefore decides which structures the model can detect: one layer sees triangles, two layers see 4- and 5-cycles. Section 3.6 of the technical note quantifies this.

In Rust

use metatron_qso::prelude::*;

fn main() -> Result<()> {
    let h2 = models::h2_sto3g();
    let ansatz = ansatz::real_amplitudes(2, 1, Entanglement::Linear);
    let vqe = Vqe::new(h2.hamiltonian.clone(), ansatz)?.with_restarts(4).run()?;
    assert!((vqe.energy + h2.nuclear_repulsion + 1.137306).abs() < 1e-6);

    let qaoa = MaxCutQaoa::new(&generators::petersen(), 2)?;
    let r = qaoa.optimize(&QaoaConfig::default())?;
    println!("Petersen, p = 2: ratio {:.4}", r.approximation_ratio.unwrap());
    Ok(())
}

Noise

Channels

A Channel is a completely positive, trace-preserving map given by Kraus operators. QSO provides the standard channels:

Qiskit parametrises depolarising noise differently. Its λ is λ = 4p/3 for one qubit and λ = 16p/15 for two qubits.

Noise models and density matrices

A NoiseModel applies channels after every gate:

• single_qubit channels follow each one-qubit gate;
• multi_qubit channels follow each multi-qubit gate, a one-qubit channel on each of its qubits, a two-qubit channel on the pair.

Optional readout errors flip measured bits. DensityMatrix.simulate runs a circuit with such a model:

import metatron_qso as q

bell = q.Circuit(2).h(0).cx(0, 1)
model = q.NoiseModel(
    single_qubit=[q.Channel.depolarizing(0.01), q.Channel.amplitude_damping(0.02)],
    multi_qubit=[q.Channel.depolarizing_2q(0.03)],
    readout=(0.02, 0.05),
)
rho = q.DensityMatrix.simulate(bell, noise=model)
print(rho.purity(), rho.fidelity_with_pure(q.simulate(bell)))
print(rho.sample(1000, seed=1, readout=model.readout))

Density matrices use the same kernels as pure states, applied to the vectorised ρ, and support up to 14 qubits. The weekly cross-validation compares them with Qiskit Aer’s density-matrix method on random noisy circuits. They agree to 3.5·10⁻¹⁵.

Shot-based estimates

estimate_expectation measures an observable the way hardware would. It groups qubit-wise commuting terms, measures each group with the given number of shots, and can add noise:

import metatron_qso as q

bell = q.Circuit(2).h(0).cx(0, 1)
h = q.PauliSum("Z0 Z1 + X0 X1")
exact = h.expectation(q.simulate(bell))                         # 2.0
noisy = q.estimate_expectation(bell, h, shots=20_000, seed=3,
                               noise=q.NoiseModel.depolarizing(0.01, 0.02))
print(exact, noisy)

T₁ and T₂

thermal_relaxation composes amplitude damping and dephasing. The test suite checks the decay against the closed forms ⟨1|ρ(t)|1⟩ = e^{−t/T₁} and |⟨0|ρ(t)|1⟩| = ½ e^{−t/T₂}:

import math
import metatron_qso as q

t1, t2, t = 50.0, 30.0, 20.0
rho = q.DensityMatrix.simulate(q.Circuit(1).h(0))
rho = rho.apply_channel(q.Channel.thermal_relaxation(t1, t2, t), [0])
m = rho.to_matrix()
assert abs(abs(m[0, 1]) - 0.5 * math.exp(-t / t2)) < 1e-12

Graphs

Building graphs

Graph is an undirected graph with real, non-zero edge weights, no loops and no multi-edges. Graphs can be built edge by edge, from an edge list, by name or from a short specification:

import metatron_qso as q

g = q.Graph.from_edges(4, [(0, 1), (1, 2), (2, 3, 0.5)])
g.add_edge(3, 0)
print(g.summary())

for spec in ["metatron", "metatron:1", "fruit-of-life", "petersen", "heawood",
             "hypercube:4", "torus:4x5", "platonic:dodecahedron",
             "random-regular:12,3,7", "erdos-renyi:20,0.2,1"]:
    h = q.Graph.from_spec(spec)
    print(f"{spec:24} {h.num_nodes:3} nodes {h.num_edges:3} edges")

Specifications and the Python constructors reject graphs with more than 2²⁰ nodes or 2²² edges before building anything. Hypercubes up to dimension 20 are the exception.

Structure and spectra

Automorphism groups are computed exactly, by colour refinement with individualisation and orbit pruning. The order is an exact integer even for 13! = 6 227 020 800:

import metatron_qso as q

assert q.Graph.metatron_cube().automorphism_group_order() == 6227020800
assert q.Graph.fruit_of_life().automorphism_group_order() == 12      # D6
assert q.Graph.platonic("dodecahedron").automorphism_group_order() == 120
assert q.Graph.heawood().girth() == 6

Metatron’s Cube

The Fruit of Life consists of 13 equal circles: one in the centre, six touching it, and six more along the same rays. Joining every pair of centres gives the 78 lines of Metatron’s Cube.

• As a graph it is K₁₃. The full symmetric group S₁₃ acts on it, and the Laplacian spectrum is {0, 13¹²}.
• The geometry is in the embedding. The lines have seven distinct lengths: 1, √3, 2, √7, 3, 2√3 and 4.
• Short lines recover the structure. metatron:<length> keeps only the lines up to a length. The contact graph (metatron:1, tangent circles) has exactly the hexagonal symmetry D₆.

import metatron_qso as q

for length, count in q.metatron_line_lengths():
    print(f"{length:.4f}  {count}")

QSO makes no other claim about Metatron’s Cube. What can be said is graph theory or physics, and every such statement is tested.

Quantum walks

Continuous-time walks

A continuous-time quantum walk evolves a state on the nodes of a graph with exp(−iHt), where H = γL (Laplacian, the default) or H = −γA (adjacency). Up to 2048 nodes QSO diagonalises H exactly, which gives probabilities at any time, limiting and time-averaged distributions, and mixing times. Beyond that size, evolve and probabilities propagate with a Krylov method.

import metatron_qso as q

walk = q.ContinuousWalk(q.Graph.metatron_cube())
print(walk.probabilities(0, 1.3))
limit = walk.limiting_distribution(0)
assert abs(limit[0] - 145 / 169) < 1e-12

The time-averaged distribution converges to the limiting distribution. The mixing time M_ε (Aharonov et al. 2001) is the time from which on the average stays within total-variation distance ε of the limit:

import metatron_qso as q

walk = q.ContinuousWalk(q.Graph.metatron_cube())
print(walk.mixing_time(0, 0.02, 20.0, steps=2000))

Spatial search

Childs–Goldstone search adds an oracle term −|w⟩⟨w| for the marked nodes w. At the critical rate γ the walk moves amplitude to the marked node. SpatialSearch.optimize tunes γ and the measurement time together:

import math
import metatron_qso as q

cube = q.Graph.metatron_cube()
search = q.SpatialSearch(cube, marked=[4], gamma=1 / 13)
assert abs(search.success_probability(math.pi * math.sqrt(13) / 2) - 1) < 1e-9

best = q.SpatialSearch.optimize(q.Graph.fruit_of_life(), [0])
print(best)            # gamma, time, success probability, expected time
print(q.classical_hitting_times(cube, 4)[0])        # 12 steps for a random walk

On complete graphs the search time grows as √N, the quadratic speed-up over classical search. The technical note checks this for N up to 1024.

Decoherence

DephasingWalk adds Lindblad dephasing at a given rate κ. With κ = 0 it reproduces the coherent walk. For strong dephasing the populations follow the classical master equation dp/dt = −(2γ²/κ)·L·p. This is the quantum Zeno regime: the walk becomes classical and slows down as κ grows.

import numpy as np
import metatron_qso as q

g = q.Graph.fruit_of_life()
for rate in [0.0, 1.0, 10.0]:
    p = q.DephasingWalk(g, rate=rate).populations(0, 3.0)
    print(rate, np.round(p[:3], 4))
print(np.round(q.classical_walk_probabilities(g, 0, 3.0)[:3], 4))

In Rust

use metatron_qso::prelude::*;

fn main() -> Result<()> {
    let cube = metatron::metatron_cube();
    let walk = ContinuousWalk::new(&cube, Generator::Laplacian, 1.0)?;
    let limit = walk.limiting_distribution(0)?;
    assert!((limit[0] - 145.0 / 169.0).abs() < 1e-12);

    let search = SpatialSearch::new(&cube, &[3], Generator::Laplacian, 1.0 / 13.0)?;
    let t = std::f64::consts::PI * 13f64.sqrt() / 2.0;
    assert!((search.success_probability(t) - 1.0).abs() < 1e-9);
    Ok(())
}

Boson sampling

Interferometers and permanents

An m-mode linear interferometer is an m × m unitary U. When n photons enter with occupation S and leave with occupation T, the output probability is a matrix permanent (Aaronson & Arkhipov 2011):

• for indistinguishable photons, P(T) = |perm(U_{T,S})|² / (Π s! Π t!);
• for distinguishable photons, P(T) = perm(|U_{T,S}|²) / Π t!.

QSO computes permanents with Glynn’s formula in Gray-code order, which costs O(n·2ⁿ) operations. The result is the same for any number of threads. For n = 24 it takes about 0.14 s on four cores, and it agrees with The Walrus to 3·10⁻¹⁴.

import numpy as np
import metatron_qso as q

assert q.permanent(np.ones((6, 6))) == 720           # perm(J) = n!
u = q.haar_unitary(8, seed=1)                        # Haar-random (Mezzadri)
assert np.allclose(u.conj().T @ u, np.eye(8), atol=1e-12)

Exact probabilities and sampling

import math
import numpy as np
import metatron_qso as q

# Hong–Ou–Mandel: two photons on a 50:50 beam splitter leave together.
hom = q.BosonSampler(q.beam_splitter(math.pi / 4), [1, 1])
assert hom.probability([1, 1]) < 1e-30
assert abs(hom.probability_distinguishable([1, 1]) - 0.5) < 1e-15

# Four photons in eight modes: full distribution and exact samples.
s = q.BosonSampler(q.haar_unitary(8, seed=1), [1, 1, 1, 1, 0, 0, 0, 0])
patterns, probs = s.distribution()
assert abs(probs.sum() - 1) < 1e-12
samples = s.sample(10_000, seed=2)           # Clifford & Clifford, exact

# Full bunching is enhanced by n! over distinguishable photons.
bunched = [0, 0, 4, 0, 0, 0, 0, 0]
ratio = s.probability(bunched) / s.probability_distinguishable(bunched)
assert abs(ratio - 24) < 1e-9

Exact sampling uses the algorithm of Clifford & Clifford (2018) and needs at most one photon per input mode. The test suite checks the samples against the exact distribution with a χ² test.

Limits

Command line

qso <command> [options]. Every command accepts --json for machine-readable output; qso <command> --help lists all options.

Graphs are given as specifications (see Graphs), for example metatron, fruit-of-life, petersen, hypercube:8 or random-regular:16,3,1.

$ qso graph fruit-of-life
Fruit of Life (contact graph)
  nodes 13, edges 18, degree 1..6, connected true, bipartite false
  diameter Some(4), triangles 6, algebraic connectivity 0.585786
  |Aut| = 12 (log10 1.079), orbits 3

$ qso qaoa --graph petersen --p 4 --layerwise
p=1  <C>=10.386751  C_max=12  ratio=0.865563  P(optimal)=0.1682  most likely 1110001011 (cut 12)
p=2  <C>=11.105320  C_max=12  ratio=0.925443  P(optimal)=0.4490  most likely 1110001011 (cut 12)
p=3  <C>=11.648892  C_max=12  ratio=0.970741  P(optimal)=0.7626  most likely 1100110110 (cut 12)
p=4  <C>=11.936517  C_max=12  ratio=0.994710  P(optimal)=0.9767  most likely 1001101101 (cut 12)

Circuits and exact gradients

This tutorial builds circuits, samples them, and computes gradients of expectation values in three ways. It also shows why the adjoint method is the default: its cost does not grow with the number of parameters.

import time

import matplotlib.pyplot as plt
import numpy as np

import metatron_qso as q

print("QSO", q.__version__)

QSO 1.2.0

A GHZ state

Qubit 0 is the least significant bit, so bit strings print qubit 2 first.

ghz = q.Circuit(3).h(0).cx(0, 1).cx(1, 2)
state = q.simulate(ghz)
print(np.round(state.amplitudes, 4))
counts = state.sample(2000, seed=7)
print(counts)
print("entanglement entropy of qubit 0:", state.entanglement_entropy([0]), "bit")

[0.7071+0.j 0.    +0.j 0.    +0.j 0.    +0.j 0.    +0.j 0.    +0.j
 0.    +0.j 0.7071+0.j]
{'000': 1044, '111': 956}
entanglement entropy of qubit 0: 0.9999999999999999 bit

fig, ax = plt.subplots(figsize=(5, 2.5))
labels = [q.bitstring(k, 3) for k in range(8)]
ax.bar(labels, [counts.get(b, 0) for b in labels])
ax.set_xlabel("outcome")
ax.set_ylabel("counts (2000 shots)")
plt.tight_layout()
plt.show()

An energy landscape and its derivative

One parameter θ, observable H = Z₀Z₁ − 0.5·X₀. The adjoint gradient is exact; finite differences agree to their truncation error.

c = q.Circuit(2)
theta = c.new_parameter()
c.ry(0, theta).cx(0, 1).rx(1, 0.5 * theta + 0.2)
h = q.PauliSum("Z0 Z1 - 0.5 X0")

xs = np.linspace(-np.pi, np.pi, 121)
values, adjoint, fd = [], [], []
for x in xs:
    v, g = q.expectation_and_gradient(c, h, [x])
    _, g_fd = q.expectation_and_gradient(c, h, [x], method="finite_difference")
    values.append(v)
    adjoint.append(g[0])
    fd.append(g_fd[0])
print("largest |adjoint − finite difference|:", np.max(np.abs(np.array(adjoint) - fd)))

fig, ax = plt.subplots(figsize=(6, 3))
ax.plot(xs, values, label="⟨H⟩(θ)")
ax.plot(xs, adjoint, label="d⟨H⟩/dθ (adjoint)")
ax.set_xlabel("θ")
ax.legend()
plt.tight_layout()
plt.show()

largest |adjoint − finite difference|: 2.3853558017705723e-11

Cost of a full gradient

For an ansatz with P parameters, parameter shift needs about 2P circuit evaluations; the adjoint method needs about three, whatever P is.

h = q.transverse_field_ising(8, 1.0, 0.7)

def seconds(f, repeats=5):
    f()
    t = time.perf_counter()
    for _ in range(repeats):
        f()
    return (time.perf_counter() - t) / repeats

layers = [1, 2, 4, 8, 12]
params, t_adj, t_shift = [], [], []
for L in layers:
    ansatz = q.efficient_su2(8, L)
    theta = np.linspace(0, 1, ansatz.num_parameters)
    params.append(ansatz.num_parameters)
    t_adj.append(seconds(lambda: q.expectation_and_gradient(ansatz, h, theta)))
    t_shift.append(seconds(lambda: q.expectation_and_gradient(ansatz, h, theta, method="parameter_shift")))

fig, ax = plt.subplots(figsize=(5, 3))
ax.loglog(params, t_adj, "o-", label="adjoint")
ax.loglog(params, t_shift, "s-", label="parameter shift")
ax.set_xlabel("number of parameters")
ax.set_ylabel("seconds per gradient")
ax.legend()
plt.tight_layout()
plt.show()
print(f"{params[-1]} parameters: parameter shift is {t_shift[-1] / t_adj[-1]:.0f}× slower")

208 parameters: parameter shift is 122× slower

OpenQASM round trip

Export binds the parameters; import gives an equivalent circuit.

text = c.to_qasm([0.4], measure=True)
print(text)
again = q.Circuit.from_qasm(text)
print("fidelity:", q.simulate(c, [0.4]).fidelity(q.simulate(again)))

OPENQASM 2.0;
include "qelib1.inc";
qreg q[2];
creg c[2];
ry(0.4) q[0];
cx q[0],q[1];
rx(0.4) q[1];
measure q -> c;

fidelity: 1.0

VQE and QAOA

Variational algorithms on exact simulations, compared with exact answers: the full-CI energy of H₂, Lanczos ground energies of spin chains, and brute-force maximum cuts.

import matplotlib.pyplot as plt
import numpy as np

import metatron_qso as q

H₂ in the STO-3G basis

Two qubits after the parity mapping with two-qubit reduction. VQE with a one-layer RealAmplitudes ansatz reaches the full-CI energy.

h2 = q.h2_sto3g()
result = q.run_vqe(h2["hamiltonian"], q.real_amplitudes(2, 1), restarts=4, seed=0)
energy = result["energy"] + h2["nuclear_repulsion"]
exact = h2["hamiltonian"].ground_state_energy() + h2["nuclear_repulsion"]
print(f"VQE   {energy:.10f} Ha")
print(f"exact {exact:.10f} Ha (full CI of the basis: -1.137306)")
print("iterations:", result["iterations"], "energies of the restarts:", np.round(result["restart_energies"], 8))

VQE   -1.1373060358 Ha
exact -1.1373060358 Ha (full CI of the basis: -1.137306)
iterations: 12 energies of the restarts: [-1.85727503 -1.85727503 -1.85727503 -1.85727503]

plt.figure(figsize=(5, 3))
plt.semilogy(np.array(result["history"]) - h2["hamiltonian"].ground_state_energy() + 1e-16)
plt.xlabel("iteration")
plt.ylabel("energy error (Ha)")
plt.tight_layout()
plt.show()

Depth matters: a transverse-field Ising chain

Six spins at the critical point J = h = 1. The error of VQE falls with the number of ansatz layers.

h = q.transverse_field_ising(6, 1.0, 1.0)
exact = h.ground_state_energy()
layers = [1, 2, 3, 4, 5, 6]
errors = []
for L in layers:
    r = q.run_vqe(h, q.efficient_su2(6, L), restarts=4, seed=1)
    errors.append(r["energy"] - exact)
    print(f"{L} layers: E = {r['energy']:.8f}, error {errors[-1]:.2e}")

plt.figure(figsize=(5, 3))
plt.semilogy(layers, np.maximum(errors, 1e-12), "o-")
plt.xlabel("layers")
plt.ylabel("E_VQE − E_0")
plt.tight_layout()
plt.show()

1 layers: E = -7.23420165, error 6.20e-02
2 layers: E = -7.29130728, error 4.92e-03
3 layers: E = -7.29213077, error 4.10e-03
4 layers: E = -7.29593020, error 3.00e-04


5 layers: E = -7.29620836, error 2.14e-05
6 layers: E = -7.29622473, error 5.08e-06

QAOA for MaxCut

On a ring, depth-p QAOA reaches exactly (2p+1)/(2p+2) of the maximum cut. Layer-wise optimisation (INTERP) carries each optimum to the next depth.

graphs = {
    "ring C10": q.Graph.cycle(10),
    "Petersen": q.Graph.petersen(),
    "Metatron's Cube": q.Graph.metatron_cube(),
}
ratios = {}
for name, g in graphs.items():
    ratios[name] = [r["approximation_ratio"] for r in q.optimize_qaoa_layerwise(g, 4, seed=0)]
    print(name, np.round(ratios[name], 5))
print("ring, theory:", [round((2 * p + 1) / (2 * p + 2), 5) for p in range(1, 5)])

plt.figure(figsize=(5, 3))
for name, r in ratios.items():
    plt.plot(range(1, 5), r, "o-", label=name)
plt.xlabel("depth p")
plt.ylabel("approximation ratio")
plt.legend()
plt.tight_layout()
plt.show()

ring C10 [0.75    0.83333 0.875   0.9    ]
Petersen [0.86556 0.9178  0.97074 0.98906]


Metatron's Cube [0.99282 0.9991  0.97005 0.99419]
ring, theory: [0.75, 0.83333, 0.875, 0.9]

The most likely bit string of the optimised state is an optimal cut:

qaoa = q.MaxCutQaoa(q.Graph.petersen(), 4)
best = qaoa.optimize(restarts=8, seed=0)
print(best["most_likely_bitstring"], "cut", best["most_likely_cut"], "of", qaoa.max_cut)
print("probability of an optimal cut:", round(best["optimal_probability"], 4))

1001101101 cut 12.0 of 12.0
probability of an optimal cut: 0.9767

Noise and decoherence

Density-matrix simulation with noise channels, compared with closed forms, and shot-based estimation as on hardware.

import matplotlib.pyplot as plt
import numpy as np

import metatron_qso as q

A Bell state under depolarising noise

Noise after every gate lowers the purity and the fidelity with the ideal state.

bell = q.Circuit(2).h(0).cx(0, 1)
ideal = q.simulate(bell)
ps = np.linspace(0, 0.3, 31)
fid, pur = [], []
for p in ps:
    rho = q.DensityMatrix.simulate(bell, noise=q.NoiseModel.depolarizing(p, p))
    fid.append(rho.fidelity_with_pure(ideal))
    pur.append(rho.purity())

plt.figure(figsize=(5, 3))
plt.plot(ps, fid, label="fidelity with |Φ⁺⟩")
plt.plot(ps, pur, label="purity")
plt.xlabel("depolarising probability p (one- and two-qubit)")
plt.legend()
plt.tight_layout()
plt.show()

T₁ and T₂

Thermal relaxation of |+⟩ over a time t: the excited population decays as e^{−t/T₁}, the coherence as e^{−t/T₂}.

t1, t2 = 50.0, 30.0
times = np.linspace(0, 150, 61)
pop, coh = [], []
for t in times:
    rho = q.DensityMatrix.simulate(q.Circuit(1).h(0))
    rho = rho.apply_channel(q.Channel.thermal_relaxation(t1, t2, t), [0])
    m = rho.to_matrix()
    pop.append(m[1, 1].real)
    coh.append(2 * abs(m[0, 1]))
print("largest deviation from the closed forms:",
      max(np.max(np.abs(np.array(pop) - 0.5 * np.exp(-times / t1))),
          np.max(np.abs(np.array(coh) - np.exp(-times / t2)))))

plt.figure(figsize=(5, 3))
plt.plot(times, pop, label="excited population")
plt.plot(times, 0.5 * np.exp(-times / t1), "k:", label="½ e^{−t/T₁}")
plt.plot(times, coh, label="coherence 2|ρ₀₁|")
plt.plot(times, np.exp(-times / t2), "k--", label="e^{−t/T₂}")
plt.xlabel("t")
plt.legend()
plt.tight_layout()
plt.show()

largest deviation from the closed forms: 2.220446049250313e-16

Shot noise

estimate_expectation measures groups of commuting Pauli terms with a finite number of shots. The error falls as 1/√shots.

c = q.efficient_su2(4, 2).bind(np.linspace(0.1, 1.5, 24))
h = q.heisenberg_xxz(4, 1.0, 0.5)
exact = h.expectation(q.simulate(c))
shots = [100, 300, 1000, 3000, 10000, 30000]
rms = []
for s in shots:
    est = [q.estimate_expectation(c, h, shots=s, seed=k) for k in range(20)]
    rms.append(np.sqrt(np.mean((np.array(est) - exact) ** 2)))

plt.figure(figsize=(5, 3))
plt.loglog(shots, rms, "o-", label="RMS error (20 seeds)")
plt.loglog(shots, rms[0] * np.sqrt(shots[0] / np.array(shots)), "k:", label="∝ 1/√shots")
plt.xlabel("shots per measurement group")
plt.legend()
plt.tight_layout()
plt.show()

Readout errors

Readout flips act on measured bits only; they do not change the state.

ghz = q.Circuit(3).h(0).cx(0, 1).cx(1, 2)
rho = q.DensityMatrix.simulate(ghz)
print("ideal:      ", rho.sample(4000, seed=1))
print("with readout:", rho.sample(4000, seed=1, readout=(0.02, 0.05)))

ideal:       {'000': 2021, '111': 1979}
with readout: {'000': 1837, '001': 41, '010': 38, '011': 89, '100': 43, '101': 110, '110': 82, '111': 1760}

Quantum walks on Metatron’s Cube

Metatron’s Cube is the complete graph K₁₃ on the 13 centres of the Fruit of Life. This tutorial looks at its symmetry, at continuous-time quantum walks on it and on the Fruit-of-Life contact graph, and at spatial search.

import math

import matplotlib.pyplot as plt
import numpy as np

import metatron_qso as q

cube = q.Graph.metatron_cube()
fol = q.Graph.fruit_of_life()

Geometry and symmetry

for length, count in q.metatron_line_lengths():
    print(f"line length {length:.4f}: {count} lines")
print("|Aut(Metatron's Cube)| =", cube.automorphism_group_order(), "= 13!")
print("|Aut(Fruit of Life)|   =", fol.automorphism_group_order(), "(the dihedral group D6)")
print("orbits of the Fruit of Life:", fol.automorphism_orbits())

line length 1.0000: 18 lines
line length 1.7321: 18 lines
line length 2.0000: 15 lines
line length 2.6458: 12 lines
line length 3.0000: 6 lines
line length 3.4641: 6 lines
line length 4.0000: 3 lines
|Aut(Metatron's Cube)| = 6227020800 = 13!
|Aut(Fruit of Life)|   = 12 (the dihedral group D6)
orbits of the Fruit of Life: [[0], [1, 2, 3, 4, 5, 6], [7, 8, 9, 10, 11, 12]]

pos = fol.positions
fig, axes = plt.subplots(1, 2, figsize=(7, 3.5))
for ax, g, title in [(axes[0], cube, "Metatron's Cube (K13)"), (axes[1], fol, "Fruit of Life contact graph")]:
    p = g.positions
    for u, v, _ in g.edges:
        ax.plot([p[u, 0], p[v, 0]], [p[u, 1], p[v, 1]], "k-", lw=0.6)
    ax.scatter(p[:, 0], p[:, 1], s=40, zorder=3)
    ax.set_title(title)
    ax.set_aspect("equal")
    ax.axis("off")
plt.tight_layout()
plt.show()

A walk from the centre

On K₁₃ the Laplacian has eigenvalues 0 and 13 (twelve-fold), so the walk from the centre returns periodically. Its time average converges to 145/169: the limiting distribution keeps a large share at the start.

walk = q.ContinuousWalk(cube)
ts = np.linspace(0, 3, 301)
p_centre = [walk.probabilities(0, t)[0] for t in ts]
avg = np.cumsum(p_centre) / np.arange(1, len(ts) + 1)
print("limiting P(centre) =", walk.limiting_distribution(0)[0], " 145/169 =", 145 / 169)

fol_walk = q.ContinuousWalk(fol)
p_fol = [fol_walk.probabilities(0, t)[0] for t in ts]
p_classical = [q.classical_walk_probabilities(fol, 0, t)[0] for t in ts]

plt.figure(figsize=(6, 3))
plt.plot(ts, p_centre, label="K13, quantum")
plt.plot(ts, avg, label="K13, running time average")
plt.axhline(145 / 169, color="k", ls=":", lw=0.8)
plt.plot(ts, p_fol, label="Fruit of Life, quantum")
plt.plot(ts, p_classical, label="Fruit of Life, classical")
plt.xlabel("t")
plt.ylabel("P(centre)")
plt.legend(fontsize=8)
plt.tight_layout()
plt.show()

limiting P(centre) = 0.8579881656804739  145/169 = 0.8579881656804734

Spatial search

With the oracle −|w⟩⟨w| and γ = 1/N, the walk on the complete graph moves all amplitude to the marked node at t = π√N/2. A classical random walk on K₁₃ needs 12 steps on average.

search = q.SpatialSearch(cube, marked=[4], gamma=1 / 13)
ts = np.linspace(0, 12, 400)
success = [search.success_probability(t) for t in ts]
t_star = math.pi * math.sqrt(13) / 2
print(f"P(t = π√13/2 = {t_star:.4f}) = {search.success_probability(t_star):.12f}")
print("classical hitting time:", q.classical_hitting_times(cube, 4)[0])

best = q.SpatialSearch.optimize(fol, [0])
print("Fruit of Life, marked centre:", best)

plt.figure(figsize=(5, 3))
plt.plot(ts, success)
plt.axvline(t_star, color="k", ls=":", lw=0.8)
plt.xlabel("t")
plt.ylabel("P(marked)")
plt.tight_layout()
plt.show()

P(t = π√13/2 = 5.6636) = 1.000000000000
classical hitting time: 11.999999999999998
Fruit of Life, marked centre: {'gamma': 0.3147337510441124, 'time': 5.46332670498062, 'success_probability': 0.5117699808266946, 'expected_time': 10.675355940485922}

From quantum to classical

Dephasing at rate κ turns the walk into a classical random walk; for large κ the populations follow dp/dt = −(2γ²/κ)·L·p (quantum Zeno regime).

ts = np.linspace(0, 6, 121)
plt.figure(figsize=(6, 3))
for rate in [0.0, 0.5, 2.0, 8.0]:
    walk = q.DephasingWalk(fol, rate=rate)
    plt.plot(ts, [walk.populations(0, t)[0] for t in ts], label=f"κ = {rate}")
plt.plot(ts, [q.classical_walk_probabilities(fol, 0, t)[0] for t in ts], "k:", label="classical")
plt.xlabel("t")
plt.ylabel("P(centre)")
plt.legend(fontsize=8)
plt.tight_layout()
plt.show()

How QSO is validated

QSO makes no claim it cannot test. This page lists the kinds of evidence behind its results and where to find them.

References for every algorithm

Every algorithm is checked against an independent reference:

These checks live in the test suites, which run on every change on Linux (x86_64 and arm64), macOS and Windows:

• Rust: unit, validation, property and robustness tests;
• Python: tests of the bindings.

Cross-validation with other frameworks

benchmarks/cross_framework compares QSO with established simulators. It runs every week, and a deviation above rounding level fails the run.

Robustness

API stability

QSO 1.x adds and never breaks:

• cargo-semver-checks compares the public Rust API with the last release;
• stubtest checks that the Python type stubs match the compiled module.

Performance

The benchmark report states the machine it ran on. A weekly run compares the timings with the median of the previous four runs. A slowdown of more than 25 % is flagged.

The technical note

Validated Quantum Algorithm Benchmarks of the Metatron QSO Simulation Framework (DOI 10.5281/zenodo.23623226) reports all of the above with the data and the programs that produce it.

API reference

Further documents in the repository: