[ondrik]

Fixes #65: F1 (SMT, from this PR) and F15 (MILP, from #106, now included here). F58 continues in #104.

Both encodings bring the output and the target to a common k (the FiveTuple rescaling), and both took the bound of that k from Synthesizer.max_k, which is computed once from the input. Past the input depth, or with custom pairs whose k is not reduced, the bound is too small:

• SMT: the rescaled vectors were left unconstrained, so wrong circuits were accepted.
• MILP: the correct circuit became infeasible, so valid circuits were missed. CBC never returned a wrong one.

SMT (F1)

Problem

Generator.add_rescaling brings the output and the target to a common k. It covers k = floor(|k_out - k_target| / 2) only for k < Synthesizer.max_k, which is computed once from the input. Past the input depth, the rescaled vectors were left unconstrained, so any circuit passed the check. This happens with incremental search, -d, restricted gate sets, custom vector pairs and pareto-incremental. Example: t;x;t;x (= ω·I) with -d 1 gave h;h;h;h (= I), and z from {h, sx} gave sx;sx (= X).

Fix

• encode_equivalence passes a bound for each pair and depth: max(k_in + 2·depth - k_target, k_target - k_in) // 2 + 1. Here 2 is the largest k increment of a layer; the MILP commit replaces it with the gate set's own (see below).
• add_rescaling asserts k < d. If a gap ever falls outside the case split, the result is now UNSAT, not an unconstrained formula. The d = 2 special case is gone.
• The pairs' k must be an integer; 4.0 is accepted, anything else raises a ValueError.

Testing

• New tests: tests/test_rescaling_bound.py has 8 cases, with Qiskit as the reference. The 5 original ones fail on master and all 8 pass here.
• Full suite: 278 tests passed before the second commit. The affected suites (custom vectors, global phase, incremental, validation, pareto) pass after it.
• Differential testing against Qiskit:
  • 441 random and hand-built cases, compared with a Qiskit brute-force search that finds the true optimum or proves none exists.
  • Coverage: 1–2 qubits, restricted gate sets, -d 1/2, every search mode including every pareto_front/ file, global phase on and off, custom pairs with unreduced k, yices2 and cvc5.
  • Results here: 0 wrong successes, 0 spurious failures, 0 crashes. The bound is tight: sx⁴ from {sx} needs exactly the last case.
  • On master, 66 of the same cases returned wrong circuits.
• Mutation check: dropping the +1, or setting the largest k increment to 1 (then MAX_K_INCREMENT), makes the edge-case tests fail.
• Performance: I compared yices2 runs on 8 benchmarks (GHZ, Bernstein–Vazirani, mqt-bench dj/bv/grover/modular_adder), in incremental and basic modes, with -v all and -v zero. Each configuration ran 3 times on each version, interleaved. Every circuit was checked against Qiskit.
  • Wall time, typical run (geometric mean): 0.66× master in incremental mode and 0.62× in basic mode. For example, modular_adder basic/all went from 30.1 s to 5.5 s and grover_2 incremental/all from 2.1 s to 0.5 s.
  • Formula size: the written .smt2 is 10–26% smaller, and the incremental stream is 0.22–0.69× of master's. Past about twice the input depth, the case split is larger than master's. This is needed for correctness.
  • Slower case: in basic mode, yices is reproducibly 1.3–3× slower on 6 of 15 instances (for example dj_2 basic/zero, 0.73 s → 2.2 s). The formula is smaller in all of them.
    • In 4 of the 6, the new (< k d) line triggers the slowdown. In other instances, removing that line makes yices slower. This looks like yices being sensitive to small formula changes, not extra work.
    • Total yices CPU over the 15 basic configurations drops from 30.4 s to 22.6 s.
  • Results: incremental gate counts are identical in all configurations. The circuits differ only where master was wrong, or in basic mode, which doesn't minimise with SMT.
  • Tests: tests/test_incremental.py and tests/test_synthesis.py take 181 s and 192 s here vs 214 s and 250 s on master. Most of that is the GHZ-3 cvc5 test.

