🛠️ address review findings on the phase-sensitive algebra

Keep redundant generators instead of stripping them
---------------------------------------------------
`independent_rows()` in `StabilizerCode.__init__` also stripped redundant rows
from `CSSCode`, whose `Hx`/`Hz` are views into the generators. A check matrix
with a dependent row came back shorter than it went in, changing the syndrome
length; redundancy is deliberate for single-shot decoding and meta-checks.

The independence was only needed so `equal_stabilizer_group` could compare
`n_rows` against the rank. It now compares group order plus one-sided
containment, so redundant generating sets are accepted and kept as given, and
the `mod2.row_basis` rewriting in `CSSCode` is no longer necessary.

`equal_stabilizer_group` also raised `InvalidPauliError` when both codes had no
generators, because `from_paulis` rejects an empty list. Two such codes both
generate the trivial group; `is_equivalent` and `get_logical_mapping` were
affected through it.

Rename `phase` to `phase_exponent` / `phase_exponents`
------------------------------------------------------
The attribute kept its name while changing meaning from a binary sign to a Z4
exponent, so `Pauli(support, 0)` silently went from `+Y` to `-iY` and the last
column of `as_vector()`/`as_matrix()` from binary to {0,1,2,3}. Renaming makes
downstream code fail loudly instead of computing the wrong thing.

Stop `PauliTableau.__getitem__` from aliasing
----------------------------------------------
It returned a `Pauli` whose symplectic support was a view into the tableau while
its phase was a copy, so mutating the result wrote back only half way. Both are
copies now, which also makes `pauli_row_echelon` structurally unable to modify
its argument rather than relying on `__mul__` allocating.

Add `PauliTableau.multiply_rows`
---------------------------------
Combining rows of a signed tableau is not an XOR, and `optimize_tableau` spelled
out the phase bookkeeping four times over fancy-indexed copies. The helper makes
the correct operation the easy one and removes the duplication.

Performance
-----------
The elimination ran O(rows * n) `Pauli.__mul__` calls in a Python loop, and
membership tests redid it twice per query. The loop is vectorized and
`StabilizerCode` reuses the echelon form it already computes when checking for
-I. For an n=761 hypergraph-product code, construction goes 2.40s -> 0.83s and
`is_stabilizer` 3.12s -> under 10ms.

Also return a `PauliRowEchelon` named tuple rather than a five-tuple, and rename
its third field to `n_global_phases` since it is a count, not a rank.

Document the conventions
------------------------
`docs/PauliConventions.md` records the i^p X^x Z^z representation, why a
positive Hermitian Pauli has p = x.z rather than 0, why XOR-ing signs is not
Pauli multiplication, and which layer enforces Hermiticity.

Note: `_stabilizer_sign_corrections` in `exact/search.py` still XOR-combines
signs across a GF(2) row reduction, which is invalid for non-CSS codes; that is
left for a separate issue.

Assisted-by: Claude Opus 5 via Claude Code

Author: lsschmid

docs/PauliConventions.md

