1

私はCHOLMODを使用して、行列Aを因数分解し、システムAx = bを解きました。ここで、Aはヘッセ行列（下に印刷）であり、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ヘッセ行列に読み込まれるdoubleを返す外部関数です。

コードのエントリポイントは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;
}
4

1 に答える 1

0

スパース行列の stype を適切に設定していなかったことがわかりました。stype は対称性を決定します (したがって、 への呼び出しのその後の動作も​​決定されますcholmod_factorize)。実際には、AA' を因数分解して解くことでした。

于 2012-10-29T14:20:50.620 に答える