Exact cover, in order of preference

Originally written in Python with NumPy and pandas during my Robotics program; ported to TypeScript here so the solver runs live in the browser. Everything below computes in your own tab.

The problem

An exact cover asks: given a collection of sets, choose a subcollection so that every element is covered by exactly one chosen set. Sudoku, polyomino tiling, and N-Queens all reduce to it. It is NP-complete — no known way to find solutions in polynomial time — yet Knuth's Algorithm X solves it elegantly by depth-first backtracking, repeatedly picking the least-covered element and branching on the sets that cover it.

Algorithm X finds all solutions, but treats them as interchangeable. Many real problems have a notion of a better solution. So I modified the search: instead of backtracking depth-first, Algorithm X* explores the tree of partial solutions best-first — a uniform-cost priority queue ordered by the cost so far, with no heuristic — and yields solutions in nondecreasing cost order. If you only need the best few, you stop early and never pay for the rest.

The application: twelve-tone chord progressions

Partition the twelve notes of the octave into four three-note chords (the twelve major and twelve minor triads), each note used once — a clean exact cover. But most such partitions sound arbitrary. I score every progression by voice-leading distance from Neo-Riemannian theory: two chords are "close" if you can move between them by shifting a single note, and a progression's cost is the sum of pairwise distances among its chords. Under that cost alone six progressions tie at 20; add a small reward for chords in root position and the tie breaks, leaving a single best progression — F, C#m, Gm, B, at cost 23.

X* returns exactly that ordering below. Toggle the cost function, reveal solutions cheapest-first, and click ▶ to hear each one.

The core of X*

The whole algorithm is a priority queue over partial solutions. Pop the cheapest path; reduce the table by the rows it has chosen; if every requirement is met, it's a solution; otherwise pick the sparsest remaining column and enqueue a child for each row that covers it, scored by the cost function.

export function* xStar(t, cost) {
  const heap = new MinHeap();
  // seed from the sparsest column …
  while (heap.size) {
    const { k: curCost, v: path } = heap.pop();
    const { remRows, colCount } = reduce(t, path);
    if (colCount.size === 0) { yield { rows: path, cost: curCost }; continue; }
    const col = chooseColumn(colCount);   // sparsest requirement, or prune
    if (col === null) continue;
    for (const row of remRows)
      if (t.rowCols.get(row).has(col)) enqueue(curCost, path, row);
  }
}

The full solver and the twelve-tone setup are two small TypeScript modules (xstar.ts, tonerow.ts). The solver is generic over any cost function and supports optional columns, so the same code handles generalized exact cover — which sets up the next step.

The same problem, reduced to a neutral-atom quantum optimizer

This section runs on Pasqal's open-source stack (QoolQit / the maximum-independent-set package) against its local emulator — Python, not the browser. Numbers below are from that emulator. The full script stays private while this is active work; the excerpts below show the parts that matter.

Exact cover has a clean quantum-native form. Map each candidate set to a graph node, weight it by how many elements it covers, and connect two nodes whenever their sets overlap. An independent set — no two chosen nodes adjacent — is then a family of pairwise disjoint sets, and a maximum-weight independent set whose weight equals the size of the universe is an exact cover. On a neutral-atom processor this is almost literal: each set becomes an atom, overlapping sets are placed within the Rydberg blockade radius so they can't both be excited, and an adiabatic ramp of the global detuning prepares the ground state of the blockaded Hamiltonian — a maximum independent set. The hardware constraint is the disjointness constraint. That mapping is only literal when the conflict graph embeds as a unit-disk graph in the plane; the solver's embedding step handles the general case, and for a dense graph that embedding, not the model, is the real source of error.

# exact-cover instance -> weighted conflict graph
g = nx.Graph()
for name, s in subsets.items():
    g.add_node(name, weight=len(s))          # weight -> per-atom detuning
for a, b in combinations(subsets, 2):
    if subsets[a] & subsets[b]:
        g.add_edge(a, b)                     # overlap -> blockade edge

# solve on the local neutral-atom emulator
cfg = SolverConfig(backend=BackendConfig(backend=BackendType.QUTIP),
                   weighting=Weighting.WEIGHTED, preprocessor=None, runs=300)
