Erreur très élevée lors du calcul de l'inverse de la matrice en utilisant l'élimination de Gauss

Question

Je travaille actuellement sur une base de code C++ qui utilise une bibliothèque de matrices pour calculer diverses choses. Une de ces choses est de calculer l'inverse d'une matrice. Il utilise gauss elimation pour y parvenir. Mais le résultat est très inexact. Si bien que multiplier la matrice inverse avec la matrice originale n'est même pas proche de la matrice identitaire.

Voici le code qui est utilisé pour calculer l'inverse, la matrice est templated sur un type numérique et les lignes et colonnes:

/// \brief Take the inverse of the matrix. 
/// \return A new matrix which is the inverse of the current one. 
matrix<T, M, M> inverse() const 
{ 
    static_assert(M == N, "Inverse matrix is only defined for square matrices."); 

    // augmented the current matrix with the identiy matrix. 
    auto augmented = this->augment(matrix<T, M, M>::get_identity()); 

    for (std::size_t i = 0; i < M; i++) 
    { 
     // divide the current row by the diagonal element. 
     auto divisor = augmented[i][i]; 
     for (std::size_t j = 0; j < 2 * M; j++) 
     { 
      augmented[i][j] /= divisor; 
     } 

     // For each element in the column of the diagonal element that is currently selected 
     // set all element in that column to 0 except the diagonal element by using the currently selected row diagonal element. 

     for (std::size_t j = 0; j < M; j++) 
     { 
      if (i == j) 
      { 
       continue; 
      } 

      auto multiplier = augmented[j][i]; 
      for (std::size_t k = 0; k < 2 * M; k++) 
      { 
       augmented[j][k] -= multiplier * augmented[i][k]; 
      } 
     } 
    } 

    // Slice of the the new identity matrix on the left side. 
    return augmented.template slice<0, M, M, M>(); 
} 

Maintenant je l'ai fait un test unitaire qui test si l'inverse est corriger en utilisant des valeurs pré-calculées. J'essaie deux matrices un 3x3 et un 4x4. J'ai utilisé ce site pour calculer l'inverse: https://matrix.reshish.com/ et ils correspondent à un certain degré. puisque le test unitaire réussit. Mais une fois que je calcule la matrice d'origine *, rien ne ressemble à une matrice d'identité. Voir le commentaire dans le code ci-dessous.

BOOST_AUTO_TEST_CASE(matrix_inverse) 
{ 
    auto m1 = matrix<double, 3, 3>({ 
     {7, 8, 9}, 
     {10, 11, 12}, 
     {13, 14, 15} 
    }); 

    auto inverse_result1 = matrix<double,3, 3>({ 
     {264917625139441.28, -529835250278885.3, 264917625139443.47}, 
     {-529835250278883.75, 1059670500557768, -529835250278884.1}, 
     {264917625139442.4, -529835250278882.94, 264917625139440.94} 
    }); 

    auto m2 = matrix<double, 4, 4>({ 
     {7, 8, 9, 23}, 
     {10, 11, 12, 81}, 
     {13, 14, 15, 11}, 
     {1, 73, 42, 65} 
    }); 

    auto inverse_result2 = matrix<double, 4, 4>({ 
     {-0.928094660194201, 0.21541262135922956, 0.4117111650485529, -0.009708737864078209}, 
     {-0.9641231796116679, 0.20979975728155775, 0.3562651699029188, 0.019417475728154842}, 
     {1.7099261731391882, -0.39396237864078376, -0.6169346682848 , -0.009708737864076772 }, 
     {-0.007812499999999244, 0.01562499999999983, -0.007812500000000278, 0} 
    }); 


    // std::cout << (m1.inverse() * m1) << std::endl; 
    // results in 
    // 0.500000000  1.000000000  -0.500000000 
    // 1.000000000  0.000000000  0.500000000 
    // 0.500000000  -1.000000000 1.000000000 
    // std::cout << (m2.inverse() * m2) << std::endl; 
    // results in 
    // 0.396541262  -0.646237864 -0.689016990 -2.162317961 
    // 1.206917476  2.292475728  1.378033981  3.324635922 
    // -0.884708738 -0.958737864 -0.032766990 -3.756067961 
    // -0.000000000 -0.000000000 -0.000000000 1.000000000 

    BOOST_REQUIRE_MESSAGE(
     m1.inverse().fuzzy_equal(inverse_result1, 0.1) == true, 
     "3x3 inverse is not the expected result." 
    ); 

    BOOST_REQUIRE_MESSAGE(
     m2.inverse().fuzzy_equal(inverse_result2, 0.1) == true, 
     "4x4 inverse is not the expected result." 
    ); 
} 

Je suis à bout de nerfs. Je ne suis en aucun cas un spécialiste des mathématiques matricielles puisque j'ai dû apprendre tout cela sur le tas mais cela me bouscule vraiment.

La classe matrice de code complet est disponible à l'adresse: https://codeshare.io/johnsmith

ligne 404 est où la fonction inverse est situé.

Toute aide est appréciée.

Answer 1

Comme déjà établi dans les commentaires la matrice d'intérêt est singulière et donc il n'y a pas d'inverse.

Très bien, vos tests ont déjà trouvé le premier problème dans le code - ce cas n'est pas géré correctement, aucune erreur n'est soulevée.

Le plus gros problème est, que ce n'est pas facile à détecter: S'il n'y avait pas d'erreurs dues aux erreurs d'arrondi, ce serait un morceau de pièce - il suffit de tester ce diviseur n'est pas 0! Mais il y a des erreurs d'arrondi dans les opérations flottantes, donc divisor sera un très petit nombre non nul.

Et il n'y a aucun moyen de dire, si cette valeur non nulle en raison d'erreurs d'arrondi ou au fait que la matrice est proche du singulier (mais pas singulier). Cependant, si la matrice est proche du singulier, elle est en mauvais état et les résultats ne peuvent donc pas être approuvés. Donc idéalement, l'algorithme devrait non seulement calculer l'inverse, mais aussi (estimer) l'état de la matrice originale, afin que l'appelant puisse réagir à une mauvaise condition.

Il est probablement judicieux d'utiliser des bibliothèques bien connues et testées pour ce genre de calcul - il y a beaucoup à considérer et ce qui peut être mal fait.