From d4e225f7489ea1fd75bb34b72136c7552097351d Mon Sep 17 00:00:00 2001 From: "Michael A. Perlin" Date: Wed, 15 Jul 2026 17:56:51 -0400 Subject: [PATCH 01/17] find a SWEL operator basis --- src/qldpc/codes/common.py | 27 +++++++++++++ src/qldpc/codes/common_test.py | 32 +++++++++++++++ src/qldpc/math.py | 74 ++++++++++++++++++++++++++++++++++ src/qldpc/math_test.py | 29 +++++++++++++ 4 files changed, 162 insertions(+) diff --git a/src/qldpc/codes/common.py b/src/qldpc/codes/common.py index 45a8f4e73..850d70f40 100644 --- a/src/qldpc/codes/common.py +++ b/src/qldpc/codes/common.py @@ -963,6 +963,10 @@ def to_css(self) -> CSSCode: ) return code + def to_swel(self) -> CSSCode: + """Convert this code into a CSSCode with a SWEL logical operator basis.""" + return self.to_css().to_swel() + def get_syndrome_subgraphs(self, *, strategy: str = "smallest_last") -> tuple[nx.DiGraph, ...]: """Sequence of subgraphs of the Tanner graph that induces a syndrome extraction sequence. @@ -2211,6 +2215,29 @@ def is_swel(self) -> bool: ) ) + def to_swel(self) -> CSSCode: + """Return a copy of this code with a SWEL logical operator basis (see is_swel). + + Raise a ValueError if this code is not self-dual, or if it has no SWEL basis. A self-dual + CSS code has a SWEL basis if and only if it has an odd-weight logical operator + (arXiv:2503.19790, Theorem 1). + """ + if not self.is_self_dual: + raise ValueError("Only self-dual codes can be converted into a SWEL basis") + # A code is SWEL when its logical operator supports are orthonormal (L @ L.T = identity), + # so a SWEL basis is an orthonormal basis for the space of logical operator supports. + supports = math.get_orthonormal_basis(self.get_logical_ops(Pauli.Z)) + if supports is None: + raise ValueError("This self-dual code has no SWEL logical operator basis") + code = CSSCode(self.code_x, self.code_z, is_subsystem_code=self._is_subsystem_code) + code._distance = self._distance + code._distance_x = self._distance_x + code._distance_z = self._distance_z + code._stabilizer_ops = self._stabilizer_ops + code._gauge_ops = self._gauge_ops + code.set_logical_ops_xz(supports, supports) + return code + @functools.cached_property def canonicalized(self) -> CSSCode: """The same code with its parity matrices in reduced row echelon form.""" diff --git a/src/qldpc/codes/common_test.py b/src/qldpc/codes/common_test.py index 6e9e598c8..61b94b175 100644 --- a/src/qldpc/codes/common_test.py +++ b/src/qldpc/codes/common_test.py @@ -604,6 +604,38 @@ def test_css_code(pytestconfig: pytest.Config) -> None: assert not codes.SurfaceCode(3).is_swel +def test_to_swel() -> None: + """Convert codes into a SWEL logical operator basis (see CSSCode.is_swel).""" + # an already-SWEL code stays SWEL, with its parity checks unchanged + code = codes.SteaneCode() + swel_code = code.to_swel() + assert swel_code.is_swel + assert np.array_equal(swel_code.matrix_x, code.matrix_x) + assert np.array_equal(swel_code.matrix_z, code.matrix_z) + + # a self-dual code that is not SWEL in its default basis: the direct sum of a SWEL code + # (SteaneCode) and a self-dual non-SWEL code (a toric code) requires Lemma 2 of arXiv:2503.19790 + toric = codes.ToricCode(2) + matrix_x = math.block_matrix([[np.asarray(code.matrix_x, dtype=int), 0], [0, toric.matrix_x]]) + matrix_z = math.block_matrix([[np.asarray(code.matrix_z, dtype=int), 0], [0, toric.matrix_z]]) + sum_code = codes.CSSCode(matrix_x, matrix_z) + assert sum_code.is_self_dual and not sum_code.is_swel + assert sum_code.to_swel().is_swel + + # a non-self-dual code cannot be converted into a SWEL basis + with pytest.raises(ValueError, match="self-dual"): + codes.SurfaceCode(3).to_swel() + + # a self-dual code whose logical operators all have even weight has no SWEL basis + with pytest.raises(ValueError, match="no SWEL"): + codes.CSSCode([[1, 1, 1, 1]], [[1, 1, 1, 1]]).to_swel() + + # QuditCode.to_swel converts to a CSSCode first + assert codes.QuditCode(code.matrix).to_swel().is_swel + with pytest.raises(ValueError, match="both X and Z support"): + codes.FiveQubitCode().to_swel() + + def test_css_ops(pytestconfig: pytest.Config) -> None: """Logical and stabilizer operator construction for CSS codes.""" seed = pytestconfig.getoption("randomly_seed") diff --git a/src/qldpc/math.py b/src/qldpc/math.py index 44563b54f..1030d51b6 100644 --- a/src/qldpc/math.py +++ b/src/qldpc/math.py @@ -112,6 +112,80 @@ def get_dual_basis(basis: galois.FieldArray, *, validate: bool = True) -> galois return dual_basis.view(type(basis)) +def get_orthonormal_basis(vectors: galois.FieldArray) -> galois.FieldArray | None: + """An orthonormal basis for the row space of a matrix over the binary field GF(2). + + Given a matrix with linearly independent rows spanning a subspace V of GF(2)^n, return a matrix + L whose rows are a basis for V with L @ L.T = identity -- that is, an orthonormal basis for V + under the standard bilinear form (the dot product). If V has no orthonormal basis, return None. + + An orthonormal basis exists if and only if the bilinear form restricted to V is nondegenerate + (no nonzero vector of V is orthogonal to all of V) and non-alternating (some vector v of V has + v @ v = 1). Over GF(2) the alternating case is exactly the case in which every vector has even + weight. + + The construction is a variant of symplectic Gram-Schmidt orthogonalization (Algorithm 1 and + Lemma 2 of arXiv:2503.19790). The first loop reduces the row space to a collection of "unit" + vectors u with u @ u = 1 and "hyperbolic pairs" (b, c) with b @ b = c @ c = 0 and b @ c = 1, all + mutually orthogonal. Lemma 2 then rewrites one unit vector and one hyperbolic pair into three + unit vectors, and is applied repeatedly to eliminate every hyperbolic pair. + """ + field = type(vectors) + dimension = vectors.shape[1] + words = [row for row in vectors] + + units: list[galois.FieldArray] = [] # vectors u with u @ u = 1 + pairs: list[tuple[galois.FieldArray, galois.FieldArray]] = [] # (b, c) with b @ c = 1 + + while words: + self_overlaps = [int(word @ word) for word in words] + if 1 in self_overlaps: + # a unit vector; orthogonalize the other words against it + pp = self_overlaps.index(1) + unit = words[pp] + units.append(unit) + words = [ + word + unit if int(word @ unit) else word + for qq, word in enumerate(words) + if qq != pp + ] + else: + # a hyperbolic pair; orthogonalize the other words against both of its vectors + first = words[0] + overlaps = [int(first @ word) for word in words] + if 1 not in overlaps: + return None # "first" is orthogonal to all of V, so the form is degenerate + pp = overlaps.index(1) + partner = words[pp] + pairs.append((first, partner)) + new_words = [] + for qq, word in enumerate(words): + if qq == 0 or qq == pp: + continue + add_first = int(word @ partner) + add_partner = int(word @ first) + if add_first: + word = word + first + if add_partner: + word = word + partner + new_words.append(word) + words = new_words + + # an orthonormal basis requires a non-alternating form: at least one unit vector must exist + if not units and pairs: + return None + + # Lemma 2: rewrite a unit vector a and a hyperbolic pair (b, c) as three unit vectors + while pairs: + vec_b, vec_c = pairs.pop() + vec_a = units.pop() + units += [vec_a + vec_b + vec_c, vec_a + vec_b, vec_a + vec_c] + + if not units: + return field.Zeros((0, dimension)) + return field(np.array(units, dtype=int)) + + def block_matrix( blocks: Sequence[Sequence[npt.NDArray[np.generic] | int | object]], ) -> npt.NDArray[np.int_]: diff --git a/src/qldpc/math_test.py b/src/qldpc/math_test.py index bbff624b1..b2800e527 100644 --- a/src/qldpc/math_test.py +++ b/src/qldpc/math_test.py @@ -74,6 +74,35 @@ def test_dual_basis(pytestconfig: pytest.Config) -> None: qldpc.math.get_dual_basis(field.Random((2, 1))) +def test_orthonormal_basis() -> None: + """Orthonormal bases over GF(2) (arXiv:2503.19790, Algorithm 1 and Lemma 2).""" + field = galois.GF(2) + + # the empty subspace has an empty orthonormal basis + basis = qldpc.math.get_orthonormal_basis(field.Zeros((0, 4))) + assert basis is not None and np.array_equal(basis, field.Zeros((0, 4))) + + # a subspace that only needs Gram-Schmidt against unit vectors + vectors = field([[1, 0, 0], [1, 1, 0], [0, 0, 1]]) + basis = qldpc.math.get_orthonormal_basis(vectors) + assert basis is not None and np.array_equal(basis @ basis.T, field.Identity(3)) + + # a subspace that requires rewriting hyperbolic pairs into unit vectors (Lemma 2) + vectors = field([[1, 1, 1, 0, 0, 0], [0, 0, 0, 1, 1, 0], [0, 0, 0, 0, 1, 1]]) + basis = qldpc.math.get_orthonormal_basis(vectors) + assert basis is not None and np.array_equal(basis @ basis.T, field.Identity(3)) + + # an alternating form (every vector has even weight) admits no orthonormal basis + vectors = field( + [[1, 1, 0, 0, 0, 0], [1, 0, 1, 0, 0, 0], [0, 1, 1, 1, 1, 0], [0, 0, 0, 1, 0, 1]] + ) + assert qldpc.math.get_orthonormal_basis(vectors) is None + + # a degenerate form (a vector orthogonal to the whole subspace) admits no orthonormal basis + vectors = field([[1, 1, 0, 0], [0, 0, 1, 1]]) + assert qldpc.math.get_orthonormal_basis(vectors) is None + + def test_block_matrix() -> None: """block_matrix assembles a nested block structure into a single NumPy array.""" eye = np.eye(2, dtype=float) From e8bc989bcf33e7851fdc6a4d6221fa724f6e346a Mon Sep 17 00:00:00 2001 From: "Michael A. Perlin" Date: Wed, 15 Jul 2026 18:45:43 -0400 Subject: [PATCH 02/17] maybe_to_swel --- src/qldpc/codes/common.py | 51 ++++++++++++++++++++++++++++------ src/qldpc/codes/common_test.py | 27 +++++++++++++----- src/qldpc/math.py | 7 +++++ src/qldpc/math_test.py | 4 +++ 4 files changed, 73 insertions(+), 16 deletions(-) diff --git a/src/qldpc/codes/common.py b/src/qldpc/codes/common.py index 850d70f40..9a770cb87 100644 --- a/src/qldpc/codes/common.py +++ b/src/qldpc/codes/common.py @@ -950,7 +950,10 @@ def maybe_to_css(self) -> QuditCode: zs = np.any(matrix_z, axis=1) if np.any(xs & zs): return self - return CSSCode(matrix_x[xs], matrix_z[zs], is_subsystem_code=self._is_subsystem_code) + code = CSSCode(matrix_x[xs], matrix_z[zs], is_subsystem_code=self._is_subsystem_code) + code._dimension = self._dimension + code._distance = self._distance + return code def to_css(self) -> CSSCode: """Try to convert this QuditCode into a CSSCode. Throw an error if we fail.""" @@ -963,8 +966,23 @@ def to_css(self) -> CSSCode: ) return code + def maybe_to_swel(self) -> QuditCode: + """Try to convert this code into a CSSCode with a SWEL logical operator basis. + + Return self if we fail. See CSSCode.maybe_to_swel and CSSCode.is_swel. + """ + code = self.maybe_to_css() + if isinstance(code, CSSCode): + swel_code = code.maybe_to_swel() + if swel_code is not code: + return swel_code + return self + def to_swel(self) -> CSSCode: - """Convert this code into a CSSCode with a SWEL logical operator basis.""" + """Convert this code into a CSSCode with a SWEL logical operator basis. + + Throw an error if we fail. See CSSCode.to_swel and CSSCode.is_swel. + """ return self.to_css().to_swel() def get_syndrome_subgraphs(self, *, strategy: str = "smallest_last") -> tuple[nx.DiGraph, ...]: @@ -2215,21 +2233,22 @@ def is_swel(self) -> bool: ) ) - def to_swel(self) -> CSSCode: - """Return a copy of this code with a SWEL logical operator basis (see is_swel). + def maybe_to_swel(self) -> CSSCode: + """Try to put this code into a SWEL logical operator basis. Return self if we fail. - Raise a ValueError if this code is not self-dual, or if it has no SWEL basis. A self-dual - CSS code has a SWEL basis if and only if it has an odd-weight logical operator - (arXiv:2503.19790, Theorem 1). + A code can be put into a SWEL basis (see is_swel) if and only if it is self-dual and has an + odd-weight logical operator (arXiv:2503.19790, Theorem 1). """ if not self.is_self_dual: - raise ValueError("Only self-dual codes can be converted into a SWEL basis") + return self # A code is SWEL when its logical operator supports are orthonormal (L @ L.T = identity), # so a SWEL basis is an orthonormal basis for the space of logical operator supports. supports = math.get_orthonormal_basis(self.get_logical_ops(Pauli.Z)) if supports is None: - raise ValueError("This self-dual code has no SWEL logical operator basis") + return self + # only the logical operator basis changes, so carry over all other cached properties code = CSSCode(self.code_x, self.code_z, is_subsystem_code=self._is_subsystem_code) + code._dimension = self._dimension code._distance = self._distance code._distance_x = self._distance_x code._distance_z = self._distance_z @@ -2238,6 +2257,20 @@ def to_swel(self) -> CSSCode: code.set_logical_ops_xz(supports, supports) return code + def to_swel(self) -> CSSCode: + """Put this code into a SWEL logical operator basis. Throw an error if we fail. + + A code can be put into a SWEL basis (see is_swel) if and only if it is self-dual and has an + odd-weight logical operator (arXiv:2503.19790, Theorem 1). + """ + code = self.maybe_to_swel() + if code is self: + raise ValueError( + "Failed to put this code into a SWEL basis;" + " it must be self-dual with an odd-weight logical operator" + ) + return code + @functools.cached_property def canonicalized(self) -> CSSCode: """The same code with its parity matrices in reduced row echelon form.""" diff --git a/src/qldpc/codes/common_test.py b/src/qldpc/codes/common_test.py index 61b94b175..04a5ffb7b 100644 --- a/src/qldpc/codes/common_test.py +++ b/src/qldpc/codes/common_test.py @@ -622,19 +622,32 @@ def test_to_swel() -> None: assert sum_code.is_self_dual and not sum_code.is_swel assert sum_code.to_swel().is_swel - # a non-self-dual code cannot be converted into a SWEL basis - with pytest.raises(ValueError, match="self-dual"): - codes.SurfaceCode(3).to_swel() - - # a self-dual code whose logical operators all have even weight has no SWEL basis - with pytest.raises(ValueError, match="no SWEL"): - codes.CSSCode([[1, 1, 1, 1]], [[1, 1, 1, 1]]).to_swel() + # maybe_to_swel returns self when there is no SWEL basis: a non-self-dual code, ... + surface_code = codes.SurfaceCode(3) + assert surface_code.maybe_to_swel() is surface_code + # ... or a self-dual code whose logical operators all have even weight + even_code = codes.CSSCode([[1, 1, 1, 1]], [[1, 1, 1, 1]]) + assert even_code.maybe_to_swel() is even_code + + # to_swel raises in both of those cases + with pytest.raises(ValueError, match="SWEL basis"): + surface_code.to_swel() + with pytest.raises(ValueError, match="SWEL basis"): + even_code.to_swel() # QuditCode.to_swel converts to a CSSCode first assert codes.QuditCode(code.matrix).to_swel().is_swel with pytest.raises(ValueError, match="both X and Z support"): codes.FiveQubitCode().to_swel() + # QuditCode.maybe_to_swel returns self unless it yields a SWEL CSSCode + maybe_swel = codes.QuditCode(code.matrix).maybe_to_swel() # CSS and SWEL-able + assert isinstance(maybe_swel, codes.CSSCode) and maybe_swel.is_swel + non_swel_code = codes.QuditCode(toric.matrix) + assert non_swel_code.maybe_to_swel() is non_swel_code # CSS but not SWEL-able + five_qubit_code = codes.FiveQubitCode() + assert five_qubit_code.maybe_to_swel() is five_qubit_code # not even CSS + def test_css_ops(pytestconfig: pytest.Config) -> None: """Logical and stabilizer operator construction for CSS codes.""" diff --git a/src/qldpc/math.py b/src/qldpc/math.py index 1030d51b6..913173b7d 100644 --- a/src/qldpc/math.py +++ b/src/qldpc/math.py @@ -129,8 +129,15 @@ def get_orthonormal_basis(vectors: galois.FieldArray) -> galois.FieldArray | Non vectors u with u @ u = 1 and "hyperbolic pairs" (b, c) with b @ b = c @ c = 0 and b @ c = 1, all mutually orthogonal. Lemma 2 then rewrites one unit vector and one hyperbolic pair into three unit vectors, and is applied repeatedly to eliminate every hyperbolic pair. + + This method is only implemented over GF(2). Over a general field the "unit" and "hyperbolic + pair" steps would require rescaling by field coefficients, and over fields of odd characteristic + an orthonormal basis need not exist even for a nondegenerate non-alternating form (existence is + then governed by whether the discriminant of the form is a square). """ field = type(vectors) + if field is not galois.GF2: + raise ValueError("get_orthonormal_basis is only implemented over GF(2)") dimension = vectors.shape[1] words = [row for row in vectors] diff --git a/src/qldpc/math_test.py b/src/qldpc/math_test.py index b2800e527..1ae5baa9d 100644 --- a/src/qldpc/math_test.py +++ b/src/qldpc/math_test.py @@ -102,6 +102,10 @@ def test_orthonormal_basis() -> None: vectors = field([[1, 1, 0, 0], [0, 0, 1, 1]]) assert qldpc.math.get_orthonormal_basis(vectors) is None + # only GF(2) is supported + with pytest.raises(ValueError, match="only implemented over GF"): + qldpc.math.get_orthonormal_basis(galois.GF(3)([[1, 0], [0, 1]])) + def test_block_matrix() -> None: """block_matrix assembles a nested block structure into a single NumPy array.""" From 529d00dd4b43db3cc0682cbc277ec6b27599dfbb Mon Sep 17 00:00:00 2001 From: "Michael A. Perlin" Date: Wed, 15 Jul 2026 20:07:53 -0400 Subject: [PATCH 03/17] more caching --- src/qldpc/codes/common.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/qldpc/codes/common.py b/src/qldpc/codes/common.py index 9a770cb87..ce62ee75c 100644 --- a/src/qldpc/codes/common.py +++ b/src/qldpc/codes/common.py @@ -1769,6 +1769,8 @@ def transform_ops(ops: galois.FieldArray) -> galois.FieldArray: code._stabilizer_ops = transform_ops(self.get_stabilizer_ops()) if self._gauge_ops is not None: code._gauge_ops = transform_ops(self.get_gauge_ops()) + code._dimension = self._dimension + code._distance = self._distance return code def conjugate(self) -> QuditCode: @@ -1817,6 +1819,7 @@ def transform_ops(ops: galois.FieldArray) -> galois.FieldArray: if self._gauge_ops is not None: code._gauge_ops = transform_ops(self.get_gauge_ops()) + code._dimension = self._dimension return code @staticmethod @@ -2993,6 +2996,10 @@ def conjugate(self) -> CSSCode: code._gauge_ops = scipy.linalg.block_diag( self.get_gauge_ops(Pauli.Z), self.get_gauge_ops(Pauli.X) ).view(self.field) + code._dimension = self._dimension + code._distance = self._distance + code._distance_x = self._distance_z + code._distance_z = self._distance_x return code def deformed( From ff4e2b154382127788df6e3b71a65cc12b0cafff Mon Sep 17 00:00:00 2001 From: "Michael A. Perlin" Date: Wed, 15 Jul 2026 20:13:08 -0400 Subject: [PATCH 04/17] minimal tests --- src/qldpc/codes/common_test.py | 47 ++++++++++------------------------ src/qldpc/math_test.py | 5 ---- 2 files changed, 14 insertions(+), 38 deletions(-) diff --git a/src/qldpc/codes/common_test.py b/src/qldpc/codes/common_test.py index 04a5ffb7b..d08f71802 100644 --- a/src/qldpc/codes/common_test.py +++ b/src/qldpc/codes/common_test.py @@ -605,48 +605,29 @@ def test_css_code(pytestconfig: pytest.Config) -> None: def test_to_swel() -> None: - """Convert codes into a SWEL logical operator basis (see CSSCode.is_swel).""" - # an already-SWEL code stays SWEL, with its parity checks unchanged - code = codes.SteaneCode() - swel_code = code.to_swel() - assert swel_code.is_swel - assert np.array_equal(swel_code.matrix_x, code.matrix_x) - assert np.array_equal(swel_code.matrix_z, code.matrix_z) - - # a self-dual code that is not SWEL in its default basis: the direct sum of a SWEL code - # (SteaneCode) and a self-dual non-SWEL code (a toric code) requires Lemma 2 of arXiv:2503.19790 - toric = codes.ToricCode(2) - matrix_x = math.block_matrix([[np.asarray(code.matrix_x, dtype=int), 0], [0, toric.matrix_x]]) - matrix_z = math.block_matrix([[np.asarray(code.matrix_z, dtype=int), 0], [0, toric.matrix_z]]) - sum_code = codes.CSSCode(matrix_x, matrix_z) - assert sum_code.is_self_dual and not sum_code.is_swel - assert sum_code.to_swel().is_swel + """Put codes into a SWEL logical operator basis (see CSSCode.is_swel).""" + # a SWEL-able CSS code can be put into a SWEL basis + steane_code = codes.SteaneCode() + assert not steane_code.is_swel + assert steane_code.to_swel().is_swel # maybe_to_swel returns self when there is no SWEL basis: a non-self-dual code, ... surface_code = codes.SurfaceCode(3) assert surface_code.maybe_to_swel() is surface_code + with pytest.raises(ValueError, match="SWEL basis"): + surface_code.to_swel() # ... or a self-dual code whose logical operators all have even weight even_code = codes.CSSCode([[1, 1, 1, 1]], [[1, 1, 1, 1]]) assert even_code.maybe_to_swel() is even_code - # to_swel raises in both of those cases - with pytest.raises(ValueError, match="SWEL basis"): - surface_code.to_swel() - with pytest.raises(ValueError, match="SWEL basis"): - even_code.to_swel() - - # QuditCode.to_swel converts to a CSSCode first - assert codes.QuditCode(code.matrix).to_swel().is_swel - with pytest.raises(ValueError, match="both X and Z support"): - codes.FiveQubitCode().to_swel() - - # QuditCode.maybe_to_swel returns self unless it yields a SWEL CSSCode - maybe_swel = codes.QuditCode(code.matrix).maybe_to_swel() # CSS and SWEL-able + # QuditCode.to_swel and maybe_to_swel convert to a CSSCode first + assert codes.QuditCode(codes.SteaneCode().matrix).to_swel().is_swel + maybe_swel = codes.QuditCode(codes.SteaneCode().matrix).maybe_to_swel() # CSS and SWEL-able assert isinstance(maybe_swel, codes.CSSCode) and maybe_swel.is_swel - non_swel_code = codes.QuditCode(toric.matrix) - assert non_swel_code.maybe_to_swel() is non_swel_code # CSS but not SWEL-able - five_qubit_code = codes.FiveQubitCode() - assert five_qubit_code.maybe_to_swel() is five_qubit_code # not even CSS + non_swel_code = codes.QuditCode(codes.ToricCode(2).matrix) # CSS but not SWEL-able + assert non_swel_code.maybe_to_swel() is non_swel_code + five_qubit_code = codes.FiveQubitCode() # not even CSS + assert five_qubit_code.maybe_to_swel() is five_qubit_code def test_css_ops(pytestconfig: pytest.Config) -> None: diff --git a/src/qldpc/math_test.py b/src/qldpc/math_test.py index 1ae5baa9d..117f2c479 100644 --- a/src/qldpc/math_test.py +++ b/src/qldpc/math_test.py @@ -82,11 +82,6 @@ def test_orthonormal_basis() -> None: basis = qldpc.math.get_orthonormal_basis(field.Zeros((0, 4))) assert basis is not None and np.array_equal(basis, field.Zeros((0, 4))) - # a subspace that only needs Gram-Schmidt against unit vectors - vectors = field([[1, 0, 0], [1, 1, 0], [0, 0, 1]]) - basis = qldpc.math.get_orthonormal_basis(vectors) - assert basis is not None and np.array_equal(basis @ basis.T, field.Identity(3)) - # a subspace that requires rewriting hyperbolic pairs into unit vectors (Lemma 2) vectors = field([[1, 1, 1, 0, 0, 0], [0, 0, 0, 1, 1, 0], [0, 0, 0, 0, 1, 1]]) basis = qldpc.math.get_orthonormal_basis(vectors) From cbe17a6a55fdd086b418f5127c1381d2d133da4a Mon Sep 17 00:00:00 2001 From: "Michael A. Perlin" Date: Wed, 15 Jul 2026 20:49:41 -0400 Subject: [PATCH 05/17] change swel methods --- src/qldpc/codes/common.py | 133 ++++++++++++++------------------- src/qldpc/codes/common_test.py | 47 +++++------- 2 files changed, 75 insertions(+), 105 deletions(-) diff --git a/src/qldpc/codes/common.py b/src/qldpc/codes/common.py index ce62ee75c..00510ee99 100644 --- a/src/qldpc/codes/common.py +++ b/src/qldpc/codes/common.py @@ -34,6 +34,7 @@ import scipy.linalg import scipy.sparse import stim +from typing_extensions import Self from qldpc import abstract, decoders, external, math from qldpc.math import IntegerArray @@ -203,9 +204,10 @@ def graph_to_matrix(graph: nx.DiGraph) -> galois.FieldArray: def dimension(self) -> int: """The number of logical (qu)dits encoded by this code.""" - def forget_distance(self) -> None: + def forget_distance(self) -> Self: """Forget the known distance of this code.""" self._distance = None + return self ################################################################################ @@ -343,7 +345,7 @@ def generator(self) -> galois.FieldArray: self._generator = self.matrix.null_space() return self._generator - def set_generator(self, generator: npt.NDArray[np.int_] | Sequence[Sequence[int]]) -> None: + def set_generator(self, generator: npt.NDArray[np.int_] | Sequence[Sequence[int]]) -> Self: """Set the generator matrix of this code.""" generator = np.asanyarray(generator).view(self.field) if np.any(self.matrix @ generator.T): @@ -357,6 +359,7 @@ def set_generator(self, generator: npt.NDArray[np.int_] | Sequence[Sequence[int] f" {required_rank})" ) self._generator = generator + return self def iter_words(self, skip_zero: bool = False) -> Iterator[galois.FieldArray]: """Iterate over the code words of this code.""" @@ -966,25 +969,6 @@ def to_css(self) -> CSSCode: ) return code - def maybe_to_swel(self) -> QuditCode: - """Try to convert this code into a CSSCode with a SWEL logical operator basis. - - Return self if we fail. See CSSCode.maybe_to_swel and CSSCode.is_swel. - """ - code = self.maybe_to_css() - if isinstance(code, CSSCode): - swel_code = code.maybe_to_swel() - if swel_code is not code: - return swel_code - return self - - def to_swel(self) -> CSSCode: - """Convert this code into a CSSCode with a SWEL logical operator basis. - - Throw an error if we fail. See CSSCode.to_swel and CSSCode.is_swel. - """ - return self.to_css().to_swel() - def get_syndrome_subgraphs(self, *, strategy: str = "smallest_last") -> tuple[nx.DiGraph, ...]: """Sequence of subgraphs of the Tanner graph that induces a syndrome extraction sequence. @@ -1406,7 +1390,7 @@ def set_logical_ops( logical_ops: npt.NDArray[np.int_] | Sequence[Sequence[int]], *, skip_validation: bool = False, - ) -> None: + ) -> Self: """Set the logical operators of this code to the provided logical operators.""" logical_ops = np.asanyarray(logical_ops).view(self.field) if not skip_validation: @@ -1422,13 +1406,14 @@ def set_logical_ops( raise ValueError("An incorrect number of logical operators was provided") self._logical_ops = logical_ops self._dimension = len(logical_ops) // 2 + return self def set_logical_ops_x( self, logicals_ops_x: npt.NDArray[np.int_] | Sequence[Sequence[int]], *, skip_validation: bool = False, - ) -> None: + ) -> Self: """Set the X-type logical operators of this code. Determine suitable Z-type logical operators automatically. This choice is unique mod @@ -1453,7 +1438,7 @@ def set_logical_ops_x( old_logicals_z = self.get_logical_ops(Pauli.Z, symplectic=True) basis_change = np.linalg.inv(math.symplectic_conjugate(old_logicals_z) @ logicals_ops_x.T) new_logicals_z = basis_change @ old_logicals_z - self.set_logical_ops( + return self.set_logical_ops( np.vstack([logicals_ops_x, new_logicals_z]), skip_validation=skip_validation ) @@ -1462,7 +1447,7 @@ def set_logical_ops_z( logicals_ops_z: npt.NDArray[np.int_] | Sequence[Sequence[int]], *, skip_validation: bool = False, - ) -> None: + ) -> Self: """Set the Z-type logical operators of this code. Determine suitable X-type logical operators automatically. This choice is unique mod @@ -1487,7 +1472,7 @@ def set_logical_ops_z( old_logicals_x = self.get_logical_ops(Pauli.X, symplectic=True) basis_change = np.linalg.inv(old_logicals_x @ math.symplectic_conjugate(logicals_ops_z).T) new_logicals_x = basis_change @ old_logicals_x - self.set_logical_ops( + return self.set_logical_ops( np.vstack([new_logicals_x, logicals_ops_z]), skip_validation=skip_validation ) @@ -2221,58 +2206,43 @@ def is_self_dual(self) -> bool: @property def is_swel(self) -> bool: - """Is this code self-dual with equivalent logicals (SWEL)? - - Specifically, SWEL codes are self-dual codes for which Hadamard-transforming every logical Z - operator recovers a dual logical X operator. - - Whether a code is SWEL depends on the choice of logical operator basis. + """Can this code be self-dual with equivalent logicals (SWEL)? + + A SWEL basis is a choice of logical operators for which Hadamard-transforming every logical + Z operator recovers a dual logical X operator, equivalently for which the logical operator + supports are orthonormal (L @ L.T = identity). This is a property of the code, not of any + particular logical basis: a self-dual code admits a SWEL basis if and only if it has an + odd-weight logical operator (arXiv:2503.19790, Theorem 1). Use set_swel_logical_ops to + actually put the code into such a basis. """ - return bool( - self.is_self_dual - and np.array_equal( - (logs_z := self.get_logical_ops(Pauli.Z)) @ logs_z.T, - np.eye(self.dimension, dtype=int), - ) + return self.is_self_dual and ( + self.dimension == 0 + or any(int(logical_op @ logical_op) for logical_op in self.get_logical_ops(Pauli.Z)) ) - def maybe_to_swel(self) -> CSSCode: - """Try to put this code into a SWEL logical operator basis. Return self if we fail. + def get_swel_logical_ops(self) -> galois.FieldArray: + """Logical operator supports of a SWEL logical basis (see is_swel). - A code can be put into a SWEL basis (see is_swel) if and only if it is self-dual and has an - odd-weight logical operator (arXiv:2503.19790, Theorem 1). + The returned matrix L has L @ L.T = identity, and each row is usable as both the X-type and + Z-type support of a logical qubit. Raise a ValueError if the code has no SWEL basis. """ - if not self.is_self_dual: - return self - # A code is SWEL when its logical operator supports are orthonormal (L @ L.T = identity), - # so a SWEL basis is an orthonormal basis for the space of logical operator supports. - supports = math.get_orthonormal_basis(self.get_logical_ops(Pauli.Z)) + supports = ( + math.get_orthonormal_basis(self.get_logical_ops(Pauli.Z)) if self.is_self_dual else None + ) if supports is None: - return self - # only the logical operator basis changes, so carry over all other cached properties - code = CSSCode(self.code_x, self.code_z, is_subsystem_code=self._is_subsystem_code) - code._dimension = self._dimension - code._distance = self._distance - code._distance_x = self._distance_x - code._distance_z = self._distance_z - code._stabilizer_ops = self._stabilizer_ops - code._gauge_ops = self._gauge_ops - code.set_logical_ops_xz(supports, supports) - return code - - def to_swel(self) -> CSSCode: - """Put this code into a SWEL logical operator basis. Throw an error if we fail. - - A code can be put into a SWEL basis (see is_swel) if and only if it is self-dual and has an - odd-weight logical operator (arXiv:2503.19790, Theorem 1). - """ - code = self.maybe_to_swel() - if code is self: raise ValueError( - "Failed to put this code into a SWEL basis;" + "This code has no SWEL logical operator basis;" " it must be self-dual with an odd-weight logical operator" ) - return code + return supports + + def set_swel_logical_ops(self) -> Self: + """Put this code into a SWEL logical operator basis (see is_swel), and return self. + + Raise a ValueError if the code has no SWEL basis. + """ + supports = self.get_swel_logical_ops() + return self.set_logical_ops_xz(supports, supports) @functools.cached_property def canonicalized(self) -> CSSCode: @@ -2553,17 +2523,17 @@ def set_logical_ops_xz( logicals_ops_z: npt.NDArray[np.int_] | Sequence[Sequence[int]], *, skip_validation: bool = False, - ) -> None: + ) -> Self: """Set the logical operators of this code to the provided logical operators.""" logical_ops = scipy.linalg.block_diag(logicals_ops_x, logicals_ops_z) - self.set_logical_ops(logical_ops, skip_validation=skip_validation) + return self.set_logical_ops(logical_ops, skip_validation=skip_validation) def set_logical_ops_x( self, logicals_ops_x: npt.NDArray[np.int_] | Sequence[Sequence[int]], *, skip_validation: bool = False, - ) -> None: + ) -> Self: """Set the X-type logical operators of this code. Determine suitable Z-type logical operators automatically. This choice is unique mod @@ -2582,14 +2552,16 @@ def set_logical_ops_x( logicals_ops_x = np.asanyarray(logicals_ops_x).view(self.field) old_logicals_z = self.get_logical_ops(Pauli.Z) new_logicals_z = np.linalg.inv(old_logicals_z @ logicals_ops_x.T) @ old_logicals_z - self.set_logical_ops_xz(logicals_ops_x, new_logicals_z, skip_validation=skip_validation) + return self.set_logical_ops_xz( + logicals_ops_x, new_logicals_z, skip_validation=skip_validation + ) def set_logical_ops_z( self, logicals_ops_z: npt.NDArray[np.int_] | Sequence[Sequence[int]], *, skip_validation: bool = False, - ) -> None: + ) -> Self: """Set the Z-type logical operators of this code. Determine suitable X-type logical operators automatically. This choice is unique mod @@ -2608,7 +2580,9 @@ def set_logical_ops_z( logicals_ops_z = np.asanyarray(logicals_ops_z).view(self.field) old_logicals_x = self.get_logical_ops(Pauli.X) new_logicals_x = np.linalg.inv(old_logicals_x @ logicals_ops_z.T) @ old_logicals_x - self.set_logical_ops_xz(new_logicals_x, logicals_ops_z, skip_validation=skip_validation) + return self.set_logical_ops_xz( + new_logicals_x, logicals_ops_z, skip_validation=skip_validation + ) def get_stabilizer_ops( self, @@ -2926,11 +2900,12 @@ def get_distance_bound_with_decoder( return min_bound - def forget_distance(self) -> None: + def forget_distance(self) -> Self: """Forget the known distance of this code.""" self._distance_x = self._distance_z = self._distance = None + return self - def reduce_logical_op(self, pauli: PauliXZ, logical_index: int, **decoder_kwargs: Any) -> None: + def reduce_logical_op(self, pauli: PauliXZ, logical_index: int, **decoder_kwargs: Any) -> Self: """Reduce the weight of a logical operator. A minimal-weight logical operator is found by enforcing that it has a trivial syndrome, and @@ -2962,8 +2937,9 @@ def reduce_logical_op(self, pauli: PauliXZ, logical_index: int, **decoder_kwargs self._logical_ops.shape = (2, self.dimension, 2, len(self)) self._logical_ops[pauli, logical_index, pauli, :] = candidate_logical_op self._logical_ops.shape = (2 * self.dimension, 2 * len(self)) + return self - def reduce_logical_ops(self, pauli: PauliXZ | None = None, **decoder_kwargs: Any) -> None: + def reduce_logical_ops(self, pauli: PauliXZ | None = None, **decoder_kwargs: Any) -> Self: """Reduce the weight of all logical operators.""" assert pauli is None or pauli in PAULIS_XZ if pauli is None: @@ -2972,6 +2948,7 @@ def reduce_logical_ops(self, pauli: PauliXZ | None = None, **decoder_kwargs: Any else: for logical_index in range(self.dimension): self.reduce_logical_op(pauli, logical_index, **decoder_kwargs) + return self def conjugated(self, qudits: slice | Sequence[int]) -> QuditCode: """Apply local Fourier transforms, swapping X-type and Z-type operators. diff --git a/src/qldpc/codes/common_test.py b/src/qldpc/codes/common_test.py index d08f71802..b44fcc94b 100644 --- a/src/qldpc/codes/common_test.py +++ b/src/qldpc/codes/common_test.py @@ -599,35 +599,28 @@ def test_css_code(pytestconfig: pytest.Config) -> None: assert nx.utils.graphs_equal(subgraphs[0], code.get_graph(Pauli.X)) assert nx.utils.graphs_equal(subgraphs[1], code.get_graph(Pauli.Z)) - # self-dual codes with equivalent logicals (SWEL) - assert codes.SteaneCode().is_swel - assert not codes.SurfaceCode(3).is_swel - -def test_to_swel() -> None: - """Put codes into a SWEL logical operator basis (see CSSCode.is_swel).""" - # a SWEL-able CSS code can be put into a SWEL basis - steane_code = codes.SteaneCode() - assert not steane_code.is_swel - assert steane_code.to_swel().is_swel - - # maybe_to_swel returns self when there is no SWEL basis: a non-self-dual code, ... - surface_code = codes.SurfaceCode(3) - assert surface_code.maybe_to_swel() is surface_code - with pytest.raises(ValueError, match="SWEL basis"): - surface_code.to_swel() - # ... or a self-dual code whose logical operators all have even weight +def test_swel() -> None: + """Identify and construct SWEL logical operator bases (see CSSCode.is_swel).""" + # is_swel is a property of the code, independent of the choice of logical basis + assert codes.SteaneCode().is_swel # self-dual with an odd-weight logical operator + assert not codes.SurfaceCode(3).is_swel # not self-dual even_code = codes.CSSCode([[1, 1, 1, 1]], [[1, 1, 1, 1]]) - assert even_code.maybe_to_swel() is even_code - - # QuditCode.to_swel and maybe_to_swel convert to a CSSCode first - assert codes.QuditCode(codes.SteaneCode().matrix).to_swel().is_swel - maybe_swel = codes.QuditCode(codes.SteaneCode().matrix).maybe_to_swel() # CSS and SWEL-able - assert isinstance(maybe_swel, codes.CSSCode) and maybe_swel.is_swel - non_swel_code = codes.QuditCode(codes.ToricCode(2).matrix) # CSS but not SWEL-able - assert non_swel_code.maybe_to_swel() is non_swel_code - five_qubit_code = codes.FiveQubitCode() # not even CSS - assert five_qubit_code.maybe_to_swel() is five_qubit_code + assert not even_code.is_swel # self-dual, but every logical operator has even weight + + # get_swel_logical_ops returns an orthonormal logical basis, or raises if none exists + supports = codes.SteaneCode().get_swel_logical_ops() + assert np.array_equal(supports @ supports.T, np.eye(len(supports), dtype=int)) + with pytest.raises(ValueError, match="no SWEL"): + codes.SurfaceCode(3).get_swel_logical_ops() # not self-dual + with pytest.raises(ValueError, match="no SWEL"): + even_code.get_swel_logical_ops() # self-dual, but no odd-weight logical operator + + # set_swel_logical_ops installs such a basis in place and returns self + code = codes.SteaneCode() + assert code.set_swel_logical_ops() is code + logs_z = code.get_logical_ops(Pauli.Z) + assert np.array_equal(logs_z @ logs_z.T, np.eye(code.dimension, dtype=int)) def test_css_ops(pytestconfig: pytest.Config) -> None: From 94e6bdd285fde6b616d1278f2c0b0d075c2b17bd Mon Sep 17 00:00:00 2001 From: "Michael A. Perlin" Date: Wed, 15 Jul 2026 22:05:51 -0400 Subject: [PATCH 06/17] guard for GF2 --- src/qldpc/codes/common.py | 8 ++++++++ src/qldpc/codes/common_test.py | 7 +++++++ 2 files changed, 15 insertions(+) diff --git a/src/qldpc/codes/common.py b/src/qldpc/codes/common.py index 00510ee99..fc281575e 100644 --- a/src/qldpc/codes/common.py +++ b/src/qldpc/codes/common.py @@ -2214,7 +2214,11 @@ def is_swel(self) -> bool: particular logical basis: a self-dual code admits a SWEL basis if and only if it has an odd-weight logical operator (arXiv:2503.19790, Theorem 1). Use set_swel_logical_ops to actually put the code into such a basis. + + SWEL code support is currently only provided for qubit codes. """ + if self.field is not galois.GF2: + raise ValueError("SWEL code support is currently only provided for qubit codes") return self.is_self_dual and ( self.dimension == 0 or any(int(logical_op @ logical_op) for logical_op in self.get_logical_ops(Pauli.Z)) @@ -2225,7 +2229,11 @@ def get_swel_logical_ops(self) -> galois.FieldArray: The returned matrix L has L @ L.T = identity, and each row is usable as both the X-type and Z-type support of a logical qubit. Raise a ValueError if the code has no SWEL basis. + + SWEL code support is currently only provided for qubit codes. """ + if self.field is not galois.GF2: + raise ValueError("SWEL code support is currently only provided for qubit codes") supports = ( math.get_orthonormal_basis(self.get_logical_ops(Pauli.Z)) if self.is_self_dual else None ) diff --git a/src/qldpc/codes/common_test.py b/src/qldpc/codes/common_test.py index b44fcc94b..01b4c44b0 100644 --- a/src/qldpc/codes/common_test.py +++ b/src/qldpc/codes/common_test.py @@ -622,6 +622,13 @@ def test_swel() -> None: logs_z = code.get_logical_ops(Pauli.Z) assert np.array_equal(logs_z @ logs_z.T, np.eye(code.dimension, dtype=int)) + # SWEL is currently only supported for qubit codes + qutrit_code = codes.CSSCode([[1, 1, 1]], [[1, 1, 1]], field=3) + with pytest.raises(ValueError, match="only provided for qubit codes"): + qutrit_code.is_swel + with pytest.raises(ValueError, match="only provided for qubit codes"): + qutrit_code.get_swel_logical_ops() + def test_css_ops(pytestconfig: pytest.Config) -> None: """Logical and stabilizer operator construction for CSS codes.""" From 415fe5988c9aefaf25598534b6492d5775f33e3b Mon Sep 17 00:00:00 2001 From: "Michael A. Perlin" Date: Thu, 16 Jul 2026 09:48:31 -0400 Subject: [PATCH 07/17] generalize beyond GF2 --- src/qldpc/codes/common.py | 29 ++-- src/qldpc/codes/common_test.py | 38 +++-- src/qldpc/math.py | 289 +++++++++++++++++++++------------ src/qldpc/math_test.py | 53 +++--- 4 files changed, 249 insertions(+), 160 deletions(-) diff --git a/src/qldpc/codes/common.py b/src/qldpc/codes/common.py index fc281575e..3162d10b1 100644 --- a/src/qldpc/codes/common.py +++ b/src/qldpc/codes/common.py @@ -2208,39 +2208,30 @@ def is_self_dual(self) -> bool: def is_swel(self) -> bool: """Can this code be self-dual with equivalent logicals (SWEL)? - A SWEL basis is a choice of logical operators for which Hadamard-transforming every logical - Z operator recovers a dual logical X operator, equivalently for which the logical operator - supports are orthonormal (L @ L.T = identity). This is a property of the code, not of any - particular logical basis: a self-dual code admits a SWEL basis if and only if it has an - odd-weight logical operator (arXiv:2503.19790, Theorem 1). Use set_swel_logical_ops to - actually put the code into such a basis. - - SWEL code support is currently only provided for qubit codes. + A SWEL basis is a choice of logical operators whose supports are orthonormal + (L @ L.T = identity), equivalently one for which a transversal Fourier transform (a + Hadamard, over qubits) maps every logical Z operator to a dual logical X operator. Whether + a SWEL basis exists is a property of the code, not of any particular logical basis (see + arXiv:2503.19790, Theorem 1, for qubits). Use set_swel_logical_ops to put the code into + such a basis. """ - if self.field is not galois.GF2: - raise ValueError("SWEL code support is currently only provided for qubit codes") return self.is_self_dual and ( - self.dimension == 0 - or any(int(logical_op @ logical_op) for logical_op in self.get_logical_ops(Pauli.Z)) + math.get_orthonormal_basis(self.get_logical_ops(Pauli.Z)) is not None ) def get_swel_logical_ops(self) -> galois.FieldArray: """Logical operator supports of a SWEL logical basis (see is_swel). The returned matrix L has L @ L.T = identity, and each row is usable as both the X-type and - Z-type support of a logical qubit. Raise a ValueError if the code has no SWEL basis. - - SWEL code support is currently only provided for qubit codes. + Z-type support of a logical qudit. Raise a ValueError if the code has no SWEL basis. """ - if self.field is not galois.GF2: - raise ValueError("SWEL code support is currently only provided for qubit codes") supports = ( math.get_orthonormal_basis(self.get_logical_ops(Pauli.Z)) if self.is_self_dual else None ) if supports is None: raise ValueError( - "This code has no SWEL logical operator basis;" - " it must be self-dual with an odd-weight logical operator" + "This code has no SWEL logical operator basis; it must be self-dual, and its" + " logical operators must admit an orthonormal basis" ) return supports diff --git a/src/qldpc/codes/common_test.py b/src/qldpc/codes/common_test.py index 01b4c44b0..0baac448f 100644 --- a/src/qldpc/codes/common_test.py +++ b/src/qldpc/codes/common_test.py @@ -603,10 +603,12 @@ def test_css_code(pytestconfig: pytest.Config) -> None: def test_swel() -> None: """Identify and construct SWEL logical operator bases (see CSSCode.is_swel).""" # is_swel is a property of the code, independent of the choice of logical basis - assert codes.SteaneCode().is_swel # self-dual with an odd-weight logical operator + assert codes.SteaneCode().is_swel # self-dual with an orthonormalizable logical basis assert not codes.SurfaceCode(3).is_swel # not self-dual even_code = codes.CSSCode([[1, 1, 1, 1]], [[1, 1, 1, 1]]) - assert not even_code.is_swel # self-dual, but every logical operator has even weight + assert not even_code.is_swel # self-dual, but its logicals admit no orthonormal basis + # a self-dual code over an odd field can also fail to be SWEL (non-square discriminant) + assert not codes.CSSCode([[1, 1, 1]], [[1, 1, 1]], field=3).is_swel # get_swel_logical_ops returns an orthonormal logical basis, or raises if none exists supports = codes.SteaneCode().get_swel_logical_ops() @@ -614,20 +616,24 @@ def test_swel() -> None: with pytest.raises(ValueError, match="no SWEL"): codes.SurfaceCode(3).get_swel_logical_ops() # not self-dual with pytest.raises(ValueError, match="no SWEL"): - even_code.get_swel_logical_ops() # self-dual, but no odd-weight logical operator - - # set_swel_logical_ops installs such a basis in place and returns self - code = codes.SteaneCode() - assert code.set_swel_logical_ops() is code - logs_z = code.get_logical_ops(Pauli.Z) - assert np.array_equal(logs_z @ logs_z.T, np.eye(code.dimension, dtype=int)) - - # SWEL is currently only supported for qubit codes - qutrit_code = codes.CSSCode([[1, 1, 1]], [[1, 1, 1]], field=3) - with pytest.raises(ValueError, match="only provided for qubit codes"): - qutrit_code.is_swel - with pytest.raises(ValueError, match="only provided for qubit codes"): - qutrit_code.get_swel_logical_ops() + even_code.get_swel_logical_ops() # self-dual, but no orthonormal logical basis + + # set_swel_logical_ops installs such a basis in place and returns self, over any finite field + # and for subsystem codes (which require a self-dual gauge group) + swel_codes = [ + codes.SteaneCode(), # a qubit code + codes.CSSCode( + [[0, 0, 1, 1, 1], [0, 1, 0, 1, 2]], [[0, 0, 1, 1, 1], [0, 1, 0, 1, 2]], field=3 + ), + codes.CSSCode([[0, 0, 1, 1]], [[0, 0, 1, 1]], field=4), + codes.CSSCode([[0, 0, 0, 0, 1], [0, 0, 0, 1, 0]], [[0, 0, 0, 0, 1], [0, 0, 0, 1, 0]]), + ] + assert swel_codes[-1].is_subsystem_code + for code in swel_codes: + assert code.is_swel + assert code.set_swel_logical_ops() is code + logs_z = code.get_logical_ops(Pauli.Z) + assert np.array_equal(logs_z @ logs_z.T, np.eye(code.dimension, dtype=int)) def test_css_ops(pytestconfig: pytest.Config) -> None: diff --git a/src/qldpc/math.py b/src/qldpc/math.py index 913173b7d..2d14018a3 100644 --- a/src/qldpc/math.py +++ b/src/qldpc/math.py @@ -35,6 +35,24 @@ DenseIntegerArrayType = TypeVar("DenseIntegerArrayType", galois.FieldArray, npt.NDArray[np.int_]) +#################################################################################################### +# general math + + +@functools.cache +def log_choose(n: int, k: int) -> float: + """Natural logarithm of (n choose k) = n! / ( k! * (n-k)! ).""" + return ( + scipy.special.gammaln(n + 1) + - scipy.special.gammaln(k + 1) + - scipy.special.gammaln(n - k + 1) + ) + + +#################################################################################################### +# manipulating Pauli strings and their symplectic vector representations + + def op_to_string(op: npt.NDArray[np.int_]) -> stim.PauliString: """Convert an integer array that represents a Pauli string into a stim.PauliString. @@ -80,14 +98,21 @@ def symplectic_weight(vectors: npt.NDArray[np.int_]) -> int: return np.count_nonzero(vectors_x | vectors_z, axis=-1).reshape(vectors.shape[:-1]) +#################################################################################################### +# matrix helper functions + + def first_nonzero_cols( matrix: npt.NDArray[np.generic] | Sequence[npt.NDArray[np.generic]], ) -> npt.NDArray[np.int_]: """Get the first nonzero column for every row in a matrix. If all columns are zero in a particular row, the column value is set to the number of columns. - If the array has more than two dimensions, return for each "row" r, the first "column" c for - which array[r, c] is not all zero. + + If the input ``matrix`` has more than two dimensions, return an ``output`` vector of length + ``matrix.shape[0]``, where ``output[r]`` is the first "column" ``c`` for which + ``np.any(matrix[r, c])`` is nonzero; (or set ``output[r] = matrix.shape[1]`` if + ``not np.any(matrix[r])``). """ _matrix = np.atleast_2d(np.asarray(matrix)) if _matrix.size == 0: @@ -99,6 +124,60 @@ def first_nonzero_cols( return first_nonzero_col_index.astype(np.int_, copy=False) +#################################################################################################### +# matrix helper functions + + +def block_matrix( + blocks: Sequence[Sequence[npt.NDArray[np.generic] | int | object]], +) -> npt.NDArray[np.int_]: + """Build a block matrix. + + Literal 0 entries are replaced by zero matrices, and literal 1 entries are replaced by an + identity matrix (padded below and to the right with zeros, if necessary). + """ + if not len(set(len(row) for row in blocks)) == 1: + raise ValueError("Inconsistent numbers of blocks in each row") + + # consistency checks + row_sizes = np.array( + [[bb.shape[0] if isinstance(bb, np.ndarray) else -1 for bb in row] for row in blocks] + ) + col_sizes = np.array( + [[bb.shape[1] if isinstance(bb, np.ndarray) else -1 for bb in row] for row in blocks] + ) + dtypes = [bb.dtype for row in blocks for bb in row if isinstance(bb, np.ndarray)] + if not all(len(set(row[row != -1])) == 1 for row in row_sizes): + raise ValueError("Inconsistent row numbers") + if not all(len(set(col[col != -1])) == 1 for col in col_sizes.T): + raise ValueError("Inconsistent column numbers") + if not len(set(dtypes)) == 1: + raise ValueError("Inconsistent block data types") + + # row numbers, column numbers, and data type + row_nums = [next(size for size in row if size != -1) for row in row_sizes] + col_nums = [next(size for size in col if size != -1) for col in col_sizes.T] + dtype = dtypes[0] + + # initialize a zero matrix and populate blocks + matrix = np.zeros((sum(row_nums), sum(col_nums)), dtype=dtype) + for rr, row in enumerate(blocks): + row_slice = slice(sum(row_nums[:rr]), sum(row_nums[: rr + 1])) + for cc, block in enumerate(row): + col_slice = slice(sum(col_nums[:cc]), sum(col_nums[: cc + 1])) + if not isinstance(block, int): + matrix[row_slice, col_slice] = block + elif block == 1: + matrix[row_slice, col_slice] = np.eye(row_nums[rr], col_nums[cc], dtype=dtype) + else: + assert block == 0, f"Unrecognized block: {block}" + return matrix + + +#################################################################################################### +# basis builders + + def get_dual_basis(basis: galois.FieldArray, *, validate: bool = True) -> galois.FieldArray: """Construct a dual basis, for which dual_basis @ basis.T = identity_matrix.""" if validate and ( @@ -113,70 +192,64 @@ def get_dual_basis(basis: galois.FieldArray, *, validate: bool = True) -> galois def get_orthonormal_basis(vectors: galois.FieldArray) -> galois.FieldArray | None: - """An orthonormal basis for the row space of a matrix over the binary field GF(2). + """An orthonormal basis for the row space of a matrix over a finite field. - Given a matrix with linearly independent rows spanning a subspace V of GF(2)^n, return a matrix + Given a matrix with linearly independent rows spanning a subspace V of GF(q)^n, return a matrix L whose rows are a basis for V with L @ L.T = identity -- that is, an orthonormal basis for V under the standard bilinear form (the dot product). If V has no orthonormal basis, return None. An orthonormal basis exists if and only if the bilinear form restricted to V is nondegenerate - (no nonzero vector of V is orthogonal to all of V) and non-alternating (some vector v of V has - v @ v = 1). Over GF(2) the alternating case is exactly the case in which every vector has even - weight. - - The construction is a variant of symplectic Gram-Schmidt orthogonalization (Algorithm 1 and - Lemma 2 of arXiv:2503.19790). The first loop reduces the row space to a collection of "unit" - vectors u with u @ u = 1 and "hyperbolic pairs" (b, c) with b @ b = c @ c = 0 and b @ c = 1, all - mutually orthogonal. Lemma 2 then rewrites one unit vector and one hyperbolic pair into three - unit vectors, and is applied repeatedly to eliminate every hyperbolic pair. - - This method is only implemented over GF(2). Over a general field the "unit" and "hyperbolic - pair" steps would require rescaling by field coefficients, and over fields of odd characteristic - an orthonormal basis need not exist even for a nondegenerate non-alternating form (existence is - then governed by whether the discriminant of the form is a square). + (no nonzero vector of V is orthogonal to all of V) and, in addition: + - over a field of characteristic 2, the form is non-alternating (some vector v has v @ v != 0); + - over a field of odd characteristic, the discriminant of the form is a square in GF(q). + + The construction is a variant of symplectic/orthogonal Gram-Schmidt orthogonalization; the + characteristic-2 case follows Algorithm 1 and Lemma 2 of arXiv:2503.19790. """ field = type(vectors) - if field is not galois.GF2: - raise ValueError("get_orthonormal_basis is only implemented over GF(2)") dimension = vectors.shape[1] words = [row for row in vectors] + units = ( + _orthonormalize_char_2(words) + if field.characteristic == 2 + else _orthonormalize_odd(words, field) + ) + if units is None: + return None + if not units: + return field.Zeros((0, dimension)) + return field(np.array(units, dtype=int)) + + +def _orthonormalize_char_2(words: list[galois.FieldArray]) -> list[galois.FieldArray] | None: + """Orthonormalize a spanning set over a field of characteristic 2, or None if impossible. + Reduce the row space to mutually orthogonal "unit" vectors u with u @ u = 1 and "hyperbolic + pairs" (b, c) with b @ b = c @ c = 0 and b @ c = 1 (Algorithm 1 of arXiv:2503.19790). Every + element of a characteristic-2 field is a square, so any vector with nonzero self-overlap can be + rescaled to a unit vector. Lemma 2 then rewrites one unit vector and one hyperbolic pair into + three unit vectors, eliminating every hyperbolic pair. + """ units: list[galois.FieldArray] = [] # vectors u with u @ u = 1 pairs: list[tuple[galois.FieldArray, galois.FieldArray]] = [] # (b, c) with b @ c = 1 - while words: - self_overlaps = [int(word @ word) for word in words] - if 1 in self_overlaps: - # a unit vector; orthogonalize the other words against it - pp = self_overlaps.index(1) - unit = words[pp] + index = next((ii for ii, word in enumerate(words) if word @ word), None) + if index is not None: + # a unit vector: rescale to self-overlap 1 and orthogonalize the other words against it + pivot = words.pop(index) + unit = pivot / np.sqrt(np.atleast_1d(pivot @ pivot))[0] + words = [word - (word @ unit) * unit for word in words] units.append(unit) - words = [ - word + unit if int(word @ unit) else word - for qq, word in enumerate(words) - if qq != pp - ] else: - # a hyperbolic pair; orthogonalize the other words against both of its vectors - first = words[0] - overlaps = [int(first @ word) for word in words] - if 1 not in overlaps: + # a hyperbolic pair: orthogonalize the other words against both of its vectors + first = words.pop(0) + index = next((ii for ii, word in enumerate(words) if first @ word), None) + if index is None: return None # "first" is orthogonal to all of V, so the form is degenerate - pp = overlaps.index(1) - partner = words[pp] + partner = words.pop(index) + partner = partner / (first @ partner) # rescale so that first @ partner == 1 + words = [word - (word @ partner) * first - (word @ first) * partner for word in words] pairs.append((first, partner)) - new_words = [] - for qq, word in enumerate(words): - if qq == 0 or qq == pp: - continue - add_first = int(word @ partner) - add_partner = int(word @ first) - if add_first: - word = word + first - if add_partner: - word = word + partner - new_words.append(word) - words = new_words # an orthonormal basis requires a non-alternating form: at least one unit vector must exist if not units and pairs: @@ -187,63 +260,79 @@ def get_orthonormal_basis(vectors: galois.FieldArray) -> galois.FieldArray | Non vec_b, vec_c = pairs.pop() vec_a = units.pop() units += [vec_a + vec_b + vec_c, vec_a + vec_b, vec_a + vec_c] + return units - if not units: - return field.Zeros((0, dimension)) - return field(np.array(units, dtype=int)) +def _orthonormalize_odd( + words: list[galois.FieldArray], field: type[galois.FieldArray] +) -> list[galois.FieldArray] | None: + """Orthonormalize a spanning set over a field of odd characteristic, or None if impossible. -def block_matrix( - blocks: Sequence[Sequence[npt.NDArray[np.generic] | int | object]], -) -> npt.NDArray[np.int_]: - """Build a block matrix. - - Literal 0 entries are replaced by zero matrices, and literal 1 entries are replaced by an - identity matrix (padded below and to the right with zeros, if necessary). + Ordinary Gram-Schmidt diagonalizes the form (in odd characteristic an alternating form is zero). + Each diagonal entry with a square self-overlap is rescaled to a unit vector; the remaining + non-square-norm vectors are paired off into unit vectors. This is possible if and only if their + number is even, i.e. if and only if the discriminant of the form is a square. """ - if not len(set(len(row) for row in blocks)) == 1: - raise ValueError("Inconsistent numbers of blocks in each row") + # diagonalize: build an orthogonal basis of vectors with nonzero self-overlap + diagonal: list[tuple[galois.FieldArray, galois.FieldArray]] = [] + while words: + index = next((ii for ii, word in enumerate(words) if word @ word), None) + if index is None: + # every self-overlap vanishes; a nonzero cross-overlap yields an anisotropic vector + index = _combine_anisotropic(words) + if index is None: + return None # the form vanishes on the remaining space, hence is degenerate + pivot = words.pop(index) + overlap = pivot @ pivot + words = [word - (word @ pivot) / overlap * pivot for word in words] + diagonal.append((pivot, overlap)) + + # rescale square-norm vectors to unit vectors; collect the non-square-norm vectors + units: list[galois.FieldArray] = [] + non_squares: list[tuple[galois.FieldArray, galois.FieldArray]] = [] + for pivot, overlap in diagonal: + if field.is_square(overlap): + units.append(pivot / np.sqrt(np.atleast_1d(overlap))[0]) + else: + non_squares.append((pivot, overlap)) - # consistency checks - row_sizes = np.array( - [[bb.shape[0] if isinstance(bb, np.ndarray) else -1 for bb in row] for row in blocks] - ) - col_sizes = np.array( - [[bb.shape[1] if isinstance(bb, np.ndarray) else -1 for bb in row] for row in blocks] - ) - dtypes = [bb.dtype for row in blocks for bb in row if isinstance(bb, np.ndarray)] - if not all(len(set(row[row != -1])) == 1 for row in row_sizes): - raise ValueError("Inconsistent row numbers") - if not all(len(set(col[col != -1])) == 1 for col in col_sizes.T): - raise ValueError("Inconsistent column numbers") - if not len(set(dtypes)) == 1: - raise ValueError("Inconsistent block data types") + # an odd number of non-square norms means a non-square discriminant: no orthonormal basis + if len(non_squares) % 2: + return None - # row numbers, column numbers, and data type - row_nums = [next(size for size in row if size != -1) for row in row_sizes] - col_nums = [next(size for size in col if size != -1) for col in col_sizes.T] - dtype = dtypes[0] + # pair non-square-norm vectors into unit vectors, using a fixed non-square element epsilon + if non_squares: + epsilon = next(elt for elt in field.elements if elt and not field.is_square(elt)) + # solve alpha**2 + beta**2 = 1 / epsilon, so that alpha u + beta w has self-overlap 1 + alpha, beta = _sum_of_two_squares(field, epsilon ** (field.order - 2)) + for (u_vec, u_norm), (w_vec, w_norm) in zip(non_squares[::2], non_squares[1::2]): + # rescale u_vec and w_vec to self-overlap epsilon, then combine into two unit vectors + u_vec = u_vec * np.sqrt(np.atleast_1d(epsilon / u_norm))[0] + w_vec = w_vec * np.sqrt(np.atleast_1d(epsilon / w_norm))[0] + units += [alpha * u_vec + beta * w_vec, -beta * u_vec + alpha * w_vec] + return units - # initialize a zero matrix and populate blocks - matrix = np.zeros((sum(row_nums), sum(col_nums)), dtype=dtype) - for rr, row in enumerate(blocks): - row_slice = slice(sum(row_nums[:rr]), sum(row_nums[: rr + 1])) - for cc, block in enumerate(row): - col_slice = slice(sum(col_nums[:cc]), sum(col_nums[: cc + 1])) - if not isinstance(block, int): - matrix[row_slice, col_slice] = block - elif block == 1: - matrix[row_slice, col_slice] = np.eye(row_nums[rr], col_nums[cc], dtype=dtype) - else: - assert block == 0, f"Unrecognized block: {block}" - return matrix +def _combine_anisotropic(words: list[galois.FieldArray]) -> int | None: + """Add one word to another to obtain a nonzero self-overlap; return its index, or None. -@functools.cache -def log_choose(n: int, k: int) -> float: - """Natural logarithm of (n choose k) = n! / ( k! * (n-k)! ).""" - return ( - scipy.special.gammaln(n + 1) - - scipy.special.gammaln(k + 1) - - scipy.special.gammaln(n - k + 1) - ) + Assumes every word has zero self-overlap. In odd characteristic, if some pair has nonzero + overlap then their sum has nonzero self-overlap; if no pair does, the form is identically zero. + """ + for ii in range(len(words)): + for jj in range(ii + 1, len(words)): + if words[ii] @ words[jj]: + words[ii] = words[ii] + words[jj] + return ii + return None + + +def _sum_of_two_squares( + field: type[galois.FieldArray], target: galois.FieldArray +) -> tuple[galois.FieldArray, galois.FieldArray]: + """Field elements (a, b) with a**2 + b**2 == target; these exist in odd characteristic.""" + for alpha in field.elements: + remainder = target - alpha * alpha + if field.is_square(remainder): + return alpha, np.sqrt(np.atleast_1d(remainder))[0] + raise AssertionError("no solution to a**2 + b**2 = target") # pragma: no cover diff --git a/src/qldpc/math_test.py b/src/qldpc/math_test.py index 117f2c479..a421ad6a2 100644 --- a/src/qldpc/math_test.py +++ b/src/qldpc/math_test.py @@ -75,31 +75,34 @@ def test_dual_basis(pytestconfig: pytest.Config) -> None: def test_orthonormal_basis() -> None: - """Orthonormal bases over GF(2) (arXiv:2503.19790, Algorithm 1 and Lemma 2).""" - field = galois.GF(2) - - # the empty subspace has an empty orthonormal basis - basis = qldpc.math.get_orthonormal_basis(field.Zeros((0, 4))) - assert basis is not None and np.array_equal(basis, field.Zeros((0, 4))) - - # a subspace that requires rewriting hyperbolic pairs into unit vectors (Lemma 2) - vectors = field([[1, 1, 1, 0, 0, 0], [0, 0, 0, 1, 1, 0], [0, 0, 0, 0, 1, 1]]) - basis = qldpc.math.get_orthonormal_basis(vectors) - assert basis is not None and np.array_equal(basis @ basis.T, field.Identity(3)) - - # an alternating form (every vector has even weight) admits no orthonormal basis - vectors = field( - [[1, 1, 0, 0, 0, 0], [1, 0, 1, 0, 0, 0], [0, 1, 1, 1, 1, 0], [0, 0, 0, 1, 0, 1]] - ) - assert qldpc.math.get_orthonormal_basis(vectors) is None - - # a degenerate form (a vector orthogonal to the whole subspace) admits no orthonormal basis - vectors = field([[1, 1, 0, 0], [0, 0, 1, 1]]) - assert qldpc.math.get_orthonormal_basis(vectors) is None - - # only GF(2) is supported - with pytest.raises(ValueError, match="only implemented over GF"): - qldpc.math.get_orthonormal_basis(galois.GF(3)([[1, 0], [0, 1]])) + """Orthonormal bases over finite fields (arXiv:2503.19790, Algorithm 1 and Lemma 2).""" + # subspaces that admit an orthonormal basis: check that L @ L.T is the identity + with_basis = [ + galois.GF(2).Zeros((0, 4)), # the empty subspace + galois.GF(2)([[1, 1, 1, 0, 0, 0], [0, 0, 0, 1, 1, 0], [0, 0, 0, 0, 1, 1]]), # Lemma 2 + galois.GF(4)([[2, 0]]), # characteristic 2, a rescaled unit vector + galois.GF(3)([[1, 0], [0, 1]]), # odd characteristic, unit norms + galois.GF(7)([[1, 1]]), # odd characteristic, a rescaled square norm + galois.GF(7)([[1, 2], [1, 3]]), # odd characteristic, two non-square norms paired off + ] + for vectors in with_basis: + field = type(vectors) + basis = qldpc.math.get_orthonormal_basis(vectors) + assert basis is not None and np.array_equal(basis @ basis.T, field.Identity(len(basis))) + + # subspaces with no orthonormal basis: get_orthonormal_basis returns None + without_basis = [ + # characteristic 2, an alternating form (every vector has even weight) + galois.GF(2)( + [[1, 1, 0, 0, 0, 0], [1, 0, 1, 0, 0, 0], [0, 1, 1, 1, 1, 0], [0, 0, 0, 1, 0, 1]] + ), + galois.GF(2)([[1, 1, 0, 0], [0, 0, 1, 1]]), # a degenerate form + galois.GF(7)([[1, 2]]), # odd characteristic, a non-square discriminant + galois.GF(7)([[1, 2, 3]]), # odd characteristic, a degenerate (isotropic) form + galois.GF(7)([[1, 2, 3], [1, 2, 4]]), # odd characteristic, isotropic then non-square + ] + for vectors in without_basis: + assert qldpc.math.get_orthonormal_basis(vectors) is None def test_block_matrix() -> None: From 1826dca35df5b97bf259c9f120c8d843f810216c Mon Sep 17 00:00:00 2001 From: "Michael A. Perlin" Date: Thu, 16 Jul 2026 10:17:00 -0400 Subject: [PATCH 08/17] minor fixes --- src/qldpc/math.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/qldpc/math.py b/src/qldpc/math.py index 2d14018a3..cc686ca64 100644 --- a/src/qldpc/math.py +++ b/src/qldpc/math.py @@ -302,7 +302,7 @@ def _orthonormalize_odd( # pair non-square-norm vectors into unit vectors, using a fixed non-square element epsilon if non_squares: - epsilon = next(elt for elt in field.elements if elt and not field.is_square(elt)) + epsilon = field.primitive_element # a generator of GF(q)* is always a non-square # solve alpha**2 + beta**2 = 1 / epsilon, so that alpha u + beta w has self-overlap 1 alpha, beta = _sum_of_two_squares(field, epsilon ** (field.order - 2)) for (u_vec, u_norm), (w_vec, w_norm) in zip(non_squares[::2], non_squares[1::2]): @@ -331,7 +331,8 @@ def _sum_of_two_squares( field: type[galois.FieldArray], target: galois.FieldArray ) -> tuple[galois.FieldArray, galois.FieldArray]: """Field elements (a, b) with a**2 + b**2 == target; these exist in odd characteristic.""" - for alpha in field.elements: + for value in range(field.order): # roughly half of all "alpha" succeed, so this exits quickly + alpha = field(value) remainder = target - alpha * alpha if field.is_square(remainder): return alpha, np.sqrt(np.atleast_1d(remainder))[0] From 0b39b71f1fad6121a28e6629c02ec2ad8c90c95d Mon Sep 17 00:00:00 2001 From: "Michael A. Perlin" Date: Fri, 17 Jul 2026 08:24:16 -0400 Subject: [PATCH 09/17] fix docstrings of QuditCode.*swel* --- src/qldpc/codes/common.py | 44 +++++++++++++++++++--------------- src/qldpc/codes/common_test.py | 4 +++- 2 files changed, 28 insertions(+), 20 deletions(-) diff --git a/src/qldpc/codes/common.py b/src/qldpc/codes/common.py index 3162d10b1..47c2b4418 100644 --- a/src/qldpc/codes/common.py +++ b/src/qldpc/codes/common.py @@ -2206,39 +2206,45 @@ def is_self_dual(self) -> bool: @property def is_swel(self) -> bool: - """Can this code be self-dual with equivalent logicals (SWEL)? - - A SWEL basis is a choice of logical operators whose supports are orthonormal - (L @ L.T = identity), equivalently one for which a transversal Fourier transform (a - Hadamard, over qubits) maps every logical Z operator to a dual logical X operator. Whether - a SWEL basis exists is a property of the code, not of any particular logical basis (see - arXiv:2503.19790, Theorem 1, for qubits). Use set_swel_logical_ops to put the code into - such a basis. + """Is this code self-dual with equivalent logicals (SWEL)? + + A SWEL code is a self-dual CSS code for which there exists a basis of logical operators + (Lx, Lz) with Lx = Lz; that is, Hadamard-transforming any single-qudit logical Pauli X (in + some logical Pauli basis) recovers the associated logical Pauli Z. + + In the case of stabilizer (non-subsystem) codes over a field with characteristic 2, (that + is, over GF(2**m), which includes qubits with m = 1), a self-dual CSS code is SWEL iff if it + has odd block length (Corollary 1 of arXiv:2503.19790), or if at least one row of Lx (in any + basis) has nonzero self-overlap (Theorem 1 of arXiv:2503.19790). + + This method first checks special cases covered by arXiv:2503.19790, and otherwise directly + checks the existence of a a logical operator basis (Lx, Lz) with Lx = Lz. """ - return self.is_self_dual and ( - math.get_orthonormal_basis(self.get_logical_ops(Pauli.Z)) is not None - ) + if not self.is_self_dual: + return False + if not self.is_subsystem_code and self.field.characteristic == 2: + return len(self) % 2 == 1 or any(op @ op for op in self.get_logical_ops(Pauli.X)) + return math.get_orthonormal_basis(self.get_logical_ops(Pauli.X)) is not None def get_swel_logical_ops(self) -> galois.FieldArray: - """Logical operator supports of a SWEL logical basis (see is_swel). + """Find a self-dual logical operator basis for this code: (Lx, Lz) = (L, L). Return L. - The returned matrix L has L @ L.T = identity, and each row is usable as both the X-type and - Z-type support of a logical qudit. Raise a ValueError if the code has no SWEL basis. + Raise a ValueError if no such basis exists (see QuditCode.is_swel). """ supports = ( - math.get_orthonormal_basis(self.get_logical_ops(Pauli.Z)) if self.is_self_dual else None + math.get_orthonormal_basis(self.get_logical_ops(Pauli.X)) if self.is_self_dual else None ) if supports is None: raise ValueError( - "This code has no SWEL logical operator basis; it must be self-dual, and its" + "This code has no self-dual logical operator basis; it must be self-dual, and its" " logical operators must admit an orthonormal basis" ) return supports def set_swel_logical_ops(self) -> Self: - """Put this code into a SWEL logical operator basis (see is_swel), and return self. + """Set the logical operators of this code to those of a self-dual logical operator basis. - Raise a ValueError if the code has no SWEL basis. + Raise a ValueError if this code is not SWEL (see QuditCode.is_swel). """ supports = self.get_swel_logical_ops() return self.set_logical_ops_xz(supports, supports) @@ -2631,7 +2637,7 @@ def get_gauge_ops( return gauge_ops.reshape(-1, 2, len(self))[:, pauli, :].view(self.field) def dual(self) -> CSSCode: - """Dual to this code, which swaps the roles of logical and gauge operators. + """The dual of this code, which swaps the roles of logical and gauge operators. See help(qldpc.codes.QuditCode.dual) for an explanation. """ diff --git a/src/qldpc/codes/common_test.py b/src/qldpc/codes/common_test.py index 0baac448f..ad7b484ce 100644 --- a/src/qldpc/codes/common_test.py +++ b/src/qldpc/codes/common_test.py @@ -603,7 +603,9 @@ def test_css_code(pytestconfig: pytest.Config) -> None: def test_swel() -> None: """Identify and construct SWEL logical operator bases (see CSSCode.is_swel).""" # is_swel is a property of the code, independent of the choice of logical basis - assert codes.SteaneCode().is_swel # self-dual with an orthonormalizable logical basis + assert ( + codes.SteaneCode().is_swel + ) # self-dual qubit stabilizer code of odd length is always SWEL assert not codes.SurfaceCode(3).is_swel # not self-dual even_code = codes.CSSCode([[1, 1, 1, 1]], [[1, 1, 1, 1]]) assert not even_code.is_swel # self-dual, but its logicals admit no orthonormal basis From 878a1f419cd93c89df4c8b058c6676d94cecfbc6 Mon Sep 17 00:00:00 2001 From: "Michael A. Perlin" Date: Fri, 17 Jul 2026 08:31:49 -0400 Subject: [PATCH 10/17] simplify SWEL tests --- src/qldpc/codes/common.py | 2 +- src/qldpc/codes/common_test.py | 48 ++++++++++------------------------ 2 files changed, 15 insertions(+), 35 deletions(-) diff --git a/src/qldpc/codes/common.py b/src/qldpc/codes/common.py index 47c2b4418..c29ea4795 100644 --- a/src/qldpc/codes/common.py +++ b/src/qldpc/codes/common.py @@ -2227,7 +2227,7 @@ def is_swel(self) -> bool: return math.get_orthonormal_basis(self.get_logical_ops(Pauli.X)) is not None def get_swel_logical_ops(self) -> galois.FieldArray: - """Find a self-dual logical operator basis for this code: (Lx, Lz) = (L, L). Return L. + """Find a self-dual basis of logical operators for this code: (Lx, Lz) = (L, L). Return L. Raise a ValueError if no such basis exists (see QuditCode.is_swel). """ diff --git a/src/qldpc/codes/common_test.py b/src/qldpc/codes/common_test.py index ad7b484ce..47fba8a3f 100644 --- a/src/qldpc/codes/common_test.py +++ b/src/qldpc/codes/common_test.py @@ -600,42 +600,22 @@ def test_css_code(pytestconfig: pytest.Config) -> None: assert nx.utils.graphs_equal(subgraphs[1], code.get_graph(Pauli.Z)) -def test_swel() -> None: +def test_swel_codes() -> None: """Identify and construct SWEL logical operator bases (see CSSCode.is_swel).""" - # is_swel is a property of the code, independent of the choice of logical basis - assert ( - codes.SteaneCode().is_swel - ) # self-dual qubit stabilizer code of odd length is always SWEL - assert not codes.SurfaceCode(3).is_swel # not self-dual + steane_code = codes.SteaneCode() + surface_code = codes.SurfaceCode(3) even_code = codes.CSSCode([[1, 1, 1, 1]], [[1, 1, 1, 1]]) - assert not even_code.is_swel # self-dual, but its logicals admit no orthonormal basis - # a self-dual code over an odd field can also fail to be SWEL (non-square discriminant) - assert not codes.CSSCode([[1, 1, 1]], [[1, 1, 1]], field=3).is_swel - - # get_swel_logical_ops returns an orthonormal logical basis, or raises if none exists - supports = codes.SteaneCode().get_swel_logical_ops() - assert np.array_equal(supports @ supports.T, np.eye(len(supports), dtype=int)) - with pytest.raises(ValueError, match="no SWEL"): - codes.SurfaceCode(3).get_swel_logical_ops() # not self-dual - with pytest.raises(ValueError, match="no SWEL"): - even_code.get_swel_logical_ops() # self-dual, but no orthonormal logical basis - - # set_swel_logical_ops installs such a basis in place and returns self, over any finite field - # and for subsystem codes (which require a self-dual gauge group) - swel_codes = [ - codes.SteaneCode(), # a qubit code - codes.CSSCode( - [[0, 0, 1, 1, 1], [0, 1, 0, 1, 2]], [[0, 0, 1, 1, 1], [0, 1, 0, 1, 2]], field=3 - ), - codes.CSSCode([[0, 0, 1, 1]], [[0, 0, 1, 1]], field=4), - codes.CSSCode([[0, 0, 0, 0, 1], [0, 0, 0, 1, 0]], [[0, 0, 0, 0, 1], [0, 0, 0, 1, 0]]), - ] - assert swel_codes[-1].is_subsystem_code - for code in swel_codes: - assert code.is_swel - assert code.set_swel_logical_ops() is code - logs_z = code.get_logical_ops(Pauli.Z) - assert np.array_equal(logs_z @ logs_z.T, np.eye(code.dimension, dtype=int)) + assert steane_code.is_swel + assert not surface_code.is_swel # not self-dual + assert not even_code.is_swel # self-dual but not SWEL + + # we automatically find a self-dual logical operator basis, if it exists + code = steane_code.set_swel_logical_ops() + assert np.array_equal(code.get_logical_ops(Pauli.X), code.get_logical_ops(Pauli.Z)) + with pytest.raises(ValueError, match="no self-dual logical operator basis"): + surface_code.get_swel_logical_ops() + with pytest.raises(ValueError, match="no self-dual logical operator basis"): + even_code.get_swel_logical_ops() def test_css_ops(pytestconfig: pytest.Config) -> None: From 9f5619a9969a279acb46b2e01b61e8fe5daa0b54 Mon Sep 17 00:00:00 2001 From: "Michael A. Perlin" Date: Fri, 17 Jul 2026 09:29:53 -0400 Subject: [PATCH 11/17] minor fixes and tests --- src/qldpc/codes/common.py | 9 ++++++--- src/qldpc/codes/common_test.py | 4 ++++ src/qldpc/math.py | 32 ++++++++++++++++++++------------ src/qldpc/math_test.py | 20 ++++++++++++-------- 4 files changed, 42 insertions(+), 23 deletions(-) diff --git a/src/qldpc/codes/common.py b/src/qldpc/codes/common.py index c29ea4795..519a01b22 100644 --- a/src/qldpc/codes/common.py +++ b/src/qldpc/codes/common.py @@ -2222,9 +2222,10 @@ def is_swel(self) -> bool: """ if not self.is_self_dual: return False + ops_x = self.get_logical_ops(Pauli.X) if not self.is_subsystem_code and self.field.characteristic == 2: - return len(self) % 2 == 1 or any(op @ op for op in self.get_logical_ops(Pauli.X)) - return math.get_orthonormal_basis(self.get_logical_ops(Pauli.X)) is not None + return len(self) % 2 == 1 or any(op @ op for op in ops_x) + return math.get_orthonormal_basis(ops_x, promise_full_rank=True) is not None def get_swel_logical_ops(self) -> galois.FieldArray: """Find a self-dual basis of logical operators for this code: (Lx, Lz) = (L, L). Return L. @@ -2232,7 +2233,9 @@ def get_swel_logical_ops(self) -> galois.FieldArray: Raise a ValueError if no such basis exists (see QuditCode.is_swel). """ supports = ( - math.get_orthonormal_basis(self.get_logical_ops(Pauli.X)) if self.is_self_dual else None + math.get_orthonormal_basis(self.get_logical_ops(Pauli.X), promise_full_rank=True) + if self.is_self_dual + else None ) if supports is None: raise ValueError( diff --git a/src/qldpc/codes/common_test.py b/src/qldpc/codes/common_test.py index 47fba8a3f..07ad813db 100644 --- a/src/qldpc/codes/common_test.py +++ b/src/qldpc/codes/common_test.py @@ -609,6 +609,10 @@ def test_swel_codes() -> None: assert not surface_code.is_swel # not self-dual assert not even_code.is_swel # self-dual but not SWEL + # a self-dual code over GF(3) reaches the general orthonormal-basis test (odd characteristic) + qutrit_matrix = [[0, 0, 1, 1, 1], [0, 1, 0, 1, 2]] + assert codes.CSSCode(qutrit_matrix, qutrit_matrix, field=3).is_swel + # we automatically find a self-dual logical operator basis, if it exists code = steane_code.set_swel_logical_ops() assert np.array_equal(code.get_logical_ops(Pauli.X), code.get_logical_ops(Pauli.Z)) diff --git a/src/qldpc/math.py b/src/qldpc/math.py index cc686ca64..c1a53de58 100644 --- a/src/qldpc/math.py +++ b/src/qldpc/math.py @@ -124,10 +124,6 @@ def first_nonzero_cols( return first_nonzero_col_index.astype(np.int_, copy=False) -#################################################################################################### -# matrix helper functions - - def block_matrix( blocks: Sequence[Sequence[npt.NDArray[np.generic] | int | object]], ) -> npt.NDArray[np.int_]: @@ -191,12 +187,17 @@ def get_dual_basis(basis: galois.FieldArray, *, validate: bool = True) -> galois return dual_basis.view(type(basis)) -def get_orthonormal_basis(vectors: galois.FieldArray) -> galois.FieldArray | None: +def get_orthonormal_basis( + vectors: galois.FieldArray, *, promise_full_rank: bool = False +) -> galois.FieldArray | None: """An orthonormal basis for the row space of a matrix over a finite field. - Given a matrix with linearly independent rows spanning a subspace V of GF(q)^n, return a matrix - L whose rows are a basis for V with L @ L.T = identity -- that is, an orthonormal basis for V - under the standard bilinear form (the dot product). If V has no orthonormal basis, return None. + Given a matrix whose rows span a subspace V of GF(q)^n, return a matrix L whose rows are a + basis for V with L @ L.T = identity -- that is, an orthonormal basis for V under the standard + bilinear form (the dot product). If V has no orthonormal basis, return None. + + The rows 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. An orthonormal basis exists if and only if the bilinear form restricted to V is nondegenerate (no nonzero vector of V is orthogonal to all of V) and, in addition: @@ -208,7 +209,11 @@ def get_orthonormal_basis(vectors: galois.FieldArray) -> galois.FieldArray | Non """ field = type(vectors) dimension = vectors.shape[1] - words = [row for row in vectors] + if not promise_full_rank: + # reduce to a basis of the row space, discarding any linearly dependent rows + vectors = vectors.row_reduce() + vectors = vectors[first_nonzero_cols(vectors) < dimension] + words = list(vectors) units = ( _orthonormalize_char_2(words) if field.characteristic == 2 @@ -268,8 +273,11 @@ def _orthonormalize_odd( ) -> list[galois.FieldArray] | None: """Orthonormalize a spanning set over a field of odd characteristic, or None if impossible. - Ordinary Gram-Schmidt diagonalizes the form (in odd characteristic an alternating form is zero). - Each diagonal entry with a square self-overlap is rescaled to a unit vector; the remaining + Gram-Schmidt diagonalizes the form, pivoting on a vector of nonzero self-overlap at each step. + If every remaining self-overlap vanishes but some cross-overlap does not, the sum of that pair + is a valid pivot (in odd characteristic (u + w) @ (u + w) = 2 * u @ w); if no cross-overlap is + nonzero either, the form vanishes on the remaining space, which is therefore degenerate. Each + diagonal entry with a square self-overlap is rescaled to a unit vector; the remaining non-square-norm vectors are paired off into unit vectors. This is possible if and only if their number is even, i.e. if and only if the discriminant of the form is a square. """ @@ -304,7 +312,7 @@ def _orthonormalize_odd( if non_squares: epsilon = field.primitive_element # a generator of GF(q)* is always a non-square # solve alpha**2 + beta**2 = 1 / epsilon, so that alpha u + beta w has self-overlap 1 - alpha, beta = _sum_of_two_squares(field, epsilon ** (field.order - 2)) + alpha, beta = _sum_of_two_squares(field, field(1) / epsilon) for (u_vec, u_norm), (w_vec, w_norm) in zip(non_squares[::2], non_squares[1::2]): # rescale u_vec and w_vec to self-overlap epsilon, then combine into two unit vectors u_vec = u_vec * np.sqrt(np.atleast_1d(epsilon / u_norm))[0] diff --git a/src/qldpc/math_test.py b/src/qldpc/math_test.py index a421ad6a2..2d54fc433 100644 --- a/src/qldpc/math_test.py +++ b/src/qldpc/math_test.py @@ -60,18 +60,15 @@ def test_nonzero_cols() -> None: assert np.array_equal(qldpc.math.first_nonzero_cols(tensor), [0]) -def test_dual_basis(pytestconfig: pytest.Config) -> None: +def test_dual_basis() -> None: """Construct dual bases.""" - np.random.seed(pytestconfig.getoption("randomly_seed")) - field = galois.GF(2) - basis = field.Random((4, 5)).row_reduce() - basis = basis[qldpc.math.first_nonzero_cols(basis) < basis.shape[1]] + basis = field([[1, 1, 0, 0, 1], [0, 1, 1, 0, 0], [1, 0, 0, 1, 1]]) dual_basis = qldpc.math.get_dual_basis(basis) assert np.array_equal(dual_basis @ basis.T, field.Identity(len(basis))) with pytest.raises(ValueError, match="wide matrices of full rank"): - qldpc.math.get_dual_basis(field.Random((2, 1))) + qldpc.math.get_dual_basis(field([[1], [0]])) def test_orthonormal_basis() -> None: @@ -85,9 +82,10 @@ def test_orthonormal_basis() -> None: galois.GF(7)([[1, 1]]), # odd characteristic, a rescaled square norm galois.GF(7)([[1, 2], [1, 3]]), # odd characteristic, two non-square norms paired off ] + # these inputs are full-rank bases crafted to exercise each branch of the orthonormalization for vectors in with_basis: field = type(vectors) - basis = qldpc.math.get_orthonormal_basis(vectors) + basis = qldpc.math.get_orthonormal_basis(vectors, promise_full_rank=True) assert basis is not None and np.array_equal(basis @ basis.T, field.Identity(len(basis))) # subspaces with no orthonormal basis: get_orthonormal_basis returns None @@ -102,7 +100,13 @@ def test_orthonormal_basis() -> None: galois.GF(7)([[1, 2, 3], [1, 2, 4]]), # odd characteristic, isotropic then non-square ] for vectors in without_basis: - assert qldpc.math.get_orthonormal_basis(vectors) is None + assert qldpc.math.get_orthonormal_basis(vectors, promise_full_rank=True) is None + + # by default (promise_full_rank=False) the rows are first reduced to a basis of their row space + vectors = galois.GF(3)([[1, 0], [0, 1], [1, 1]]) # three dependent rows spanning all of GF(3)^2 + basis = qldpc.math.get_orthonormal_basis(vectors) + assert basis is not None and basis.shape == (2, 2) + assert np.array_equal(basis @ basis.T, galois.GF(3).Identity(2)) def test_block_matrix() -> None: From c171c7b10ed871e8839e610d6bf02ca3db54f35c Mon Sep 17 00:00:00 2001 From: "Michael A. Perlin" Date: Fri, 17 Jul 2026 09:37:25 -0400 Subject: [PATCH 12/17] even more cleanup --- src/qldpc/math.py | 26 ++++++++++++++++---------- 1 file changed, 16 insertions(+), 10 deletions(-) diff --git a/src/qldpc/math.py b/src/qldpc/math.py index c1a53de58..3a4822916 100644 --- a/src/qldpc/math.py +++ b/src/qldpc/math.py @@ -210,9 +210,7 @@ def get_orthonormal_basis( field = type(vectors) dimension = vectors.shape[1] if not promise_full_rank: - # reduce to a basis of the row space, discarding any linearly dependent rows - vectors = vectors.row_reduce() - vectors = vectors[first_nonzero_cols(vectors) < dimension] + vectors = vectors.row_space() # reduce to a basis, discarding linearly dependent rows words = list(vectors) units = ( _orthonormalize_char_2(words) @@ -242,7 +240,7 @@ def _orthonormalize_char_2(words: list[galois.FieldArray]) -> list[galois.FieldA if index is not None: # a unit vector: rescale to self-overlap 1 and orthogonalize the other words against it pivot = words.pop(index) - unit = pivot / np.sqrt(np.atleast_1d(pivot @ pivot))[0] + unit = pivot / _sqrt(pivot @ pivot) words = [word - (word @ unit) * unit for word in words] units.append(unit) else: @@ -300,7 +298,7 @@ def _orthonormalize_odd( non_squares: list[tuple[galois.FieldArray, galois.FieldArray]] = [] for pivot, overlap in diagonal: if field.is_square(overlap): - units.append(pivot / np.sqrt(np.atleast_1d(overlap))[0]) + units.append(pivot / _sqrt(overlap)) else: non_squares.append((pivot, overlap)) @@ -315,8 +313,8 @@ def _orthonormalize_odd( alpha, beta = _sum_of_two_squares(field, field(1) / epsilon) for (u_vec, u_norm), (w_vec, w_norm) in zip(non_squares[::2], non_squares[1::2]): # rescale u_vec and w_vec to self-overlap epsilon, then combine into two unit vectors - u_vec = u_vec * np.sqrt(np.atleast_1d(epsilon / u_norm))[0] - w_vec = w_vec * np.sqrt(np.atleast_1d(epsilon / w_norm))[0] + u_vec = u_vec * _sqrt(epsilon / u_norm) + w_vec = w_vec * _sqrt(epsilon / w_norm) units += [alpha * u_vec + beta * w_vec, -beta * u_vec + alpha * w_vec] return units @@ -339,9 +337,17 @@ def _sum_of_two_squares( field: type[galois.FieldArray], target: galois.FieldArray ) -> tuple[galois.FieldArray, galois.FieldArray]: """Field elements (a, b) with a**2 + b**2 == target; these exist in odd characteristic.""" - for value in range(field.order): # roughly half of all "alpha" succeed, so this exits quickly - alpha = field(value) + for alpha in field.elements: # roughly half of all "alpha" succeed, so this exits quickly remainder = target - alpha * alpha if field.is_square(remainder): - return alpha, np.sqrt(np.atleast_1d(remainder))[0] + return alpha, _sqrt(remainder) raise AssertionError("no solution to a**2 + b**2 = target") # pragma: no cover + + +def _sqrt(value: galois.FieldArray) -> galois.FieldArray: + """A square root of a (0-dimensional) field element that is known to be a square. + + Delegates to np.sqrt, which computes finite-field square roots but must be given an array + rather than a 0-dimensional scalar (which it rejects over some extension fields). + """ + return np.sqrt(np.atleast_1d(value))[0] From 6cfc9b6f98976435389fd7c9f18a988ac36a87a6 Mon Sep 17 00:00:00 2001 From: "Michael A. Perlin" Date: Fri, 17 Jul 2026 09:45:11 -0400 Subject: [PATCH 13/17] row_space --- src/qldpc/abstract/rings.py | 5 ++--- src/qldpc/circuits/encoding.py | 2 +- src/qldpc/codes/common.py | 16 +++++++--------- 3 files changed, 10 insertions(+), 13 deletions(-) diff --git a/src/qldpc/abstract/rings.py b/src/qldpc/abstract/rings.py index 0869eb1f1..8a56f4cd0 100644 --- a/src/qldpc/abstract/rings.py +++ b/src/qldpc/abstract/rings.py @@ -1082,10 +1082,9 @@ def _get_block_howell_form(matrix: galois.FieldArray, *, right: bool = False) -> if right and size > 1: matrix = matrix.transpose(0, 1, 3, 2) - # row-reduce as an expanded 2-D matrix and remove all-zero rows + # row-reduce as an expanded 2-D matrix, keeping a basis of the row space (no all-zero rows) shape = (num_block_rows * size, num_block_cols * size) - matrix = matrix.transpose(0, 2, 1, 3).reshape(shape).view(field).row_reduce() - matrix = matrix[qldpc.math.first_nonzero_cols(matrix) < matrix.shape[1]].view(field) + matrix = matrix.transpose(0, 2, 1, 3).reshape(shape).view(field).row_space() if size > 1: # insert zero rows to shift pivots down so that they always lie on the diagonal of a block diff --git a/src/qldpc/circuits/encoding.py b/src/qldpc/circuits/encoding.py index 6dd1c885c..8d64dfcd7 100644 --- a/src/qldpc/circuits/encoding.py +++ b/src/qldpc/circuits/encoding.py @@ -226,7 +226,7 @@ def get_state_stabilizers( # identify stabilizers that are supported entirely on the data qubits of the code stabilizers = [] - for row in matrix.row_reduce(): + for row in matrix.row_space(): state_xs = row[cols_state_x] state_zs = row[cols_state_z] other_xs = row[cols_other_x] diff --git a/src/qldpc/codes/common.py b/src/qldpc/codes/common.py index 519a01b22..21fd647e1 100644 --- a/src/qldpc/codes/common.py +++ b/src/qldpc/codes/common.py @@ -267,10 +267,9 @@ def canonicalized(self) -> ClassicalCode: """The same code with its parity matrix in reduced row echelon form.""" if self._is_canonicalized: # pragma: no cover return self - matrix_rref = self.matrix.row_reduce() - matrix_rref = matrix_rref[np.any(matrix_rref, axis=1), :] - code = ClassicalCode(matrix_rref, self.field) - code._dimension = len(self) - len(matrix_rref) + matrix = self.matrix.row_space() + code = ClassicalCode(matrix, self.field) + code._dimension = len(self) - len(matrix) code._distance = self._distance code._is_canonicalized = True return code @@ -894,11 +893,10 @@ def canonicalized(self) -> QuditCode: """The same code with its parity matrix in reduced row echelon form.""" if self._is_canonicalized: # pragma: no cover return self - matrix_rref = self.matrix.row_reduce() - matrix_rref = matrix_rref[np.any(matrix_rref, axis=1), :] - code = QuditCode(matrix_rref, self.field, is_subsystem_code=self._is_subsystem_code) + matrix = self.matrix.row_space() + code = QuditCode(matrix, self.field, is_subsystem_code=self._is_subsystem_code) if not self._is_subsystem_code: - code._dimension = len(code) - len(matrix_rref) + code._dimension = len(code) - len(matrix) code._distance = self._distance code._stabilizer_ops = self._stabilizer_ops code._gauge_ops = self._gauge_ops @@ -1241,7 +1239,7 @@ def get_standard_form_data( parity check matrix (modulo elementary row operations and array reshaping) is matrix[:, :, np.argsort(qudit_locs)]. As a sanity check, the following test should pass: matrix_2d = matrix[:, :, np.argsort(qudit_locs)].reshape(-1, 2 * len(self)) - assert np.array_equal(matrix_2d.row_reduce(), self.canonicalized.matrix) + assert np.array_equal(matrix_2d.row_space(), self.canonicalized.matrix) Finally, this method also returns slices (index sets) for all row and column sectors, which enables selecting blocks of the parity check matrix with, say, matrix[rows_sx, 1, cols_lz]. From 7f8cf14ce5dd578023dbdde39f99379ff43c286b Mon Sep 17 00:00:00 2001 From: "Michael A. Perlin" Date: Fri, 17 Jul 2026 09:49:15 -0400 Subject: [PATCH 14/17] typos --- src/qldpc/codes/common.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/qldpc/codes/common.py b/src/qldpc/codes/common.py index 21fd647e1..7fdef879a 100644 --- a/src/qldpc/codes/common.py +++ b/src/qldpc/codes/common.py @@ -2210,13 +2210,13 @@ def is_swel(self) -> bool: (Lx, Lz) with Lx = Lz; that is, Hadamard-transforming any single-qudit logical Pauli X (in some logical Pauli basis) recovers the associated logical Pauli Z. - In the case of stabilizer (non-subsystem) codes over a field with characteristic 2, (that - is, over GF(2**m), which includes qubits with m = 1), a self-dual CSS code is SWEL iff if it - has odd block length (Corollary 1 of arXiv:2503.19790), or if at least one row of Lx (in any + In the case of stabilizer (non-subsystem) codes over a field with characteristic 2 (that + is, over GF(2**m), which includes qubits with m = 1), a self-dual CSS code is SWEL iff it + has odd block length (Corollary 1 of arXiv:2503.19790), or at least one row of Lx (in any basis) has nonzero self-overlap (Theorem 1 of arXiv:2503.19790). This method first checks special cases covered by arXiv:2503.19790, and otherwise directly - checks the existence of a a logical operator basis (Lx, Lz) with Lx = Lz. + checks the existence of a logical operator basis (Lx, Lz) with Lx = Lz. """ if not self.is_self_dual: return False From 598d82ac616430d3e85ec948cdef7c2d5ba2dc3d Mon Sep 17 00:00:00 2001 From: "Michael A. Perlin" Date: Fri, 17 Jul 2026 10:04:53 -0400 Subject: [PATCH 15/17] some wording cleanup --- src/qldpc/math.py | 30 ++++++++++++++---------------- 1 file changed, 14 insertions(+), 16 deletions(-) diff --git a/src/qldpc/math.py b/src/qldpc/math.py index 3a4822916..73e05cf78 100644 --- a/src/qldpc/math.py +++ b/src/qldpc/math.py @@ -188,7 +188,7 @@ def get_dual_basis(basis: galois.FieldArray, *, validate: bool = True) -> galois def get_orthonormal_basis( - vectors: galois.FieldArray, *, promise_full_rank: bool = False + matrix: galois.FieldArray, *, promise_full_rank: bool = False ) -> galois.FieldArray | None: """An orthonormal basis for the row space of a matrix over a finite field. @@ -196,7 +196,7 @@ def get_orthonormal_basis( basis for V with L @ L.T = identity -- that is, an orthonormal basis for V under the standard bilinear form (the dot product). If V has no orthonormal basis, return None. - The rows may be linearly dependent; they are first reduced to a basis of V. Pass + The rows of the matrix 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. An orthonormal basis exists if and only if the bilinear form restricted to V is nondegenerate @@ -207,11 +207,11 @@ def get_orthonormal_basis( The construction is a variant of symplectic/orthogonal Gram-Schmidt orthogonalization; the characteristic-2 case follows Algorithm 1 and Lemma 2 of arXiv:2503.19790. """ - field = type(vectors) - dimension = vectors.shape[1] + field = type(matrix) + dimension = matrix.shape[1] if not promise_full_rank: - vectors = vectors.row_space() # reduce to a basis, discarding linearly dependent rows - words = list(vectors) + matrix = matrix.row_space() # reduce to a basis, discarding linearly dependent rows + words = list(matrix) units = ( _orthonormalize_char_2(words) if field.characteristic == 2 @@ -219,9 +219,7 @@ def get_orthonormal_basis( ) if units is None: return None - if not units: - return field.Zeros((0, dimension)) - return field(np.array(units, dtype=int)) + return field(units) if units else field.Zeros((0, dimension)) def _orthonormalize_char_2(words: list[galois.FieldArray]) -> list[galois.FieldArray] | None: @@ -230,8 +228,8 @@ def _orthonormalize_char_2(words: list[galois.FieldArray]) -> list[galois.FieldA Reduce the row space to mutually orthogonal "unit" vectors u with u @ u = 1 and "hyperbolic pairs" (b, c) with b @ b = c @ c = 0 and b @ c = 1 (Algorithm 1 of arXiv:2503.19790). Every element of a characteristic-2 field is a square, so any vector with nonzero self-overlap can be - rescaled to a unit vector. Lemma 2 then rewrites one unit vector and one hyperbolic pair into - three unit vectors, eliminating every hyperbolic pair. + rescaled to a unit vector. Lemma 2 of arXiv:2503.19790 then rewrites one unit vector and one + hyperbolic pair into three unit vectors, eliminating every hyperbolic pair. """ units: list[galois.FieldArray] = [] # vectors u with u @ u = 1 pairs: list[tuple[galois.FieldArray, galois.FieldArray]] = [] # (b, c) with b @ c = 1 @@ -273,11 +271,11 @@ def _orthonormalize_odd( Gram-Schmidt diagonalizes the form, pivoting on a vector of nonzero self-overlap at each step. If every remaining self-overlap vanishes but some cross-overlap does not, the sum of that pair - is a valid pivot (in odd characteristic (u + w) @ (u + w) = 2 * u @ w); if no cross-overlap is - nonzero either, the form vanishes on the remaining space, which is therefore degenerate. Each - diagonal entry with a square self-overlap is rescaled to a unit vector; the remaining - non-square-norm vectors are paired off into unit vectors. This is possible if and only if their - number is even, i.e. if and only if the discriminant of the form is a square. + is a valid pivot: (u + w) @ (u + w) = 2 * u @ w. If no cross-overlap is nonzero either, the + form vanishes on the remaining space, which is therefore degenerate. Each diagonal entry with a + square self-overlap is rescaled to a unit vector; the remaining non-square-norm vectors are + paired off into unit vectors. This is possible if and only if their number is even, i.e. if and + only if the discriminant of the form is a square. """ # diagonalize: build an orthogonal basis of vectors with nonzero self-overlap diagonal: list[tuple[galois.FieldArray, galois.FieldArray]] = [] From 38b05971aaeb26d300d9b29a6745326abff39df0 Mon Sep 17 00:00:00 2001 From: "Michael A. Perlin" Date: Fri, 17 Jul 2026 10:09:27 -0400 Subject: [PATCH 16/17] iff --- src/qldpc/math.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/qldpc/math.py b/src/qldpc/math.py index 73e05cf78..6205623e4 100644 --- a/src/qldpc/math.py +++ b/src/qldpc/math.py @@ -274,8 +274,8 @@ def _orthonormalize_odd( is a valid pivot: (u + w) @ (u + w) = 2 * u @ w. If no cross-overlap is nonzero either, the form vanishes on the remaining space, which is therefore degenerate. Each diagonal entry with a square self-overlap is rescaled to a unit vector; the remaining non-square-norm vectors are - paired off into unit vectors. This is possible if and only if their number is even, i.e. if and - only if the discriminant of the form is a square. + paired off into unit vectors. This is possible if and only if the number of non-square-norm + vectors is even. """ # diagonalize: build an orthogonal basis of vectors with nonzero self-overlap diagonal: list[tuple[galois.FieldArray, galois.FieldArray]] = [] From 14a1be11754db395014f1ef97e74e2c1f580d1b9 Mon Sep 17 00:00:00 2001 From: "Michael A. Perlin" Date: Fri, 17 Jul 2026 10:18:31 -0400 Subject: [PATCH 17/17] better explanations --- src/qldpc/math.py | 65 +++++++++++++++++++++++++++-------------------- 1 file changed, 37 insertions(+), 28 deletions(-) diff --git a/src/qldpc/math.py b/src/qldpc/math.py index 6205623e4..2e5f75cb8 100644 --- a/src/qldpc/math.py +++ b/src/qldpc/math.py @@ -193,19 +193,21 @@ def get_orthonormal_basis( """An orthonormal basis for the row space of a matrix over a finite field. Given a matrix whose rows span a subspace V of GF(q)^n, return a matrix L whose rows are a - basis for V with L @ L.T = identity -- that is, an orthonormal basis for V under the standard - bilinear form (the dot product). If V has no orthonormal basis, return None. + basis for V with L @ L.T = identity -- that is, a basis of V whose vectors each have + self-overlap v @ v = 1 and are pairwise orthogonal under the dot product. If V has no such + basis, return None. The rows of the matrix 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. - An orthonormal basis exists if and only if the bilinear form restricted to V is nondegenerate - (no nonzero vector of V is orthogonal to all of V) and, in addition: - - over a field of characteristic 2, the form is non-alternating (some vector v has v @ v != 0); - - over a field of odd characteristic, the discriminant of the form is a square in GF(q). + An orthonormal basis for V exists if and only if (a) V is nondegenerate, meaning no nonzero + vector of V is orthogonal to all of V, and (b), according to the field's characteristic: + - over a field of characteristic 2, some vector of V has nonzero self-overlap v @ v; + - over a field of odd characteristic, an even number of the vectors in any mutually orthogonal + basis of V have a self-overlap with no square root in GF(q). - The construction is a variant of symplectic/orthogonal Gram-Schmidt orthogonalization; the - characteristic-2 case follows Algorithm 1 and Lemma 2 of arXiv:2503.19790. + The construction is a variant of Gram-Schmidt orthogonalization; the characteristic-2 case + follows Algorithm 1 and Lemma 2 of arXiv:2503.19790. """ field = type(matrix) dimension = matrix.shape[1] @@ -223,7 +225,7 @@ def get_orthonormal_basis( def _orthonormalize_char_2(words: list[galois.FieldArray]) -> list[galois.FieldArray] | None: - """Orthonormalize a spanning set over a field of characteristic 2, or None if impossible. + """Try to orthonormalize linearly independent vectors over a field of characteristic 2. Reduce the row space to mutually orthogonal "unit" vectors u with u @ u = 1 and "hyperbolic pairs" (b, c) with b @ b = c @ c = 0 and b @ c = 1 (Algorithm 1 of arXiv:2503.19790). Every @@ -236,23 +238,23 @@ def _orthonormalize_char_2(words: list[galois.FieldArray]) -> list[galois.FieldA while words: index = next((ii for ii, word in enumerate(words) if word @ word), None) if index is not None: - # a unit vector: rescale to self-overlap 1 and orthogonalize the other words against it + # extract a unit vector and orthogonalize other words against it pivot = words.pop(index) unit = pivot / _sqrt(pivot @ pivot) words = [word - (word @ unit) * unit for word in words] units.append(unit) else: - # a hyperbolic pair: orthogonalize the other words against both of its vectors + # extract a hyperbolic pair and orthogonalize other words against both of its vectors first = words.pop(0) index = next((ii for ii, word in enumerate(words) if first @ word), None) if index is None: - return None # "first" is orthogonal to all of V, so the form is degenerate + return None # "first" is orthogonal to all of V, so no orthonormal basis exists partner = words.pop(index) partner = partner / (first @ partner) # rescale so that first @ partner == 1 words = [word - (word @ partner) * first - (word @ first) * partner for word in words] pairs.append((first, partner)) - # an orthonormal basis requires a non-alternating form: at least one unit vector must exist + # an orthonormal basis needs at least one unit vector; pure hyperbolic pairs have none if not units and pairs: return None @@ -267,25 +269,32 @@ def _orthonormalize_char_2(words: list[galois.FieldArray]) -> list[galois.FieldA def _orthonormalize_odd( words: list[galois.FieldArray], field: type[galois.FieldArray] ) -> list[galois.FieldArray] | None: - """Orthonormalize a spanning set over a field of odd characteristic, or None if impossible. - - Gram-Schmidt diagonalizes the form, pivoting on a vector of nonzero self-overlap at each step. - If every remaining self-overlap vanishes but some cross-overlap does not, the sum of that pair - is a valid pivot: (u + w) @ (u + w) = 2 * u @ w. If no cross-overlap is nonzero either, the - form vanishes on the remaining space, which is therefore degenerate. Each diagonal entry with a - square self-overlap is rescaled to a unit vector; the remaining non-square-norm vectors are - paired off into unit vectors. This is possible if and only if the number of non-square-norm - vectors is even. + """Try to orthonormalize linearly independent vectors over a field of odd characteristic. + + First, find an orthogonal basis in which each vector has nonzero self-overlap (with Gram-Schmidt + orthogonalization): at each step, pick a vector with nonzero self-overlap v @ v and subtract + its projection from the others. Out of the remaining null vectors (with zero self-overlap), if + a pair (u, v) have nonzero overlap u @ v, then their sum has nonzero self-overlap (in a field + with odd characteristic), (u + v) @ (u + v) = 2 * u @ v, so replace u <- u + v, and + orthogonalize remaining vectors against u. If no remaining null vectors have nonzero overlap, + no orthonormal basis exists. + + A vector with nonzero self-overlap can be rescaled to a unit vector (self-overlap 1) exactly + when its self-overlap has a square root in the field, by taking v -> v / sqrt(v @ v). Vectors + with non-square self-overlaps cannot be rescaled on their own, but a pair of them may be + combined into two unit vectors, allowing an orthonormal basis when the number of such vectors + is even. """ # diagonalize: build an orthogonal basis of vectors with nonzero self-overlap diagonal: list[tuple[galois.FieldArray, galois.FieldArray]] = [] while words: index = next((ii for ii, word in enumerate(words) if word @ word), None) if index is None: - # every self-overlap vanishes; a nonzero cross-overlap yields an anisotropic vector - index = _combine_anisotropic(words) + # every self-overlap vanishes; sum an overlapping pair to get a nonzero self-overlap + index = _combine_null_vectors(words) if index is None: - return None # the form vanishes on the remaining space, hence is degenerate + return None # every overlap vanishes on the remaining space; no basis exists + # extract a vector with nonzero overlap, and otrhogonalize the remaining vectors against it pivot = words.pop(index) overlap = pivot @ pivot words = [word - (word @ pivot) / overlap * pivot for word in words] @@ -300,7 +309,7 @@ def _orthonormalize_odd( else: non_squares.append((pivot, overlap)) - # an odd number of non-square norms means a non-square discriminant: no orthonormal basis + # non-square self-overlaps pair off in twos, so an odd count means no orthonormal basis if len(non_squares) % 2: return None @@ -317,11 +326,11 @@ def _orthonormalize_odd( return units -def _combine_anisotropic(words: list[galois.FieldArray]) -> int | None: +def _combine_null_vectors(words: list[galois.FieldArray]) -> int | None: """Add one word to another to obtain a nonzero self-overlap; return its index, or None. Assumes every word has zero self-overlap. In odd characteristic, if some pair has nonzero - overlap then their sum has nonzero self-overlap; if no pair does, the form is identically zero. + overlap then their sum has nonzero self-overlap; if no pair does, every overlap is zero. """ for ii in range(len(words)): for jj in range(ii + 1, len(words)):