solutions = MISSolver(MISInstance(g), cfg).solve()

On Knuth's textbook instance (universe of 7, six subsets) the emulator places six atoms, takes 300 projective measurements, and returns the unique exact cover {B, D, F} in 100% of shots. On a noiseless emulator with six atoms that clean result is the expected one, not a finding — it confirms the encoding.

Pointing the backend at Pasqal Cloud — REMOTE_EMUFREE in place of the local emulator — submits the identical program to Pasqal's hosted EMU_FREE emulator. There, over 200 shots, the exact cover {B, D, F} comes back in 84% of shots, with the rest landing on the smaller independent sets ({B, E}, {A, D}) the array sampled along the way. That sampled distribution — rather than a clean argmax — is closer to what a device returns, though EMU_FREE is still an emulator, and the max-weight cover dominates it.

# same solver, same graph — only the backend changes
backend = BackendConfig(
    backend=BackendType.REMOTE_EMUFREE,     # Pasqal Cloud's EMU_FREE emulator
    username=os.environ["PASQAL_USERNAME"],
    password=os.environ["PASQAL_PASSWORD"],
    project_id=os.environ["PASQAL_PROJECT_ID"],
)
# for the QPU: backend=BackendType.REMOTE_QPU, device=NamedDevice("FRESNEL")
solutions = MISSolver(MISInstance(g), SolverConfig(backend=backend, runs=200)).solve()

Run the same reduction on the twelve-tone instance above and a nice structure falls out. Each triad covers exactly three notes and 4 × 3 = 12, so any independent set of size four already tiles the whole octave: here maximum independent set and exact cover are the same object. Enumerated classically, all eighteen covers are exactly two major plus two minor triads. Weighting the atoms by tonal centrality — roots near C on the circle of fifths — biases the array toward the most in-key cover, the neutral-atom analogue of X*'s preference ordering. At twenty-four atoms the QuTiP backend is impractical, so I ran that instance on the emu-sv state-vector emulator instead — and it turned out to be the most instructive run of the project (below).

What the real hardware demanded

Getting the same program to validate against the FRESNEL device specification — not just an emulator's — took two changes, each a small piece of applied physics. First, weighted mode encodes node weights as per-atom detunings via a DMM (detuning-map modulator) channel, which FRESNEL's spec doesn't expose; since this instance's maximum independent set is already its exact cover, running unweighted loses nothing. Second, the solver's default pulse can violate the channel's minimum average amplitude — the atom layout is randomized, so some embeddings produce a too-weak pulse and the sequence is rejected at validation. Keeping the default's bump shape but clamping its peak up to a floor above the channel minimum — so the average clears it on every layout — passes every time (15/15, versus ~11/15 for the stock pulse). Emulating that FRESNEL-valid sequence locally returns {B, D, F} in 93% of shots. I have not yet run it on the QPU itself; the gap between "validates against the spec" and "runs on the machine" is exactly where the interesting engineering lives.

The two decisions, as an excerpt from the shaper:

class HardwarePulseShaper(DefaultPulseShaper):
    def pulse(self, config, register, backend, instance):
        p  = self._calculate_parameters(register=register, backend=backend, instance=instance)
        ch = backend.device().channels["rydberg_global"]
        # keep the default's bump shape, but clamp the peak so the pulse's *average*
        # clears the channel's minimum-average-amplitude floor
        omega = min(max(p.maximum_amplitude, 3 * ch.min_avg_amp), 0.99999 * ch.max_amp)
        amplitude = InterpolatedWaveform(p.duration_ns, [1e-9, omega, 1e-9])
        detuning  = InterpolatedWaveform(p.duration_ns, [-p.final_detuning, 0, p.final_detuning])
        return Pulse(amplitude, detuning, 0)

    def detuning(self, *args, **kwargs):
        return []   # global-only: no per-atom detuning, so no DMM channel required

What the twenty-four-atom run taught me

Running the twelve-tone instance on emu-sv (500 shots, no kernelization, the FRESNEL-ready pulse) takes 33 minutes on a laptop and does find exact covers: five of the ten returned solutions are covers — every one two major plus two minor, as the theory says — for 26.6% of shots. The stock pulse gets 15.6%. But the single most likely outcome (69% of shots) is a size-three set like Am G#m Gm: maximal but not maximum, because its three leftover notes form no triad. That is the classic independent-set trap, and it was a signal that something upstream was wrong.

