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

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)

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]