C++，Lapack Cholesky 分解实现不准确的结果

Question

我正在尝试用 C++ 实现 Cholesky 分解，它以前在 lapack dpotrf_ 中实现过。

Cholesky 分解：R' * R = A

代码：

#include <iostream>
#include <armadillo>




long my_chol(
    arma::mat &R,
    const arma::mat A,
    long lda
    )
{

    arma::arma_debug_check( (A.is_square() == false), "chol(): given matrix must be square sized" );
    arma::arma_debug_check( (lda<std::max(1l,(long)A.n_rows)), "chol(): LDA must be equal to or greater than max(1,N)" );

    double sum;
    long i, j, k;
    long n = A.n_rows;
    if(lda==(long)A.n_rows)
        R=A;
    else
    {
        R.zeros(lda,A.n_rows);
        R.submat(0,lda-1,0,A.n_rows-1)=A;
    }

    for(i=0; i<n; ++i)
        for(j=i+1; j<n; ++j)
            R(j,i)=0;


    for( i=0; i<n; ++i )
    {
        /* j == i */
        sum = R(i,i);

        for( k=(i-1); k>=0; --k )
            sum -= R(k,i)*R(k,i);

        if ( sum > 0.0 )
            R(i,i) = sqrt( sum );
        else
        {
            R(0) = sum; /* tunnel negative diagonal element to caller */
            return (long)i+1;
        }

        for( j=(i+1); j<n; ++j )
        {
            sum = R(i,j);

            for( k=(i-1); k>=0; --k )
                sum -= R(k,i) * R(k,i);

            R(i,j) = sum / R(i,i);
        }
    }
    return 0;
}

我使用以下代码测试此功能：

int main()
{
    arma::mat A={{10, 3, 5},{3, 60, 7},{5, 7, 9}};
    arma::mat B;
    my_chol(B,A,3);
    std::cout<<"---------------------------\n";
    std::cout<<"A:\n";
    A.print();
    std::cout<<"B:\n";
    B.print();
    std::cout<<"---------------------------\n";
    return 0;
}

这是结果：

---------------------------
A:
   10.0000    3.0000    5.0000
    3.0000   60.0000    7.0000
    5.0000    7.0000    9.0000
B:
   3.1623   0.9487   1.5811
        0   7.6877   0.7935
        0        0   2.4229
---------------------------

但是在八度音程中测试相同的矩阵给出了不同的结果：

A=[10,3,5;3,60,7;5,7,9];
chol(A)

   3.16228   0.94868   1.58114
   0.00000   7.68765   0.71543
   0.00000   0.00000   2.44707

尽管结果看起来很接近，但它们之间存在细微差别。

R23和R33略有不同。我检查了结果。 octave 的结果是正确的，而我的结果不是：

R=[3.1623,0.9487,1.5811;0,7.6877,0.7935;0,0,2.4229];
R'*R
ans =

   10.0001    3.0001    4.9999
    3.0001   60.0008    7.6002
    4.9999    7.6002    9.0000

为什么我的代码给出了错误的结果？

Answer 1

最里面的循环应该是

sum -= R(k,i) * R(k,j);

而不是

sum -= R(k,i) * R(k,i);

C++，Lapack Cholesky 分解实现不准确的结果

C++, Lapack Cholesky decomposition implementation inaccurate result

c++

octave

lapack

armadillo