Skip to content

Commit

Permalink
Enhancements to GF(2) Linear Algebra: in-place row echelon with pivot…
Browse files Browse the repository at this point in the history 
…s, nullspace, and basis for row space
  • Loading branch information
@Fe-r-oz
Fe-r-oz committed Dec 4, 2024
1 parent 5fcf353 commit 7b3238b
Show file tree
Hide file tree
Showing 2 changed files with 214 additions and 0 deletions.
64 changes: 64 additions & 0 deletions src/QuantumClifford.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@ export
affectedqubits, #TODO move to QuantumInterface?
# GF2
stab_to_gf2, gf2_gausselim!, gf2_isinvertible, gf2_invert, gf2_H_to_G,
gf2_row_echelon_with_pivots!, gf2_nullspace, gf2_rowspace_basis,
# Canonicalization
canonicalize!, canonicalize_rref!, canonicalize_gott!,
# Linear Algebra
Expand Down Expand Up @@ -1168,6 +1169,69 @@ function gf2_H_to_G(H)
G[:,invperm(sindx)]
end

"""Performs in-place Gaussian elimination on a binary matrix and returns
its *row echelon form*,*rank*, the *transformation matrix*, and the *pivot
columns*. The transformation matrix that converts the original matrix into
the row echelon form. The `full` parameter controls the extent of elimination:
if `true`, only rows below the pivot are affected; if `false`, both above and
below the pivot are eliminated."""
function gf2_row_echelon_with_pivots!(M::AbstractMatrix{Int}; full=false)
r, c = size(M)
N = Matrix{Int}(LinearAlgebra.I, r, r)
p = 1
pivots = Int[]
for col in 1:c
@inbounds for row in p:r
if M[row, col] == 1
if row != p
M[[row, p], :] .= M[[p, row], :]
N[[row, p], :] .= N[[p, row], :]
end
break
end
end
if M[p, col] == 1
if !full
elim_range = p+1:r
else
elim_range = 1:r
end
@simd for j in elim_range
@inbounds if j != p && M[j, col] == 1
M[j, :] .= (M[j, :] .+ M[p, :]) .% 2
N[j, :] .= (N[j, :] .+ N[p, :]) .% 2
end
end
p += 1
push!(pivots, col)
end
if p > r
break
end
end
rank = p - 1
return M, rank, N, pivots
end

