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'.