問題タブ [numerical-stability]

私はPython用のウィルソンのスペクトル密度因数分解アルゴリズム[1]の実装を書き込もうとしています。このアルゴリズムは、[QxQ]行列関数を反復的に平方根に因数分解します（これは、スペクトル密度行列のニュートンラプソン平方根ファインダーの拡張のようなものです）。

問題は、私の実装がサイズ45x45以下の行列に対してのみ収束することです。したがって、20回の反復後、行列間の二乗差の合計は約2.45e-13になります。ただし、サイズ46x46の入力を行うと、100回程度の反復まで収束しません。47x47以上の場合、行列は収束しません。エラーは約100回の反復で100から1000の間で変動し、その後非常に急速に大きくなり始めます。

このようなものをデバッグするにはどうすればよいですか？気が狂うような特定のポイントはないようで、行列が大きすぎて実際に手動で計算を試みることはできません。誰かがこのような奇妙な数値のバグを見つけるためのヒント/チュートリアル/ヒューリスティックを持っていますか？

私はこれまでこのようなことを扱ったことがありませんが、あなたの何人かが持っていることを願っています...

ありがとう、-ダン

[1]GTウィルソン。「マトリックススペクトル密度の因数分解」。SIAMJ.Appl。数学（第23巻、第4号、1972年12月）

次のように用語文書行列(TD)から計算される共起行列(C)の因子分析を実行しようとしています: C=TD*TD'

理論的には C は正の半正定値である必要がありますが、そうではなく、因数分解アルゴリズムはこのため動作しません。速度上の理由からアルゴリズムを変更できません)。

私はそれを調べて、それは数値安定性の問題かもしれません: 半正定行列を生成するための単純なアルゴリズム - 答え2.

ここで進める良い方法は何ですか?

多くの数値アルゴリズムは、32/64 ビット浮動小数点で実行される傾向があります。

しかし、精度の低い (そして消費電力の少ない) コプロセッサにアクセスできたらどうでしょうか? では、数値アルゴリズムでどのように利用できるでしょうか?

これらの問題に対処する優れた本/記事を知っている人はいますか?

ありがとう！

私は Python でプログラムされたベイジアン分類器を持っています。問題は、機能の確率を乗算すると、2.5e-320 などの非常に小さな float 値が得られ、突然 0.0 に変わることです。MAX 値 (より大きい値) を返すクラスに基づいて「最適な」クラスを見つける必要があるため、0.0 は明らかに役に立ちません。

これに対処する最善の方法は何ですか？数値の指数部分 (-320) を見つけ、それが低すぎる場合は、その値に 1e20 またはそのような値を掛けることを考えました。しかし、もっと良い方法があるのではないでしょうか?

ln（1-x）を計算するアルゴリズムを探しています。xは小さい場合が多い（<0.01）が、大きい場合もあります。アルゴリズムは正確である必要があり、遅すぎないようにする必要があります。精度が低下する可能性があるため、ln（x）にライブラリを使用したくありません。

私はJavaでFuzzyC-Meansアルゴリズムのバージョンを実装しようとしています。また、一度だけ計算できるすべてのものを一度だけ計算することによって、いくつかの最適化を実行しようとしています。

これは反復アルゴリズムであり、マトリックスの更新に関して、ピクセルxクラスターメンバーシップマトリックスU（行の値の合計は1.0である必要があります）、これは最適化する更新ルールです。

ここで、xは行列の要素X（ピクセルx特徴）であり、vは行列V （クラスターx特徴）に属します。Andは、からまでmの範囲のパラメーターであり、クラスターの数です。使用される距離はユークリッドノルムです。1.1infinityc

この式を平凡な方法で実装する必要がある場合は、次のようにします。

Xこのようにして、いくつかの最適化がすでに行われ、との間のすべての可能な二乗距離を事前に計算し、Vそれらを行列に格納しましたが、 2回Dの要素を循環させて、2つのネストされたループを生成するため、それだけでは不十分です。V式を見ると、分数の分子は合計に依存しないため、分子と分母を個別に計算でき、分母はピクセルごとに1回だけ計算できます。だから私はこのような解決策に到達しました：

D距離を計算しているときに、分母を行列の追加の列に事前計算した場合：

