私は[[1,2]、[3,4]] mod 7のような行列の逆行列をPythonで取りたいと思います。私はnumpyを見てきましたが(行列反転はしますが、行列の逆行列はありません)、オンラインでいくつかの数学理論パッケージを見ましたが、この比較的一般的な手順を行うようなものは何もありませんでした(少なくとも、Pythonでモジュール式の行列反転を行う最も簡単な方法は?
ところで、上記の行列の逆数は、[[5,1]、[5,3]](mod 7)です。私はPythonでそれをやって欲しいです。
私は[[1,2]、[3,4]] mod 7のような行列の逆行列をPythonで取りたいと思います。私はnumpyを見てきましたが(行列反転はしますが、行列の逆行列はありません)、オンラインでいくつかの数学理論パッケージを見ましたが、この比較的一般的な手順を行うようなものは何もありませんでした(少なくとも、Pythonでモジュール式の行列反転を行う最も簡単な方法は?
ところで、上記の行列の逆数は、[[5,1]、[5,3]](mod 7)です。私はPythonでそれをやって欲しいです。
わかりました。気になる人には、私自身の問題を解決しました。私はしばらくかかったが、これはうまくいくと思う。これはおそらく、最もエレガントではない、といくつかのより多くのエラー処理を含める必要がありますが、それは動作します:
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
まだ完全ではありません。私はちょうどint(linalg.det(A))が常にあなたに正しい決定要因を与えるとは限りません。うーん、numpyの行列式アルゴリズムの大ファンではない。私が扱っている行列(今はちょうど小さな3x3行列)に対して、行列式は整数でなければなりません!なぜnumpyのアルゴリズムが間違っているのですか? – John
私は現在int(round(linalg.det(A)))を使用しています。キモい。しかし、私はそれが動作すると思います。 – John
共有してくれてありがとう、私は後のスタジアムで使うことができるように保存しました:)。 – bastijn
丸め誤差が問題でないときに動作ハックトリック:
によって国防省項目別でください。ここで私のコードは私が自分の目的で書いたガウス消去法を使っています(丸め誤差が私の問題でした)。 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:コメント者が指摘するように、このアルゴリズムは失敗することがあります。修正するのはやや難しいですが、今日は時間がありません。当時私のケースではランダム行列のために働いていました(モジュロは大きな素数の積です)。基本的には、最初の非ゼロエントリはモジュラスに対して相対的に素数ではないかもしれません。異なる行を検索してスワップすることができるので、プライムケースは簡単です。非プライム場合、私はそれがすべて有数のエントリは、あなたがそれらを結合する必要がありますので互いに素ではないということができると思い
コードをありがとう。私は今、ソリューションが必要なことを長い間残しています(去年の授業を受講しました)が、あなたの努力に感謝します。私はあなたの提案が本当に好きです - あなたの最初の提案に特に良い数学的な推論、。あなたのガウス消去ソリューションでは、あなたが提供したコードは、私が与えたコードとほぼ同じ量の作業ですが、よりエレガントなケースを作ることができます。いずれにせよ、素晴らしい仕事!質問に答える時間をとってくれてありがとう。 – John
私のコードでは丸めの問題はまったく発生しません。しかし、最初の "ハッキー"な提案はします。しかし問題はない! – WuTheFWasThat
この有益なモジュールをありがとう! – SAbbasizadeh
それは
ようセージ( www.sagemath.org)を用いて計算することができます。行列(IntegerModRing(7)、[[1,2]、[3,4]])。逆()
セージをインストールするには巨大であり、あなたは痛みでそれに付属しているのpythonのバージョンを使用する必要がありますが。
あなたはセージを見ましたか? – Alejandro
あなた自身の小さなコードを書くことになったら。私たちの多くが興味を持っていると思うので、ここでそれを共有することを検討してください:)。 – bastijn
モジュラ行列の逆変換はsympyに組み込まれています(この質問は尋ねられたので新しいかもしれません)。モジュール化された行の削減はかなり簡単です。 http://stackoverflow.com/a/37015283/2747370 –