当稀疏矩阵变得太密集时，CHOLMOD 超节点分解失败

Question

我在 SuiteSparse 中使用 CHOLMOD 来分解 N 与 N 相对稀疏的大带对角矩阵，即它只包含几个非零的对角线。矩阵的稀疏度由协方差长度参数l 设置。 l 越大，非零非对角线元素的数量就越多。

当l 变大并且许多元素不为零时，超节点 CHOLMOD 分解突然开始失败，并显示错误消息“CHOLMOD 警告：矩阵不是正定的”。通过使用单独的 Python 实现检查数学，我知道矩阵应该是正定的。此外，当我从

Common->.supernodal = CHOLMOD_SUPERNODAL; 

到

Common->supernodal = CHOLMOD_SIMPLICIAL;

然后分解将成功。对于下面的以下代码示例，只要l < 2.5，超节点分解就会成功。如果我增加l >= 3.0，我会得到矩阵不是正定的错误。然而，如果我随后选择 CHOMOD_SIMPLICIAL，分解将成功。谁能帮我确定为什么 CHOLMOD_SUPERNODAL 在一定的稀疏度/密度之上突然失败？谢谢！

//file is cov.c
//Compile with $ gcc -Wall -o cov -Icholmod cov.c -lm -lcholmod -lamd -lcolamd -lblas -llapack -lsuitesparseconfig
#include <stdio.h>
#include <math.h>
#include "cholmod.h"

#define PI 3.14159265

float k_3_2 (float r, float a, float l, float r0)
{
    return (0.5 + 0.5 * cos(PI * r/r0)) * pow(a, 2.) * (1 + sqrt(3) * r/l) * exp(-sqrt(3) * r/l) ;
}

// function to initialize an array of increasing wavelengths
void linspace (double *wl, int N, double start, double end)
{
    double j; //double index
    double Ndist = (double) N;
    double increment = (end - start)/(Ndist -1.);

    int i;
    for (i = 0; i < N; i++) {
        j = (double) i;
    wl[i] = start + j * increment;
    }
}

// create and return a sparse matrix using a wavelength array and parameters
// for a covariance kernel.
cholmod_sparse *create_sparse(double *wl, int N, double a, double l, cholmod_common *c)
{
    double r0 = 6.0 * l; //Beyond r0, all entries will be 0
    //Pairwise calculate all of the r distances
    int i = 0, j = 0;
    double r;

    //First loop to determine the number non-zero elements
    int M = 0; //number of non-zero elements
    for (i = 0; i < N; i++)
    {
        for (j = 0; j < N; j++)
        {
        r = fabs(wl[i] - wl[j]);
        if (r < r0) //Is the separation below our cutoff?
            M++;
        }
    }

    /* Initialize a cholmod_triplet matrix, which we will subsequently fill with
     * values. This matrix is NxN sparse with M total non-zero elements. 1 means we
     * want a square and symmetric matrix. */
    cholmod_triplet *T = cholmod_allocate_triplet(N, N, M, 1, CHOLMOD_REAL, c);

    if (T == NULL || T->stype == 0)         /* T must be symmetric */
    {
      cholmod_free_triplet (&T, c) ;
      cholmod_finish (c) ;
      return (0) ;
    }

    //Do the loop again, this time to fill in the matrix
    int * Ti = T->i;
    int * Tj = T->j;
    double * Tx = T->x;
    int k = 0;

    //This time, only fill in the lower entries (and diagonal). 
    for (i = 0; i < N; i++)
    {
        for (j = 0; j <= i; j++)
        {
        r = fabs(wl[i] - wl[j]);

            if (r < r0) //If the distance is below our cutoff, initialize
            {
                Ti[k] = i;
                Tj[k] = j;
                Tx[k] = k_3_2(r, a, l, r0);
                k++;
            }
        }
    }
    T->nnz = k;

    //The conversion will transpose the entries and add to the upper half.
    cholmod_sparse *A = cholmod_triplet_to_sparse(T, k, c);
    cholmod_free_triplet(&T, c);
    return A;

}

int main(void)
{
    //Create a sample array of wavelengths for testing purposes. 
    int N = 3000; 
    double wl[N];
    linspace(wl, N, 5100., 5200.); //initialize with wavelength values

    cholmod_common c ; //c is actually a struct, not a pointer to it.
    cholmod_start (&c) ; // start CHOLMOD
    c.print = 5;

    //c.supernodal = CHOLMOD_SIMPLICIAL;
    c.supernodal = CHOLMOD_SUPERNODAL;

    float l = 2.5;
    cholmod_sparse *A = create_sparse(wl, N, 1.0, l, &c);

    cholmod_factor *L ;
    L = cholmod_analyze (A, &c) ;           

    cholmod_factorize (A, L, &c) ;          
    printf("L->minor = %d\n", (int) L->minor);

    cholmod_dense *b, *x, *r;

    //Create a vector with the same number of rows as A
    b = cholmod_ones (A->nrow, 1, A->xtype, &c) ;   // b = ones(n,1) 
    x = cholmod_solve (CHOLMOD_A, L, b, &c) ;       // solve Ax=b    
    //the reason these are length two is because they can be complex
    double alpha [2] = {1,0}, beta [2] = {0,0} ;        // basic scalars 

    r = cholmod_copy_dense (b, &c) ;            // r = b 
    cholmod_sdmult (A, 0, alpha, beta, x, r, &c) ;      // r = Ax 

    cholmod_print_dense(r, "r", &c); //This should be equal to b

    cholmod_free_sparse(&A, &c); // free all of the variables
    cholmod_free_factor(&L, &c);
    cholmod_free_dense(&b, &c);
    cholmod_free_dense(&x, &c);
    cholmod_free_dense(&r, &c);
    cholmod_finish (&c) ;  // finish CHOLMOD 

    return (0) ;
}

Answer 1

在处理您的代码时，您似乎正在创建相当病态的矩阵并要求 CHOLMOD 处理它们。当我告诉它对失败的情况进行简单分解时，我得到大约 1e-7 的倒数条件数（使用cholmod_rcond 计算）。当l = 2.5 时，我在2e-6 附近得到一个倒数条件数。请注意，float 的机器 epsilon 就在附近...

如果我在您的k_3_2 声明中将所有floats 替换为doubles，它似乎不再失败。 （我没有看过你的矩阵里面有什么，所以除了说pow(a, 2.) 是一种糟糕的平方方法之外，我不会再发表评论。）

我不知道为什么你会因为显然不是那么糟糕的矩阵而失败。我不完全确定 CHOLMOD 的超节点因式分解是如何实现的细节，但我相信它需要 BLAS 的dpotrf 进行小密集分解和dsyrk 来处理块的其余部分。快速的dpotrf 或dsyrk 可能会让您感到悲伤。此外，CHOLMOD 会为矩阵的一些结构上为零的元素创建非零值，以便它可以使用 BLAS，这可能会让您有点失望。