Skip to content

Commit 6065fab

Browse files
Fe-r-ozFe-r-oz
authored
Enhancements to GF(2) Linear Algebra: in-place row echelon with pivots, nullspace, and basis for row space (#445)
1 parent c6bc6ee commit 6065fab

File tree

3 files changed

+232
-0
lines changed

3 files changed

+232
-0
lines changed

‎CHANGELOG.md

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,10 @@
55

66
# News
77

8+
## v0.9.16 - dev
9+
10+
- Enhancements to `GF(2)` Linear Algebra: unexported, experimental `gf2_row_echelon_with_pivots!`, `gf2_nullspace`, `gf2_rowspace_basis`.
11+
812
## v0.9.15 - 2024-12-22
913

1014
- `pftrajectories` now supports fast multiqubit measurements with `PauliMeasurement` in addition to the already supported single qubit measurements `sMX/Z/Y` and workarounds like `naive_syndrome_circuit`.

‎src/QuantumClifford.jl

Lines changed: 63 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1168,6 +1168,69 @@ function gf2_H_to_G(H)
11681168
G[:,invperm(sindx)]
11691169
end
11701170

1171+
"""Performs in-place Gaussian elimination on a binary matrix and returns
1172+
its *row echelon form*,*rank*, the *transformation matrix*, and the *pivot
1173+
columns*. The transformation matrix that converts the original matrix into
1174+
the row echelon form. The `full` parameter controls the extent of elimination:
1175+
if `true`, only rows below the pivot are affected; if `false`, both above and
1176+
below the pivot are eliminated."""
1177+
function gf2_row_echelon_with_pivots!(M::AbstractMatrix{Int}; full=false)
1178+
r, c = size(M)
1179+
N = Matrix{Int}(LinearAlgebra.I, r, r)
1180+
p = 1
1181+
pivots = Int[]
1182+
for col in 1:c
1183+
@inbounds for row in p:r
1184+
if M[row, col] == 1
1185+
if row != p
1186+
M[[row, p], :] .= M[[p, row], :]
1187+
N[[row, p], :] .= N[[p, row], :]
1188+
end
1189+
break
1190+
end
1191+
end
1192+
if M[p, col] == 1
1193+
if !full
1194+
elim_range = p+1:r
1195+
else
1196+
elim_range = 1:r
1197+
end
1198+
@simd for j in elim_range
1199+
@inbounds if j != p && M[j, col] == 1
1200+
M[j, :] .= (M[j, :] .+ M[p, :]) .% 2
1201+
N[j, :] .= (N[j, :] .+ N[p, :]) .% 2
1202+
end
1203+
end
1204+
p += 1
1205+
push!(pivots, col)
1206+
end
1207+
if p > r
1208+
break
1209+
end
1210+
end
1211+
rank = p - 1
1212+
return M, rank, N, pivots
1213+
end
1214+
1215+
"""The nullspace of a binary matrix."""
1216+
function gf2_nullspace(H::AbstractMatrix{Int})
1217+
m = size(H',1)
1218+
_, matrix_rank, transformation_matrix, _ = gf2_row_echelon_with_pivots!(copy(H)')
1219+
if m == matrix_rank
1220+
# By the rank-nullity theorem, if rank(M) = m, then nullity(M) = 0
1221+
return zeros(Bool, 1, m)
1222+
end
1223+
# Extract the nullspace from the transformation matrix
1224+
return transformation_matrix[matrix_rank+1:end, :]
1225+
end
1226+
1227+
"""The basis for the row space of the binary matrix."""
1228+
function gf2_rowspace_basis(H::AbstractMatrix{Int})
1229+
pivots = gf2_row_echelon_with_pivots!(copy(H)')[4]
1230+
# Extract the rows corresponding to the pivot columns
1231+
H[pivots,:]
1232+
end
1233+
11711234
##############################
11721235
# Error classes
11731236
##############################

‎test/test_row_echelon_with_pivots.jl

Lines changed: 165 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,165 @@
1+
@testitem "Row echelon with pivots" begin
2+
using Random
3+
using Nemo
4+
using Nemo: echelon_form, matrix, GF
5+
using QuantumClifford
6+
using QuantumClifford: gf2_row_echelon_with_pivots!, gf2_nullspace, gf2_rowspace_basis
7+
test_sizes = [1,2,10,63,64,65,127,128,129]
8+
9+
@testset "GF(2) row echelon form with transformation matrix, pivots etc." begin
10+
for n in test_sizes
11+
for rep in 1:10
12+
gf2_matrices = [rand(Bool, size, size) for size in test_sizes]
13+
for (i, mat) in enumerate(gf2_matrices)
14+
naive_echelon_form, _, transformation, _ = gf2_row_echelon_with_pivots!(Matrix{Int}(mat), full=true) # in-place
15+
# Check the correctness of the transformation matrix
16+
@test (transformation*mat) .%2 == naive_echelon_form
17+
# Check the correctness of Gaussian elimination
18+
@test naive_echelon_form == gf2_gausselim!(mat)
19+
# Consistency check with Nemo.jl's echelon_form
20+
nemo_mat = matrix(GF(2), Matrix{Int}(mat))
21+
@test echelon_form(nemo_mat) == matrix(GF(2), naive_echelon_form)
22+
end
23+
end
24+
end
25+
end
26+
27+
function is_in_nullspace(A, x)
28+
# Ensure x is the correct orientation
29+
if size(x, 1) != size(A, 2)
30+
x = transpose(x)
31+
end
32+
# Perform modulo 2 arithmetic: A * x must be zero mod 2
33+
if size(x, 2) == 1 # x is a single column vector
34+
result = A * x
35+
return all(result .% 2 .== 0) # Check if A * x = 0 mod 2
36+
else # x is a matrix, check each column vector
37+
for i in 1:size(x, 2)
38+
result = A * x[:, i] # Multiply A with the i-th column of x
39+
if !all(result .% 2 .== 0) # Check if A * column = 0 mod 2
40+
return false
41+
end
42+
end
43+
return true # All columns are in the null space mod 2
44+
end
45+
end
46+
47+
@testset "GF(2) nullspace of the binary matrix" begin
48+
for n in test_sizes
49+
for rep in 1:10
50+
gf2_matrices = [rand(Bool, size, size) for size in test_sizes]
51+
for (i, matrix) in enumerate(gf2_matrices)
52+
imat = Matrix{Int}(matrix)
53+
ns = gf2_nullspace(imat)
54+
@test is_in_nullspace(imat, ns)
55+
end
56+
end
57+
end
58+
end
59+
60+
@testset "Consistency check with ldpc" begin
61+
# sanity checks for comparison to https://github.com/quantumgizmos/ldpc
62+
# results compared with 'from ldpc.mod2 import nullspace, row_basis, row_echelon'
63+
# Consistency check 1
64+
H = [1 1 1; 1 1 1; 0 1 0]
65+
echelon_form, rank, transformation, pivots = gf2_row_echelon_with_pivots!(copy(H)) # in-place
66+
@test echelon_form == [1 1 1; 0 1 0; 0 0 0]
67+
@test rank == 2
68+
@test transformation == [1 0 0; 0 0 1; 1 1 0]
69+
@test pivots == [1, 2] # in python, it's [0, 1] due to zero-indexing
70+
@test mod.((transformation*copy(H)), 2) == echelon_form
71+
@test gf2_nullspace(copy(H)) == [1 0 1]
72+
@test gf2_rowspace_basis(copy(H)) == [1 1 1; 0 1 0]
73+
# Consistency check 2
74+
H = [0 0 0 1 1 1 1;
75+
0 1 1 0 0 1 1;
76+
1 0 1 0 1 0 1]
77+
echelon_form, rank, transformation, pivots = gf2_row_echelon_with_pivots!(copy(H)) # in-place
78+
@test echelon_form == [1 0 1 0 1 0 1;
79+
0 1 1 0 0 1 1;
80+
0 0 0 1 1 1 1]
81+
@test rank == 3
82+
@test transformation == [0 0 1;
83+
0 1 0;
84+
1 0 0]
85+
@test pivots == [1, 2, 4] # in python, it's [0, 1, 3] due to zero-indexing
86+
@test mod.((transformation*copy(H)), 2) == echelon_form
87+
@test gf2_nullspace(copy(H)) == [1 1 1 0 0 0 0;
88+
0 1 1 1 1 0 0;
89+
0 1 0 1 0 1 0;
90+
0 0 1 1 0 0 1]
91+
@test gf2_rowspace_basis(copy(H)) == [0 0 0 1 1 1 1;
92+
0 1 1 0 0 1 1;
93+
1 0 1 0 1 0 1]
94+
# Consistency check 3
95+
H = [1 1 0; 0 1 1; 1 0 1]
96+
echelon_form, rank, transformation, pivots = gf2_row_echelon_with_pivots!(copy(H)) # in-place
97+
@test echelon_form == [1 1 0;
98+
0 1 1;
99+
0 0 0]
100+
@test rank == 2
101+
@test transformation == [1 0 0;
102+
0 1 0;
103+
1 1 1]
104+
@test pivots == [1,2 ] # in python, it's [0, 1] due to zero-indexing
105+
@test mod.((transformation*copy(H)), 2) == echelon_form
106+
@test gf2_nullspace(copy(H)) == [1 1 1]
107+
@test gf2_rowspace_basis(copy(H)) == [1 1 0;
108+
0 1 1]
109+
# Consistency check 4
110+
H = [1 1 0; 0 1 0; 0 0 1]
111+
echelon_form, rank, transformation, pivots = gf2_row_echelon_with_pivots!(copy(H)) # in-place
112+
@test echelon_form == [1 1 0;
113+
0 1 0;
114+
0 0 1]
115+
@test rank == 3
116+
@test transformation == [1 0 0;
117+
0 1 0;
118+
0 0 1]
119+
@test pivots == [1, 2, 3] # in python, it's [0, 1, 2] due to zero-indexing
120+
@test mod.((transformation*copy(H)), 2) == echelon_form
121+
@test gf2_nullspace(copy(H)) == [0 0 0]
122+
@test gf2_rowspace_basis(copy(H)) == [1 1 0;
123+
0 1 0;
124+
0 0 1]
125+
# Consistency check 5
126+
H = [1 1 0; 0 1 0; 0 0 1; 0 1 1]
127+
echelon_form, rank, transformation, pivots = gf2_row_echelon_with_pivots!(copy(H)) # in-place
128+
@test echelon_form == [1 1 0;
129+
0 1 0;
130+
0 0 1;
131+
0 0 0]
132+
@test rank == 3
133+
@test transformation == [1 0 0 0;
134+
0 1 0 0;
135+
0 0 1 0;
136+
0 1 1 1]
137+
@test pivots == [1, 2, 3] # in python, it's [0, 1, 2] due to zero-indexing
138+
@test mod.((transformation*copy(H)), 2) == echelon_form
139+
@test gf2_nullspace(copy(H)) == [0 0 0]
140+
@test gf2_rowspace_basis(copy(H)) == [1 1 0;
141+
0 1 0;
142+
0 0 1]
143+
# Consistency check 6
144+
H = [0 0 0 1 1 1 1;
145+
0 1 1 0 0 1 1;
146+
1 0 1 0 1 0 1]
147+
echelon_form, rank, transformation, pivots = gf2_row_echelon_with_pivots!(copy(H)) # in-place
148+
@test echelon_form == [1 0 1 0 1 0 1;
149+
0 1 1 0 0 1 1;
150+
0 0 0 1 1 1 1]
151+
@test rank == 3
152+
@test transformation == [0 0 1;
153+
0 1 0;
154+
1 0 0]
155+
@test pivots == [1, 2, 4] # in python, it's [0, 1, 3] due to zero-indexing
156+
@test mod.((transformation*copy(H)), 2) == echelon_form
157+
@test gf2_nullspace(copy(H)) == [1 1 1 0 0 0 0;
158+
0 1 1 1 1 0 0;
159+
0 1 0 1 0 1 0;
160+
0 0 1 1 0 0 1]
161+
@test gf2_rowspace_basis(copy(H)) == [0 0 0 1 1 1 1;
162+
0 1 1 0 0 1 1;
163+
1 0 1 0 1 0 1]
164+
end
165+
end

0 commit comments

Comments
 (0)