"""The nullspace of a binary matrix."""
function gf2_nullspace(H::AbstractMatrix{Int})
m = size(H',1)
_, matrix_rank, transformation_matrix, _ = gf2_row_echelon_with_pivots!(copy(H)')
if m == matrix_rank
# By the rank-nullity theorem, if rank(M) = m, then nullity(M) = 0
return zeros(Bool, 1, m)
end
# Extract the nullspace from the transformation matrix
return transformation_matrix[matrix_rank+1:end, :]
end

"""The basis for the row space of the binary matrix."""
function gf2_rowspace_basis(H::AbstractMatrix{Int})
pivots = gf2_row_echelon_with_pivots!(copy(H)')[4]
# Extract the rows corresponding to the pivot columns
H[pivots,:]
end

##############################
# Error classes
##############################
Expand Down
150 changes: 150 additions & 0 deletions test/test_gf2.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,4 +24,154 @@
end
end
end
@testset "GF(2) row echelon form with transformation matrix, pivots etc." begin
for n in test_sizes
for rep in 1:100
gf2_matrices = [rand(Bool, size, size) for size in test_sizes]
for (i, matrix) in enumerate(gf2_matrices)
echelon_form, _, transformation, _ = gf2_row_echelon_with_pivots!(Matrix{Int}(matrix), full=true) # in-place
# Check the correctness of the transformation matrix
@test (transformation*matrix) .%2 == echelon_form
# Check the correctness of Gaussian elimination
@test echelon_form == gf2_gausselim!(matrix)
end
end
end
end
function is_in_nullspace(A, x)
# Ensure x is the correct orientation
if size(x, 1) != size(A, 2)
x = transpose(x)
end
# Perform modulo 2 arithmetic: A * x must be zero mod 2
if size(x, 2) == 1 # x is a single column vector
result = A * x
return all(result .% 2 .== 0) # Check if A * x = 0 mod 2
else # x is a matrix, check each column vector
for i in 1:size(x, 2)
result = A * x[:, i] # Multiply A with the i-th column of x
if !all(result .% 2 .== 0) # Check if A * column = 0 mod 2
return false
end
end
return true # All columns are in the null space mod 2
end
end
@testset "GF(2) nullspace of the binary matrix" begin
for n in test_sizes
for rep in 1:100
gf2_matrices = [rand(Bool, size, size) for size in test_sizes]
for (i, matrix) in enumerate(gf2_matrices)
imat = Matrix{Int}(matrix)
ns = gf2_nullspace(imat)
@test is_in_nullspace(imat, ns)
end
end
end
end
@testset "Consistency check with ldpc" begin
# sanity checks for comparison to https://github.com/quantumgizmos/ldpc
# results compared with 'from ldpc.mod2 import nullspace, row_basis, row_echelon'
# Consistency check 1
H = [1 1 1; 1 1 1; 0 1 0]
echelon_form, rank, transformation, pivots = gf2_row_echelon_with_pivots!(copy(H)) # in-place
@test echelon_form == [1 1 1; 0 1 0; 0 0 0]
@test rank == 2
@test transformation == [1 0 0; 0 0 1; 1 1 0]
@test pivots == [1, 2] # in python, it's [0, 1] due to zero-indexing
@test mod.((transformation*copy(H)), 2) == echelon_form
@test gf2_nullspace(copy(H)) == [1 0 1]
@test gf2_rowspace_basis(copy(H)) == [1 1 1; 0 1 0]
# Consistency check 2
H = [0 0 0 1 1 1 1;
0 1 1 0 0 1 1;
1 0 1 0 1 0 1]
echelon_form, rank, transformation, pivots = gf2_row_echelon_with_pivots!(copy(H)) # in-place
@test echelon_form == [1 0 1 0 1 0 1;
0 1 1 0 0 1 1;
0 0 0 1 1 1 1]
@test rank == 3
@test transformation == [0 0 1;
0 1 0;
1 0 0]
@test pivots == [1, 2, 4] # in python, it's [0, 1, 3] due to zero-indexing
@test mod.((transformation*copy(H)), 2) == echelon_form
@test gf2_nullspace(copy(H)) == [1 1 1 0 0 0 0;
0 1 1 1 1 0 0;
0 1 0 1 0 1 0;
0 0 1 1 0 0 1]
@test gf2_rowspace_basis(copy(H)) == [0 0 0 1 1 1 1;
0 1 1 0 0 1 1;
1 0 1 0 1 0 1]
# Consistency check 3
H = [1 1 0; 0 1 1; 1 0 1]
echelon_form, rank, transformation, pivots = gf2_row_echelon_with_pivots!(copy(H)) # in-place
@test echelon_form == [1 1 0;
0 1 1;
0 0 0]
@test rank == 2
@test transformation == [1 0 0;
0 1 0;
1 1 1]
@test pivots == [1,2 ] # in python, it's [0, 1] due to zero-indexing
@test mod.((transformation*copy(H)), 2) == echelon_form
@test gf2_nullspace(copy(H)) == [1 1 1]
@test gf2_rowspace_basis(copy(H)) == [1 1 0;
0 1 1]
# Consistency check 4
H = [1 1 0; 0 1 0; 0 0 1]
echelon_form, rank, transformation, pivots = gf2_row_echelon_with_pivots!(copy(H)) # in-place
@test echelon_form == [1 1 0;
0 1 0;
0 0 1]
@test rank == 3
@test transformation == [1 0 0;
0 1 0;
0 0 1]
@test pivots == [1, 2, 3] # in python, it's [0, 1, 2] due to zero-indexing
@test mod.((transformation*copy(H)), 2) == echelon_form
@test gf2_nullspace(copy(H)) == [0 0 0]
@test gf2_rowspace_basis(copy(H)) == [1 1 0;
0 1 0;
0 0 1]
# Consistency check 5
H = [1 1 0; 0 1 0; 0 0 1; 0 1 1]
echelon_form, rank, transformation, pivots = gf2_row_echelon_with_pivots!(copy(H)) # in-place
@test echelon_form == [1 1 0;
0 1 0;
0 0 1;
0 0 0]
@test rank == 3
@test transformation == [1 0 0 0;
0 1 0 0;
0 0 1 0;
0 1 1 1]
@test pivots == [1, 2, 3] # in python, it's [0, 1, 2] due to zero-indexing
@test mod.((transformation*copy(H)), 2) == echelon_form
@test gf2_nullspace(copy(H)) == [0 0 0]
@test gf2_rowspace_basis(copy(H)) == [1 1 0;
0 1 0;
0 0 1]
# Consistency check 6
H = [0 0 0 1 1 1 1;
0 1 1 0 0 1 1;
1 0 1 0 1 0 1]
echelon_form, rank, transformation, pivots = gf2_row_echelon_with_pivots!(copy(H)) # in-place
@test echelon_form == [1 0 1 0 1 0 1;
0 1 1 0 0 1 1;
0 0 0 1 1 1 1]
@test rank == 3
@test transformation == [0 0 1;
0 1 0;
1 0 0]
@test pivots == [1, 2, 4] # in python, it's [0, 1, 3] due to zero-indexing
@test mod.((transformation*copy(H)), 2) == echelon_form
@test gf2_nullspace(copy(H)) == [1 1 1 0 0 0 0;
0 1 1 1 1 0 0;
0 1 0 1 0 1 0;
0 0 1 1 0 0 1]
@test gf2_rowspace_basis(copy(H)) == [0 0 0 1 1 1 1;
0 1 1 0 0 1 1;
1 0 1 0 1 0 1]
end
end

0 comments on commit 7b3238b

Please sign in to comment.