Inverser matrice 4x4 - Solution numérique la plus stable nécessaire

Question

Je veux inverser une matrice 4x4. Mes numéros sont stockés au format virgule fixe (le 1.15.16 pour être exact).

Avec l'arithmétique en virgule flottante, je construis généralement la matrice adjointe et me divise par le déterminant (par exemple, la force brute de la solution). Cela a fonctionné pour moi jusqu’à présent, mais lorsqu’il s’agit de nombres à points fixes, j’ai une perte de précision inacceptable en raison de toutes les multiplications utilisées.

Remarque: dans l'arithmétique en virgule fixe, je rejette toujours certains des résultats les moins significatifs des résultats immédiats.

Alors - Quelle est la manière la plus stable et numérique d’inverser une matrice? Les performances ne me dérangent pas beaucoup, mais aller simplement en virgule flottante reviendrait à ralentir mon architecture cible.

Answer 1

Je pense que la réponse à cela dépend de la forme exacte de la matrice. Une méthode de décomposition standard (LU, QR, Cholesky, etc.) avec pivot (un élément essentiel) est assez bonne sur un point fixe, en particulier pour une petite matrice 4x4. Voir le livre 'Recettes numériques' de Press et al. pour une description de ces méthodes.

Cet article donne quelques algorithmes utiles, mais est derrière un paywall malheureusement. Ils recommandent une décomposition (pivotée) de Cholesky avec des fonctionnalités supplémentaires trop compliquées pour être listées ici.

Answer 2

Méta-réponse: Est-ce vraiment une matrice 4x4 générale? Si votre matrice a une forme spéciale, il existe des formules directes pour inverser qui seraient rapides et garder le compte à rebours de votre opération.

Par exemple, s'il s'agit d'une transformation de coordonnées homogène standard à partir de graphiques, telle que:

[ux vx wx tx]
[uy vy wy ty]
[uz vz wz tz]
[ 0  0  0  1]

(en supposant une composition de matrices de rotation, d’échelle et de translation)

Ensuite, il existe une formule directe facile à déduire , qui est

[ux uy uz -dot(u,t)]
[vx vy vz -dot(v,t)]
[wx wy wz -dot(w,t)]
[ 0  0  0     1    ]

(matrices ASCII volées à partir de la page liée.)

Vous ne pouvez probablement pas battre cela pour la perte de précision en point fixe.

Si votre matrice provient d'un domaine où vous savez qu'elle est plus structurée, la réponse sera probablement simple.

Answer 3

Je voudrais appuyer la question posée par Jason S: êtes-vous certain de devoir inverser votre matrice? Ce n'est presque jamais nécessaire. Non seulement cela, c'est souvent une mauvaise idée. Si vous avez besoin de résoudre Ax = b, il est plus stable numériquement de résoudre le système directement que de multiplier b par A inverse.

Même si vous devez résoudre Ax = b à plusieurs reprises pour de nombreuses valeurs de b, ce n'est pas une bonne idée d'inverser A. Vous pouvez factoriser A (par exemple, la factorisation LU ou la factorisation de Cholesky) et sauvegardez les facteurs pour ne pas refaire ce travail à chaque fois, mais résolvez toujours le système à chaque fois en utilisant la factorisation.

Answer 4

Permettez-moi de poser une question différente: avez-vous vraiment besoin d'inverser la matrice (appelez-la M), ou avez-vous besoin d'utiliser l'inverse de la matrice pour résoudre d'autres équations? (par exemple, Mx = b pour M connu, b) Il existe souvent d'autres moyens de le faire sans avoir explicitement besoin de calculer l'inverse. Ou si la matrice M est fonction du temps & amp; cela change lentement, vous pouvez alors calculer l'inverse complet une fois, & amp; il existe des moyens itératifs de le mettre à jour.

Answer 5

Pour minimiser les erreurs de troncature et autres dommages, utilisez "pivotant". - voir le chapitre sur l'inversion des matrices dans Recettes numériques. Ils ont la meilleure explication que j'ai trouvée jusqu'à présent.

Answer 6

Vous pouvez envisager de doubler à 1,31 avant d'utiliser votre algorithme normal. Cela doublera le nombre de multiplications, mais vous effectuez une inversion de matrice et tout ce que vous ferez sera assez lié au multiplicateur de votre processeur.

Pour tous ceux qui souhaitent trouver les équations d'un inverseur 4x4, vous pouvez utiliser un logiciel mathématique symbolique pour les résoudre à votre place. La TI-89 le fera même, même si cela prendra plusieurs minutes.

Si vous nous donnez une idée de ce que la matrice inverse fait pour vous et de la manière dont elle s'intègre au reste de votre traitement, nous pourrons peut-être vous suggérer des alternatives.

-Adam

Answer 7

Une élimination simple, ancienne et gaussienne, fonctionnerait bien.

Cela dépend des bibliothèques / classes / structures que vous utilisez. Vous pouvez consulter la GSL .

Answer 8

Si la matrice représente une transformation affine (c'est souvent le cas avec les matrices 4x4 tant que vous n'introduisez pas de composant de mise à l'échelle), l'inverse est simplement la transposition de la partie de rotation supérieure 3x3 avec la dernière colonne inversée. Si vous avez besoin d’une solution généralisée, il est évident que l’élimination gaussienne est probablement la plus simple.