上三角矩阵递归求逆：128维矩阵验证失败问题排查

问题背景 #

我实现了一个递归函数来计算上三角矩阵的逆，采用分治策略将矩阵划分为Top、Bottom和Corner三个部分，遵循分块矩阵求逆的数学方法。仅支持维度为2的幂的矩阵，测试2、4、8、16、32、64维矩阵时都能通过断言（A*A_inv接近单位矩阵），但到128维时断言失败，test_matrix的部分非对角元不为0，明显超出了可接受的误差范围。

核心实现代码 #

用C++的std::vector<std::vector<double>>存储矩阵，核心递归逻辑如下：

Matrix inverse_matrix(Matrix A){
    Matrix A_inv;
    if(A.dimension == 2){
        A_inv = simple_inverse(A);
    } else{
        // 拆分矩阵为Top、Corner、Bottom（省略具体拆分逻辑）
        Matrix Top = get_top_block(A);
        Matrix Corner = get_corner_block(A);
        Matrix Bottom = get_bottom_block(A);

        Matrix Top_inv = inverse_matrix(Top);
        // 验证Top*Top_inv为单位矩阵
        Matrix Bottom_inv = inverse_matrix(Bottom);
        // 验证Bottom*Bottom_inv为单位矩阵

        Matrix Corner_inv = multiply_matrix(Top_inv, Corner);
        Corner_inv = multiply_matrix(Corner_inv, Bottom_inv);
        Corner_inv = negate(Corner_inv); // 矩阵元素取反

        // 将各块组装到A_inv中（省略组装逻辑）
        assemble_inverse(A_inv, Top_inv, Corner_inv, Bottom_inv);
    }
    return A_inv;
}

小维度时，test_matrix的非对角元是极小值（如-9.58122e-14），递归过程中所有子矩阵的逆验证都通过，推测数学逻辑正确，问题出在double类型的精度特性上。

问题根源分析 #

你遇到的是典型的浮点数累积误差问题：

• double类型的有效精度约为15-17位十进制数，每次矩阵乘法和递归求逆都会引入微小的舍入误差。
• 128维矩阵需要7层递归（128=2^7），每一层的误差会不断叠加，尤其是Corner_inv的计算涉及三次矩阵乘法（Top_inv * Corner * Bottom_inv），每一次乘法都会放大误差，最终在顶层矩阵的非对角元上表现为超出断言阈值的非零值。
• 子矩阵的逆验证看似没问题，是因为子矩阵维度小，累积误差还没达到可观测的程度，但顶层是多层误差的叠加结果。

解决方案 #

1. 调整断言的误差阈值 #

不要严格判断矩阵是否等于单位矩阵，而是判断每个元素与单位矩阵对应元素的差值是否在可接受范围内：

bool is_identity(const Matrix& mat, double eps = 1e-8) {
    int n = mat.size();
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            double expected = (i == j) ? 1.0 : 0.0;
            if (std::abs(mat[i][j] - expected) > eps) {
                return false;
            }
        }
    }
    return true;
}

可以根据维度动态调整eps：比如128维时放宽到1e-6或1e-7，维度越大，累积误差的容忍度可以适当提高。

2. 优化矩阵乘法的精度 #

• 针对上三角矩阵的特性优化乘法：只计算上三角部分的元素，避免对零元素的无效操作，减少舍入次数。
• 使用Kahan求和算法降低点积计算的累加误差：

double kahan_dot(const std::vector<double>& a, const std::vector<double>& b) {
    double sum = 0.0;
    double error = 0.0;
    for (size_t i = 0; i < a.size(); ++i) {
        double y = a[i] * b[i] - error;
        double t = sum + y;
        error = (t - sum) - y;
        sum = t;
    }
    return sum;
}

用这个函数代替普通的点积计算，能有效减少累加过程中的舍入误差。

3. 改用更高精度的浮点数类型 #

如果精度要求严格，可以替换double为long double（多数平台为80位精度，有效位数约18-19位），将矩阵存储改为std::vector<std::vector<long double>>，能显著降低累积误差，但会带来一定的性能开销。

4. 再次验证分块逻辑 #

虽然子矩阵的逆验证通过，但仍需检查矩阵拆分和组装的代码，确保Top、Corner、Bottom的划分完全正确（比如128维矩阵拆分后，Top和Bottom都是64x64，Corner是64x64），避免索引越界或元素复制错误——这类问题可能在大维度时才暴露，容易和精度误差混淆。

总结 #

大维度下的断言失败本质是浮点数累积误差的正常表现，并非数学逻辑错误。通过调整误差阈值、优化乘法精度或改用更高精度类型，就能解决这个问题。

内容的提问来源于stack exchange，提问作者Tapojyoti Mandal