@@ -0,0 +1,121 @@
+# Pauli Conventions
+
+This page pins down how Pauli operators, their phases, and stabilizer groups are
+represented in `mqt.qecc.codes.core.pauli`. The conventions matter: a phase that
+is read as if it were a sign bit produces silently wrong results rather than an
+error.
+
+## The representation
+
+A `Pauli` on $n$ qubits is stored as a binary symplectic support $(x \mid z)$
+together with a **phase exponent** $p \in \{0,1,2,3\}$:
+
+$$
+P = i^{p}\, X^{x} Z^{z}
+\qquad\text{where}\qquad
+X^{x} Z^{z} = \bigotimes_{j=1}^{n} X^{x_j} Z^{z_j}
+$$
+
+The exponent is stored in `Pauli.phase_exponent`; for a `PauliTableau` the
+per-row exponents are in `PauliTableau.phase_exponents`.
+
+The single-qubit letters follow from this with $Y = iXZ$:
+
+| letter | $(x_j \mid z_j)$ | contribution to $p$ |
+| ------ | ---------------- | ------------------- |
+| `I`    | $(0 \mid 0)$     | 0                   |
+| `X`    | $(1 \mid 0)$     | 0                   |
+| `Z`    | $(0 \mid 1)$     | 0                   |
+| `Y`    | $(1 \mid 1)$     | 1                   |
+
+So a Hermitian Pauli with a `+` sign has $p = x \cdot z$, which is the number of
+`Y` letters modulo four — **not** zero. `Pauli(support)` with no explicit
+exponent picks exactly this canonical positive Hermitian choice.
+
+## Exponents are not sign bits
+
+`phase_exponent` carries four values, not two. The two are related by
+
+$$
+P = (-1)^{r}\, i^{\,x \cdot z}\, X^{x} Z^{z},
+\qquad
+r = \frac{p - x \cdot z}{2} \bmod 2
+$$
+
+Use the explicit converters rather than reaching for the raw exponent:
+
+- `Pauli.sign()` / `PauliTableau.signs()` return the binary sign $r$, and raise
+  `InvalidPauliError` on a non-Hermitian operator, which has no real sign.
+- `Pauli.from_symplectic_and_sign(support, sign)` and
+  `PauliTableau.phase_from_signs(matrix, signs)` go the other way.
+- `Pauli.is_hermitian()` / `PauliTableau.is_hermitian()` test whether a sign
+  exists at all.
+
+## Multiplication needs a correction term
+
+Because $Z^{z_1} X^{x_2} = (-1)^{z_1 \cdot x_2} X^{x_2} Z^{z_1}$, the product of
+two Paulis is
+
+$$
+P_1 P_2 = i^{\,p_1 + p_2 + 2 (z_1 \cdot x_2)}
+X^{x_1 \oplus x_2} Z^{z_1 \oplus z_2}
+$$
+
+The extra $2(z_1 \cdot x_2)$ is why
+**XOR-ing symplectic rows and XOR-ing their signs is not Pauli multiplication**.
+Concretely:
+
+$$(X \otimes X)(Z \otimes Z) = -\,Y \otimes Y$$
+
+while XOR-ing the two `+` signs would predict $+\,Y \otimes Y$. Note the
+correction depends on the number of qubits: $(XXXX)(ZZZZ) = +YYYY$.
+
+Consequences for anyone combining rows of a signed tableau:
+
+- Use `PauliTableau.multiply_rows(target, source)`, never a raw XOR on
+  `tableau.tableau.data` followed by an XOR on the phases.
+- Use `pauli_row_echelon`, not `mod2.row_echelon`, whenever phases must survive
+  the reduction. A plain mod-2 reduction is only safe on a CSS tableau, where
+  the pivoting never combines an X-type row with a Z-type row and the correction
+  term vanishes.
+- `PauliTableau.independent_rows()` selects rows by support only and is
+  explicitly phase-insensitive; do not use it to decide anything about signs.
+
+## Which layer enforces what
+
+The two layers deliberately allow different things:
+
+- **`Pauli` / `PauliTableau` represent the full $n$-qubit Pauli group
+  $\mathfrak{P}_n$.** Non-Hermitian elements such as `+iX` are legal and
+  necessary: row reduction genuinely produces them as intermediates, since
+  $X \cdot Z = -iY$.
+- **`StabilizerCode` enforces the stabilizer conditions.** Its constructor
+  rejects generators that do not commute, are not Hermitian, or together
+  generate $-I$. Those checks — not the Pauli layer — are what guarantee a valid
+  code.
+
+Generators need **not** be independent. A redundant generating set is accepted
+and kept as given, so `CSSCode` preserves the check matrices you pass in,
+including redundant rows that matter for single-shot decoding and meta-checks.
+Group-level comparisons (`equal_stabilizer_group`, `stabilizer_equivalent`,
+`is_stabilizer`) compare the generated groups and are unaffected by redundancy.
+
+## Subgroups and rank
+
+`pauli_row_echelon` returns a `PauliRowEchelon`, whose `rank` is $\log_2 |G|$
+for the generated subgroup $G$. This includes the central scalars, so it can
+exceed the number of pivot columns:
+
+| generators    | generated subgroup                                | order | `rank` |
+| ------------- | ------------------------------------------------- | ----- | ------ |
+| `["XX","ZZ"]` | $\{I, XX, ZZ, -YY\}$ — no scalars beyond $I$      | 4     | 2      |
+| `["+iX"]`     | $\{I, iX, -I, -iX\}$ — one pivot, four scalars    | 4     | 2      |
+| `["X","Z"]`   | $\{\pm I, \pm X, \pm Z, \pm iY\}$ — anticommuting | 8     | 3      |
+
+The middle row has a single pivot column yet rank 2: the generator squares to
+$-I$, so the subgroup contains scalars the support alone cannot account for. The
+last row picks up $-I$ from the anticommutator.
+
+To test many Paulis against one subgroup, compute the echelon once and call
+`pauli_in_reduced_subgroup`; `PauliTableau.is_in_subgroup` redoes the
+elimination on every call.

docs/index.md

@@ -30,6 +30,7 @@ self
 :maxdepth: 1
 
 installation
+PauliConventions
 LightsOutDecoder
 StatePrep
 CatStates

