Improve performances of Gauss-Jordan (XOR-SAT) Algorithm?

In this question I was looking for an algorithm to solve what seems to be a XOR-SAT problem.

Let's consider this equation system :

A ^ B ^ C = 0
A ^ B ^ D = 1
E ^ F ^ G = 2
D ^ E ^ F = 3
G ^ F ^ H = 4
A ^ F ^ G = 5


where ^ is a bit-wise XOR operation. I can decompose this system for each bit of the right-hand value (here 5 -> 3-bit number):

Bit 0

A ^ B ^ C = 0
A ^ B ^ D = 1
E ^ F ^ G = 0
D ^ E ^ F = 1
G ^ F ^ H = 0
A ^ F ^ G = 1


Bit 1

A ^ B ^ C = 0
A ^ B ^ D = 0
E ^ F ^ G = 1
D ^ E ^ F = 1
G ^ F ^ H = 0
A ^ F ^ G = 0


Bit 2

A ^ B ^ C = 0
A ^ B ^ D = 0
E ^ F ^ G = 0
D ^ E ^ F = 0
G ^ F ^ H = 1
A ^ F ^ G = 1


Algorithm

The system can be written such as:

system = [(0, 1, 2),
(0, 1, 3),
(4, 5, 6),
(3, 4, 5),
(6, 5, 7),
(0, 5, 6)]


The algorithm is the following:

nvars = max(max(i) for i in system) + 1
matrix = [[int(i in system[k]) for i in range(nvars)] for k in range(len(system))]
solutions = []
bits = range(ceil(log2(len(system))))
for bit in bits:
aug_matrix = np.array([k + [int(bool(i & (1 << bit)))]
for i, k in enumerate(matrix)], dtype=int)

reduced_matrix = gauss_jordan_gf2(aug_matrix.copy())

# Extract the solution by solving the system
def extract_solution(matrix):
rows, cols = matrix.shape
solution = np.zeros(cols - 1, dtype=int)
for i in range(rows):
pivot_col = np.where(matrix[i, :-1] == 1)[0]
if len(pivot_col) > 0:
pivot = pivot_col[0]
solution[pivot] = matrix[i, -1]
return solution

solution = extract_solution(reduced_matrix)

# Verify the bit solution
for i, equation in enumerate(system):
got = reduce(xor, [solution[k] for k in equation])
assert got == int(bool(i & (1 << bit)))
solutions.append(np.array(solution))

# Combine solutions
s = np.sum([s * exp for s, exp in zip(solutions, (2**b for b in bits))], axis=0)
for i, equation in enumerate(system):
assert reduce(xor, [s[k] for k in equation])== i


It works well and gives this result: array([1, 1, 0, 1, 6, 4, 0, 0])

However, for a much larger problem (200 system) given below, it takes about 3 seconds. Is there any way I can improve it?

Not that I had to modify my gauss_jordan_gf2 to accept larger matrices:

def gauss_jordan_gf2(matrix):
rows, cols = matrix.shape
for r in range(rows):
return matrix
i = r
i += 1
if i == rows:
i = r
return matrix
matrix[[i, r]] = matrix[[r, i]]
matrix[r] = matrix[r] / lv
for i in range(rows):
if i != r:
matrix[i] = matrix[i] ^ (matrix[r] * lv)
return matrix

