固定维数（N=9）、对称、半正定的密集线性系统的快速求解

Posted 2023-02-22

【中文标题】固定维数（N=9）、对称、半正定的密集线性系统的快速求解

【英文标题】：Fast solution of dense linear system of fixed dimension (N=9), symmetric, positive-semidefinite

【发布时间】：2012-11-13 21:53:03

【问题描述】：

对于固定维度 (N=9) 的密集线性系统（矩阵是对称的，半正定的）的快速求解，您会推荐哪种算法？

高斯消元法 LU 分解 Cholesky 分解 等等？

类型为 32 位和 64 位浮点数。

此类系统将被求解数百万次，因此算法在维度 (n=9) 方面应该相当快。

附：推荐算法的健壮 C++ 实现示例。

1) “解决百万次”是什么意思？相同的系数矩阵具有一百万个不同的右手项，还是一百万个不同的矩阵？

数百万个不同的矩阵。

2) 正 _semi_defined 表示矩阵可以是奇异的（机器精度）。你想如何处理这个案子？只是提出一个错误，还是尝试返回一些合理的答案？

引发错误是可以的。

【问题讨论】：

您的问题可能更适合scicomp.stackexchange.com

我肯定会建议在 scicomp 上重新发布问题：但是，您应该指定一些附加信息。 1）“解决百万次”是什么意思？相同的系数矩阵具有一百万个不同的右手项，还是一百万个不同的矩阵？ 2）正_semi_defined意味着矩阵可以是奇异的（机器精度）。你想如何处理这个案子？只是提出一个错误，或者尝试返回一些合理的答案？

【参考方案1】：

矩阵是对称的，半正定的，Cholesky 分解 严格优于 LU 分解。 （无论矩阵大小如何，大约比 LU 快两倍。来源：Trefethen 和 Bau 的“数值线性代数”）

它实际上也是小型密集矩阵的标准（来源：我攻读计算数学博士学位）除非系统变得足够大，否则迭代方法的效率低于直接方法（快速的经验法则没有任何意义，但是总是很高兴：在任何现代计算机上，任何小于 100*100 的矩阵绝对是一个需要直接方法而不是迭代方法的小矩阵）

现在，我不建议自己做。那里有大量经过彻底测试的优秀库。但如果非要我给你推荐一个，那就是Eigen：

无需安装（只有标头库，所以没有要链接的库，只有#include） 强大而高效（他们在主页上有很多基准测试，结果很好） 易于使用且有据可查

顺便说一句，here in the documentation，您可以在一个简洁的表格中了解他们的 7 个直接线性求解器的各种优缺点。似乎在你的情况下，LDLT（Cholesky 的变体）获胜

【讨论】：

“但如果我必须向您推荐一个，那就是 Eigen” - 是的，我考虑过 - 它非常酷，并且具有内置的固定大小的向量和矩阵。但不幸的是，文档说它不适用于较旧的编译器（如 GCC 3.x、VS2005 等）。虽然这样的限制对于非常高性能的库来说是相当合理的，但我从事的项目应该支持这样的旧编译器。 “LDLT（Cholesky 的一种变体）获胜” - 是的，我检查了它，与 LLT 相比，它不需要 sqrt。我主要担心的是——也许 LU 更适合像 9x9 这样的小矩阵？

哦，该死的，我从来没有考虑过支持旧库...嗯，我不知道有什么库可以满足您的需求，但我知道当矩阵是 s.p.d 时，Cholesky 是大约比 LU 快两倍（来源：Trefethen 和 Bau 的“数值线性代数”）它与大小无关。也许你必须自己实现它......

“也许你必须自己实现它......” - 我期待 :) 这就是为什么我的主要问题是关于算法，而不是关于库。

@Fezvez - 以上对于更大和稀疏的矩阵是否也正确？ （更大的意思是 n:~700-~300 有 ~1/16 的矩阵有非零值）？谢谢。

对不起，我在线性代数这方面有点生疏了。我希望你能找到答案（不过，最好的方法是自己进行基准测试）

【参考方案2】：

一般来说，最好使用现有的库，而不是自己动手的方法，因为要实现快速、稳定的数值实现，需要处理许多繁琐的细节。

这里有一些可以帮助您入门：

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

犰狳： http://arma.sourceforge.net/

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

