Port synth_circuit_from_stabilizers to Rust - #16610
Conversation
Port the Gaussian-elimination stabilizer-to-circuit synthesis to Rust. The Python wrapper keeps the input parsing (PauliList construction, phase and commutation validation) and passes the stabilizers' z/x/phase arrays to the Rust function, which returns the synthesized CircuitData. Instead of re-simulating the circuit built so far with Clifford(circuit) for every stabilizer as the Python implementation did, the Rust implementation maintains the Clifford tableau incrementally while emitting gates. This required a new Clifford::conjugate_pauli method in qiskit-quantum-info computing the Schrodinger-frame conjugation C P C-dag (the existing evolve_pauli computes the Heisenberg-frame C-dag P C), built on the same row-product phase bookkeeping. Error messages and the synthesized circuits are identical to the previous Python implementation. Towards Qiskit#12243 (part of Qiskit#12212). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
|
Thank you for opening a new pull request. Before your PR can be merged it will first need to pass continuous integration tests and be reviewed. Sometimes the review process can be slow, so please be patient. While you're waiting, please feel free to review other open PRs. While only a subset of people are authorized to approve pull requests for merging, everyone is encouraged to review open pull requests. Doing reviews helps reduce the burden on the core team and helps make the project's code better for everyone. One or more of the following people are relevant to this code:
|
|
|
jakelishman
left a comment
There was a problem hiding this comment.
Thanks for your interest. Please note: this is not a fully review including a validation of the mathematics. No maintainer will spend time to do that until we are sure that there is a human submitting the PR who also understands the mathematics, can explain it, and has reviewed it and verified it themselves (all without an LLM - it is fine to use as a tool, but it is not helpful to us if you couldn't do it yourself).
Also note: LLM-generated code should be attributed inline, similar to a licence disclosure for ported code.
The PR comment says "Towards #12243". Which part of that issue does it not solve?
| /// Tracks the circuit built so far both as a gate sequence and as a Clifford | ||
| /// tableau, so that stabilizers can be evolved through it. | ||
| struct TrackedCircuit { | ||
| gates: CliffordGatesVec, |
There was a problem hiding this comment.
What is the benefit of doing this rather than building a CircuitData directly?
There was a problem hiding this comment.
Two reasons, though I'm happy to change it if you'd prefer CircuitData directly:
- When
invert=Falsethe emitted circuit has to be inverted at the end (reverse order, S to Sdg). On the plainVecthat's a reversed iterator adapter feedingCircuitData::from_standard_gates; CircuitDatahas no cheap reverse, so I'd end up building a second one anyway. - It matches how the other Clifford-synthesis routines in this crate work (
greedy_synthesis.rs,bm_synthesis.rsboth build aCliffordGatesVecand convert at the boundary), and the Rust unit tests replay the gate list back through the tableau to check the diagonalization invariant, which is straightforward on the vec.
(1) is the main one. If you'd still rather build CircuitData directly, the invert=False path could emit into a fresh CircuitData in reverse at the end, let me know and I'll restructure.
| StandardGate::H => self.clifford.append_h(qubit), | ||
| StandardGate::S => self.clifford.append_s(qubit), | ||
| StandardGate::X => self.clifford.append_x(qubit), | ||
| _ => unreachable!(), |
There was a problem hiding this comment.
How is this provably unreachable?
There was a problem hiding this comment.
I guess it isn't; it's only a module-private invariant (push1/push2 are private and only called with these gates from within this file), which the compiler can't see, so unreachable!() was the wrong tool then . I'll restructure TrackedCircuit into one method per gate (h(), s(), x(), cx(), swap()), which removes both matches and both unreachable!()s entirely and makes the emitted gate set structural, that also localizes the S to Sdg inversion assumption you flagged further down.
| let in_z: Vec<bool> = z.to_vec(); | ||
| let in_x: Vec<bool> = x.to_vec(); | ||
| self.clifford.conjugate_pauli(&in_z, &in_x) |
There was a problem hiding this comment.
It should not be necessary to always allocate a new Vec; if nothing else we could have arranged for all the X and Z matrices to be contiguous by copying them into a new Vec on entry to the outer function, then we can guarantee that we can always take a slice.
There was a problem hiding this comment.
Agreed. The to_vec() was defensiveness against as_slice() returning None for a non-contiguous row view, but you're right that contiguity is better guaranteed once at entry. I'll have the outer function put the matrices in standard layout up front (a borrow, not a copy, in the common case the numpy arrays coming from astype are already C-contiguous), after which evolve can pass row slices straight through with no per-call allocation.
| /// stabilizers, using Gaussian elimination with Clifford gates. Based on the | ||
| /// stim implementation; see the Python-space | ||
| /// ``qiskit.synthesis.synth_circuit_from_stabilizers`` for details. |
There was a problem hiding this comment.
What Python-space implementation? This PR clears it all out.
This was a pre-existing problem, but please change "based on stim implementation" to a correct Apache 2 licence attribution of the Stim code - ports are still derivative works.
There was a problem hiding this comment.
Sorry, that sentence was meant to point at the public API function qiskit.synthesis.synth_circuit_from_stabilizers (which survives this PR as the wrapper, with the docstring and references), not at the removed implementation. But as written it reads exactly like the latter, so it's misleading either way; I'll reword it to cite Stim directly.
Agreed on the attribution. I'll replace "based on stim implementation" with a proper Apache-2.0 derivative-work notice in this file's header
| allow_redundant: bool, | ||
| allow_underconstrained: bool, | ||
| invert: bool, | ||
| ) -> Result<(usize, CliffordGatesVec), String> { |
There was a problem hiding this comment.
Result<T, String> is rarely correct. We should use a proper error type here - we can either return QiskitError directly, or use a thiserror-derived struct to defer creation of the format string until the message is displayed to the user.
| circuit.push1(StandardGate::H, pivot); | ||
| circuit.push1(StandardGate::S, j); | ||
| } | ||
| (false, false) => unreachable!(), |
There was a problem hiding this comment.
This would be better handled by not having the if !(curr_x[j] || curr_z[j]) continue; continue on lines 139--140, and putting a no-op on line 163, rather than an unexplained unreachable.
There was a problem hiding this comment.
That's cleaner for sure, I'll drop the !(curr_x[j] || curr_z[j]) half of the continue and make (false, false) an explicit no-op arm, so the match covers all four cases on its own with no unreachable!().
| // Invert the circuit: all emitted gates except S are self-inverse. | ||
| circuit | ||
| .gates | ||
| .into_iter() | ||
| .rev() | ||
| .map(|(gate, params, qubits)| { | ||
| let gate = match gate { | ||
| StandardGate::S => StandardGate::Sdg, | ||
| other => other, | ||
| }; | ||
| (gate, params, qubits) | ||
| }) | ||
| .collect() |
There was a problem hiding this comment.
The comment should at least by near line 198, but more importantly: line 202 is a very long-range assumption. If the function assumes that only S, H, X, CX and Swap are present, it's clearer to spell that out.
There was a problem hiding this comment.
You're right on both counts. I'll move the comment onto the match and spell the assumption out, something like:
let gate = match gate {
StandardGate::S => StandardGate::Sdg,
// The synthesis above emits only H, S, X, CX and Swap; all but S are self-inverse.
StandardGate::H | StandardGate::X | StandardGate::CX | StandardGate::Swap => gate,
_ => unreachable!("synthesis emits only H, S, X, CX and Swap"),▎
};
and with TrackedCircuit restructured into per-gate methods (your comment above), "only these five gates" becomes structurally true rather than an unstated assumption. If you'd rather keep other => other with the comment than list the variants, happy to do that instead.
| /// together with its group phase (0 or 2) in ``phase``. ``labels`` are the | ||
| /// stabilizer string representations, used only in error messages. |
There was a problem hiding this comment.
It seems rather unnecessary to allocate an entire list[str] in Python space and then extract it into an allocated Vec<String> in Rust, when it's not even supposed to get used, because it's only for errors. If it's absolutely necessary, can it not be recalculated if there's an error?
There was a problem hiding this comment.
Good point, it can go entirely and probably makes more sense. The label in the old Python messages was str(pauli) of the offending stabilizer, which is fully determined by the (z, x, phase) row the Rust code already has (sign prefix plus I/X/Y/Z per bit pair, highest qubit leftmost). So I'll drop the argument on both sides and reconstruct the label in Rust only on the error path. That removes the Python list[str] and the Vec<String>, keeps the messages byte-identical to before, and slots into the error-enum change from your other comment.
| /// Conjugate a dense Pauli, given by its per-qubit ``z`` and ``x`` bits, by | ||
| /// the Clifford ``C``, returning `C P C†` as ``(sign, z, x)`` dense bits. | ||
| /// | ||
| /// This is the Schrödinger-picture counterpart of [Clifford::evolve_pauli] | ||
| /// (which computes `C† P C`): the result is the ordered product of the | ||
| /// tableau rows selected by the input Pauli's bits. | ||
| pub fn conjugate_pauli(&self, z: &[bool], x: &[bool]) -> (bool, Vec<bool>, Vec<bool>) { |
There was a problem hiding this comment.
If this is so related to evolve_pauli, how come it's written in a completely different manner and has a different name?
There was a problem hiding this comment.
If this is so related to
evolve_pauli, how come it's written in a completely different manner and has a different name?
The relationship is "the opposite-direction conjugation", and that direction change turns out to change the natural algorithm completely:
For C P C† the tableau already is the answer: its rows are by definition the images of the generators under conjugation by C, so the result is just the phase-tracked product of the rows selected by P's bits, accessed through &self.
C† P C can't be read off those same rows (that would need the adjoint's tableau), which is why evolve_pauli instead goes through the PPR basis-change appends (_append_initial_part_ppr/_append_final_part_ppr) to reduce the dense Pauli to a single-qubit Z, which is why it takes &mut self. Reusing it here would have meant constructing the adjoint tableau per call, including recomputing its phase column which is more work than the direct row product.
They do share the part that's genuinely common; conjugate_pauli bottoms out in the same compute_phase_product_pauli helper that evolve_pauli uses via evolve_single_qubit_pauli.
On the name, I don't have a strong attachment, both operations are conjugations, so conjugate_pauli vs evolve_pauli admittedly doesn't convey the actual difference. The Python-space precedent is Pauli.evolve(..., frame="s"/"h"), so maybe evolve_pauli_schrodinger? Happy to rename to whatever you prefer, and to fold the explanation above into the docstring either way.
There was a problem hiding this comment.
see also the function evolve_single_qubit_pauli_dense in #16137.
although I think it's also evolving a Pauli by Heisenberg and not Schrodinger frame.
|
I know that the python code is pre-existing (and I have even reviewed it at some point), but I have a high-level conceptual question. If we had a complete stabilizer list, couldn't we first recover the corresponding destabilizers by solving a linear algebra problem (essentially inverting a matrix or performing Gaussian elimination), and then apply one of our existing Clifford synthesis algorithms (greedy, layered, etc.)? The greedy synthesis algorithm, for example, already performs a form of symplectic Gaussian elimination, and in addition has heuristics to reduce the number of CX gates in the synthesized circuit. Is the main motivation for this algorithm that we may not have a complete stabilizer list, or that we deliberately want to avoid computing the destabilizers? In other words, I'm trying to understand how this algorithm conceptually relates to our existing Clifford synthesis algorithms, and whether one can be viewed as a reduction of the other. As a possibly interesting experiment, we could generate random Cliffords, extract their complete stabilizer sets, and compare this algorithm against greedy synthesis in terms of both runtime and resulting CX count. |
|
All of these answers appear to me to be entirely LLM generated despite an explicit call out from me in the top comment that this was unacceptable. I will not participate further in this PR. |
I will admit that I did perhaps rely a little too heavily on LLMs to generate responses to the comments made on the code. In further work I will ensure that all reasoning and logic will come from myself only and I will only use LLMs for code generation and debugging. I will review the I do apologise for the complacency with with I have relied on LLMs for this PR but will ensure from here on that I have the mathematical reasoning and code implementation can be derived without said tools. I hope that the trust can be built again. |
| features_synthesis: | ||
| - | | ||
| The function :func:`.synth_circuit_from_stabilizers` was ported to Rust, leading to a | ||
| significant increase in performance. The Rust implementation maintains the Clifford |
There was a problem hiding this comment.
is there some benchmark that shows that there is a significant increase in performance?
could you present the details in a comment to this PR?
Port the Gaussian-elimination stabilizer-to-circuit synthesis to Rust. The Python wrapper keeps the input parsing (PauliList construction, phase and commutation validation) and passes the stabilizers' z/x/phase arrays to the Rust function, which returns the synthesized CircuitData.
Instead of re-simulating the circuit built so far with Clifford(circuit) for every stabilizer as the Python implementation did, the Rust implementation maintains the Clifford tableau incrementally while emitting gates. This required a new Clifford::conjugate_pauli method in qiskit-quantum-info computing the Schrodinger-frame conjugation C P C-dag (the existing evolve_pauli computes the Heisenberg-frame C-dag P C), built on the same row-product phase bookkeeping.
Error messages and the synthesized circuits are identical to the previous Python implementation.
Towards #12243 (part of #12212).
AI/LLM disclosure