qLDPCOrg
diff --git a/‎src/qldpc/codes/common.py‎
Lines changed: 44 additions & 44 deletions b/‎src/qldpc/codes/common.py‎
Lines changed: 44 additions & 44 deletions
diff --git a/‎src/qldpc/codes/common_test.py‎
Lines changed: 120 additions & 0 deletions b/‎src/qldpc/codes/common_test.py‎
Lines changed: 120 additions & 0 deletions
diff --git a/‎src/qldpc/math.py‎
Lines changed: 61 additions & 0 deletions b/‎src/qldpc/math.py‎
Lines changed: 61 additions & 0 deletions
@@ -915,7 +915,7 @@ def __str__(self) -> str:
     def is_subsystem_code(self) -> bool:
         """Is this code a subsystem code?
 
-        That is, do all parity checks commute?
+        That is, do some parity checks fail to commute with each other?
         """
         if self._is_subsystem_code is None:
             self._is_subsystem_code = bool(
@@ -1172,17 +1172,17 @@ def get_logical_ops(
         If this method is passed a pauli operator (Pauli.X or Pauli.Z), it returns only the logical
         operators of that type.
 
-        Due to the way that logical operators are constructed in this method, logical Z-type
-        operators only address physical qudits by physical Z-type operators, while logical X-type
-        operators address at least one physical qudit with a physical X-type operator, and may
-        additionally address physical qudits with physical Z-type operators.
+        For a non-subsystem code, logical X-type operators have X-support and may additionally have
+        Z-support, while logical Z-type operators have only Z-support and no X-support.  These
+        logical operators are constructed with a method similar to that in Section 4.1 of
+        Gottesman's thesis (arXiv:9705052): fix the values of the logical operator matrix in the GL
+        sector of the parity check matrix when written in standard form (see
+        QuditCode.get_standard_form_data), and fill in the remaining entries of the logical operator
+        matrix as required by commutation constraints.
 
-        Logical operators are constructed with the method similar to that in Section 4.1 of
-        Gottesman's thesis (arXiv:9705052), generalized for subsystem qudit codes.  The basic
-        strategy is to fix the values of the logical operator matrix in the GL sector of the parity
-        check matrix when written in standard form (see QuditCode.get_standard_form_data), and then
-        fill in the remaining entries of the logical operator matrix as required by parity check
-        constraints.
+        For a subsystem code, logical operators are constructed as a symplectic basis of the gauge
+        group's centralizer, whose symplectic radical is the stabilizer group.  These logical
+        operators make no guarantees on their physical support.
 
         The symplectic argument is provided for compatibility with CSSCode.get_logical_ops, and must
         always be True for a non-CSS code.
@@ -1198,12 +1198,24 @@ def get_logical_ops(
         if not (self._logical_ops is None or recompute):
             return self._logical_ops
 
+        # For a subsystem code, extract the logical operators as a symplectic basis of the gauge
+        # group's centralizer C(G): the operators commuting with every gauge generator.  The
+        # symplectic radical of C(G) is the stabilizer group, so a symplectic (hyperbolic) basis of
+        # C(G) is exactly the k dual pairs of logical operators.
+        if self.is_subsystem_code:
+            centralizer = math.symplectic_conjugate(self.canonicalized.matrix).null_space()
+            logical_ops, _radical = math.symplectic_gram_schmidt(
+                centralizer, promise_full_rank=True
+            )
+            self._logical_ops = logical_ops
+            return self._logical_ops
+
         # construct the standard-form parity check matrix
         (
             matrix,
             qudit_locs,
-            (rows_sx, rows_gx, rows_sz, rows_gz),
-            (cols_sx, cols_gx, cols_lx, cols_sz, _cols_gz, cols_lz),
+            (rows_sx, _rows_gx, rows_sz, _rows_gz),
+            (cols_sx, _cols_gx, cols_lx, cols_sz, _cols_gz, cols_lz),
         ) = self.get_standard_form_data()
         matrix_x = matrix[:, 0, :]
         matrix_z = matrix[:, 1, :]
@@ -1213,34 +1225,8 @@ def get_logical_ops(
         logicals_zz = self.field.Zeros((len(self), self.dimension))
 
         # "seed" the logical operators in the GL sector
-        if not self.is_subsystem_code:
-            logicals_xx[cols_lz] = self.field.Identity(self.dimension)
-            logicals_zz[cols_lx] = self.field.Identity(self.dimension)
-
-        else:
-            cols_gl = np.sort(_join_slices(cols_gx, cols_lx))  # indices for all GL columns
-            """
-            Focusing on the gauge-qudit rows (i.e., constraints) of the parity check matrix, define
-                A = matrix_z[rows_gz, cols_gl],
-                B = matrix_x[rows_gx, cols_gl],
-            and denote the logical operator components in the GL sector by
-                U = logicals_xx[cols_gl],
-                V = logicals_zz[cols_gl].
-            These components need to satisfy the system of matrix equations
-                (1) A @ U.T = 0,
-                (2) B @ V.T = 0,
-                (3) U.T @ V = I.
-            Without loss of generality, we can satisfy (1) and (2) by setting
-                U = null_space(A).T
-                V = null_space(B).T @ M,
-            where the matrix M is determined by subsituting U and V back into (3),
-                U.T @ W @ M = I.
-            """
-            mat_U = matrix_z[rows_gz, cols_gl].view(self.field).null_space().T
-            mat_W = matrix_x[rows_gx, cols_gl].view(self.field).null_space().T
-            mat_M = np.linalg.inv(mat_U.T @ mat_W)
-            logicals_xx[cols_gl] = mat_U
-            logicals_zz[cols_gl] = mat_W @ mat_M
+        logicals_xx[cols_lz] = self.field.Identity(self.dimension)
+        logicals_zz[cols_lx] = self.field.Identity(self.dimension)
 
         # fill in remaining entries by enforcing parity check constraints
         logicals_xx[cols_sz] = -matrix_z[rows_sz] @ logicals_xx
@@ -1249,7 +1235,6 @@ def get_logical_ops(
         # Z support of X-type logicals, as column vectors
         logicals_xz = self.field.Zeros((len(self), self.dimension))
         logicals_xz[cols_lx] = self.field.Identity(self.dimension)
-        logicals_xz[cols_gx] = -matrix_x[rows_gx] @ logicals_xz + matrix_z[rows_gx] @ logicals_xx
         logicals_xz[cols_sx] = -matrix_x[rows_sx] @ logicals_xz + matrix_z[rows_sx] @ logicals_xx
 
         # full X and Z logicals as row vectors
@@ -2406,7 +2391,7 @@ def from_qecdb_id(code_id: str) -> CSSCode:
     def is_subsystem_code(self) -> bool:
         """Is this code a subsystem code?
 
-        That is, do all parity checks commute?
+        That is, do some parity checks fail to commute with each other?
         """
         if self._is_subsystem_code is None:
             self._is_subsystem_code = bool(np.any(self.matrix_x @ self.matrix_z.T))
@@ -2597,7 +2582,22 @@ def get_logical_ops(
             logicals_z[cols_lx] = self.field.Identity(self.dimension)
 
         else:
-            # see QuditCode.get_logical_ops for an explanation of what's happening here
+            # Restrict to the gauge-qudit rows (constraints) and define
+            #   A = matrix_z[rows_gz, cols_gl] and
+            #   B = matrix_x[rows_gx, cols_gl].
+            # The GL-sector components of the logical operators,
+            #   U = logicals_x[cols_gl] and
+            #   V = logicals_z[cols_gl],
+            # must satisfy
+            #   (1) A @ U = 0,
+            #   (2) B @ V = 0, and
+            #   (3) U.T @ V = I.
+            # Setting
+            #   U = null_space(A).T and
+            #   V = null_space(B).T @ M for some M
+            # satisfies (1) and (2), and setting
+            #   M = inv(U.T @ W) with W = null_space(B).T
+            # satisfies (3), since U.T @ V = U.T @ W @ M = I.
             cols_gl = np.sort(_join_slices(cols_gx, cols_lx))
             mat_U = matrix_z[rows_gz, cols_gl].view(self.field).null_space().T
             mat_W = matrix_x[rows_gx, cols_gl].view(self.field).null_space().T
 
@@ -562,6 +562,126 @@ def test_qudit_ops(pytestconfig: pytest.Config) -> None:
     assert np.array_equal(code.get_stabilizer_ops(canonicalized=True), stabilizer_ops[:-1])
 
 
+def test_qudit_subsystem_logical_ops() -> None:
+    """Non-CSS subsystem codes get a valid symplectic basis of logical operators.
+
+    The base QuditCode.get_logical_ops extracts the logical operators of a subsystem code as a
+    symplectic basis of the gauge group's centralizer.  This test checks that basis is valid --
+    correct commutation relations, and commuting with the gauge group -- for codes whose parity
+    checks mix X-type and Z-type support.
+    """
+
+    def assert_valid_basis(code: codes.QuditCode) -> None:
+        assert code.is_subsystem_code
+        logical_ops = code.get_logical_ops()
+        assert len(logical_ops) == 2 * code.dimension
+        # a valid basis has the full symplectic Gram matrix [[0, I], [-I, 0]]: X-type logicals
+        # mutually commute, Z-type logicals mutually commute, and each logical anticommutes only
+        # with its dual
+        assert np.array_equal(
+            logical_ops @ math.symplectic_conjugate(logical_ops).T,
+            get_symplectic_form(code.dimension, code.field),
+        )
+        # logical operators commute with every gauge generator (they lie in the centralizer)
+        assert not np.any(code.matrix @ math.symplectic_conjugate(logical_ops).T)
+        # gauge operators (dual().get_logical_ops()) are likewise a valid symplectic basis
+        gauge_ops = code.get_gauge_ops()
+        assert np.array_equal(
+            gauge_ops @ math.symplectic_conjugate(gauge_ops).T,
+            get_symplectic_form(code.gauge_dimension, code.field),
+        )
+        # the basis round-trips through the commutation-relation validation of set_logical_ops
+        code.set_logical_ops(logical_ops)
+
+    # conjugating some qudits mixes a Bacon-Shor code's X/Z support into a non-CSS subsystem code.
+    # The shipped reproduction must yield a valid basis without error, including via get_gauge_ops
+    # and get_distance, which build the logical operators of the dual code.
+    shipped = codes.BaconShorCode(3).conjugated([1, 3, 5])
+    assert_valid_basis(shipped)
+    assert shipped.get_distance() == 3
+
+    # property test over several fields, for one and two logical qudits.  Away from GF(2), build
+    # the non-CSS codes with _local_fourier, as conjugated() is only symplectic over GF(2).
+    for field in [galois.GF(2), galois.GF(3), galois.GF(4)]:
+        bacon_shor = codes.BaconShorCode(3, field=field.order).matrix
+        single = _local_fourier(bacon_shor, [1, 3, 5])
+        double = _local_fourier(_direct_sum(bacon_shor, bacon_shor), [1, 3, 5, 10, 14])
+        for matrix, expected_dimension in [(single, 1), (double, 2)]:
+            code = codes.QuditCode(matrix)
+            # the parity checks genuinely mix X-type and Z-type support (the code is not CSS)
+            has_x = np.any(code.matrix[:, : len(code)], axis=1)
+            has_z = np.any(code.matrix[:, len(code) :], axis=1)
+            assert np.any(has_x & has_z)
+            assert code.dimension == expected_dimension
+            assert_valid_basis(code)
+
+    # a subsystem code with no logical qudits yields an empty basis rather than an error
+    five_qubit = codes.FiveQubitCode()
+    gauged = codes.QuditCode(np.vstack([five_qubit.matrix, five_qubit.get_logical_ops()]))
+    assert gauged.is_subsystem_code and gauged.dimension == 0
+    assert gauged.get_logical_ops().shape == (0, 2 * len(gauged))
+
+
+def _local_fourier(matrix: galois.FieldArray, qudits: Sequence[int]) -> galois.FieldArray:
+    """Apply the local symplectic Fourier map ``(x_i, z_i) -> (z_i, -x_i)`` on the given qudits.
+
+    This preserves every symplectic inner product, so it maps a valid code to a valid code over any
+    field.  Applied to some qudits of a CSS code, it can mix X-type and Z-type support to build a
+    non-CSS code -- the regime that exercises the base QuditCode.get_logical_ops construction (CSS
+    codes use their own override).
+    """
+    field = type(matrix)
+    num_qudits = matrix.shape[1] // 2
+    x_bits, z_bits = matrix[:, :num_qudits].copy(), matrix[:, num_qudits:].copy()
+    x_bits[:, qudits], z_bits[:, qudits] = z_bits[:, qudits].copy(), -x_bits[:, qudits]
+    return np.hstack([x_bits, z_bits]).view(field)
+
+
+def _direct_sum(matrix_a: galois.FieldArray, matrix_b: galois.FieldArray) -> galois.FieldArray:
+    """Combine two symplectic parity check matrices into a code acting on disjoint qudits.
+
+    The result encodes the logical qudits of both summands, giving a simple way to build subsystem
+    codes with more than one logical qudit.
+    """
+    field = type(matrix_a)
+    num_a, num_b = matrix_a.shape[1] // 2, matrix_b.shape[1] // 2
+    num_qudits = num_a + num_b
+    result = field.Zeros((len(matrix_a) + len(matrix_b), 2 * num_qudits))
+    result[: len(matrix_a), :num_a] = matrix_a[:, :num_a]
+    result[: len(matrix_a), num_qudits : num_qudits + num_a] = matrix_a[:, num_a:]
+    result[len(matrix_a) :, num_a:num_qudits] = matrix_b[:, :num_b]
+    result[len(matrix_a) :, num_qudits + num_a :] = matrix_b[:, num_b:]
+    return result
+
+
+def test_get_standard_form_data_subsystem() -> None:
+    """QuditCode.get_standard_form_data reduces a subsystem code's checks to standard form.
+
+    get_logical_ops does not route subsystem codes through get_standard_form_data, so this covers
+    its subsystem branch directly against the identity-block structure that the method documents.
+    """
+    code = codes.QuditCode(codes.BaconShorCode(3).matrix)
+    assert code.is_subsystem_code
+    matrix, qudit_locs, row_sectors, col_sectors = code.get_standard_form_data()
+    rows_sx, rows_gx, rows_sz, rows_gz = row_sectors
+    cols_sx, cols_gx, _cols_lx, cols_sz, cols_gz, _cols_lz = col_sectors
+
+    # each stabilizer/gauge pivot block is an identity matrix, as documented
+    for rows, pauli_index, cols in [
+        (rows_sx, 0, cols_sx),
+        (rows_gx, 0, cols_gx),
+        (rows_sz, 1, cols_sz),
+        (rows_gz, 1, cols_gz),
+    ]:
+        block = matrix[rows, pauli_index, cols]
+        assert block.shape[0] == block.shape[1]
+        assert np.array_equal(block, code.field.Identity(len(block)))
+
+    # undoing the qudit permutation recovers the canonicalized parity check matrix
+    reordered = matrix[:, :, np.argsort(qudit_locs)].reshape(-1, 2 * len(code)).view(code.field)
+    assert np.array_equal(reordered.row_space(), code.canonicalized.matrix)
+
+
 def test_set_logical_ops_single_type_support() -> None:
     """QuditCode.set_logical_ops_x/z accept width-n single-type support."""
     css_code = codes.SteaneCode()
 
@@ -375,3 +375,64 @@ def _sqrt(value: galois.FieldArray) -> galois.FieldArray:
     than a 0-dimensional scalar (which it rejects over some extension fields).
     """
     return np.sqrt(np.atleast_1d(value))[0]
+
+
+def symplectic_gram_schmidt(
+    vectors: galois.FieldArray, *, promise_full_rank: bool = False
+) -> tuple[galois.FieldArray, galois.FieldArray]:
+    """Reduce vectors to symplectic hyperbolic pairs and a symplectic radical.
+
+    The rows of ``vectors`` span a subspace V of ``GF(q)^(2n)`` equipped with the symplectic inner
+    product ``⟨a, b⟩_s = a @ symplectic_conjugate(b)`` (see symplectic_conjugate).  Return a pair
+    ``(hyperbolic, radical)``:
+
+    - ``hyperbolic`` has shape ``(2m, 2n)`` and holds ``m`` mutually orthogonal hyperbolic pairs.
+      Its rows are ordered ``[b_0, ..., b_{m-1}, c_0, ..., c_{m-1}]``, so that
+      ``hyperbolic @ symplectic_conjugate(hyperbolic).T`` is the block matrix ``[[0, I], [-I, 0]]``:
+      ``⟨b_i, c_j⟩_s = δ_ij`` and all other products vanish.
+    - ``radical`` spans the symplectic radical of V -- the vectors of V that are orthogonal to all
+      of V.  Its rows are isotropic and orthogonal to every row of ``hyperbolic`` and ``radical``.
+
+    Together the rows of ``hyperbolic`` and ``radical`` form a basis for V.  The rows of ``vectors``
+    may be linearly dependent; they are first reduced to a basis of V.  Pass promise_full_rank=True
+    to skip this reduction when the rows are already independent; passing it for dependent rows
+    leaves the dependent directions in ``radical`` as spurious (possibly zero) rows, though the
+    ``hyperbolic`` pairs stay correct.
+
+    Because the symplectic form is alternating, ``⟨v, v⟩_s = 0`` for every vector in every
+    characteristic, so -- unlike get_orthonormal_basis -- there is no unit-vector case: the
+    construction peels off one hyperbolic pair at a time and collects the leftover radical.
+    """
+    field = type(vectors)
+    dimension = vectors.shape[1]
+    if not promise_full_rank:
+        vectors = vectors.row_space()  # reduce to a basis, discarding linearly dependent rows
+    words = list(vectors)
+
+    firsts: list[galois.FieldArray] = []  # the b_j
+    partners: list[galois.FieldArray] = []  # the c_j, with ⟨b_j, c_j⟩_s = 1
+    radical: list[galois.FieldArray] = []
+    while words:
+        first = words.pop(0)
+        index = next(
+            (ii for ii, word in enumerate(words) if first @ symplectic_conjugate(word) != 0), None
+        )
+        if index is None:
+            # "first" is orthogonal to every remaining word (and, by prior projections, to the
+            # extracted pairs and radical), so it belongs to the symplectic radical
+            radical.append(first)
+            continue
+        partner = words.pop(index)
+        # rescale so that ⟨first, partner⟩_s = 1
+        partner = partner / (first @ symplectic_conjugate(partner))
+        # project the remaining words to be symplectically orthogonal to both "first" and "partner"
+        conj_first = symplectic_conjugate(first)
+        conj_partner = symplectic_conjugate(partner)
+        words = [
+            word - (word @ conj_partner) * first + (word @ conj_first) * partner for word in words
+        ]
+        firsts.append(first)
+        partners.append(partner)
+
+    hyperbolic = field(firsts + partners) if firsts else field.Zeros((0, dimension))
+    return hyperbolic, field(radical) if radical else field.Zeros((0, dimension))