Inverti matrice 4x4 - Necessaria soluzione numerica più stabile

Question

Voglio invertire una matrice 4x4. I miei numeri sono memorizzati in formato a virgola fissa (1.15.16 per l'esattezza).

Con l'aritmetica in virgola mobile di solito costruisco solo la matrice aggiunta e la divido per il determinante (ad esempio forza bruta la soluzione). Finora ha funzionato per me, ma quando ho a che fare con numeri a virgola fissa ottengo una inaccettabile perdita di precisione dovuta a tutte le moltiplicazioni utilizzate.

Nota: nell'aritmetica a virgola fissa butto sempre via alcuni dei bit meno significativi di risultati immediati.

Quindi - Qual è il modo più numerico stabile per invertire una matrice? Non mi preoccupo molto delle prestazioni, ma semplicemente andare in virgola mobile sarebbe rallentare la mia architettura di destinazione.

Answer 1

Penso che la risposta a questo dipenda dalla forma esatta della matrice. Un metodo di decomposizione standard (LU, QR, Cholesky ecc.) Con pivot (un elemento essenziale) è abbastanza buono su un punto fisso, specialmente per una matrice 4x4 piccola. Vedi il libro "Ricette numeriche" di Press et al. per una descrizione di questi metodi.

Questo documento fornisce alcuni utili algoritmi, ma sfortunatamente è dietro un paywall. Raccomandano una decomposizione (imperniata) di Cholesky con alcune funzionalità aggiuntive troppo complicate per essere elencate qui.

Answer 2

Meta-risposta: è davvero una matrice 4x4 generale? Se la tua matrice ha una forma speciale, allora ci sono formule dirette per l'inversione che sarebbero veloci e mantengono il conto alla rovescia dell'operazione.

Ad esempio, se si tratta di una trasformazione di coordinate omogenea standard dalla grafica, come:

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

(assumendo una composizione di rotazione, scala, matrici di traduzione)

poi c'è una formula diretta facilmente derivabile , che è

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

(matrici ASCII rubate dalla pagina collegata.)

Probabilmente non puoi batterlo per perdita di precisione in punto fisso.

Se la tua matrice proviene da un dominio in cui sai che ha più struttura, allora probabilmente ci sarà una risposta facile.

Answer 3

Vorrei rispondere alla seconda domanda posta da Jason S: sei sicuro di dover invertire la tua matrice? Questo non è quasi mai necessario. Non solo, è spesso una cattiva idea. Se è necessario risolvere Ax = b, è più numericamente stabile risolvere direttamente il sistema che moltiplicare b per A inverso.

Anche se devi risolvere Ax = b più e più volte per molti valori di b, non è comunque una buona idea invertire A. Puoi fattorizzare A (ad esempio fattorizzazione LU o fattorizzazione di Cholesky) e salva i fattori in modo da non ripetere il lavoro ogni volta, ma risolveresti comunque il sistema ogni volta utilizzando la fattorizzazione.

Answer 4

Lasciami fare una domanda diversa: devi assolutamente invertire la matrice (chiamala M) o devi usare l'inverso della matrice per risolvere altre equazioni? (ad esempio Mx = b per M nota, b) Spesso ci sono altri modi per farlo senza che sia necessario esplicitamente calcolare l'inverso. O se la matrice M è una funzione del tempo & amp; cambia lentamente quindi puoi calcolare l'intero inverso una volta, & amp; ci sono modi iterativi per aggiornarlo.

Answer 5

Per ridurre al minimo gli errori di troncamento e altri problemi, utilizzare "pivoting" - vedere il capitolo sull'inversione delle matrici in Ricette numeriche. Hanno la migliore spiegazione che ho trovato finora.

Answer 6

Potresti considerare di raddoppiare a 1,31 prima di eseguire il tuo normale algoritmo. Raddoppierà il numero di moltiplicazioni, ma stai facendo un'inversione di matrice e qualsiasi cosa tu faccia sarà abbastanza legata al moltiplicatore nel tuo processore.

Per chiunque sia interessato a trovare le equazioni per un'inversione 4x4, è possibile utilizzare un pacchetto matematico simbolico per risolverle. La TI-89 lo farà anche, anche se ci vorranno diversi minuti.

Se ci dai un'idea di ciò che l'inversione di matrice fa per te e di come si adatta al resto della tua elaborazione, potremmo essere in grado di suggerire alternative.

-Adam

Answer 7

La semplice eliminazione gaussiana funzionerebbe bene.

Dipende dalle librerie / classi / strutture che stai utilizzando. Puoi dare un'occhiata al GSL .

Answer 8

Se la matrice rappresenta una trasformazione affine (molte volte questo è il caso delle matrici 4x4 fintanto che non si introduce un componente di ridimensionamento) l'inverso è semplicemente la trasposizione della parte di rotazione 3x3 superiore con l'ultima colonna negata. Ovviamente se hai bisogno di una soluzione generalizzata, esaminare l'eliminazione gaussiana è probabilmente il più semplice.