7

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

  • 高斯消元
  • LU分解
  • 乔列斯基分解
  • ETC？

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

这样的系统将被求解数百万次，因此算法在维度（n = 9）方面应该相当快。

对所提出算法的鲁棒C++ 实现的PS 示例表示赞赏。

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

数百万个不同的矩阵。

2）正_semi_defined意味着矩阵可以是奇异的（机器精度）。你想如何处理这个案子？只是提出一个错误，或者尝试返回一些明智的答案？

引发错误是可以的。

4

6 回答 6

7

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

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

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

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

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

于 2012-11-14T07:14:06.533 回答
6

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

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

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

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

搜索周围，你会发现很多其他的。

于 2012-11-13T23:06:27.243 回答
3

我会推荐 LU 分解，特别是如果“解决数百万次”真的意味着“解决一次并应用于数百万个向量”。您将创建 LU 分解，保存它，并根据需要对尽可能多的 rhs 向量应用前向替换。

如果您使用旋转，它在四舍五入时会更稳定。

于 2012-11-13T23:55:08.030 回答
2

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

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

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

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

于 2012-11-14T08:28:32.890 回答
2

像上面的其他人一样，我推荐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;
}
/*  ----------------------------------------------------------------
*/
于 2012-11-14T12:17:36.193 回答
0

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

于 2014-09-09T07:00:46.443 回答