当前位置: 动力学知识库 > 问答 > 编程问答 >

c - CHOLMOD sparse matrix cholesky decomposition: incorrect factor?

问题描述:

I have been using CHOLMOD to factorise the matrix A and solve the system Ax = b, for A being the Hessian matrix (printed below) and b = [1, 1, 1] created by the cholmod_ones function.

Unfortunately, the solution for x is incorrect (should be [1.5, 2.0, 1.5]) and to confirm I then multiplied A and x back together and don't get [1, 1, 1]. I don't quite understand what I am doing wrong.

Additionally, I've looked at the factor and the values of the matrix elements don't make sense either.

Output

Hessian:

2.000 -1.000 0.000

-1.000 2.000 -1.000

0.000 -1.000 2.000

Solution:

2.500 0.000 0.000

3.500 0.000 0.000

2.500 0.000 0.000

B vector:

1.500 0.000 0.000

2.000 0.000 0.000

1.500 0.000 0.000

Code

iterate_hessian() is an external function that returns doubles which are read into the CHOLMOD hessian matrix.

The entry point for the code is cholesky_determinant which is called with an argument which gives the dimension of the (square) matrix.

#include <cholmod.h>

#include <string.h>

// Function prototype that gives the next value of the Hessian

double iterate_hessian();

cholmod_sparse *cholmod_hessian(double *hessian, size_t dimension, cholmod_common *common) {

// This function assigns the Hessian matrix from OPTIM to a dense matrix for CHOLMOD to use.

// Allocate a dense cholmod matrix of appropriate size

cholmod_triplet *triplet_hessian;

triplet_hessian = cholmod_allocate_triplet(dimension, dimension, dimension*dimension, 0, CHOLMOD_REAL, common);

// Loop through values of hessian and assign their row/column index and values to triplet_hessian.

size_t loop;

for (loop = 0; loop < (dimension * dimension); loop++) {

if (hessian[loop] == 0) {

continue;

}

((int*)triplet_hessian->i)[triplet_hessian->nnz] = loop / dimension;

((int*)triplet_hessian->j)[triplet_hessian->nnz] = loop % dimension;

((double*)triplet_hessian->x)[triplet_hessian->nnz] = hessian[loop];

triplet_hessian->nnz++;

}

// Convert the triplet to a sparse matrix and return.

cholmod_sparse *sparse_hessian;

sparse_hessian = cholmod_triplet_to_sparse(triplet_hessian, (dimension * dimension), common);

return sparse_hessian;

}

void print_matrix(cholmod_dense *matrix, size_t dimension) {

// matrix->x is a void pointer, so first copy it to a double pointer

// of an appropriate size

double *y = malloc(sizeof(matrix->x));

y = matrix->x;

// Loop variables

size_t i, j;

// Row

for(i = 0; i < dimension; i++) {

// Column

for(j = 0; j < dimension; j++) {

printf("% 8.3f ", y[i + j * dimension]);

}

printf("\n");

}

}

cholmod_dense *factorized(cholmod_sparse *sparse_hessian, cholmod_common *common) {

cholmod_factor *factor;

factor = cholmod_analyze(sparse_hessian, common);

cholmod_factorize(sparse_hessian, factor, common);

cholmod_dense *b, *x;

b = cholmod_ones(sparse_hessian->nrow, 1, sparse_hessian->xtype, common);

x = cholmod_solve(CHOLMOD_LDLt, factor, b, common);

cholmod_free_factor(&factor, common);

// Return the solution, x

return x;

}

double cholesky_determinant(int *dimension) {

// Declare variables

double determinant;

cholmod_sparse *A;

cholmod_dense *B, *Y;

cholmod_common common;

// Start CHOLMOD

cholmod_start (&common);

// Allocate storage for the hessian (we want to copy it)

double *hessian = malloc(*dimension * *dimension * sizeof(hessian));

// Get the hessian from OPTIM

int i = 0;

for (i = 0; i < (*dimension * *dimension); i++) {

hessian[i] = iterate_hessian();

}

A = cholmod_hessian(hessian, *dimension, &common);

printf("Hessian:\n");

print_matrix(cholmod_sparse_to_dense(A, &common), *dimension);

B = factorized(A, &common);

printf("Solution:\n");

print_matrix(B, *dimension);

double alpha[] = {1, 0};

double beta[] = {0, 0};

Y = cholmod_allocate_dense(*dimension, 1, *dimension, CHOLMOD_REAL, &common);

cholmod_sdmult(A, 0, alpha, beta, B, Y, &common);

printf("B vector:\n");

print_matrix(Y, *dimension);

determinant = 0.0;

// Free up memory and finish CHOLMOD

cholmod_free_sparse (&A, &common);

cholmod_free_dense (&B, &common);

cholmod_finish (&common);

return determinant;

}

网友答案:

It turns out that I hadn't set the stype for my sparse matrix properly. The stype determines the symmetry (and thus the subsequent behaviour of calls to cholmod_factorize). It was in fact factorising and solving for AA'.

分享给朋友:
您可能感兴趣的文章:
随机阅读: