7

私はn実数変数を持っています（わからない、本当に気にしないでください）、それらを呼びましょうX[n]。私はm >> nそれらの間に関係もあります。それらR[m]を次の形式で呼びましょう。

X[i] = alpha*X[j]は、alphaゼロ以外の正の実数でiありj、別個のものですが、(i, j)ペアは必ずしも一意ではありません（つまり、異なるアルファ係数を持つ同じ変数間に2つの関係が存在する可能性があります）

私がやろうとしているのは、alpha最小二乗の意味で過剰決定系を解決する一連のパラメーターを見つけることです。理想的な解決策は、各方程式パラメーターとその選択された値の間の差の2乗和を最小化することですが、次の近似に満足しています。

m個の方程式をn個の未知数の過剰決定系に変換すると、疑似逆行列ベースの数値ソルバーは明白な解（すべてゼロ）を与えます。したがって、現在私が行っているのは、別の方程式をミックスに追加しx[0] = 1（実際には任意の定数で可能です）、ムーア・ペンローズ疑似逆行列を使用して、生成されたシステムを最小二乗の意味で解きます。(x[0] - 1)^2これはの合計との二乗和を最小化しようとしx[i] - alpha*x[j]ますが、これは私の問題に対する適切で数値的に安定した近似であることがわかります。次に例を示します。

a = 1
a = 2*b
b = 3*c
a = 5*c

オクターブ：

A = [
  1  0  0;
  1 -2  0;
  0  1 -3;
  1  0 -5;
]

B = [1; 0; 0; 0]

C = pinv(A) * B or better yet:
C = pinv(A)(:,1)

これにより、、、の値が得られます。これによりa、 次の（合理的な）関係が得られます。bc[0.99383; 0.51235; 0.19136]

a = 1.9398*b
b = 2.6774*c
a = 5.1935*c

したがって、今はこれをC / C ++ / Javaで実装する必要があり、次の質問があります。

問題を解決するためのより高速な方法はありますか、それとも過剰決定系を生成して疑似逆行列を計算することで正しい方向に進んでいますか？

私の現在のソリューションでは、特異値分解と3つの行列乗算が必要です。これは、5000または10000になる可能性があることを少し考慮していmます。疑似逆行列を計算するより高速な方法はありますか（実際には、最初の列のみが必要であり、行列のスパース性を考慮して、Bがゼロである場合の行列全体（最初の行を除く）（各行にはゼロ以外の値が2つ含まれ、一方は常に1で、もう一方は常に負です）

このためにどの数学ライブラリを使用することをお勧めしますか？LAPACKは大丈夫ですか？

数値的に安定していて漸近的に高速であるという条件で、他の提案も受け入れます（たとえばk*n^2、k大きくなる可能性があります）。

4

2 に答える 2

4

あなたの問題は不適切です。問題をn個の変数の関数（差の最小二乗）として扱う場合、関数には正確に1つのグローバル最小超平面があります。

変数の1つを非ゼロに修正するか、他の方法で関数ドメインを縮小しない限り、そのグローバル最小値には常にゼロが含まれます。

ソリューション超平面のパラメーター化が必要な場合は、Moore-Penrose Pseudo-Inverse http://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_pseudoinverseから取得して、すべての取得に関するセクションを確認できます。ソリューション。

（私は技術的に間違った方法で「超平面」という言葉を使用していることに注意してください。私はあなたのパラメータ空間の「平らな」無制限のサブセットを意味します...線、平面、1つまたは複数のベクトルによってパラメータ化できるもの。何らかの理由で、そのようなオブジェクトの一般的な名詞が見つかりません）

于 2011-08-19T09:32:36.813 に答える
3

SVDアプローチは数値的に非常に安定していますが、それほど高速ではありません。SVDを使用する場合、LAPACKは使用するのに適したライブラリです。1回限りの計算であれば、おそらく十分に高速です。

大幅に高速なアルゴリズムが必要な場合は、安定性を犠牲にする必要があります。1つの可能性は、QR分解を使用することです。詳細を確認するには、これを読む必要がありますが、理由の一部は次のようになります。AP = QR（Pは置換行列、Qは直交行列、Rは三角行列）がAの経済QR分解である場合、方程式AX=BはQRP^{-1} X=Bになります。解はX=PR ^ {-1} Q ^ T Bです。次のオクターブコードは、コードと同じAとBを使用してこれを示しています。

[Q,R,P] = qr(A,0)
C(P) = R \ (Q' * B)

これの良いところは、スパースQR分解を実行することでAのスパース性を利用できることです。qr関数のOctaveヘルプにいくつかの説明がありますが、すぐには機能しませんでした。

さらに高速な（ただし安定性がさらに低い）のは、通常の方程式を使用することです。AX = Bの場合、A ^ TAX = A ^TBです。行列A^TAは（うまくいけば）フルランクの正方行列なので、次のように使用できます。線形方程式のソルバー。Octaveコード：

C = (A' * A) \ (A' * B)

この場合も、このアプローチではスパース性を利用できます。スパース線形システムを解くための多くの方法とライブラリがあります。人気があるのはUMFPACKのようです。

後で追加：このフィールドについて、定量化するのに十分な知識がありません。これについては本全体が書かれています。おそらく、QRは約3倍または5倍速いSVDであり、通常の方程式は再び2倍速くなります。数値安定性への影響は行列Aに依存します。スパースアルゴリズムははるかに高速です（たとえばm倍）が、計算コストと数値安定性は問題に大きく依存し、よく理解されていない場合があります。

あなたのユースケースでは、SVDを使用してソリューションを計算し、所要時間を確認し、それが許容できる場合はそれを使用することをお勧めします（n=1000およびm=10000の場合は約1分になると思います） ）。さらに詳しく調べたい場合は、QRと通常の方程式も試して、どれだけ速く、どれだけ正確かを確認してください。それらがSVDとほぼ同じソリューションを提供する場合、あなたはそれらがあなたの目的のために十分に正確であるとかなり確信することができます。これらがすべて遅すぎて、ある程度の時間を費やしても構わないと思っている場合にのみ、スパースアルゴリズムを調べてください。

于 2011-08-20T02:56:21.377 に答える