Самый простой способ выполнения модульной инверсии матрицы с помощью Python?

Ответ 1

Хакинг-трюк, который работает при ошибках округления, не является проблемой:

• найти регулярные обратные (могут иметь нецелые записи) и детерминант (целое число), оба реализованы в numpy
• умножить инверсию на определитель и округлить до целых (hacky)
• теперь умножаем все на детерминантный мультипликативный обратный (по модулю ваш модуль, код ниже)
• введите mod mod по модулю

Менее опасным способом является реализация гауссовой элиминации. Здесь мой код с использованием исключения Гаусса, который я написал для моих целей (ошибки округления были для меня проблемой). 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

EDIT: как отмечают комментаторы, в некоторых случаях этот алгоритм терпит неудачу. Это немного нетривиально для исправления, и сейчас у меня нет времени. Тогда он работал для случайных матриц в моем случае (модули были продуктами больших простых чисел). В принципе, первая ненулевая запись может быть не взаимно простой с модулем. Простой случай легко, так как вы можете искать другую строку и своп. В случае, не относящемся к делу, я думаю, что все ведущие записи не являются относительно простыми, поэтому вы должны их комбинировать

Ответ 2

Хорошо... для тех, кто заботится, я решил свою проблему. Мне потребовалось некоторое время, но я думаю, что это работает. Он, вероятно, не самый элегантный и должен включать в себя еще некоторую обработку ошибок, но он работает:

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

Ответ 3

Его можно вычислить с помощью Sage (www.sagemath.org) как

Матрица (IntegerModRing (7), [[1, 2], [3,4]]). inverse()

Хотя Sage огромный для установки, и вы должны использовать версию python, которая поставляется с ней, что является болью.

Ответ 4

К сожалению, numpy не имеет модульных арифметических реализаций. Вы всегда можете кодировать предложенный алгоритм с помощью сокращения строк или детерминант, как показано здесь. Модульный обратный, по-видимому, весьма полезен для криптографии.

Ответ 5

Этот маленький фрагмент кода, похоже, делает это: ссылка

Обратите внимание на приведенный ниже комментарий для небольшого улучшения. Насколько я вижу, кажется, что нужно сделать правильную линейную алгебру. Я никогда не встречал никаких опций в обычных пакетах, поэтому, возможно, самый простой способ получить фрагмент кода из Интернета (есть намного больше).