Matrix inverse calculation of upper triangular matrix gives error for large matrix dimensions

Question

I have a recursive function to calculate the inverse of an upper triangular matrix. I have divided the matrix into Top, Bottom and Corner sections and then followed the methodology as laid down in https://math.stackexchange.com/a/2333418. Here is a pseudocode form:

//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)
}

For simplicity I have implemented the code such that only power of 2 dimension matrices are supported.

My problem is that I have tested this code for matrix dimensions of 2, 4, 8, 16, 32 and 64. All of these pass all of the assertion checks as shown in code. But for matrix dimension of 128 I get failure is the assertion in main(). And when I check, I observer that the test_matrix is not Identity matrix. Some non-diagonal elements are not equal to 0. I am wondering what could be the reason for this:-

1. I am using C++ std::vector<std::vector<double>> for Matrix representation.
2. Since the data type is double the non-diagonal elements of test_matrix for cases 2, 4, 8, ..., 64 do have some value but very small. For example, -9.58122e-14
3. All my matrices at any recursion stage are square matrix
4. I am performing checks that Top*Top_inv = Identity and Bottom*Bottom_inv = Identity.
5. Finally, for dimensions 2, 4, ..., 64 I generated random numbers(b/w 1 and 10) to create my upper triangular matrix. Since, these cases passed, I guess my mathematical implementation is correct.

I feel like there is some aspect of C++ datatypes about double which I am unaware of that could be causing the error. Otherwise the sudden error from 64->128 doesn't make sense.

Answer 1

Could you please elaborate on how the matrix == identity operation is implemented? My guess is that the problem might be resumed to the floating point comparison.

The matrix inversion can be O(n^3) in the worst case. This means that, as the matrix size increases, the amount of computations involved also increase. Real numbers cannot be perfectly represented even when using 64 bit floating point, they are always an approximation.

For operations such as matrix inversion this can cause problems of numerical error propagation, due to the loss of precision on the accumulated multiply adds operations.

For this, there has been discussions already in the StackOverflow: How should I do floating point comparison?

EDIT: Other thing to consider if the full matrix is actually invertible. Perhaps the Top and/or Bottom matrices are invertible, but the full matrix (when composing with the Corner matrix) is not.