そうすることで、「banal」計算と2番目の計算の間に数値の違いが生じ、U（最初の反復ではなく、すぐに）の数値が異なります。問題は、非常に小さい数を高い値にべき乗することです（の要素はU0.0から1.0の範囲でありexp、の場合m = 1.1、はです10）は非常に小さい値になりますが、分子と分母を除算して結果をべき乗すると、数値的に良くなるために。問題は、それがはるかに多くの操作を伴うことです。

アップデート

私が得るいくつかの値ITERATION 0：

Dこれは、最適化されていない行列の最初の行です。

384.6632 44482.727 17379.088 1245.4205

これは、最適化された方法でマトリックスの最初の行ですD（最後の値は事前に計算された分母であることに注意してください）。

384.6657 44482.7215 17379.0847 1245.4225 1.4098E-26

Uこれは、最適化されていない最初の行です。

0.99999213 2.3382613E-21 2.8218658E-17 7.900302E-6

これは、U最適化された最初の行です。

0.9999921 2.338395E-21 2.822035E-17 7.900674E-6

ITERATION 1：

Dこれは、最適化されていない行列の最初の行です。

414.3861 44469.39 17300.092 1197.7633

これは、最適化された方法でマトリックスの最初の行ですD（最後の値は事前に計算された分母であることに注意してください）。

414.3880 44469.38 17300.090 1197.7657 2.0796E-26

Uこれは、最適化されていない最初の行です。

0.99997544 4.9366603E-21 6.216704E-17 2.4565863E-5

これは、U最適化された最初の行です。

0.3220644 1.5900239E-21 2.0023086E-17 7.912171E-6

最後の値のセットは、伝播されたエラー（私はまだいくつかの間違いをしていることを願っています）のためにそれらが非常に異なっており、それらの値の合計が1.0でなければならないという制約さえ違反していることを示しています。

私は何か間違ったことをしていますか？コードを最適化し、数値的に安定させるための可能な解決策はありますか？任意の提案や批判をいただければ幸いです。

私はアルゴリズムの過程でいくつかの線形プログラミングの演習を行っています。これを行う際に、分数を使用して多くの操作を手動で解決しています。これを行うことで、人間は数値の不安定性に悩まされていないことに気付きました。値を分数表現で保持し、最終的に式の値を (おそらく電卓を使用して) 評価します。

これを自動的に行う技術はありますか？

ある種のシンボリック計算を実現し、数値を内部的に単純化し、最終的に式の評価中にのみ値を生成するものを考えています。

非常に小さい、ほぼゼロの数値を生成しているように見える数学 (C++) があります (trig 関数呼び出しが私の本当の問題であると思われます) が、これらのケースを検出して、もっと詳しく。

現在、以下を試していますが、正しいですか？

第二に、数学の性質は本質的に三角関数です（別名、多くのラジアン/度変換および//呼び出しsinなどを使用します）。数学的なエラーを回避するには、どのような変換を行うことができますか?costan

明らかに、乗算には対数変換を使用できます-他に何がありますか?

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

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

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

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

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

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

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

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

ありがとう。

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

2D のベジエ曲線 (P0、P1、P2、P3) のポイントが与えられた場合、特定の x 座標の y 座標を見つけたいと思います。この問題は、次の制限により明確に定義されています。

• P0 = (0,0)、P3 = (1,1)
• P1 = (t, 1-t) for t for 0, 1
• P2 = 1 - P1 (x と y)

上記のすべての制限をベジエ曲線の式CubicBezier.htmlに入れて、答えを計算する次の関数があります。私は Newton-Raphson を使用して、必要なポイントのパラメーターを計算しています。(定義された許容範囲内で) 終了するまでループを終了させないため、これが機能することがわかっています。

この関数を使用して、画像にコントラスト フィルターを適用しています。この 0.5 では同じ画像が返され、0.0 ではコントラストが最大に減少し、1.0 では最大に増加します。

編集次の関数は修正され、完全に機能するようになりました。

p = 0.5 に設定すると、直線が得られるはずですが、linspace に対してこれを行って点をプロットすると、0.5 と 1.0 の間で曲がります。なぜこれが起こっているのか、誰にもわかりますか?何かできることがあれば教えてください。