首页
学习
活动
专区
圈层
工具
发布
社区首页 >问答首页 >使用Python执行模矩阵求逆的最简单方法?

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

Stack Overflow用户
提问于 2010-11-27 02:28:03
回答 5查看 15K关注 0票数 22

我想取一个矩阵的模逆,就像Python中的[1,2,3,4] mod7。我看过numpy (它做矩阵求逆,但不做模矩阵求逆),我在网上看到了一些数论软件包,但似乎没有一个做这个相对常见的过程(至少对我来说是相对常见的)。

顺便说一句,上面矩阵的逆矩阵是[5,1],[5,3]。不过,我想让Python为我做这件事。

EN

回答 5

Stack Overflow用户

发布于 2011-07-15 06:38:56

当四舍五入误差不是问题时,这是一个黑客技巧:

  • 找到正则逆(可能有非整数条目)和行列式(整数),两者都在numpy
  • 中实现将逆乘以行列式,并四舍五入为整数(Hacky)
  • 现在将一切乘以行列式的乘法逆(模数,下面的代码)
  • 按模数

进行入口式mod <>f29

一种不那么老套的方法是实际实现高斯消除。这是我使用高斯消去法的代码,这是我自己写的(舍入误差对我来说是个问题)。Q是模数,它不一定是素数。

代码语言:javascript
复制
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

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

票数 12
EN

Stack Overflow用户

发布于 2010-11-28 02:10:29

Okay...for那些关心的人,我解决了我自己的问题。我花了一段时间,但我认为这是可行的。它可能不是最优雅的,应该包括一些更多的错误处理,但它是有效的:

代码语言:javascript
复制
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
票数 11
EN

Stack Overflow用户

发布于 2014-10-22 07:27:55

可以使用Sage (www.sagemath.org)将其计算为

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

尽管Sage安装起来非常庞大,而且您必须使用随附的python版本,但这是一件痛苦的事情。

票数 5
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/4287721

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档