Solve Integer Matrix Equation Ax = B (mod q) for x

Question

I have a known Python numpy array A whose shape is n * m (m > n), a known Python numpy array B whose shape is n * 1, and a positive integer q that is larger than 2. All values in A and B are integers that are in the interval [0, q). I wonder how to compute x to satisfy Ax = B (mod q) using Python. There may be many solutions of x but it is adequate to find out one of the solutions. Thanks a lot.

For example, A, B, and q are set as follows.

from numpy import array
A = array([[2, 6, 2, 0], [4, 0, 2, 1], [3, 6, 5, 2]])
B = array([[1], [6], [0]])
q = 7

I have tried the following codes but none of them are correct.

1. lstsq (usually output answers with floats)

Code:

from numpy.linalg import lstsq
return lstsq(A, B, rcond = None)[0].astype("int") % q

Output:

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

However, Ax is not equalled to B.

2. solve

Code:

from numpy.linalg import solve
return solve(A, B)

Output:

File "test.py", line 12, in sol2
    return solve(A, B)
  File "D:\Program Files\Python\lib\site-packages\numpy\linalg\linalg.py", line 396, in solve
    _assert_stacked_square(a)
  File "D:\Program Files\Python\lib\site-packages\numpy\linalg\linalg.py", line 213, in _assert_stacked_square
    raise LinAlgError('Last 2 dimensions of the array must be square')
numpy.linalg.LinAlgError: Last 2 dimensions of the array must be square

3. inv_mod from sympy

Code:

from numpy import asarray
from sympy import Matrix
A_inv = Matrix(A).inv_mod(q)
return asarray(A_inv.dot(B)).astype("int") % q

Output:

File "test.py", line 17, in sol3
    A_inv = Matrix(A).inv_mod(q)
  File "D:\Program Files\Python\lib\site-packages\sympy\matrices\matrices.py", line 2153, in inv_mod
    return _inv_mod(self, m)
  File "D:\Program Files\Python\lib\site-packages\sympy\matrices\inverse.py", line 169, in _inv_mod
    raise NonSquareMatrixError()
sympy.matrices.common.NonSquareMatrixError

4. Using (A * A.T).inv_mod(q)

Code:

from numpy import asarray, dot
from sympy import Matrix
return dot(dot(A.T, asarray(Matrix(dot(A, A.T)).inv_mod(q)).astype("int")), B) % q

Output:

File "test.py", line 23, in sol4
    x = dot(dot(A.T, asarray(Matrix(dot(A, A.T)).inv_mod(q)).astype("int")), B) % q
  File "D:\Program Files\Python\lib\site-packages\sympy\matrices\matrices.py", line 2153, in inv_mod
    return _inv_mod(self, m)
  File "D:\Program Files\Python\lib\site-packages\sympy\matrices\inverse.py", line 178, in _inv_mod
    raise NonInvertibleMatrixError('Matrix is not invertible (mod %d)' % m)
sympy.matrices.common.NonInvertibleMatrixError: Matrix is not invertible (mod 7)

Answer 1

For what you want - choosing one (arbitrary) integer solution - ILP does fine:

import numpy as np
from scipy.optimize import milp, Bounds, LinearConstraint

A = np.array((
    (2, 6, 2, 0),
    (4, 0, 2, 1),
    (3, 6, 5, 2),
))
B = np.array((1, 6, 0))
q = 7
m, n = A.shape

'''
Ax = B (mod q)
Ax = B + qf
Variables: x (n), f (m)
'''
constraint = LinearConstraint(
    A=np.hstack((
        A,
        np.diag(np.full(shape=m, fill_value=-q)),
    )),
    lb=B, ub=B,
)

result = milp(
    c=np.zeros(n + m),
    integrality=np.ones(n + m),
    bounds=Bounds(lb=np.zeros(n + m)),
    constraints=constraint,
)
assert result.success
print(result.message)
x, f = np.split(result.x.astype(int), (n,))
print(f'A{x} = {B} + {q}{f} = {A@x} = {B + q*f}')

Optimization terminated successfully. (HiGHS Status 7: Optimal)
A[0 4 6 1] = [1 6 0] + 7[5 1 8] = [36 13 56] = [36 13 56]

Answer 2

Assuming that q is a constraint for the solutions, and that you only want to have integer solutions in x, this problem is not solvable using "normal" optimization. You will need integer linear programming for that. Scipy provides a function for linear programming. With this you can use the integrality parameter, to specify whether solutions should be in the integer space.

Unconstrained Optimization

If q is a constraint for the inputs A and B only, and you're fine with real number solutions, lstsq does exactly what you need. The reason why you're getting garbage results is that you're cutting off everything after the decimal separator by casting to int with astype("int").

import numpy as np
from numpy.linalg import lstsq

if __name__ == '__main__':
    A = np.array([[2, 6, 2, 0], [4, 0, 2, 1], [3, 6, 5, 2]])
    b = np.array([1, 6, 0])

    q = 7

    solution = lstsq(A, b, rcond=None)[0]
    print("real number solution:", solution)
    reconstructed = A @ solution
    # after rounding to 10 significant digits, this should be equal to b
    rounded_reconstructed = np.round(reconstructed, 10)
    print(rounded_reconstructed, b)

    assert all(np.isclose(reconstructed, b)), "unequal after rounding"

You can see, that the solution is a very close real approximation. In this case, to about 15 significant digits. And that's all that lstsq will give you. If you expect certain constraints to be held, you have to make that explicit, NumPy won't guess.

Constrained Optimization

One way to turn a constrained optimization problem into an unconstrained one, that lstseq et al. could work with, is the Lagrangian. Here is a pretty decent question, that deals with lagrange multipliers. Another useful thing to look into for constrained optimization problems is the simplex algorithm.

This will not help you if you are looking for integer solutions. But if you want those, it gets a lot harder very quickly. This is "real math".^[1]

General Advice

Put more effort into your questions. There's a guide on asking good questions.

For example: It helps to know what you're working on, to understand the problem and why you need something. Here, it's pretty hard to know what kind of solution you're looking for. It also helps to get some examples not only for inputs, but also the outputs that you would expect for them.

---

[1]: actually it's discrete math =D