CHOLMOD supernodal factorization fails when sparse matrix becomes too dense

Question

I am using CHOLMOD in SuiteSparse to factor an N by N large band-diagonal matrix that is relatively sparse, ie, it contains only a few diagonals that are non-zero. The sparsity of the matrix is set by a covariance length parameter l. The larger l the greater number of off-diagonal elements that are non-zero.

When l becomes large and many elements are non-zero, the supernodal CHOLMOD factorization suddenly begins to fail with the error message that "CHOLMOD warning: matrix not positive definite." From checking the math with a separate Python implementation, I know that the matrix should be positive definite. Moreover, when I change from

Common->.supernodal = CHOLMOD_SUPERNODAL;

to

Common->supernodal = CHOLMOD_SIMPLICIAL;

the factorization will then succeed. For the following code example below, the supernodal factorization will succeed as long as l < 2.5. If I increase l >= 3.0, I get the error that the matrix is not positive definite. Yet, if I then choose CHOMOD_SIMPLICIAL, the factorization will succeed. Can anyone help me identify why CHOLMOD_SUPERNODAL suddenly fails above a certain sparsity/density? Thank you!

//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) ;
}

tmyklebu tmyklebu · Accepted Answer · 2014-04-01T01:21:41

Playing with your code, it looks like you're creating rather ill-conditioned matrices and asking CHOLMOD to deal with them. When I tell it to do a simplicial factorisation on the cases that fail, I get reciprocal condition numbers (computed using cholmod_rcond) around 1e-7. When l = 2.5, I get a reciprocal condition number around 2e-6. Note that machine epsilon for a float is somewhere around here...

If I replace all the floats with doubles in your declaration of k_3_2, it no longer appears to fail. (I haven't looked at what your matrix has inside of it, so I'm not going to comment further beyond saying that pow(a, 2.) is a poor way to square something.)

I don't know why you're getting failure for matrices that are apparently not so poorly-conditioned. I'm not completely sure of the details of how CHOLMOD's supernodal factorisation is implemented, but I believe it calls out to BLAS's dpotrf to do the small dense factorisations and dsyrk to deal with the rest of the block. It's possible that a fast dpotrf or dsyrk could cause you grief. Also, CHOLMOD creates nonzeros for some structurally zero elements of the matrix just so that it can use BLAS, and this can throw you off a little bit.

CHOLMOD supernodal factorization fails when sparse matrix becomes too dense

1 Answers