How to solve a linear rectangular matrix equation over mod 2?

I am trying to solve a linear matrix equation of the AX = B using python, where A,B,X are binary matrices i.e. over GF(2) :

A = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0],
   [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1],
   [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1],
   [0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
   [0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
   [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0],
   [1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
   [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
   [0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0],
   [0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0],
   [0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0]])

I know that A has rank 11, and shape (11,12). i.e. A is a wide rectangular matrix with more variables than equations. So there is likely to be a free variable and a parameterized solution.

The RHS, B in the equation has the form,

B = np.array([[1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
   [0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
   [0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
   [0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0],
   [0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0],
   [0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0],
   [0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0],
   [0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0],
   [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0],
   [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0],
   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1]])

Which I can also calculate has rank 11.

I want find linear combinations of rows in A that give B. Ideally row-swaps are included.

• The usual methods in python (np.linalg.solve(A,B)) can't be used straightforwardly since A and B are not square matrices.

• I tried the numpy solver in galois next. Since these are finite field matrices, I use first galois to wrap these into mod 2 objects,

  import galois GF = galois.GF(2) A1 = GF(A)

• But the non-square matrix problem remains: Padding with redundant rows makes the matrix singular, and the linear system unsolvable using np.linalg.solve or scipy.linalg.solve.

• Next I obtained the row-reduced echelon form, calculated using A1.row_reduce() on the galois object

A1.row_reduce() = GF([[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
[0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1],
[0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1],
[0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1],
[0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1],
[0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1],
[0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1],
[0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1]], order=2)

(From the above, it's easy to see consecutive rows sum to give RHS, ie row1+row2 of row-reduced gives B), but the galois rref function doesn't help me arrive at an X or sequence of explicit operations to transform A to B. I know a solution exists because of the rank equality and row-reduced form, I just don't know how to find X.

tldr; how does one solve for a wide rectangular (i.e. non-square) system of equations of mod 2 matrices, using python?

Edit: Mathematica LinearSolve[A,B, Modulus->2] and Sage A.solve_right(B) where A and B are defined as GF(2) matrices both helped get a solution. I'm hoping to find a way to solve this from within Python.

By sheer chance, I made such a tool ages ago specifically for this purpose, called solve_gf2, available here: https://github.com/nneonneo/pwn-stuff/blob/master/math/gf2.py. It takes as input a matrix A of any size and a vector b, and solves for x in Ax=b over GF(2). It yields all possible solutions; generally, there will be 2^k solutions where k is the number of free variables.

To extend this to solving for a matrix X in AX = B, it suffices to run the solver once per column of B. You can combine any set of solutions to obtain the desired matrix X.

Here's how you would use this library:

from gf2 import solve_gf2
import numpy as np

A = np.array(
  [[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0],
   [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1],
   [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1],
   [0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
   [0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
   [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0],
   [1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
   [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
   [0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0],
   [0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0],
   [0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0]])

B = np.array(
  [[1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
   [0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
   [0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
   [0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0],
   [0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0],
   [0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0],
   [0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0],
   [0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0],
   [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0],
   [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0],
   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1]])

Am = list(map(list, A))
for i in range(B.shape[1]):
    b = list(B[:, i])
    for x in solve_gf2(Am, b):
        print(x, (A @ x) % 2)
    print()

and here's the sample output:

[1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1] [1 0 0 0 0 0 0 0 0 0 0]
[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0] [1 0 0 0 0 0 0 0 0 0 0]

[0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0] [1 1 0 0 0 0 0 0 0 0 0]
[1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1] [1 1 0 0 0 0 0 0 0 0 0]

[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0] [0 1 1 0 0 0 0 0 0 0 0]
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1] [0 1 1 0 0 0 0 0 0 0 0]

[0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1] [0 0 1 1 0 0 0 0 0 0 0]
[1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0] [0 0 1 1 0 0 0 0 0 0 0]

[1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] [0 0 0 1 1 0 0 0 0 0 0]
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] [0 0 0 1 1 0 0 0 0 0 0]

[1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1] [0 0 0 0 1 1 0 0 0 0 0]
[0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0] [0 0 0 0 1 1 0 0 0 0 0]

[0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0] [0 0 0 0 0 1 1 0 0 0 0]
[1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1] [0 0 0 0 0 1 1 0 0 0 0]

[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] [0 0 0 0 0 0 1 1 0 0 0]
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] [0 0 0 0 0 0 1 1 0 0 0]

[0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1] [0 0 0 0 0 0 0 1 1 0 0]
[1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0] [0 0 0 0 0 0 0 1 1 0 0]

[1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1] [0 0 0 0 0 0 0 0 1 1 0]
[0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0] [0 0 0 0 0 0 0 0 1 1 0]

[1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0] [0 0 0 0 0 0 0 0 0 1 1]
[0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1] [0 0 0 0 0 0 0 0 0 1 1]

[0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1] [0 0 0 0 0 0 0 0 0 0 1]
[1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0] [0 0 0 0 0 0 0 0 0 0 1]

Note: this library expects that the inputs are pure Python lists, not NumPy arrays; it may silently fail (producing nonsensical results) if passed NumPy arrays.

Since A is an underdetermined system, you can try using numpy's implementation of least squares.

X = np.linalg.lstsq(A, B)[0]
AX = A @ X

tol = 1e-9 # clean up some very small floats
AX = np.where(AX < tol, 0, AX)

print(AX)

# [[1. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
#  [0. 1. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
#  [0. 0. 1. 1. 0. 0. 0. 0. 0. 0. 0. 0.]
#  [0. 0. 0. 1. 1. 0. 0. 0. 0. 0. 0. 0.]
#  [0. 0. 0. 0. 1. 1. 0. 0. 0. 0. 0. 0.]
#  [0. 0. 0. 0. 0. 1. 1. 0. 0. 0. 0. 0.]
#  [0. 0. 0. 0. 0. 0. 1. 1. 0. 0. 0. 0.]
#  [0. 0. 0. 0. 0. 0. 0. 1. 1. 0. 0. 0.]
#  [0. 0. 0. 0. 0. 0. 0. 0. 1. 1. 0. 0.]
#  [0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 1. 0.]
#  [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 1.]]

However, the X you get this way is not a "nice" one, with row swaps.

Since it is an under-determined system, it may have infinitely many solutions.