计算复杂 Cholesky 分解时的大错误

计算科学 矩阵分解
2021-12-26 11:12:18

我正在使用自己的例程来计算复杂的正定对称矩阵的 Cholesky 分解。

我的代码如下所示:

void CholeskyDecomposition(Complex **L, Complex **A, const int dim)
{
    Complex sum, t;
    Complex sum1;

    for (int k = 0; k < dim; k++)
    {
        // First, calculate the diagonal element of a column
        sum1 = A[k][k];
        for (int j = 0; j < k; j++)
        {
            Complex product;
            Cconjmul(product, L[j][k], L[j][k]);
            Csub(sum1, sum1, product);
        }
        Csqrt(&sum1, sum1);
        L[k][k] = sum1;

        // Now, calculate the elements right the diagonal element
        for (int i = k + 1; i < dim; i++)
        {
            sum = A[k][i];
            for (int j = 0; j < k; j++)
            {
                Cconjmul(t, L[j][k], L[j][i]);
                Csub(sum, sum, t);
            }
            Cdiv(L[k][i], sum, sum1);
        }
    }
}

我针对 MATLAB Cholesky-Factorization 测试了这段代码,结果是

A =

  13.6393 + 0.0000i   1.8844 + 3.4319i  -4.8788 - 4.0195i
   1.8844 - 3.4319i   5.7265 + 0.0000i   1.1568 - 0.4937i
  -4.8788 + 4.0195i   1.1568 + 0.4937i   6.6512 + 0.0000i

K>> L = chol(A)

L =

   3.6931 + 0.0000i   0.5102 + 0.9293i  -1.3210 - 1.0884i
   0.0000 + 0.0000i   2.1454 + 0.0000i   1.3248 - 0.5435i
   0.0000 + 0.0000i   0.0000 + 0.0000i   1.2927 + 0.0000i

我的 C-Implementation 的结果似乎并非完全错误。第一行与 MATLAB 解决方案相同。第二行和第三行在最后一列包含错误:

MATLAB         <-->       C
1.3248-0.5435i <--> 1.3248+0.08324i
1.2927+0.0000i <--> 1.3998+0.00000i

谁能帮我确定错误来源?

1个回答

抱歉打扰各位了。错误很简单:

// Now, calculate the elements right the diagonal element
for (int i = k + 1; i < dim; i++)
{
    sum = A[k][i];
    for (int j = 0; j < k; j++)
    {
        Cconjmul(t, L[j][k], L[j][i]);
        Csub(sum, sum, t);
    }
    Cdiv(L[k][i], sum, sum1);
}

只需交换复共轭乘法的 L。必须:

Cconjmul(t, L[j][i], L[j][k]);