Eigen(C++)与Fortran 90极小值矩阵运算结果差异咨询

使用C++ Eigen库与Fortran PowerStation 4.0（调用MSIMSL模块）计算极小值矩阵的逆矩阵与原矩阵的乘积时，结果存在细微差异。该矩阵行列式极低（接近0），请问此差异是否与行列式接近零有关？期望两类计算结果相近，求相关技术建议。

C++测试代码 #

#include <Eigen/Dense>
#include <iostream>

int main(){

    MatrixXd test_matrix
    {
        {1.3e-17, -3.2e-17, -4.1e-17, 0.2e-17},
        {-2.3e-17, 0.2e-17, 1.1e-17, -2.2e-17},
        {3.3e-17, -3.2e-17, -3.1e-17, -7.2e-17},
        {-4.3e-17, 0.2e-17, 6.1e-17, 0.2e-17},
    };
    
    std::cout << "Test Matrix: " << std::endl;
    std::cout << test_matrix << std::endl;
    std::cout << " " << std::endl;
    
    std::cout << "Inverse of Test Matrix: " << std::endl;
    std::cout << test_matrix.inverse() << std::endl;
    std::cout << " " << std::endl;

    std::cout << "Inverse * Original: " << std::endl;
    std::cout << test_matrix * test_matrix.inverse() << std::endl;
    
    std::cout << "Determinant of Test Matrix: " << std::endl;
    std::cout << test_matrix.determinant() << std::endl;
    
    return 0;
}

Fortran测试代码 #

USE MSIMSL
INTEGER M
DOUBLE PRECISION A1(4,4),A2(4,4),B(4,4)
    OPEN(unit=1, file='output')

1 FORMAT(16(E20.5,1X))


M=4
A1(1,1)=1.3e-17
A1(1,2)=-3.2e-17
A1(1,3)=-4.1e-17
A1(1,4)=0.2e-17

A1(2,1)=-2.3e-17
A1(2,2)=0.2e-17
A1(2,3)=1.1e-17
A1(2,4)=-2.2e-17

A1(3,1)=3.3e-17
A1(3,2)=-3.2e-17
A1(3,3)=-3.1e-17
A1(3,4)=-7.2e-17

A1(4,1)=-4.3e-17
A1(4,2)=0.2e-17
A1(4,3)=6.1e-17
A1(4,4)=0.2e-17

C       returns inverse of matrix into A2
CALL DLINRG(M,A1,M,A2,M)
C       B = A1 * A2
        CALL DMRRRR(M,M,A1,M,M,M,A2,M,M,M,B,M)

        DO 777 I=1,4
        WRITE(1,1)(A1(I,K),K=1,4)
777     CONTINUE
        WRITE(1,*) '       '

        DO 888 I=1,4
        WRITE(1,1)(A2(I,K),K=1,4)
888     CONTINUE
        WRITE(1,*) '       '

        DO 999 I=1,4
        WRITE(1,1)(B(I,K),K=1,4)
999     CONTINUE
        WRITE(1,*) '       '

        END

C++输出结果 #

Test Matrix:
 1.3e-17 -3.2e-17 -4.1e-17    2e-18
-2.3e-17    2e-18  1.1e-17 -2.2e-17
 3.3e-17 -3.2e-17 -3.1e-17 -7.2e-17
-4.3e-17    2e-18  6.1e-17    2e-18

Inverse of Test Matrix:
-1.46582e+16  -4.1296e+16  1.23181e+16  3.85461e+15
-2.41195e+16  1.97719e+16 -7.36473e+15 -2.35196e+16
-9.81172e+15 -2.92629e+16  9.21779e+15  1.97601e+16
 8.22593e+15 -1.51155e+16 -8.93865e+15  3.71206e+15

Inverse * Original:
           1  2.42861e-17            0  1.14492e-16
 1.38778e-16            1 -2.77556e-17  2.77556e-17
           0  4.44089e-16            1            0
 1.38778e-17  2.42861e-17            0            1

Determinant of Test Matrix:
-3.3674e-66

Fortran输出结果 #

A1:
.13000E-16  -.32000E-16  -.41000E-16   .20000E-17
-.23000E-16   .20000E-17   .11000E-16  -.22000E-16
.33000E-16  -.32000E-16  -.31000E-16  -.72000E-16
-.43000E-16   .20000E-17   .61000E-16   .20000E-17
A2 = Inverse of A1:
-.14658E+17  -.41296E+17   .12318E+17   .38546E+16
-.24119E+17   .19772E+17  -.73647E+16  -.23520E+17
-.98117E+16  -.29263E+17   .92178E+16   .19760E+17
.82259E+16  -.15116E+17  -.89386E+16   .37121E+16
A2 * A1
.10000E+01   .21511E-15   .31225E-16  -.13878E-15
.55511E-16   .10000E+01  -.27756E-16   .13878E-16
.00000E+00  -.22204E-15   .10000E+01   .11102E-15
.33307E-15  -.69389E-17  -.19082E-15   .10000E+01

差异与行列式接近0的关联 #

是，二者高度相关。行列式接近0意味着矩阵是病态矩阵（条件数极大），这类矩阵的逆对输入的微小误差极其敏感：

1. 浮点数计算本身存在精度限制（双精度约15-17位有效数字），矩阵元素本身的微小舍入误差，在求逆过程中会被放大几个数量级；
2. Eigen和MSIMSL的求逆算法实现细节不同（比如LU分解的选主策略、数值稳定化技巧），误差放大的程度和方向也会有差异，最终导致逆矩阵与原矩阵的乘积结果出现细微偏差。

技术建议 #

1. 避免直接求逆 #

如果实际需求是求解线性方程组（而非单纯计算逆矩阵），直接使用A.x = b的求解接口，而非先求逆再乘：

• Eigen中用test_matrix.colPivHouseholderQr().solve(b)或test_matrix.fullPivLu().solve(b)，这类方法的数值稳定性远高于直接求逆；
• Fortran MSIMSL中选择专门的线性方程组求解子程序（如DLINSL），而非先调用DLINRG求逆再做矩阵乘法。

2. 矩阵缩放处理 #

将原矩阵的元素放大到合理数量级（比如乘以1e17，让元素值接近1），计算完成后再对结果做反向缩放：

• 缩放后矩阵的行列式会从~1e-66变为~1，条件数大幅降低，数值稳定性显著提升；
• 注意缩放因子的选择要确保元素不会超出浮点数的有效范围。

3. 使用更高精度的计算 #

如果必须保留原矩阵的数量级，可以尝试扩展精度计算：

• Eigen支持long double类型（部分编译器支持），或借助第三方高精度库；
• Fortran中使用QUAD PRECISION（若编译器支持），替换DOUBLE PRECISION。

4. 统一数值计算策略 #

• 强制Eigen使用与MSIMSL一致的分解算法：比如Eigen的FullPivLU分解（全选主元LU），通常是数值最稳定的通用求逆方法，可显式调用test_matrix.fullPivLu().inverse()；
• 检查MSIMSLDLINRG的文档，确认其使用的分解策略，尽量对齐两边的算法选择。

5. 评估矩阵的条件数 #

通过计算矩阵的条件数（如2-范数条件数），量化矩阵的病态程度：

• Eigen中用test_matrix.conditionNumber()；
• 若条件数远大于1e16（双精度浮点数的精度阈值），说明矩阵的逆计算本质上就不可靠，此时应重新审视问题本身是否需要求逆。

内容的提问来源于stack exchange，提问作者Yerlan Amir