5

Cで作業している数値ソルバーでは、2x2行列を反転する必要があり、右側で別の行列が乗算されます。

C = B . inv(A)

私は逆2x2行列の次の定義を使用しています。

a = A[0][0];
b = A[0][1];
c = A[1][0];
d = A[1][1];
invA[0][0] = d/(a*d-b*c);
invA[0][1] = -b/(a*d-b*c);
invA[1][0] = -c/(a*d-b*c);
invA[1][1] = a/(a*d-b*c);

私のソルバーの最初の数回の反復では、これは正しい答えを与えるように見えますが、いくつかのステップの後、物事は成長し始め、最終的に爆発します。

さて、SciPyを使用した実装と比較すると、同じ数学が爆発しないことがわかりました。私が見つけることができる唯一の違いは、SciPyコードがを使用していることですscipy.linalg.inv()。これは内部でLAPACKを使用して反転を実行します。

の呼び出しを上記の計算に置き換えるとinv()、Pythonバージョンが爆発するので、これが問題であると確信しています。計算のわずかな違いが忍び寄っています。これは数値的な問題であると私に信じさせます。反転操作にとってはまったく驚くべきことではありません。

数値の問題が問題にならないことを期待して、倍精度浮動小数点数（64ビット）を使用していますが、そうではないようです。

しかし、LAPACKのようなライブラリを呼び出さなくても、Cコードでこれを解決したいと思います。これは、純粋なCに移植する理由は、ターゲットシステムで実行するためです。また、ブラックボックスを呼び出すだけでなく、問題を理解したいと思います。最終的には、可能であれば単精度でも実行したいと思います。

だから、私の質問は、そのような小さな行列の場合、Aの逆行列を計算するための数値的により安定した方法がありますか？

ありがとう。

編集：現在、を解くことによって反転を回避できるかどうかを理解しようとしていCます。

4

5 に答える 5

6

行列を反転しないでください。ほとんどの場合、逆行列を使用して達成することは、行列を逆行列にしなくても、より高速かつ正確に実行できます。逆行列は本質的に不安定であり、それを浮動小数点数と混合すると問題が発生します。

と言うのは、 CC = B . inv(A)を解きたいと言うのと同じです。これは、各とを 2 つの列に分割することで実現できます。解いてC を生成します。AC = BBCA C1 = B1A C2 = B2

于 2011-12-31T22:25:29.720 に答える
5

あなたのコードは問題ありません。ただし、 4 つの減算のいずれかから精度が失われるリスクがあります。

matfunc.pyで使用されているような、より高度な手法の使用を検討してください。そのコードは、Householder リフレクションで実装されたQR 分解を使用して反転を実行します。結果は、反復改良によってさらに改善されます。

于 2011-12-31T20:50:05.020 に答える
4

行列式の計算は安定していません。より良い方法は、ここで明示的に簡単に解決できる、部分的なピボットで Gauss-Jordan を使用することです。

2x2 システムを解く

連立方程式を解いてみましょう (c, f = 1, 0 を使用し、次に c, f = 0, 1 を使用して逆数を取得します)。

a * x + b * y = c
d * x + e * y = f

擬似コードでは、これは読み取ります

if a == 0 and d == 0 then "singular"

if abs(a) >= abs(d):
    alpha <- d / a
    beta <- e - b * alpha
    if beta == 0 then "singular"
    gamma <- f - c * alpha
    y <- gamma / beta
    x <- (c - b * y) / a
else
    swap((a, b, c), (d, e, f))
    restart

これは、行列式 + comatrix (beta安定した方法で計算された行列式 * 何らかの定数) よりも安定しています。完全なピボットの同等物を計算することができます (つまり、x と y を潜在的に交換して、最初の除算がa、b、d、e の中で最大の数になるようにします)。状況によっては、これがより安定する場合がありますa。a、しかし、上記の方法は私にとってはうまく機能しています。

これは、LU 分解を実行することと同じです (この LU 分解を保存したい場合は、ガンマ、ベータ、a、b、c を保存します)。

QR 分解の計算は、明示的に行うこともできます (また、正しく行うと非常に安定します) が、遅くなります (また、平方根を取る必要があります)。選択はあなた次第です。

精度の向上

より高い精度が必要な場合 (上記の方法は安定していますが、固有値の比率に比例して丸め誤差が発生します)、「修正のために解く」ことができます。

実際、上記の方法でを解いA * x = bたとします。xここでを計算するA * xと、 と完全に等しくないことがわかりbます。わずかな誤差があります。

A * x - b = db

で を解くとdx、次のようになりA * dx = dbます。

A * (x - dx) = b + db - db - ddb = b - ddb

ここddbで、 は の数値解法によって生じる誤差でA * dx = db、通常は よりもはるかに小さいdb(dbは よりもはるかに小さいためb)。

上記の手順を繰り返すことができますが、通常、完全なマシン精度を復元するには 1 つの手順が必要です。

于 2012-04-17T09:09:45.933 に答える
1

ヤコビ法を使用します。これは、A の主対角のみを "反転" する反復法です。これは、行列全体を反転するよりも非常に簡単で、数値的に不安定になる傾向が少なくなります。

于 2012-01-02T20:15:49.027 に答える