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

Question

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.

FWIW: the solution consists of a bunch of 12-element vectors. There are only 2^12 binary 12-element vectors, which makes this actually brute-forceable if you want - iterate through all 4096 vectors x, and see which ones give you Ax = b for each column b of B. — nneonneo, Commented Jul 10 at 5:14
True, but I'll be solving equations like this with arbitrary dimension as well :) I also gave a shot to sage and mathematica which turn out to have solvers for this - those were great. But a tool like the one you shared allows one to stay in python. — clearski, Commented Jul 10 at 12:28

nneonneo · Accepted Answer · 2024-07-10 12:56:51Z

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.

Thanks! This is exactly what I needed. A question- the two rows in each line of the output look like 1s complement of each other. Why are there two rows - are both the solution? — clearski, Commented Jul 10 at 8:14
Yes, both rows are solutions. On the left is the solution vector x; on the right is the corresponding b column calculated by using Ax=b. You can see that in both cases you get the same b. — nneonneo, Commented Jul 10 at 12:16
I see! And for the final matrix solution, you could pick either option for each row? Giving 2^r possible solutions if there are r rows with 2 solutions each? — clearski, Commented Jul 10 at 12:30
Yes, exactly. This is what I meant by You can combine any set of solutions to obtain the desired matrix X. — nneonneo, Commented Jul 10 at 13:35

Vorac · Accepted Answer · 2024-07-19 12:18:08Z

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.

The question asks about matrices in GF(2); using least-squares on such problems won't produce meaningful results. For example, note that 1 + 1 = 0 in GF(2). — nneonneo, Commented Jul 10 at 5:09

Collectives™ on Stack Overflow

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

2 Answers 2

Not the answer you're looking for? Browse other questions tagged
python
numpy
or ask your own question.

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

Not the answer you're looking for? Browse other questions tagged pythonnumpy or ask your own question.

Related

Not the answer you're looking for? Browse other questions tagged
python
numpy
or ask your own question.