The embedding was. Across five random layouts, the weakest true conflict had an interaction of essentially zero, and 131–140 of the 144 conflict edges sat below the Rabi peak — un-blockaded. The register was not encoding the constraint graph. Kernelization is a no-op on this graph, and the optimized embedder only gets to 110–117. So where did the valid solutions come from? The solver's default postprocessor strips every raw shot to an independent set and greedily augments it. The decisive test is to feed that same cleanup no quantum input at all:

On this instance the quantum prior is worse than random through the same classical cleanup. That is not an emulator bug and not a pulse problem; it is the reduction being applied to a graph the hardware cannot represent. Analog neutral-atom MIS natively solves unit-disk graphs — two atoms conflict when they are within the blockade radius, in a plane. The twelve-tone conflict graph is twelve overlapping six-cliques at density 0.52, about the worst shape for a plane. Knuth's six-node instance embeds cleanly, which is why it gives 93%.

Two things came out of that. First, a pre-flight gate: embed, then count the true edges whose interaction falls below the planned Rabi peak and the non-edges above it, and refuse to submit anything — emulator minutes or QPU credits — when the register does not encode the graph. Knuth passes at 0 of 10; twelve-tone is refused at 96 of 144 on its best draw. The metric is a dozen lines:

p     = shaper._calculate_parameters(register=register, backend=backend, instance=instance)
omega = shaper.pulse(config, register, backend, instance).amplitude.samples.max()
conn  = np.asarray(p.connected)      # interaction on every true edge
disc  = np.asarray(p.disconnected)   # interaction on every non-edge
report = dict(
    edges_below_omega     = int((conn < omega).sum()),   # conflicts the drive overrides
    non_edges_above_omega = int((disc > omega).sum()),   # spurious conflicts
    unit_disk_valid       = bool(conn.min() > disc.max()),
)
if report["edges_below_omega"] or report["non_edges_above_omega"]:
    raise RuntimeError("register does not encode this graph — not submitting")