system = [[0, 1, 2],
[3, 4, 5],
[6, 7, 8],
[9, 10, 11],
[12, 13, 14],
[15, 16, 17],
[18, 10, 19],
[20, 21, 22],
[23, 24, 25],
[26, 27, 28],
[29, 30, 31],
[32, 33, 34],
[35, 36, 37],
[38, 39, 40],
[41, 42, 17],
[43, 44, 45],
[46, 47, 48],
[35, 49, 50],
[51, 52, 53],
[54, 55, 50],
[32, 56, 57],
[58, 7, 59],
[60, 61, 62],
[63, 64, 65],
[66, 30, 67],
[43, 68, 69],
[70, 71, 72],
[73, 74, 75],
[54, 1, 76],
[9, 10, 77],
[78, 21, 5],
[60, 79, 80],
[23, 81, 82],
[6, 74, 75],
[83, 71, 84],
[18, 56, 34],
[85, 55, 31],
[86, 87, 11],
[88, 89, 14],
[90, 91, 37],
[92, 87, 93],
[94, 95, 96],
[97, 10, 65],
[98, 89, 22],
[99, 100, 101],
[94, 102, 103],
[104, 30, 105],
[106, 107, 59],
[78, 4, 93],
[43, 68, 108],
[94, 44, 109],
[110, 111, 57],
[70, 112, 113],
[51, 114, 115],
[110, 116, 117],
[78, 118, 119],
[120, 121, 122],
[0, 123, 82],
[15, 124, 125],
[90, 126, 2],
[127, 128, 129],
[130, 131, 115],
[20, 132, 80],
[6, 133, 31],
[92, 87, 31],
[83, 91, 28],
[134, 135, 14],
[88, 21, 93],
[136, 74, 137],
[130, 71, 53],
[104, 30, 138],
[139, 10, 140],
[141, 102, 142],
[38, 143, 103],
[139, 16, 144],
[145, 52, 146],
[41, 147, 137],
[46, 148, 149],
[43, 27, 150],
[151, 152, 142],
[70, 112, 153],
[78, 118, 62],
[154, 49, 155],
[43, 118, 117],
[0, 156, 157],
[158, 36, 159],
[160, 47, 67],
[161, 162, 163],
[164, 165, 69],
[54, 116, 166],
[83, 167, 77],
[23, 168, 113],
[169, 74, 170],
[171, 172, 115],
[92, 1, 173],
[130, 174, 175],
[26, 121, 176],
[29, 177, 75],
[134, 68, 22],
[178, 128, 179],
[180, 181, 129],
[130, 182, 34],
[58, 121, 93],
[85, 47, 48],
[151, 68, 48],
[51, 182, 153],
[58, 183, 122],
[73, 184, 163],
[145, 39, 166],
[154, 49, 185],
[66, 61, 62],
[51, 49, 115],
[160, 52, 186],
[187, 33, 67],
[20, 56, 122],
[188, 165, 125],
[189, 111, 2],
[120, 167, 138],
[90, 4, 2],
[127, 143, 129],
[88, 91, 5],
[23, 112, 34],
[29, 181, 190],
[191, 21, 173],
[192, 193, 176],
[104, 194, 138],
[178, 36, 28],
[192, 64, 105],
[180, 184, 37],
[86, 95, 96],
[139, 124, 146],
[141, 49, 65],
[58, 148, 22],
[85, 56, 67],
[70, 126, 153],
[43, 195, 103],
[46, 135, 196],
[51, 49, 196],
[35, 197, 125],
[32, 183, 113],
[78, 71, 8],
[158, 68, 119],
[198, 10, 67],
[189, 111, 57],
[161, 1, 82],
[127, 107, 199],
[26, 200, 170],
[171, 182, 115],
[120, 52, 31],
[201, 95, 11],
[88, 81, 173],
[92, 68, 76],
[3, 202, 22],
[201, 172, 53],
[88, 203, 179],
[26, 55, 105],
[188, 116, 40],
[85, 123, 176],
[180, 89, 45],
[88, 128, 48],
[29, 194, 50],
[92, 102, 142],
[141, 30, 153],
[35, 56, 101],
[139, 100, 53],
[106, 91, 28],
[97, 56, 204],
[43, 197, 40],
[99, 126, 93],
[73, 24, 69],
[78, 182, 122],
[198, 81, 59],
[12, 148, 205],
[60, 123, 5],
[164, 33, 186],
[9, 162, 80],
[198, 156, 173],
[188, 16, 69],
[206, 168, 176],
[164, 197, 75],
[6, 131, 170],
[26, 202, 5],
[120, 7, 101],
[164, 100, 207],
[18, 183, 11],
[130, 133, 14],
[208, 132, 96],
[209, 143, 210],
[46, 165, 28],
[104, 100, 207],
[32, 61, 108],
[54, 148, 93],
[110, 121, 45],
[145, 184, 65],
[211, 52, 157],
[98, 121, 149],
[70, 33, 57],
[151, 79, 196],
[73, 39, 125],
[161, 116, 69]]

• – D.W.
May 15 at 22:59
• @D.W. I searched a bit, but lu or solve don't work on the GF(2) Galois field. I somehow already work with numpy. yesterday
• Some possibilities: Does bitbucket.org/malb/m4ri/src/master help? I read that it is used by SAGE. It looks like it might have an efficient implementation of solving linear systems over GF(2). Alternatively, is github.com/mhostetter/galois faster than what you've written? I've never used it before and I don't know whether it will be significantly faster.
– D.W.
yesterday
• Wonderful, galois does it with galois.GF2(augmented_matrix).row_reduce(). yesterday

Thanks to D.W.'s comment, I discovered galois and now have a much better understanding of what the Gauss-Jordan elimination, specifically galois.GF2(m).row_reduce(), does.

The last column in the augmented matrix contains the constants. These constants are not considered during the Gaussian reduction. However, we know the steps during this reduction could involve either:

1. Swapping rows
2. XOR-ing rows

Using the built-in row_reduce() is tempting because it is well optimized. However, in our case, it would require performing the Gaussian reduction for each bit. Despite this, the RREF matrix without the constants will be the same for all bits, leading to wasted time.

