固定次元 (N=9) (行列は対称、半正定値) の密な線形システムの高速なソリューションに推奨するアルゴリズムはどれですか?
型は 32 ビットと 64 ビットの浮動小数点です。
このようなシステムは何百万回も解決されるため、アルゴリズムは次元 (n=9) に関してかなり高速である必要があります。
提案されたアルゴリズムの堅牢なC++ 実装の PS の例は高く評価されます。
1) 「何百万回も解決した」とはどういう意味ですか? 100 万の異なる右辺項を持つ同じ係数行列、または 100 万の異なる行列ですか?
何百万もの異なるマトリックス。
2) 正の _semi_defined は、行列が (マシンの精度で) 特異になる可能性があることを意味します。この場合、どのように対処したいですか？エラーを発生させるだけですか、それとも賢明な回答を返そうとしますか?
上げ間違いOK。
行列は対称で正半有限であるため、コレスキー分解はLU分解よりも厳密に優れています。（行列のサイズに関係なく、LUの約2倍高速です。出典：TrefethenとBauによる「NumericalLinearAlgebra」）
これは、事実上、小さな密な行列の標準でもあります（出典：計算数学で博士号を取得します）システムが十分に大きくならない限り、反復法は直接法よりも効率が低くなります（簡単な経験則では意味がありませんが、それは常に素晴らしいことです）持っている：最近のコンピューターでは、100 * 100より小さい行列は間違いなく小さな行列であり、反復法ではなく直接法が必要です）
今、私はそれを自分で行うことをお勧めしません。徹底的にテストされた優れたライブラリがたくさんあります。しかし、私があなたに1つを推薦しなければならないなら、それはEigenでしょう：
ちなみに、ここのドキュメントには、7つの直接線形ソルバーのさまざまな長所と短所が簡潔な表にまとめられています。あなたの場合、LDLT（コレスキーのバリエーション）が勝つようです
一般に、高速で安定した数値実装を追求するには面倒な詳細が数多くあるため、独自のアプローチではなく、既存のライブラリを使用することをお勧めします。
開始するためのいくつかを次に示します。
Eigen ライブラリ (個人的な好み):
http://eigen.tuxfamily.org/dox/QuickRefPage.html#QuickRef_Headers
アルマジロ: http://arma.sourceforge.net/
周りを検索すると、他にもたくさん見つかります。
特に「何百万回も解いた」が本当に「一度解いて何百万ものベクトルに適用した」ことを意味する場合は、LU 分解をお勧めします。LU 分解を作成して保存し、希望する数の rhs ベクトルに対して前方後方置換を適用します。
ピボットを使用すると、ラウンドオフに直面してもより安定します。
対称半正定行列の LU はあまり意味がありません。不必要な操作を実行して、入力データの優れたプロパティを破棄します。
LLT と LDLT のどちらを選択するかは、行列の条件数と、エッジ ケースをどのように処理するかによって異なります。LDLT は、統計的に有意な精度の向上を証明できる場合、または堅牢性がアプリにとって最も重要である場合にのみ使用してください。
(行列のサンプルがなければ適切なアドバイスをするのは難しいですが、このような小さな次数 N=9 では、小さな対角項を D の下部に向かってピボットする必要はないと思います。したがって、古典から始めます。適切に選択された許容誤差に対して diag 項が小さくなった場合は、Cholesky を実行し、単純に因数分解を中止します)。
Cholesky は非常に簡単にコーディングできます。非常に高速なコードを作成しようとする場合は、自分で実装することをお勧めします。
上記の他の人たちと同様に、コレスキーをお勧めします。加算、減算、メモリ アクセスの回数が増えるということは、LDLt が cholesky よりも遅いということです。
実際、コレスキーにはさまざまなバリエーションがあり、どのバリエーションが最も高速かは、行列に選択する表現によって異なります。私は通常、Fortran スタイルの表現を使用します。つまり、行列 M は double* M で、M(i,j) は m[i+dim*j] です。このため、上三角コレスキーが (少し) 最も速いと思います。つまり、U'*U = M で上三角 U を探します。
固定された小さい次元の場合、ループを使用しないバージョンを作成することを検討する価値があります。これを行うための比較的簡単な方法は、それを行うプログラムを作成することです。私が思い出したように、テンプレートとして一般的なケースを処理するルーチンを使用して、特定の固定次元バージョンを作成するプログラムを作成するのに朝しかかかりませんでした。節約はかなりのものになる可能性があります。たとえば、私の一般的なバージョンでは 100 万回の 9x9 因数分解を実行するのに 0.47 秒かかりますが、ループレス バージョンでは 0.17 秒かかります。これらのタイミングは、2.6 GHz の PC でシングル スレッドで実行されます。
これが主要な作業ではないことを示すために、そのようなプログラムのソースを以下に示します。コメントとして因数分解の一般的なバージョンが含まれています。行列が特異に近くない状況でこのコードを使用しましたが、そこでは問題なく動作すると思います。ただし、より繊細な作業には粗雑すぎる可能性があります。
/* ----------------------------------------------------------------
** to write fixed dimension ut cholesky routines
** ----------------------------------------------------------------
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <strings.h>
/* ----------------------------------------------------------------
*/
#if 0
static inline double vec_dot_1_1( int dim, const double* x, const double* y)
{
double d = 0.0;
while( --dim >= 0)
{ d += *x++ * *y++;
}
return d;
}
/* ----------------------------------------------------------------
** ut cholesky: solve U'*U = P for ut U in P (only ut of P accessed)
** ----------------------------------------------------------------
*/
int mat_ut_cholesky( int dim, double* P)
{
int i, j;
double d;
double* Ucoli;
for( Ucoli=P, i=0; i<dim; ++i, Ucoli+=dim)
{ /* U[i,i] = P[i,i] - Sum{ k<i | U[k,i]*U[k,i]} */
d = Ucoli[i] - vec_dot_1_1( i, Ucoli, Ucoli);
if ( d < 0.0)
{ return 0;
}
Ucoli[i] = sqrt( d);
d = 1.0/Ucoli[i];
for( j=i+1; j<dim; ++j)
{ /* U[i,j] = (P[i,j] - Sum{ k<i | U[k,i]*U[k,j]})/U[i,i] */
P[i+j*dim] = d*(P[i+j*dim] - vec_dot_1_1( i, Ucoli, P+j*dim));
}
}
return 1;
}
/* ----------------------------------------------------------------
*/
#endif
/* ----------------------------------------------------------------
**
** ----------------------------------------------------------------
*/
static void write_ut_inner_step( int dim, int i, int off)
{
int j, k, l;
printf( "\td = 1.0/P[%d];\n", i+off);
for( j=i+1; j<dim; ++j)
{ k = i+j*dim;
printf( "\tP[%d] = d * ", k);
if ( i)
{ printf( "(P[%d]", k);
printf( " - (P[%d]*P[%d]", off, j*dim);
for( l=1; l<i; ++l)
{ printf( " + P[%d]*P[%d]", l+off, l+j*dim);
}
printf( "));");
}
else
{ printf( "P[%d];", k);
}
printf( "\n");
}
}
static void write_dot( int n, int off)
{
int i;
printf( "P[%d]*P[%d]", off, off);
for( i=1; i<n; ++i)
{ printf( "+P[%d]*P[%d]", off+i, off+i);
}
}
static void write_ut_outer_step( int dim, int i, int off)
{
printf( "\td = P[%d]", off+i);
if ( i)
{ printf( " - (");
write_dot( i, off);
printf( ")");
}
printf( ";\n");
printf( "\tif ( d <= 0.0)\n");
printf( "\t{\treturn 0;\n");
printf( "\t}\n");
printf( "\tP[%d] = sqrt( d);\n", i+off);
if ( i < dim-1)
{ write_ut_inner_step( dim, i, off);
}
}
static void write_ut_chol( int dim)
{
int i;
int off=0;
printf( "int\tut_chol_%.2d( double* P)\n", dim);
printf( "{\n");
printf( "double\td;\n");
for( i=0; i<dim; ++i)
{ write_ut_outer_step( dim, i, off);
printf( "\n");
off += dim;
}
printf( "\treturn 1;\n");
printf( "}\n");
}
/* ----------------------------------------------------------------
*/
/* ----------------------------------------------------------------
**
** ----------------------------------------------------------------
*/
static int read_args( int* dim, int argc, char** argv)
{
while( argc)
{ if ( strcmp( *argv, "-h") == 0)
{ return 0;
}
else if (strcmp( *argv, "-d") == 0)
{ --argc; ++argv;
*dim = atoi( (--argc, *argv++));
}
else
{ break;
}
}
return 1;
}
int main( int argc, char** argv)
{
int dim = 9;
if( read_args( &dim, --argc, ++argv))
{ write_ut_chol( dim);
}
else
{ fprintf( stderr, "usage: wchol (-d dim)? -- writes to stdout\n");
}
return EXIT_SUCCESS;
}
/* ----------------------------------------------------------------
*/
使いやすいため、比較のためだけに固有値ソルバーを使用できます。特定のユースケースでは、特定のソルバーがより高速である可能性がありますが、別のソルバーの方が優れているはずです。そのために、選択のためだけに各アルゴリズムの実行時間を測定できます。その後、目的のオプションを実装できます (または、ニーズに最適な既存のオプションを見つけます)。