使用Python执行模矩阵求逆的最简单方法？

Question

我想在Python中采用矩阵的模逆，例如 [[1,2],[3,4]] mod 7 。我看过 numpy （它进行矩阵求逆，但不进行模矩阵求逆），并且在网上看到了一些数论包，但似乎没有任何东西可以执行这个相对常见的过程（至少，对我来说似乎相对常见）。

顺便说一句，上述矩阵的逆矩阵是 [[5,1],[5,3]] (mod 7)。不过我希望 Python 能为我做这件事。

Answer 1

好吧...对于那些关心的人，我解决了我自己的问题。我花了一段时间，但我认为这有效。它可能不是最优雅的，应该包括更多的错误处理，但它可以工作：

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 2

当舍入误差不是问题时有效的黑客技巧：

• 找到常规逆（可能有非整数条目）和行列式（整数），两者都在 numpy 中实现，
• 将逆乘以行列式，然后舍入为整数（hacky）
• 现在将所有内容乘以行列式的乘法逆元（以您的模数为模，下面的代码）
• 用您的模数执行entrywise mod

一个不太hacky的方法是实际实现高斯消除。这是我使用高斯消去法的代码，这是我为了自己的目的而编写的（舍入错误对我来说是一个问题）。 q 是模数，不一定是素数。

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

编辑：正如评论者指出的那样，在某些情况下该算法会失败。修复这个问题有点不简单，而且我现在没有时间。当时，在我的例子中，它适用于随机矩阵（模数是大素数的乘积）。基本上，第一个非零条目可能与模数不互质。主要情况很简单，因为您可以搜索不同的行并交换。在非素数情况下，我认为可能所有领先条目都不是相对素数，因此您必须将它们组合起来

Answer 3

可以使用 Sage (www.sagemath.org) 计算：

<块引用>

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

虽然Sage安装起来很大，而且你必须使用它附带的python版本，这很痛苦。

Answer 4

'sympy' 包 Matrix 类函数 'sqMatrix.inv_mod(mod)' 计算小模数和任意大模数的模矩阵逆。通过将 sympy 与 numpy 结合起来，计算二维 numpy 数组的模逆变得很容易（请参见下面的代码片段）：

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 5

不幸的是 numpy 没有模块化算术实现。您始终可以使用行缩减或行列式对建议的算法进行编码，如所示 在这里。模逆似乎对于密码学非常有用。

Answer 6

您可以使用此处的答案计算矩阵的共轭https://stackoverflow.com/a/75566371/4454877，然后除以行列式（或者更确切地说，乘以它的倒数！）并在最后取元素模。我不确定这是否受到舍入问题的影响。

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

主要缺点是对于较大的矩阵来说它似乎相当慢。