用Python进行模矩阵求逆的最简单方法是什么？

Question

我想在Python中计算一个矩阵的模逆，比如[[1,2],[3,4]]模7。我查过numpy（它可以做矩阵的逆，但不能做模矩阵的逆），还在网上看到一些数论的包，但没有找到一个能做这个相对常见的操作的工具（至少在我看来是相对常见的）。

顺便说一下，上面那个矩阵的逆是[[5,1],[5,3]]（模7）。不过我希望Python能帮我算出来。

Answer 1

可以通过Sage这个工具来计算，Sage的官网是www.sagemath.org。

你可以用这段代码来计算：Matrix(IntegerModRing(7), [[1, 2], [3,4]]).inverse()

不过，Sage的安装比较复杂，而且你必须使用它自带的Python版本，这点让人觉得麻烦。

Answer 2

这里有一个小技巧，当四舍五入误差不是问题的时候可以用：

• 先找到常规的逆矩阵（可能会有非整数的值）和行列式（这是一个整数），这两个都可以用numpy来实现。
• 把逆矩阵乘以行列式，然后四舍五入成整数（这有点小技巧的味道）。
• 接着，把所有的结果乘以行列式的乘法逆元（要注意模运算，下面有代码）。
• 最后，对每个元素进行模运算。

还有一种不那么“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

好吧……对于在意这个的人，我自己解决了问题。花了我一段时间，但我觉得这个方法可行。可能不是最优雅的写法，也应该加一些错误处理，但总的来说是能用的：

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