MILP (F15, from #106 by @jkhav)

Problem

add_milp_rescaling took its k bound and both big-Ms from Synthesizer.max_k. Past the input depth, and with custom pairs whose k is not reduced, the correct circuit was infeasible:

• X from {sx, z}: no sx;sx
• Z from {h, sx}: no h;sx;sx;h
• custom pair |0⟩ → |1⟩ with the target written as a=4, k=4: no x
• custom input with k=6 (a=1): no x

Fix

• Same bound as the SMT side. max_gap = max(k_in + inc·depth − k_target, k_target − k_in) bounds |k1 − k2|. The rescaling k is in [0, max_gap // 2], the parity big-M is 2·max_gap + 1, and the value big-M is B + R, with B the vectors' entry bound and R = 2B·2^(max_gap // 2) the bound of the Rescaled vectors (no longer 2**max_k).
• Stale PuLP bounds. In PuLP the bound belongs to the variable, and a redeclared name gets back the first variable:
  • k of the rescaling is redeclared at each depth, so it has no ub now; the one-hot bounds it.
  • The target's Rescaled2 vector was named after the target only, so every depth got back the first depth's variables and bound. Since incremental and pareto-incremental never try the empty circuit #67, the first depth is 0, with the smallest gap, and CBC again missed h;sx;sx;h for Z from {h, sx}. It is named per depth now (last commit, added while consolidating).
• Intermediate k bound. Inter's k bound was 2·(d+1), ignoring the input's own k. It is k_in + inc·d now. This is the k half of MILP/Gurobi state bounds ignore the input state #69 F4.
• Per-gate increments. inc is the largest k increment in the gate set, from a k_increments table in gates.py. encode_layer reads the same table, so a new gate that raises k must be added there (KeyError otherwise) instead of silently breaking a hardcoded 2. Gate sets without h/sx get a gap of 0 and smaller big-Ms. It replaces MAX_K_INCREMENT on the SMT side too.

Why it is sound

With Δ = k1 − k2, K = ⌊|Δ|/2⌋ and q = |Δ| mod 2, the equality is v1·2^K·√2^q = v2 (or swapped), and √2·(a,b,c,d) = (b−d, a+c, b+d, c−a), as encoded.

• Each layer holds exactly one gate (ExactlyOne), so |Δ| ≤ max_gap and K ≤ max_gap // 2. A larger gap is infeasible, never unconstrained.
• Parity rows: the relaxed side reaches 2|Δ| ≤ 2·max_gap. The active side forces 2k + q = |Δ|; lb=0 rules out the wrong sleq.
• Value rows: an inactive row is at most c·|expr| + R, with c ≤ 2^(max_gap // 2) and |expr| ≤ 2B. The row multipliers in constrained_rescaling (2·tmp, 2c·tmp, 3c·tmp, with tmp = B + R) cover each case.

This holds in exact arithmetic. MILP solvers accept values within a tolerance, so very large big-Ms are a numerical risk; that is #104.

Testing

• New tests in tests/test_rescaling_bound.py: the four MILP cases above (all fail without the MILP commit), and max_k_increment for four gate sets.
• The rescaling on its own vs exact values: a script builds only add_milp_rescaling with fixed vectors and compares CBC's feasibility with complex-float ground truth from the definition. 1,848 solves (630 equal pairs, 870 not), with tight, off-by-one-parity and loose max_gap, plus too-small bounds that must be infeasible: 0 mismatches. As a check on the checker, a parity big-M of 2·max_gap − 1 gives 58 mismatches and a value big-M of (B+R)/4 gives 11.
• Still to do: benchmark the MILP times with the new big-Ms.

Rebase

Rebased onto master (501bc62). The MILP commit conflicted with #79 (Pauli basis removed, encode_layer dedented); its k_increments[gate] and k_bound changes were reapplied to master's code. Full suite: 653 passed.

Not in this PR

• MILP FiveTuple equivalence relies on big-M rows that solver tolerances can relax #104: numerical robustness of the big-M rows (F58), including PuLP's Gurobi backend running with Gurobi's default IntFeasTol of 1e-5.
• MILP/Gurobi state bounds ignore the input state #69, MILP big-M constants are too small, so CBC results are not gate-count optimal #68: custom inputs with large entries (e.g. a=4, k=4) are still infeasible in MILP, because of the coefficient bound of the intermediate states (MILP/Gurobi state bounds ignore the input state #69 F4) and the per-layer big-M (MILP big-M constants are too small, so CBC results are not gate-count optimal #68 F5), not the rescaling.
• Approximate synthesis (Approximate synthesis (-ap) accepts wrong circuits and rejects exact ones #66, F2): the -ap path has the same stale bound in FiveTuple.to_real. z from {h, sx} with approx=True still returns sx;sx.

🤖 Generated with Claude Code

https://claude.ai/code/session_0177AhE9SXm3iJ1KP8GTSk6o

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.

📢 Thoughts on this report? Let us know!