# 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.0000.000 -1.000 2.000Solution:2.500 0.000 0.0003.500 0.000 0.0002.500 0.000 0.000B vector:1.500 0.000 0.0002.000 0.000 0.0001.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 Hessiandouble 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 sizecholmod_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 sizedouble *y = malloc(sizeof(matrix->x));y = matrix->x;// Loop variablessize_t i, j;// Rowfor(i = 0; i < dimension; i++) {// Columnfor(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, xreturn x;}double cholesky_determinant(int *dimension) {// Declare variablesdouble determinant;cholmod_sparse *A;cholmod_dense *B, *Y;cholmod_common common;// Start CHOLMODcholmod_start (&common);// Allocate storage for the hessian (we want to copy it)double *hessian = malloc(*dimension * *dimension * sizeof(hessian));// Get the hessian from OPTIMint 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 CHOLMODcholmod_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'.