src/mqt/qecc/circuit_synthesis/encoding.py

@@ -18,8 +18,7 @@
 from mqt.qecc import mod2
 
 from ..codes import CSSCode
-from ..codes.core.pauli import CheckMatrix, Pauli, PauliTableau, complete_stabilizer_tableau_with_destabilizers
-from ..codes.core.symplectic import SymplecticVector
+from ..codes.core.pauli import CheckMatrix, PauliTableau, complete_stabilizer_tableau_with_destabilizers
 from .circuits import CliffordIsometry, CNOTCircuit
 from .exact import (
     Objective,
@@ -370,17 +369,17 @@ def synthesize_encoding_circuit(
         additional_z_rows,
     ))
     gens_phase: npt.NDArray[np.int8] = np.hstack((
-        code.generators.phase,
-        code.x_logicals.phase[additional_x_checks],
-        code.z_logicals.phase[additional_z_checks],
+        code.generators.phase_exponents,
+        code.x_logicals.phase_exponents[additional_x_checks],
+        code.z_logicals.phase_exponents[additional_z_checks],
     ))
     log_mat: npt.NDArray[np.int8] = np.vstack((
         np.delete(code.x_logicals.tableau.data, additional_checks, axis=0),
         np.delete(code.z_logicals.tableau.data, additional_checks, axis=0),
     ))
     log_phase: npt.NDArray[np.int8] = np.hstack((
-        np.delete(code.x_logicals.phase, additional_checks, axis=0),
-        np.delete(code.z_logicals.phase, additional_checks, axis=0),
+        np.delete(code.x_logicals.phase_exponents, additional_checks, axis=0),
+        np.delete(code.z_logicals.phase_exponents, additional_checks, axis=0),
     ))
 
     return encoder_from_stabilizers_and_logicals(
@@ -465,61 +464,33 @@ def optimize_tableau(tableau: PauliTableau, stab_rows: list[int]) -> PauliTablea
                 if i == j:
                     continue
                 tab = tableau.copy()
-                mat = tab.tableau.data
-                destabs = mat[:half][stab_rows]
-                stabs = mat[half:][stab_rows]
 
                 stab_i_abs, stab_j_abs = half + stab_rows[i], half + stab_rows[j]
                 destab_i_abs, destab_j_abs = stab_rows[i], stab_rows[j]
 
-                new_stab_i = Pauli(SymplecticVector(stabs[i].copy()), tab.phase[stab_i_abs]) * Pauli(
-                    SymplecticVector(stabs[j].copy()), tab.phase[stab_j_abs]
-                )
-                new_destab_j = Pauli(SymplecticVector(destabs[j].copy()), tab.phase[destab_j_abs]) * Pauli(
-                    SymplecticVector(destabs[i].copy()), tab.phase[destab_i_abs]
-                )
-
-                stabs[i] ^= stabs[j]
-                destabs[j] ^= destabs[i]
-                mat[:half][stab_rows] = destabs
-                mat[half:][stab_rows] = stabs
-                tab.phase[stab_i_abs] = new_stab_i.phase
-                tab.phase[destab_j_abs] = new_destab_j.phase
-
-                new_score, _ = score_symplectic(PauliTableau(mat, tab.phase.copy()))
+                # Multiplying a stabilizer onto another must be mirrored on the
+                # destabilizers in the opposite direction to preserve the symplectic form.
+                tab.multiply_rows(stab_i_abs, stab_j_abs)
+                tab.multiply_rows(destab_j_abs, destab_i_abs)
+
+                new_score, _ = score_symplectic(tab)
                 if lexicographical_compare_np(new_score, best[1]):
                     best = (tab, new_score)
                     improved = True
             for j in range(len(logical_rows)):
                 tab = tableau.copy()
-                mat = tab.tableau.data
-                destabs = mat[:half][stab_rows]
-                stabs = mat[half:][stab_rows]
 
-                other_log = mat[logical_rows[(j + k) % (2 * k)]]
                 other_log_abs = logical_rows[(j + k) % (2 * k)]
                 destab_i_abs = stab_rows[i]
                 log_j_abs = logical_rows[j]
                 stab_i_abs = half + stab_rows[i]
 
-                new_destab_i = Pauli(SymplecticVector(destabs[i].copy()), tab.phase[destab_i_abs]) * Pauli(
-                    SymplecticVector(other_log.copy()), tab.phase[other_log_abs]
-                )
-
-                destabs[i] ^= other_log
-                logj = mat[logical_rows[j]]
-
-                new_logj = Pauli(SymplecticVector(logj.copy()), tab.phase[log_j_abs]) * Pauli(
-                    SymplecticVector(stabs[i].copy()), tab.phase[stab_i_abs]
-                )
-
-                logj ^= stabs[i]
-                mat[:half][stab_rows] = destabs
-                mat[logical_rows[j]] = logj
-                tab.phase[destab_i_abs] = new_destab_i.phase
-                tab.phase[log_j_abs] = new_logj.phase
+                # The stabilizer row is read by the second multiplication, so it must
+                # not be one of the rows written by the first.
+                tab.multiply_rows(destab_i_abs, other_log_abs)
+                tab.multiply_rows(log_j_abs, stab_i_abs)
 
-                new_score, _ = score_symplectic(PauliTableau(mat, tab.phase.copy()))
+                new_score, _ = score_symplectic(tab)
                 if lexicographical_compare_np(new_score, best[1]):
                     best = (tab, new_score)
                     improved = True
@@ -547,11 +518,11 @@ def combine_stabilizer_and_logical_tableau(stabilizers: PauliTableau, logicals:
     # Combine stabilizers and logicals into a single tableau
     x_logicals = logicals.tableau.data[: logicals.num_rows() // 2]
     z_logicals = logicals.tableau.data[logicals.num_rows() // 2 :]
-    x_logicals_phase = logicals.phase[: logicals.num_rows() // 2]
-    z_logicals_phase = logicals.phase[logicals.num_rows() // 2 :]
+    x_logicals_phase = logicals.phase_exponents[: logicals.num_rows() // 2]
+    z_logicals_phase = logicals.phase_exponents[logicals.num_rows() // 2 :]
     combined_matrix = np.vstack([x_logicals, z_logicals, stabilizers.tableau.data])
 
-    combined_phase = np.hstack([x_logicals_phase, z_logicals_phase, stabilizers.phase])
+    combined_phase = np.hstack([x_logicals_phase, z_logicals_phase, stabilizers.phase_exponents])
     combined_tableau = PauliTableau(combined_matrix, combined_phase)
 
     # Complete with destabilizers for the stabilizers only

src/mqt/qecc/circuit_synthesis/exact/search.py

@@ -506,9 +506,9 @@ def _combine_stabilizers_and_logicals(
     ])
 
     combined_phase = np.concatenate([
-        x_logicals.phase,
-        z_logicals.phase,
-        stabilizers.phase,
+        x_logicals.phase_exponents,
+        z_logicals.phase_exponents,
+        stabilizers.phase_exponents,
     ])
 
     return StabilizerTableau(combined_matrix, combined_phase)

src/mqt/qecc/circuit_synthesis/exact/verification.py

@@ -49,13 +49,13 @@ def verify_clifford_unitary(circuit: CliffordIsometry, target: StabilizerTableau
 
     act_x: npt.NDArray[np.int8] = actual.tableau.data[:n]
     act_z: npt.NDArray[np.int8] = actual.tableau.data[n:]
-    act_xp: npt.NDArray[np.int8] = actual.phase[:n]
-    act_zp: npt.NDArray[np.int8] = actual.phase[n:]
+    act_xp: npt.NDArray[np.int8] = actual.phase_exponents[:n]
+    act_zp: npt.NDArray[np.int8] = actual.phase_exponents[n:]
 
     tgt_x: npt.NDArray[np.int8] = target.tableau.data[:n]
     tgt_z: npt.NDArray[np.int8] = target.tableau.data[n:]
-    tgt_xp: npt.NDArray[np.int8] = target.phase[:n]
-    tgt_zp: npt.NDArray[np.int8] = target.phase[n:]
+    tgt_xp: npt.NDArray[np.int8] = target.phase_exponents[:n]
+    tgt_zp: npt.NDArray[np.int8] = target.phase_exponents[n:]
 
     # Find permutation π such that actual.X[q] == target.X[π(q)] for each q.
     perm: dict[int, int] = {}
@@ -117,9 +117,9 @@ def verify_clifford_isometry(
     if num_rows != expected_rows:
         return False
 
-    x_logicals = StabilizerTableau(target.tableau.data[:k, :], target.phase[:k])
-    z_logicals = StabilizerTableau(target.tableau.data[k : 2 * k, :], target.phase[k : 2 * k])
-    stabilizers = StabilizerTableau(target.tableau.data[2 * k :, :], target.phase[2 * k :])
+    x_logicals = StabilizerTableau(target.tableau.data[:k, :], target.phase_exponents[:k])
+    z_logicals = StabilizerTableau(target.tableau.data[k : 2 * k, :], target.phase_exponents[k : 2 * k])
+    stabilizers = StabilizerTableau(target.tableau.data[2 * k :, :], target.phase_exponents[2 * k :])
 
     circuit_code = circuit.get_code()
     target_code = StabilizerCode(stabilizers, x_logicals=x_logicals, z_logicals=z_logicals)

src/mqt/qecc/circuit_synthesis/operations.py

@@ -125,7 +125,7 @@ def apply_stabilizer_tableau_inplace(self, tableau: PauliTableau) -> None:
         _apply_transvection_numba(
             mat[:, self.i], mat[:, self.i + n], mat[:, self.j], mat[:, self.j + n], signs, *self.v
         )
-        tableau.phase = PauliTableau.phase_from_signs(mat, signs)
+        tableau.phase_exponents = PauliTableau.phase_from_signs(mat, signs)
 
     def apply_stabilizer_tableau(self, tableau: PauliTableau, inplace: bool = False) -> PauliTableau:
         """Apply the transvection operation to a stabilizer tableau."""
@@ -138,7 +138,7 @@ def apply_stabilizer_tableau(self, tableau: PauliTableau, inplace: bool = False)
         _apply_transvection_numba(
             mat[:, self.i], mat[:, self.i + n], mat[:, self.j], mat[:, self.j + n], signs, *self.v
         )
-        out.phase = PauliTableau.phase_from_signs(mat, signs)
+        out.phase_exponents = PauliTableau.phase_from_signs(mat, signs)
 
         return out
 

src/mqt/qecc/circuit_synthesis/transvection.py

@@ -138,7 +138,7 @@ def _score_transvections(
     scored: list[tuple[TableauOperation, int | tuple[int, ...]]] = []
 
     original_state = tableau.tableau.data.copy()
-    original_phase = tableau.phase.copy()
+    original_phase = tableau.phase_exponents.copy()
     for op in operations:
         op.apply_stabilizer_tableau_inplace(tableau)
         h_vec, _ = score_symplectic(tableau)
@@ -148,7 +148,7 @@ def _score_transvections(
             scored.append((op, score_value))
 
         tableau.tableau.data[:] = original_state
-        tableau.phase[:] = original_phase
+        tableau.phase_exponents[:] = original_phase
     scored.sort(key=operator.itemgetter(1))
     return scored
 

src/mqt/qecc/codes/constructions/concatenation.py

@@ -82,7 +82,7 @@ def _outer_pauli_to_physical(self, p: Pauli) -> Pauli:
             msg = "The Pauli operator must have the same number of qubits as the outer code."
             raise InvalidStabilizerCodeError(msg)
         concatenated = SymplecticVector.zeros(self.n)
-        phase = p.phase
+        phase = p.phase_exponent
         offset: int = 0
         for i in range(self.outer_code.n):
             c = self.inner_codes[i]
@@ -100,7 +100,7 @@ def _outer_pauli_to_physical(self, p: Pauli) -> Pauli:
                 continue
             concatenated[offset:new_offset] = block.x_part()
             concatenated[offset + self.n : new_offset + self.n] = block.z_part()
-            phase = (phase + block.phase) % 4
+            phase = (phase + block.phase_exponent) % 4
             offset = new_offset
         return Pauli(concatenated, phase)
 

src/mqt/qecc/codes/core/css_code.py

@@ -53,17 +53,13 @@ def __init__(
             self._check_valid_check_matrices(Hx, Hz)
             if Hx is not None:
                 inferred_n = Hx.shape[1]
-                hx = np.asarray(mod2.row_basis(Hx), dtype=np.int8)
-                hz = (
-                    np.asarray(mod2.row_basis(Hz), dtype=np.int8)
-                    if Hz is not None
-                    else np.zeros((0, inferred_n), dtype=np.int8)
-                )
+                hx = np.asarray(Hx, dtype=np.int8)
+                hz = np.asarray(Hz, dtype=np.int8) if Hz is not None else np.zeros((0, inferred_n), dtype=np.int8)
             else:
                 assert Hz is not None
                 inferred_n = Hz.shape[1]
                 hx = np.zeros((0, inferred_n), dtype=np.int8)
-                hz = np.asarray(mod2.row_basis(Hz), dtype=np.int8)
+                hz = np.asarray(Hz, dtype=np.int8)
 
         num_qubits = hx.shape[1]
         if n is not None and n != num_qubits:
@@ -92,7 +88,7 @@ def __init__(
         x_logicals = PauliTableau.from_check_matrix(CheckMatrix(lx, "X"))
         z_logicals = PauliTableau.from_check_matrix(CheckMatrix(lz, "Z"))
 
-        self._num_x_checks = mod2.rank(hx)
+        self._num_x_checks = hx.shape[0]
         super().__init__(
             generators,
             distance,