【讨论】：

虽然我同意库在大多数情况下是最佳选择，但对于如此小的矩阵，与成熟库相关的开销可能会过多。

我已经考虑过 Eigen。感谢犰狳点，我会检查它。

Armadillo 是一个不错的选择，因为它使用经过充分测试的 LAPACK 进行 LU 分解、Cholesky 分解等。您还可以使用高度优化的插入式替换轻松更改底层 LAPACK 库，例如作为英特尔的MKL libraries。

【参考方案3】：

我会推荐 LU 分解，特别是如果“解决数百万次”真的意味着“解决一次并应用于数百万个向量”。您将创建 LU 分解，保存它，然后对尽可能多的 r.h.s 应用前后替换。随心所欲的矢量。

如果你使用旋转，它在面对舍入时更稳定。

【讨论】：

不，抱歉，我的意思是数以百万计的不同 Ax=b 系统（其中 A 也不同）。无论如何+1

【参考方案4】：

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

在 LLT 或 LDLT 之间进行选择实际上取决于矩阵的条件数，以及您打算如何处理边缘情况。仅当您可以证明准确性在统计上显着提高，或者稳健性对您的应用至关重要时，才应使用 LDLT。

（如果没有您的矩阵样本，很难给出合理的建议，但我怀疑对于如此小的阶数 N=9，确实没有必要将小对角线项转向 D 的底部。所以我会从经典的 Cholesky 开始，如果 diag 项相对于某些明智选择的容差变得很小，则简单地中止分解。）

Cholesky 的代码非常简单，如果您力求编写真正快速的代码，最好自己实现。

【讨论】：

我担心“经典 Cholesky”LLT 涉及 N 个平方根。也许对大 N 来说不是问题，但对 9 来说 - 这是非常可疑的。

您可以在不进行旋转的情况下使用 LDLT：在因式分解过程中，您以 (N-1)(N-2)/2 乘法为代价节省了 N 个平方根：因此它大约是 9 个平方根对 28 个乘法.这应该与主要成本相比，大约是 (N^3-N)/6 madds（madd=multiply and add），所以 720 madds。对反向解的影响大约是 N 次乘法。如果您担心这些数字，则别无选择实施您自己的非常激进的子例程。对于此类问题，您必须在非常低的级别进行优化，同时考虑不同缓存级别的内存带宽和延迟、CPU 亲和性……

【参考方案5】：

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

实际上，cholesky 有许多变体，哪一种最快取决于您为矩阵选择的表示形式。我通常使用 fortran 样式表示，即矩阵 M 是 double* M ， M(i,j) 是 m[i+dim*j];为此，我认为上三角 cholesky 是（有点）最快的，即寻找 U'*U = M 的上三角 U。

对于固定的、小尺寸的，绝对值得考虑编写一个不使用循环的版本。一个相对简单的方法是编写一个程序来完成它。我记得，使用处理一般情况的例程作为模板，只用了一个上午就编写了一个可以编写特定固定尺寸版本的程序。节省的费用可能相当可观。例如，我的通用版本需要 0.47 秒来完成一百万个 9x9 因式分解，而无环版本需要 0.17 秒——这些时间在 2.6GHz 电脑上运行单线程。

为了表明这不是一项主要任务，我在下面提供了此类程序的源代码。它包括作为注释的分解的一般版本。我在矩阵不接近单数的情况下使用了此代码，并且我认为它在那里可以正常工作；然而，对于更精细的工作来说，它可能太粗糙了。

/*  ----------------------------------------------------------------
**  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;

/*  ----------------------------------------------------------------
*/

【讨论】：

+1 表示贡献的代码。编译器可以完成某种程度的循环展开，但代价是有时寄存器分配不佳。内部 Chol 循环包含 720 乘法和加法：所以完全展开是可行的，正如您的代码所示，但我不确定它是否会比仔细编码的循环版本更好......

gcc 的循环展开并没有给我留下深刻的印象。如果我采用可变维度代码，删除 dim 参数并将其替换为 9，生成的函数在大约 0.43 微秒内执行，而可变维度一的执行时间为 0.47，而“手动编码”版本为 0.17

也许您可以提供一个示例线性系统来求解以及如何将其作为矩阵输入到此代码中。

【参考方案6】：

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