It is therefore better to re-implement the Gauss-Jordan elimination algorithm in GF(2) using vectorized functions from numpy and directly apply each step to the constant term vector. There is no need to decompose it into bits, thanks to bit-wise XOR operations.

Here is the optimized version, which is about 30 times faster than the original:

def gauss_jordan_elimination_gf2(A, b):
""" Perform Gauss-Jordan elimination on a matrix in GF(2).
Return a solution of the system of equations in GF(2).
"""
A = A.copy()
m, n = A.shape
h, k = 0, 0

# Compute RREF, apply each row operation to b
while h < m and k < n:
# Find the k-th pivot
i_max = np.argmax(A[h:m, k]) + h

if A[i_max, k] == 0:
# No pivot in this column, pass to the next
k += 1
continue

# Swap rows
A[[h, i_max]] = A[[i_max, h]]
b[[h, i_max]] = b[[i_max, h]]

# Vectorized row reduction
non_pivot_rows = np.where((np.arange(m) != h) & (A[:, k] == 1))[0]
A[non_pivot_rows] ^= A[h]
b[non_pivot_rows] ^= b[h]

h += 1
k += 1

solution = np.zeros(A.shape[1], dtype=int)
solution[np.argmax(A, axis=1)] = b

return solution


Full Benchmark

1. OP's solution
2. Optimized solution
3. Galois row_reduce ran for each bit

The OP's long system (200 entries) is used.

import numpy as np
import galois
import timeit
from math import ceil, log2
from functools import reduce
from operator import xor

def create_matrix(system):
""" Create a matrix from a system of equations.
Each row of the system describe which columns are set to 1.
The number of columns is the maximum value found in the system.
>>> system = np.array([(0, 1, 2),
(4, 5, 6)])
array([[1, 1, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 1, 1]], dtype=int32)
"""
rows = system.shape[0]
cols = system.max() + 1
matrix = np.zeros((rows, cols), dtype=np.int32)
matrix[np.arange(rows)[:, None], system] = 1
return matrix

def augmented_matrix(A, b):
""" Build the augmented matrix [A|b].
>>> A = array([[ 1,  2,  3,  4],
[ 5,  6,  7,  8],
[ 9, 10, 11, 12]])
>>> b = array([100, 200, 300])
>>> augmented_matrix(A, b)
array([[  1,   2,   3,   4, 100],
[  5,   6,   7,   8, 200],
[  9,  10,  11,  12, 300]])
"""
return np.hstack([A, b.reshape(-1, 1)])

def decompose_bits(n):
""" Decompose a range of integers into bits.
>>> decompose_bits(4)
array([[0, 1, 0, 1],
[0, 0, 1, 1]])
"""
bits = int(np.ceil(np.log2(n)))
matrix = np.arange(n)[:, None] >> np.arange(bits) & 1
return matrix

def verify_solution(system, solution):
""" Check if the solution is correct. """
v = solution[system]
s = np.bitwise_xor.reduce(v, axis=1)
return np.array_equal(s, np.arange(len(s)))

A = galois.GF(2)(create_matrix(system))

ret = None

def experiment1():
global ret

# OP's solution

ret = s
return s

def experiment2():
global ret
b = np.arange(system.shape[0])
ret = gauss_jordan_elimination_gf2(A.copy(), b)
return ret

def experiment3():
global ret
m = create_matrix(system)
b = decompose_bits(system.shape[0]).T

ma = galois.GF2(np.ones(np.shape(m), dtype=int))
constants = galois.GF2(np.zeros((m.shape[0], b.shape[0]), dtype=int))

# Solve the system of equations (Gauss-Jordan elimination)
for i in range(b.shape[0]):
am = galois.GF2(augmented_matrix(m, b[i]))
rref = am.row_reduce()
ma = np.bitwise_and(ma, rref[:,:-1])
constants[:,i] = rref[:, -1]

# Aggregate constants bit-wise
s = np.zeros((b.shape[0], m.shape[1]), dtype=np.int32)
shifts = np.arange(s.shape[0])
c = np.bitwise_or.reduce(constants << shifts, axis=1)

positions = np.argmax(ma, axis=1)

sol = np.zeros(m.shape[1], dtype=int)
sol[positions] = c

ret = sol
return sol

time1 = timeit.timeit(experiment1, number=10) / 10
assert (verify_solution(system, ret))
print(time1)

time2 = timeit.timeit(experiment2, number=10) / 10
assert (verify_solution(system, ret))
print(time2)

time3 = timeit.timeit(experiment3, number=10) / 10
assert (verify_solution(system, ret))
print(time3)

print(time1 / time2)
print(time3 / time2)

• OP's solution: 0.69s
• Optimized version: 0.019s
• Galois: 0.20s

Here we are 35x faster, or only 10x faster using galois.