上三角矩阵的矩阵逆计算给出了大矩阵维度的误差

Posted 2023-02-22

【中文标题】上三角矩阵的矩阵逆计算给出了大矩阵维度的误差

【英文标题】：Matrix inverse calculation of upper triangular matrix gives error for large matrix dimensions

【发布时间】：2020-02-08 22:50:46

【问题描述】：

我有一个递归函数来计算上三角矩阵的逆矩阵。我将矩阵分为顶部、底部和角落部分，然后遵循https://math.stackexchange.com/a/2333418 中规定的方法。这是一个伪代码形式：

//A diagram structure of the Matrix
Matrix = [Top Corner]
         [0   Bottom]

Matrix multiply_matrix(Matrix A, Matrix B)
    Simple Code to multiply two matrices and return a Matrix


Matrix simple_inverse(Matrix A)
    Simple Code to get inverse of a 2x2 Matrix


Matrix inverse_matrix(Matrix A)
   //Creating an empty A_inv matrix of dimension equal to A
   Matrix A_inv;
   if(A.dimension == 2)
       A_inv = simple_inverse(A) 
    
   else
       Top_inv = inverse_matrix(Top);
       (Code to check Top*Top_inv == Identity Matrix)

       Bottom_inv = inverse_matrix(Bottom);
       (Code to check Bottom*Bottom_inv == Identity Matrix)

       Corner_inv = multiply_matrix(Top_inv, Corner);
       Corner_inv = multiply_matrix(Corner_inv, Bottom_inv);
       Corner_inv = negate(Corner_inv);    //Just a function for negation of the matrix elements

       //Code to copy Top_inv, Bottom_inv and Corner_inv to A_inv
       ...
       ...
   
   return A_inv;


int main()
    matrix A = An upper triangular matrix with random integers between 1 and 9;
    A_inv = inverse_matrix(A);
    test_matrix = multiply_matrix(A, A_inv);
    (Code to Raise error if test_matrix != Identity matrix)

为简单起见，我已经实现了只支持二维矩阵的幂的代码。

我的问题是我已经针对 2、4、8、16、32 和 64 的矩阵维度测试了这段代码。所有这些都通过了代码中所示的所有断言检查。 但是对于 128 的矩阵维度，我得到的失败是 main() 中的断言。当我检查时，我观察到 test_matrix 不是身份矩阵。一些非对角元素不等于 0。 我想知道这可能是什么原因：-

我正在使用 C++ std::vector<std::vector<double>> 进行矩阵表示。 由于数据类型是double，对于情况 2、4、8、...、64，test_matrix 的非对角元素确实有一些价值，但非常小。例如，-9.58122e-14 我在任何递归阶段的所有矩阵都是方阵 我正在检查Top*Top_inv = Identity 和Bottom*Bottom_inv = Identity。 最后，对于维度 2、4、...、64，我生成了随机数（b/w 1 和 10）来创建我的上三角矩阵。既然这些案例都通过了，我想我的数学实现是正确的。

我觉得 C++ 数据类型的某些方面是关于 double 的，我不知道这可能会导致错误。否则 64->128 的突然错误没有意义。

【问题讨论】：

当增加大小时，随机矩阵越来越病态，这增加了最终的路由错误，即使是double。在这里，您可以通过计算最大对角元素与最小对角元素之间的比率（绝对值）来评估它。

【参考方案1】：

能否详细说明matrix == identity操作是如何实现的？ 我的猜测是，问题可能会恢复到浮点比较。

在最坏的情况下，矩阵求逆可能是 O(n^3)。这意味着，随着矩阵大小的增加，所涉及的计算量也会增加。即使使用 64 位浮点数也无法完美表示实数，它们始终是近似值。

对于诸如矩阵求逆之类的运算，这可能会导致数值误差传播的问题，因为累积的乘加运算会损失精度。

对此，*** 中已经有讨论：How should I do floating point comparison?

编辑：如果整个矩阵实际上是可逆的，则需要考虑其他事情。 也许顶部和/或底部矩阵是可逆的，但完整的矩阵（与 Corner 矩阵组合时）不是。

【讨论】：

我认为 matrix==identity 计算没有错误。因为我有代码来检查“test_matrix”的元素。 test_matrix 应该几乎等同于身份矩阵，但这在我的情况下是失败的。有趣的是，只有矩阵“角”区域中的元素具有非零值。因此，“数值误差传播”可能是问题所在，但我的“顶部”和“底部”矩阵是正确的，所以不确定为什么“角”矩阵计算会出错。

那我可能离这里很远，但是，如果顶部和底部矩阵不是奇异的，但整个矩阵是奇异的或条件不好怎么办？

@Caceres 我认为您关于“数字错误传播”的想法是正确的。为了进行实验，我将数据类型从“双”更改为“浮点”。现在我也遇到了 32 和 64 矩阵维度的错误。我的猜测是，在矩阵乘法中，C[i][j] += A[i][k]*B[k][j]; 是发生错误传播的地方。但是我不知道如何解决这个问题。

注意：当所有对角线元素都不为空时，三角形矩阵在理论上是可转换的。这里的问题似乎是大型随机矩阵的病态

@Damien 是的，你是对的。我的矩阵是病态的。我认为我对此很好，因为我了解错误的根本原因是什么。谢谢！！