3

大きな正方行列 (100-200 行/列) を頻繁に反転する必要があるコードを python で書いています。

私はマシンの精度の限界に達しているので、mpmath任意精度の行列反転を行うために使用しようとし始めましたが、 を使用しても非常に遅いgmpyです。

精度 30 (10 進数) でサイズ 20、30、60 のランダム行列を反転すると、~ 0.19、0.60、および 4.61 秒かかりますが、同じ操作にmathematicaは 0.0084、0.015、および 0.055 秒かかります。

これは、arch Linux マシンでpython3and mpmath 0.17(gmpy のバージョンは不明) を使用しています。mpmath が非常に遅い理由はわかりませんが、これに対して mathematica が管理する速度に近づくオープンソース ライブラリはありますか (1/2 の速度でも良いでしょう)。

任意の精度は必要ありません。おそらく 128 ビットで十分でしょう。また、mpmath がどのように非常に遅くなる可能性があるのか​​ もわかりません。非常に異なる行列反転アルゴリズムを使用している必要があります。具体的には、私が使用してM**-1います。

より高速なアルゴリズムを使用したり、高速化したりする方法はありますか。

4

3 に答える 3

2

残念ながら、mpmath の Linera 代数はかなり遅いです。この問題をより適切に解決する多くのライブラリがあります (たとえばSage )。とはいえ、Stuart の提案のフォローアップとして、固定小数点演算を使用して、ライブラリをインストールせずに Python でかなり高速な高精度の行列乗算を実行するのはかなり簡単です。入力と出力に mpmath 行列を使用するバージョンを次に示します。

def fixmul(A, B, prec):
    m = A.rows; p = B.rows; n = B.cols;
    A = [[A[i,j].to_fixed(prec) for j in range(p)] for i in range(m)]
    B = [[B[i,j].to_fixed(prec) for j in range(n)] for i in range(p)]
    C = [([0] * n) for r in range(m)]
    for i in range(m):
        for j in range(n):
            s = 0
            for k in range(p):
                s += A[i][k] * B[k][j]
            C[i][j] = s
    return mp.matrix(C) * mpf(2)**(-2*prec)

256 ビットの精度で、これは 2 つの 200x200 行列を mpmath よりも 16 倍高速に乗算します。この方法で逆行列ルーチンを直接記述することも難しくありません。もちろん、行列のエントリが非常に大きいか非常に小さい場合は、最初に再スケーリングする必要があります。より確実な解決策は、gmpyの浮動小数点型を使用して独自の行列関数を作成することです。これは本質的に高速である必要があります。

于 2013-03-11T06:42:58.490 に答える
1

倍精度は最終結果の精度の問題ではなく、逆の中間結果で問題を引き起こしている特定の行列では問題ないと思います。その場合、通常の numpy (倍精度) の逆数の結果を単なる適切な近似値として扱い、これをニュートン法を数回反復して逆数を解くための開始点として使用します。

Aを反転しようとしている行列とし、Xを逆行列の推定値とします。ニュートンの方法の反復は、単純に次のもので構成されます。

X = X*(2I - AX)

大きな行列の場合、上記の数回の反復を計算する労力は、逆行列を見つける労力と比較してほとんど重要ではなく、最終結果の精度が大幅に向上する可能性があります。これを試してみてください。

ところで、Iは上記の式の単位行列です。

編集して、浮動小数点型の精度をテストするコードを追加します。

このコードを使用して、float 型の精度をテストします。

x = float128('1.0')
wun = x
two = wun + wun
cnt = 1
while True:
   x = x/two
   y = wun + x
   if y<=wun: break
   cnt +=1

print 'The effective number of mantissa bits is', cnt
于 2013-03-10T16:21:26.417 に答える