Easiest way to perform modular matrix inversion with Python?

Question

I'd like to take the modular inverse of a matrix like [[1,2],[3,4]] mod 7 in Python. I've looked at numpy (which does matrix inversion but not modular matrix inversion) and I saw a few number theory packages online, but nothing that seems to do this relatively common procedure (at least, it seems relatively common to me).

By the way, the inverse of the above matrix is [[5,1],[5,3]] (mod 7). I'd like Python to do it for me though.

Answer 1

You can calculate the adjugate of the matrix using the answer here https://stackoverflow.com/a/75566371/4454877, then divide by the determinant (or rather, multiply by its inverse!) and take the elementwise modulus at the end. I'm unsure if this is affected by rounding issues.

def modular_matrix_inverse(mat: np.matrix, modulus: int):
    return (adjugate(mat) * pow(int(np.linalg.det(mat)), -1, modulus)) % modulus

The main downside is that it seems rather slow for larger matrices.

Answer 2

'sympy' package Matrix class function 'sqMatrix.inv_mod(mod)' computes modulo matrix inverse for small and arbitrarily large modulus. By combining sympy with numpy, it becomes easy to compute modulo inverse of 2-D numpy arrays (see the code snippet below):

import numpy
from sympy import Matrix

def matInvMod (vmnp, mod):
    nr = vmnp.shape[0]
    nc = vmnp.shape[1]
    if (nr!= nc):
        print "Error: Non square matrix! exiting"
        exit()
    vmsym = Matrix(vmnp)
    vmsymInv = vmsym.inv_mod(mod)
    vmnpInv = numpy.array(vmsymInv)
    print "vmnpInv: ", vmnpInv, "\n"
    k = nr
    vmtest = [[1 for i in range(k)] for j in range(k)]  # just a 2-d list
    vmtestInv = vmsym*vmsymInv
    for i in range(k):
      for j in range(k):
         #print i, j, vmtrx2[i,j] % mod
         vmtest[i][j] = vmtestInv[i,j] % mod
    print "test vmk*vkinv % mod \n:", vmtest
    return vmnpInv

if __name__ == '__main__':
    #p = 271
    p = 115792089210356248762697446949407573530086143415290314195533631308867097853951
    vm = numpy.array([[1,1,1,1], [1, 2, 4, 8], [1, 4, 16, 64], [1, 5, 25, 125]])   
    #vminv = modMatInv(vm, p)
    vminv = matInvMod(vm, p)
    print vminv
    vmtestnp = vm.dot(vminv)%p     # test mtrx inversion
    print vmtestnp

Answer 3

A hackish trick which works when rounding errors aren't an issue:

• find the regular inverse (may have non-integer entries), and the determinant (an integer), both implemented in numpy
• multiply the inverse by the determinant, and round to integers (hacky)
• now multiply everything by the determinant's multiplicative inverse (modulo your modulus, code below)
• do entrywise mod by your modulus

A less hackish way is to actually implement gaussian elimination. Here's my code using Gaussian elimination, which I wrote for my own purposes (rounding errors were an issue for me). q is the modulus, which is not necessarily prime.

def generalizedEuclidianAlgorithm(a, b):
    if b > a:
        return generalizedEuclidianAlgorithm(b,a);
    elif b == 0:
        return (1, 0);
    else:
        (x, y) = generalizedEuclidianAlgorithm(b, a % b);
        return (y, x - (a / b) * y)

def inversemodp(a, p):
    a = a % p
    if (a == 0):
        print "a is 0 mod p"
        return None
    if a > 1 and p % a == 0:
        return None
    (x,y) = generalizedEuclidianAlgorithm(p, a % p);
    inv = y % p
    assert (inv * a) % p == 1
    return inv

def identitymatrix(n):
    return [[long(x == y) for x in range(0, n)] for y in range(0, n)]

def inversematrix(matrix, q):
    n = len(matrix)
    A = np.matrix([[ matrix[j, i] for i in range(0,n)] for j in range(0, n)], dtype = long)
    Ainv = np.matrix(identitymatrix(n), dtype = long)
    for i in range(0, n):
        factor = inversemodp(A[i,i], q)
        if factor is None:
             raise ValueError("TODO: deal with this case")
        A[i] = A[i] * factor % q
        Ainv[i] = Ainv[i] * factor % q
        for j in range(0, n):
            if (i != j):
                factor = A[j, i]
                A[j] = (A[j] - factor * A[i]) % q
                Ainv[j] = (Ainv[j] - factor * Ainv[i]) % q
    return Ainv

EDIT: as commenters point out, there are some cases this algorithm fails. It's slightly nontrivial to fix, and I don't have time nowadays. Back then it worked for random matrices in my case (the moduli were products of large primes). Basically, the first non-zero entry might not be relatively prime to the modulus. The prime case is easy since you can search for a different row and swap. In the non-prime case, I think it could be that all leading entries aren't relatively prime so you have to combine them

Answer 4

It can be calculated using Sage (www.sagemath.org) as

Matrix(IntegerModRing(7), [[1, 2], [3,4]]).inverse()

Although Sage is huge to install and you have to use the version of python that comes with it which is a pain.

Answer 5

Okay...for those who care, I solved my own problem. It took me a while, but I think this works. It's probably not the most elegant, and should include some more error handling, but it works:

import numpy
import math
from numpy import matrix
from numpy import linalg

def modMatInv(A,p):       # Finds the inverse of matrix A mod p
  n=len(A)
  A=matrix(A)
  adj=numpy.zeros(shape=(n,n))
  for i in range(0,n):
    for j in range(0,n):
      adj[i][j]=((-1)**(i+j)*int(round(linalg.det(minor(A,j,i)))))%p
  return (modInv(int(round(linalg.det(A))),p)*adj)%p

def modInv(a,p):          # Finds the inverse of a mod p, if it exists
  for i in range(1,p):
    if (i*a)%p==1:
      return i
  raise ValueError(str(a)+" has no inverse mod "+str(p))

def minor(A,i,j):    # Return matrix A with the ith row and jth column deleted
  A=numpy.array(A)
  minor=numpy.zeros(shape=(len(A)-1,len(A)-1))
  p=0
  for s in range(0,len(minor)):
    if p==i:
      p=p+1
    q=0
    for t in range(0,len(minor)):
      if q==j:
        q=q+1
      minor[s][t]=A[p][q]
      q=q+1
    p=p+1
  return minor

Answer 6

Unfortunately numpy does not have modular arithmetic implementations. You can always code up the proposed algorithm using row reduction or determinants as demonstrated here. A modular inverse seems to be quite useful for cryptography.