从 Matlab 到 C 的 Cholesky 分解函数

Question

我尝试将执行乔列斯基分解的 Matlab 循环转换为 C。仅第一列就得到了正确的结果。

Matlab代码:

M = [5 1.2 0.3 -0.6; 1.2 6 -0.4 0.9; 0.3 -0.4 8 1.7; -0.6 0.9 1.7 10];
n = length( M );
L = zeros( n, n );
for i=1:n
    L(i, i) = sqrt( M(i, i) - L(i, :)*L(i, :)' );

    for j=(i + 1):n
        L(j, i) = ( M(j, i) - L(i, :)*L(j, :)' )/L(i, i);
    end
end

我的C代码:

int  main()
{

    double M[4*4]={5, 1.2, 0.3, -0.6, 1.2, 6, -0.4, 0.9, 0.3, -0.4, 8, 1.7, -0.6, 0.9, 1.7, 10};
    int n=4;
    double L[4*4]={0};
    double sum;
    for (int i=0;i<n;i++){
        sum=0;
        for (int j=0;j<n;j++){
            sum+=L[i*n+j]*L[j*n+i];
            }
    L[i*n+i]=sqrt(M[i*n+i]-sum);

    for (int j=i+1;j<n;j++){
        sum=0;
        for(int k=0;k<n;k++){
            sum+=L[i*n+k]*L[k*n+j];
           }
        L[j*n+i]=(M[j*n+i]-sum)/L[i*n+i];

      }

}
 return 0;

}

Matlab 的结果是:

L =

    2.2361         0         0         0
    0.5367    2.3900         0         0
    0.1342   -0.1975    2.8183         0
   -0.2683    0.4368    0.6466    3.0527

但在 C 语言中结果是:

L_c =

    2.2361         0         0         0
    0.5367    2.4495         0         0
    0.1342   -0.1633    2.8284         0
   -0.2683    0.3674    0.6010    3.1623

哪里出错了？

Answer 1

检查线路

sum+=L[i*n+j]*L[j*n+i];

和

sum+=L[i*n+k]*L[k*n+j];

我认为应该是这样

sum+=L[i*n+j]*L[i*n+j];

和

sum+=L[i*n+k]*L[j*n+k];

分别。 Matlab代码中的转置只是改变 vector 的方向，但不会改变参数或其他东西。

我建议您创建一个计算点积的函数，以便您可以重用代码。