定维(N=9）对称半正定稠密线性系统的快速解

矩阵是对称的，半正定的，Cholesky分解严格优于LU分解。(无论矩阵大小，大约是LU的两倍。来源：Trefethen和Bau的"数值线性代数")

事实上，它也是小稠密矩阵的标准(来源：我是计算数学博士)迭代方法的效率不如直接方法，除非系统变得足够大(快速经验法则，这毫无意义，但这总是很好的：在任何现代计算机上，任何小于100*100的矩阵都绝对是一个小矩阵，需要直接方法，而不是迭代方法)

现在，我不建议你自己做。有很多好的图书馆都经过了彻底的测试。但如果我必须向你推荐一个，那就是Eigen：

• 不需要安装(只有头的库，所以没有要链接的库，只有#include<>)
• 稳健高效(他们在主页上有很多基准测试，结果很好)
• 易于使用且记录良好

顺便说一句，在文档中，您可以在一个简洁的表格中了解他们的7个直接线性求解器的各种优点和缺点。似乎在你的情况下，LDLT(Cholesky的变体)赢得了

通常，最好使用现有的库，而不是自己滚动的方法，因为要实现快速、稳定的数字实现，需要注意许多繁琐的细节。

以下是一些让你开始的方法：

特征库(我个人偏好)：
http://eigen.tuxfamily.org/dox/QuickRefPage.html#QuickRef_Headers

Armadillo：http://arma.sourceforge.net/

四处搜索，你会发现很多其他的。

我建议LU分解，特别是如果"求解数百万次"实际上意味着"求解一次并应用于数百万个向量"。您将创建LU分解，保存它，并根据需要对任意多个r.h.s.向量应用前向-后向替换。

如果你使用旋转，它在面对圆角时会更稳定。

对称半定矩阵的LU没有多大意义：您破坏了输入数据的一个良好属性，执行了不必要的操作。

LLT和LDLT之间的选择实际上取决于矩阵的条件数，以及如何处理边缘情况。只有当您能够证明准确性在统计上有显著提高，或者稳健性对您的应用程序至关重要时，才应使用LDLT。

(如果没有矩阵的样本，很难给出合理的建议，但我怀疑，对于这样一个小阶N=9，将小对角项向D的底部旋转真的没有必要。因此，我将从经典的Cholesky开始，如果对角项相对于一些合理选择的容差变得太小，则简单地中止因子分解。)

Cholesky的代码非常简单，如果你努力获得一个真正快速的代码，最好自己实现它。

和上面的其他人一样，我推荐cholesky。我发现加法、减法和内存访问次数的增加意味着LDLt比cholesky慢。

事实上，cholesky有很多变体，哪种变体最快取决于您为矩阵选择的表示形式。我通常使用fortran风格的表示，即矩阵M是双*M，其中M(I，j)是M[I+dim*j]；为此，我认为上三角cholesky是(有点)最快的，也就是说，我们寻找U’*U=M的上三角U。

对于固定的、小的维度，写一个不使用循环的版本绝对值得考虑。一个相对简单的方法是编写一个程序来完成这项工作。据我回忆，使用一个处理一般情况的例程作为模板，只花了一个上午就编写了一个程序，可以编写特定的固定维度版本。节省的费用可能相当可观。例如，我的通用版本需要0.47秒来进行一百万次9x9分解，而无环版本需要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( "inttut_chol_%.2d( double* P)n", dim);
printf( "{n");
printf( "doubletd;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 stdoutn");
}
return EXIT_SUCCESS;
}
/*  ----------------------------------------------------------------
*/

由于其易用性，您可以使用特征解算器进行比较。对于特定的用例，一个特定的解算器可能更快，尽管另一个应该更好。为此，您可以测量每个算法的运行时间，仅用于选择。之后，您可以实现所需的选项(或找到最适合您需求的现有选项)。