固定寸法の高密度線形システム（N = 9）、対称、正定 - 半定形の高速ソリューション

Question

固定寸法の高密度線形システム（N = 9）の高速ソリューションに推奨するアルゴリズムはどれですか（マトリックスは対称、 semidefinite）？

• ガウスの消去
• LU分解
• コレスキー分解
• など？

タイプは32ビットと64ビットの浮動小数点です。

このようなシステムは何百万回も解決されるので、アルゴリズムは次元（n = 9）に関してかなり高速でなければなりません。

P.S. の堅牢な例提案されたアルゴリズムのためのC++の実装は評価されています。

1）「何百万回も解決した」とはどういう意味ですか？何百万もの異なる右辺用語を持つ同じ係数行列、または何百万もの異なる行列？

百万の異なるマトリックス。

2）正の_semi_definiteは、行列が（機械精度に）特異である可能性があることを意味します。どのようにこのケースに対処したいですか？ちょうどエラーを起こすか、または賢明な答えを返そうとしますか？

エラーが発生してもOKです。

Answer 1

行列は、対称、正半正定値、コレスキー分解はLU分解よりも厳密に優れています。 （行列のサイズに関係なく、LUより約2倍高速です。出典：TrefethenとBauによる「Numerical Linear Algebra」

小規模の密度行列の標準でもあります（出典：私は計算上の数学で博士号を取得しています）。反復法は、 100×100より小さい行列は、反復的なものではなく、直接的な方法を必要とする小さな行列です）

今、私はあなた自身でそれをすることをお勧めしません。徹底的にテストされた数多くの良いライブラリがあります。しかし、私はあなたに1つをお勧めする必要がある場合、それはEigen次のようになります（

• （ヘッダのみのライブラリ、唯一<の#include、リンクするので、ノーライブラリ>）必要ありませんインストール
• 堅牢かつ効率的に彼らが持っていますちなみに

を使用し、十分に文書化するには、メインページ上のベンチマークの多くは、結果いいです）

簡単、here in the documentation、あなたはAでの7つの直接線形ソルバーの様々な長所と短所を持っていますいい、簡潔なテーブル。あなたのケースでは、LDLT（Choleskyのバリエーション）が勝ったようです。

Answer 2

一般的に、高速で安定した数値実装を追求するには、多くの退屈な細部があるので、独自のアプローチではなく、既存のライブラリを使用するのが最適です。 http://arma.sourceforge.net/

検索の周りに、あなたは他の人をたくさん見つけることができます：
http://eigen.tuxfamily.org/dox/QuickRefPage.html#QuickRef_Headers

アルマジロ：

固有のライブラリ（私の個人的な好み）：

ここでは、始めるのは少数です。

Answer 3

私はLU分解をお勧めします。特に、「何百万回も解決された」というのは、「一度解決して数百万のベクトルに適用する」ことを意味します。 LU分解を作成して保存し、多くのr.h.sに対してフォワードバック置換を適用します。あなたが望むようにベクトル。

ピボットを使用すると、丸めの面でより安定します。

Answer 4

対称半限定マトリックスのLUはあまり意味がありません。入力データの素敵なプロパティを破壊して不要な操作を実行します。

LLTまたはLDLTの選択は、実際にはマトリックスの条件数とエッジケースの扱い方によって異なります。 LDLTは、統計的に有意な精度の向上が見られた場合、またはアプリケーションにとって堅牢性が最も重要な場合にのみ使用してください。

（あなたの行列のサンプルがなければ、助言を与えるのは難しいですが、N = 9のような小さな順序で、小さな対角項をDの底部にピボットすることは本当に必要ではないと思います。古典的なコレスキーから始まり、慎重に選択された許容誤差に関してdiagの項が小さくなるならば、単に分解を中止します。）Choleskyは非常にシンプルでコードが簡単です。自分で実装してください。

Answer 5

他のものと同様に、私はコレスキーをお勧めします。私は、加算、減算、およびメモリアクセスの数が増えると、LDLtがコレスキーよりも遅いことが分かりました。

実際には、コレスキーにはいくつかのバリエーションがあります。どちらが最速になるかは、使用するマトリックスの表現によって異なります。私は一般的に行列MであるM（i、j）がm [i + dim * j]であるdouble * Mです。これについては、上部の三角形のコレスキーが最も速いと考えています。つまり、U '* U = Mの上三角形のUを探します。

固定小規模な次元については、ループを使用しません。これを行うための比較的簡単な方法は、それを行うプログラムを書くことです。私が思い出しているように、一般的なケースをテンプレートとして扱うルーチンを使用すると、特定の固定ディメンションバージョンを書き込むプログラムを書くのは朝のみでした。節約はかなりの可能性があります。たとえば、私の一般的なバージョンでは9倍の9倍分解を行うのに0.47秒かかりますが、無能なバージョンでは0.17秒かかります。これらのタイミングは2.6GHzの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; 
} 
/* ---------------------------------------------------------------- 
*/ 

Answer 6

使いやすさのため、比較のためにEigenソルバを使用することができます。特定のユースケースでは、特定のソルバーは速くなるかもしれませんが、別のソルバーはより良いと思われます。そのためには、選択したアルゴリズムだけのランタイムを測定することができます。その後、希望のオプションを実装することができます（または、ニーズに最適な既存のオプションを見つけることができます）。