Second, a set of go/no-go rules for which exact-cover problems fit this hardware today, tested against three problems of my own. The universe has to live in two dimensions and each subset has to be a compact patch (twelve-tone fails here: each triad spans most of the octave). There can be at most about six candidate subsets per location, because co-located candidates form a clique and a clique has to fit inside one blockade disk (an edge-matching tiler of mine has forty variants per cell; constraint propagation from two fixed neighbours brings that to six or fewer, which is the honest place to hand a classical search's frontier to the machine). And the pairwise relation has to be overlap-if-and-only-if-proximity — the tiler fails this one, since "keys match" is neither. Equal-size patches placed on candidate centres satisfy all three by construction, which is why coverage planning is the canonical neutral-atom problem and chord progressions are not.

Exploring the boundary: cutting the torus

If the twelve-tone instance fails because it is dense and global, the next question is whether it is dense in every coordinate system. It isn't. Triads live on the Tonnetz — a lattice of pitches with fifths along one axis and major thirds along the other — where every major triad is an upward triangle and every minor triad a downward one, and two triads share a note exactly when their triangles touch. On the plane that is a local, geometric relation. What makes twelve-tone dense is that equal temperament wraps the lattice into a torus, and a torus will not lie flat in the atoms' plane.

Cutting the torus open at the best of its twelve seams and placing each triad once as a triangle: 95 of the 144 conflicts become planar, with a 5.4× energy margin between touching and non-touching centroids — a more robust embedding than the domino board. The other 49 are seam edges, conflicts that exist only through the wrap-around, and 21 of the 24 triads touch one. They are not dense in any basis; they are topological, and there is no cut that removes them. Relaxing them and filtering classically doesn't work either: without the seam the largest packings have six triangles rather than four, so the relaxed ground state and the real covers live in different energy sectors — zero of the four relaxed maxima are covers.

The clever move is to relax the tuning system instead. In just intonation the lattice is honestly planar — no wrap, no seam — and a chord cover over a finite patch of it is the same musical problem with one physical assumption removed. Planting an instance (the patch is whatever a random set of disjoint triangles happens to cover, so a cover is guaranteed), placing the atoms at the triangles' actual centroids, and running the gate: it passes at exactly the predicted 5.36×.

class TonnetzEmbedder(BaseEmbedder):
    """Atoms at the triads' centroids, scaled so the nearest pair clears min_atom_distance."""
    def embed(self, instance, config, backend):
        pts = [self.centroid[instance.index_to_node_label[i]] for i in range(len(instance.graph))]
        cx, cy = mean(p[0] for p in pts), mean(p[1] for p in pts)
        return Register({f"q{i}": ((x - cx) * self.scale, (y - cy) * self.scale)
                         for i, (x, y) in enumerate(pts)})

At fifteen atoms the atoms beat the classical cleanup for the first time in this project: same solver, same postprocessor, same shot count, and the unique cover sits at the top of the distribution instead of under a pile of traps. At twenty-four atoms — the same count as twelve-tone — the cover is still the most likely outcome, but the enrichment over random is gone. That was worth understanding rather than reporting.

Where the boundary actually is

Two probes on the fifteen-atom instance. First, the raw shots with no classical cleanup at all: 78.6% are valid packings — every one of the 48 conflicts respected — but most are small, sizes two and three, and only 2.4% are full covers. The blockade is doing its job; the sweep is not reaching the ground state. Second, the same geometry on a virtual device with no sequence-length cap: 6 µs gives 54% covers, 12 µs 84%, 24 µs 98%, and 48 µs gives 100%. Lowering the drive amplitude, by contrast, only hurt.

So the limit here is not shape — the instance passes every geometric test — but an energy-time budget. Six candidate triads around a pitch force the nearest atoms to the device's minimum spacing, the weakest constraint is then enforced by about 0.7 rad/µs, and multiplied by the device model's 6 µs maximum sequence that is a dimensionless budget of roughly four; the virtual device reaches full fidelity at about two hundred. The budget is set by three numbers on the spec sheet — the interaction coefficient, the minimum atom spacing, and the maximum sequence duration — and it tightens as the instance grows, which is exactly the difference between fifteen atoms and twenty-four. It is a fourth go/no-go test, the one the embedding gate cannot see, and the next thing to add to it once I have FRESNEL's own values.

That is the honest shape of the boundary as I've mapped it: the reduction is native to the hardware for planar, overlap-only packings with a small menu per location; the atoms enforce those constraints faithfully; and what decides whether they also find the answer is whether the device's sweep can be slow enough, or its interactions strong enough, for the instance's weakest constraint. Each of those is a measurement, not a guess, and each moved the project closer to the machine rather than further from it.

One more check, for completeness: the same Knuth instance as a textbook QUBO — minimize the squared over- and under-coverage of every element — with the atoms' positions fitted so their 1/r⁶ interactions match the couplings, and a global-detuning sweep on the QuTiP emulator. It returns the cover in 93% of shots, identical to the reduction; a per-atom-detuning version of the same sweep changes nothing. That is the expected answer, since an independent set on a blockade graph is the QUBO the hardware implements natively. The instructive part is the fit: ten graded couplings plus six that must vanish exceed what twelve planar coordinates can satisfy, and the optimizer quietly sacrifices one conflict — the answer survives only because that constraint happens to be redundant for this instance. The embedding wall is already visible at six atoms; the toy just doesn't lean on it.

Next steps

• Run it on the QPU. Two FRESNEL-valid sequences are ready and pass the gate — the six-atom Knuth instance and the fifteen-atom just-intonation Tonnetz. The remaining step is hardware access; a modest shot count first.
• Quantify the energy-time budget for FRESNEL. Its interaction coefficient, minimum spacing, and maximum sequence duration give the budget directly; add it to the gate as the fourth test, and use it to size Tonnetz patches the device can actually finish.
• Send the findings upstream. The amplitude-floor clamp and the missing embedding-validity warning are both a few lines in the solver library; drafts are written.
• Move emulation past sixteen atoms to the tensor-network backend. emu-mps is installed and should beat an hour at modest bond dimension — the 24-atom runs took 33 and 56 minutes on the state-vector emulator.