矩阵形成中的奇怪行为（C++，犰狳）

Posted 2023-02-16

【中文标题】矩阵形成中的奇怪行为（C++，犰狳）

【英文标题】：Strange behavior in matrix formation (C++, Armadillo)

【发布时间】：2021-04-18 20:38:41

【问题描述】：

只要能量变量（双精度类型）没有收敛到某个阈值以下，我就有一个 while 循环。计算此能量所需的变量之一是双精度的犰狳矩阵，名为 f_mo。在while循环中，这个f_mo迭代更新，所以我在每个循环开始时计算f_mo为：

    arma::mat f_mo = h_core_mo;//h_core_mo is an Armadillo matrix of doubles
    for(size_t p = 0; p < n_mo; p++) //n_mo is of type size_t
    for(size_t q = 0; q < n_mo; q++) 
        double sum = 0.0;
            for(size_t i = 0; i < n_occ; i++) //n_occ is of type size_t
                //f_mo(p,q) += 2.0*g_mat_full_qp1_qp1_mo(p*n_mo + q, i*n_mo + i)-g_mat_full_qp1_qp1_mo(p*n_mo+i,i*n_mo+q); //all g_mat_ are Armadillo matrices of doubles
                sum += 2.0*g_mat_full_qp1_qp1_mo(p*n_mo + q, i*n_mo + i)-g_mat_full_qp1_qp1_mo(p*n_mo+i,i*n_mo+q);
            
            for(size_t i2 = 0; i2 < n_occ2; i2++) //n_occ2 is of type size_t
                //f_mo(p,q) -= 1.0*g_mat_full_qp1_qp2_mo(p*n_mo + q, i2*n_mo2 + i2);
                sum -= 1.0*g_mat_full_qp1_qp2_mo(p*n_mo + q, i2*n_mo2 + i2);
            
        
        f_mo(p,q) +=sum;
    

但是假设我直接将总和（我在 f_mo(p,q) 的末尾添加）替换为 f_mo(p,q) 的添加（注释掉的代码）。输出 f_mo 矩阵与机器精度相同。有关代码的任何内容都不应更改。循环中唯一受影响的变量是 sum 和 f_mo。然而，代码会收敛到不同的能量，并且会在非常不同的 while 循环迭代次数中收敛。我不知道差异的原因。当我运行相同的代码 2、3、4、5 次时，每次都会得到相同的结果。当我在没有优化的情况下重新编译时，我遇到了同样的问题。当我在另一台计算机上运行（控制环境）时，尽管 f_mo 相同，但我再次在 # of while 循环中发现差异，但每种方法的迭代总数（sum += 和 f_mo(p,q) += ) 不同。

值得注意的是，代码输出不同的点始终是g_mat_full_qp1_qp2_mo，稍后在while循环中重新计算。但是，进入 g_mat_full_qp1_qp2_mo 计算的每个变量在两个代码之间都是相同的。这使我认为 C++ 有一些我不理解的更深刻的东西。我欢迎任何关于如何继续调试此行为的想法（我几乎可以肯定这不是一个典型的错误，并且我已经控制了环境和优化）

【问题讨论】：

这可能是一个精度错误。 sum 从零开始，因此将小幅度数字添加到 0 将保持完整精度。因为f_mo(p,q) 中有一个值，所以添加小幅度数字可能会导致精度损失，因为较低精度的位会丢失。见this question。

@1201ProgramAlarm 谢谢你的提示。我应该补充一点，能量值不同的第一次迭代，它们在 10^-10 的数量级上不同（并且在上一次迭代中它们匹配到 10^-50）。如果我的解释有误，请纠正我，但因此在我看来这不太可能是由精度错误引起的。

没有您的代码和minimal reproducible example，只能推测原因。

【参考方案1】：

我将假设这是一个 Hartree-Fock 或某种其他类型的电子结构计算，您将两电子积分添加到核心哈密顿量中，并应用一些领域知识。

部分假设是两电子积分的单个元素非常小，特别是与核心哈密顿量相比。因此，正如 1201ProgramAlarm 在他们的评论中提到的那样，添加的顺序很重要。 You will get a more accurate result if you add smaller numbers together first to avoid loosing precision when adding two numbers many orders of magintude apart.。因为您会迭代此过程，直到 Fock 矩阵 f_mo 紧密收敛，您最终会收敛到相同的值。

为了以更准确的顺序将数字相加，并希望更快地收敛，大多数电子结构程序都有一个单独的例程来计算二电子积分，然后将它们添加到核心哈密顿量，这就是您在示例代码中逐个元素地做的事情。

Presentation on numerical computing.