気付いていないかもしれませんが、あなたの質問はおそらくプログラミングの質問というよりも線形代数の質問ですが、プログラミングの側面があります。
この回答では、実際の問題は、図に示した N = 8 のスケッチよりもはるかに大きくなる可能性があると想定しています。数学的な議論が最初に来ます。サンプルコードは次のとおりです。
線形代数の問題の性質は、一般的な線形代数の手法によって実際に処理できるか、そうでないかのいずれかであると思われます。LU 因数分解や特異値分解などの一般的な線形代数手法は、1960 年代から徹底的に開発され、確固たる理論的基盤を持っています。問題は、そのような一般的な手法は計算処理で O(N*N*N) になる傾向があることです。つまり、ベクトルを 10 倍長くすると、処理に 1000 倍のコンピューター時間が必要になります。実用的な最大サイズを超えると、問題は計算上扱いにくくなります (少なくとも一般的な線形代数の手法では)。
問題が非常に大きい場合
幸いなことに、O(N*N*N) 問題を O(N*N) 問題に縮小するための当惑するほどさまざまな専用手法が存在します。これらの技術は主に 1960 年代以降に出現し、活発な研究が行われている分野です。残念ながら、どの手法をどのように適用するかを知るには、問題を深く理解し、行列数学をかなりよく理解している必要があります。O(N*N*N) では、一般的な手法が知られています。O(N*N) では、一般的な手法はなく、限られた種類のシステムで動作する多くの専用の手法のみがあります。十分に長く作業すれば、適切な制限を見つけることができますが、行列を汎用ソルバーにダンプするほど単純ではありません。O(N*N) ではありません。
このテーマについて私が読んだ中で最も優れた本は、Henk A. van der Vorst のIterative Krylov Methods for Large Linear Systems です。 幸いなことに、本の価格は妥当です。もちろん、ファン デル フォールストに取り組む前に、一般的な線形代数の手法をしっかりと理解する必要があります。Joel N. Franklin の 1968 年の短い本Matrix Theoryをお勧めします。Dover から 1993 年に安価なペーパーバックで入手できます。(Van der Vorst は、Kapp と Brown の巧妙で単純で、しばしば効果的な順序付き複数相互作用の方法について言及していません。これは、私が繰り返し有用であると感じているので、必要に応じて調べられるように、ここで言及させてください。開示: 私はブラウンを個人的に知っており、この理由で部分的である可能性があります; それでも、彼とカップ.
理由はよくわかりませんが、最近の手法のほとんどは (Kapp や Brown のものではありませんが) 実数で機能することを好み、複素数を特殊なケースとして扱うか、完全に無視しているようです。数値が複雑かどうかについては言及していませんが、複雑な場合は、オプションが多少制限される可能性があります。
フランクリンとファン デル フォルストの中間レベルとして、後者に取り組む時間や興味がない場合は、少なくとも調べて、共役勾配の方法に精通する必要があります。
問題が適度に大きい場合
問題が適度に大きい場合 (たとえば、N < 20,000 程度) の場合、作業ははるかに簡単になります。この場合、van der Vorst のことは忘れてください。彼は必要ありません。答えがミリ秒、秒、分、または時間のいずれであるかに応じて (どれについては言及していませんが、N の実際の制限にのみ影響します)、O(N*N*N) のパフォーマンスを許容できます。フランクリンの 1960 年代からの一般的なテクニックは非常に堅牢です。
高速、効率的、完全かつ正確な LAPACK ライブラリとその低レベル ヘルパー BLAS は、この分野の標準です。それらを使用する必要があります。
LAPACK、BLAS、および C++
私の同僚と私はそれぞれ、LAPACK と BLAS へのさまざまな C および C++ インターフェイスを試しました。さまざまな理由から、役立つように努めていますが、完全に満足できるものはありません。私たちのそれぞれは、最終的にインターフェイスを完全にスキップし、LAPACK と BLAS を直接使用することにしました。これは私があなたに勧める傾向があるものです。
LAPACK と BLAS は、C++ からではなく、FORTRAN 77 から呼び出されることを意図していました。それでも、C++ からそれらを呼び出すことはそれほど難しくなく、それを行うために C++ インターフェイスは必要ありません。実際、C++ インターフェイスは必要ないでしょう。少なくとも、他の誰かが作成した一般的なインターフェイスは必要ありません。LAPACK と BLAS は、1 つにまとめない方が幸せになれます。(FORTRAN 77 では、リンク機能は有効ではありましたが制限されていたため、C スタイルのヘッダー ファイルがなかったことを思い出してください。)
LAPACK と BLAS を直接使用するために最初に知っておかなければならないことはこれです。これを理解するまでは悲惨なことになるでしょう: 行列は列優先で格納されます。 つまり、行列の各行ではなく、各列がメモリ内で 1 つの単位として順番に表されます。したがって、行列要素は、その左右の要素の間ではなく、すぐ上と下の要素の間に格納されます。実際、LAPACK と BLAS はこの方法で行列を格納するのが賢明です。なぜなら、C を発明した Kernighan と Ritchie は線形代数にあまり興味を持っていなかったからです。C は優れた言語ですが、C のデフォルトの行列格納規則はおそらく最初から間違いでした。数学者の観点から見ると、LAPACK と BLAS は行列を列優先で格納する必要があります。また、線形代数コードを作成している場合は、行列もこの方法で保存する必要があります。C のデフォルトの行優先の規則を無視します。これは、行列数学の規則とは何の関係もありません。
完全な例を次に示します。
#include <iostream>
extern "C" {
void dgesv_(
const int *N,
const int *NRHS,
double *A,
const int *LDA,
int *IPIV,
double *B,
const int *LDB,
int *INFO
);
}
int main()
{
const int N = 2;
// Let A = | 3.0 6.1 |
// |-1.2 1.7 |
double A[N*N] = {3.0, -1.2, 6.1, 1.7};
// Let b = | 5.5 |
// | 0.4 |
double x[N] = {5.5, 0.4};
// (Why is b named x? Answer: because b and x share
// the same storage, because Lapack will write x over b.)
// Solve the linear system A*x = b.
const int NRHS = 1;
const int LDA = N;
int IPIV [N];
const int LDB = N;
int INFO;
dgesv_(&N, &NRHS, A, &LDA, IPIV, x, &LDB, &INFO);
// Uncomment the next line if you wish to see
// the error code Lapack returns.
//std::cout << "INFO == " << INFO << "\n";
// Output x.
for (int i = 0; i < N; ++i) std::cout << x[i] << "\n";
return 0;
}
[関数の名前が dgesv_() なのはなぜですか? 回答: 少なくとも DGESV という名前の FORTRAN 77 ではそうではありません。ただし、FORTRAN 77 では大文字と小文字が区別されないため、これを C のリンク規則に変換すると、dgesv_() が得られます。少なくとも、それは私の同僚と私が試したすべての Linux、BSD、または OSX マシンで取得したものです。Debian または Ubuntu では、シェル コマンド「readelf --symbols /usr/lib/liblapack.so | grep -i DGESV」がこれを検出します。Microsoft プラットフォームでは、C リンクされたシンボルに別の名前が付いている場合があります。これが自分に関係があるかどうかを調査する必要があります。]
LAPACK と BLAS は、その機能が非常に優れています。それらを使用する必要があります。
データを C スタイルの配列 (例のように) に保存するか、C++ スタイルの std::valarrays に保存するかはあなた次第です。ストレージが列優先である限り、LAPACK と BLAS は気にしません。
係数の最小化
係数で何をしたいのかを正確に伝えるのに十分な情報が提供されていませんが、必要なのは N 行 (2*N) 列の劣決定システムの解であると思われます。そのようなことはフランクリンの本の範囲を超えていますが、フランクリンの資料をすでに吸収している場合は、疑似逆数に関する私自身のメモが役立つかもしれません.
明らかに、迅速な解決策を探しているなら、私にはありません。あなたが望むことを正確に実行する既製のソフトウェアパッケージがあるかもしれませんが、私の経験では、これはあなたに不利であることを示唆しています. 事前に用意されたソフトウェアがクラックするには、実際に発生する何百もの異なる種類の行列問題が多すぎます。ただし、LAPACK と BLAS は、Franklin、van der Vorst、および現在読んでいる回答とともに、特定の問題を解決するために必要なすべてのツールを提供します。
幸運を。