CHOLMOD 희소 매트릭스 콜레 스키 분해 : 잘못된 요소?

Question

나는 행렬 A를 인수 분해하고 시스템 A = b를 풀기 위해 CHOLMOD를 사용했다. A는 cholmod_ones 함수로 생성 된 헤 시안 행렬 (아래에 인쇄 됨)과 b = [1, 1, 1]이다.

불행히도 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 독일인 행렬로 판독 복식 반환 외부 함수이다.

코드의 진입 점은 (제곱) 행렬의 차원을 제공하는 인수와 함께 호출되는 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 '에 대한 인수 분해 및 해결이었습니다.