我使用CHOLMOD在suiteparse中对一个相对稀疏的大带宽对角矩阵进行因子分解,即它只包含几个非零的对角线矩阵的稀疏性由协方差长度参数设置N
越大,非零的非对角元素的数量就越多。
当N
变大且许多元素为非零时,超级节点CHOLMOD分解突然开始失败,并出现错误消息“CHOLMOD warning:matrix not positive definity”。通过使用单独的Python实现检查数学,我知道矩阵应该是正定的而且,当我从
Common->.supernodal = CHOLMOD_SUPERNODAL;
到
Common->supernodal = CHOLMOD_SIMPLICIAL;
然后因子分解将成功对于下面的代码示例,只要
l
,超级节点分解就会成功如果我增加l
,就会得到矩阵不是正定的误差然而,如果我选择CHOMOD_SIMPLICIAL,则分解将成功有谁能帮我弄清楚为什么CHOLMOD_超级结在一定的稀疏度/密度上突然失效谢谢您!//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) ;
}
最佳答案
玩你的代码,看起来你正在创建相当病态的矩阵,并要求CHOLMOD处理它们当我告诉它对失败的情况进行简单的因子分解时,我得到了大约1e-7的倒数条件数(使用cholmod_rcond
计算)当l = 2.5
时,我得到一个关于2e-6
的倒数条件数注意afloat
的机器epsilon就在这附近。。。
如果我在您的float
声明中将所有的double
s替换为k_3_2
s,它将不再失败(我还没有看到矩阵里面有什么,所以除了说pow(a, 2.)
是一个很糟糕的平方方法之外,我不会做进一步的评论。)
我不知道你为什么会因为矩阵的条件不那么差而失败我不完全清楚CHOLMOD的超节分解是如何实现的,但我相信它需要BLAS的dpotrf
来做小的密集分解,而dsyrk
来处理块的其余部分很有可能一个快速的dpotrf
或dsyrk
会引起你的悲伤另外,CHOLMOD为矩阵的一些结构上为零的元素创建了非零,这样它就可以使用BLAS,这会让您有点不快。