CHOLMOD 稀疏矩阵 Cholesky 分解：不正确的因子？

Question

我一直在使用 CHOLMOD 分解矩阵 A 并求解系统 Ax = b，因为 A 是 Hessian 矩阵（打印在下面），b = [1, 1, 1] 由 cholmod_ones 函数创建。

不幸的是，x 的解是不正确的（应该是 [1.5, 2.0, 1.5]），为了确认我然后将 A 和 x 相乘并没有得到 [1, 1, 1]。我不太明白我做错了什么。

此外，我查看了因子，矩阵元素的值也没有意义。

输出

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 

代码

iterate_hessian() 是一个外部函数，它返回读入 CHOLMOD hessian 矩阵的双精度数。

代码的入口点是cholesky_determinant，它使用一个参数调用，该参数给出（方形）矩阵的维数。

#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;
}

Answer 1

事实证明，我没有正确设置稀疏矩阵的 stype。 stype 确定对称性（以及随后调用cholmod_factorize 的行为）。它实际上是对 AA' 进行因式分解和求解。