用Python执行模块化矩阵求逆的最简单方法? [英] Easiest way to perform modular matrix inversion with Python?

问题描述

我想采用Python [[1,2]，[3,4]] mod 7矩阵的模逆.我已经看过numpy(它执行矩阵求逆，但不进行模块化矩阵求逆)，并且在网上看到了一些数字理论软件包，但是似乎没有什么东西可以做这种相对通用的过程(至少，对我来说似乎是相对通用的). /p>

顺便说一下，上述矩阵的逆是[[5,1]，[5,3]](mod 7).我希望Python可以帮我做.

解决方案

当舍入错误不成问题时可以使用的骇客技巧:

• 查找正则逆(可能具有非整数条目)和行列式(整数)，均以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

正如评论者所指出的，在某些情况下此算法失败.修复起来有些微不足道，如今我没有时间.当时，它适用于我的情况下的随机矩阵(模是大质数的乘积).基本上，第一个非零条目可能不是模数的质数.最基本的情况很容易，因为您可以搜索其他行并交换.在非主要情况下，我认为可能是所有前导项都不是相对主要的，因此您必须将它们组合起来

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.

解决方案

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

这篇关于用Python执行模块化矩阵求逆的最简单方法?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！