From 7b3238ba4057b5888e2c843e1f0f0b5bee257a77 Mon Sep 17 00:00:00 2001 From: Fe-r-oz Date: Thu, 5 Dec 2024 04:07:01 +0500 Subject: [PATCH] Enhancements to GF(2) Linear Algebra: in-place row echelon with pivots, nullspace, and basis for row space --- src/QuantumClifford.jl | 64 ++++++++++++++++++ test/test_gf2.jl | 150 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 214 insertions(+) diff --git a/src/QuantumClifford.jl b/src/QuantumClifford.jl index 6b90af7a8..f5815e97c 100644 --- a/src/QuantumClifford.jl +++ b/src/QuantumClifford.jl @@ -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 @@ -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 ############################## diff --git a/test/test_gf2.jl b/test/test_gf2.jl index 45aa6701b..7ebfcc4a6 100644 --- a/test/test_gf2.jl +++ b/test/test_gf2.jl @@ -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