Invertir matriz 4x4: la solución numérica más estable que se necesita

Question

Quiero invertir una matriz 4x4. Mis números se almacenan en formato de punto fijo (1.15.16 para ser exactos).

Con la aritmética de punto flotante, por lo general, solo construyo la matriz adjunta y la divido por el determinante (por ejemplo, fuerza bruta de la solución). Eso me funcionó hasta ahora, pero al tratar con números de puntos fijos obtengo una pérdida de precisión inaceptable debido a todas las multiplicaciones utilizadas.

Nota: en aritmética de punto fijo, siempre elimino algunos de los bits menos significativos de resultados inmediatos.

Entonces, ¿cuál es la forma más numérica y estable de invertir una matriz? No me importa mucho el rendimiento, pero simplemente ir a punto flotante sería ralentizar mi arquitectura de destino.

Answer 1

Creo que la respuesta a esto depende de la forma exacta de la matriz. Un método de descomposición estándar (LU, QR, Cholesky, etc.) con pivote (un elemento esencial) es bastante bueno en el punto fijo, especialmente para una matriz pequeña de 4x4. Consulte el libro 'Recetas numéricas' de Press et al. para una descripción de estos métodos.

Este documento da algunos algoritmos útiles, pero está detrás de un muro de pago desafortunadamente. Recomiendan una descomposición de Cholesky (pivotada) con algunas características adicionales demasiado complicadas para enumerar aquí.

Answer 2

Meta-respuesta: ¿Es realmente una matriz general de 4x4? Si su matriz tiene una forma especial, entonces hay fórmulas directas para invertir que serían rápidas y mantendrán su cuenta regresiva de la operación.

Por ejemplo, si se trata de una transformación de coordenadas homogéneas estándar de los gráficos, como:

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

(asumiendo una composición de rotación, escala, matrices de traducción)

luego hay un fórmula directa fácilmente derivable , que es

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

(Matrices ASCII robadas de la página enlazada.)

Probablemente no puedas superar eso por pérdida de precisión en un punto fijo.

Si su matriz proviene de algún dominio donde sabe que tiene más estructura, entonces es probable que haya una respuesta fácil.

Answer 3

Me gustaría secundar la pregunta que Jason S planteó: ¿estás seguro de que necesitas invertir tu matriz? Esto casi nunca es necesario. No solo eso, a menudo es una mala idea. Si necesita resolver Ax = b, es más estable numéricamente resolver el sistema directamente que multiplicar b por A inverso.

Incluso si tiene que resolver Ax = b una y otra vez para muchos valores de b, todavía no es una buena idea invertir A. Puede factor A (por ejemplo, factorización LU o factorización Cholesky) y guarda los factores para que no estés rehiciendo ese trabajo cada vez, pero aún así resolverías el sistema cada vez que usas la factorización.

Answer 4

Permítame hacer una pregunta diferente: ¿definitivamente necesita invertir la matriz (llámela M), o necesita usar la matriz inversa para resolver otras ecuaciones? (por ejemplo, Mx = b para M conocido, b) A menudo hay otras formas de hacer esto sin necesidad explícita de calcular el inverso. O si la matriz M es una función del tiempo & amp; cambia lentamente y luego puede calcular el inverso completo una vez, & amp; hay formas iterativas de actualizarlo.

Answer 5

Para minimizar los errores de truncamiento y otros daños, use " pivotando " - vea el capítulo sobre la inversión de matrices en Recetas Numéricas. Tienen la mejor explicación que he encontrado hasta ahora.

Answer 6

Puede considerar duplicar a 1.31 antes de hacer su algoritmo normal. Doblará el número de multiplicaciones, pero estás haciendo una inversión de matriz y todo lo que hagas estará bastante ligado al multiplicador de tu procesador.

Para cualquier persona interesada en encontrar las ecuaciones para un 4x4 invertido, puede usar un paquete de matemáticas simbólicas para resolverlas por usted. La TI-89 lo hará incluso, aunque llevará varios minutos.

Si nos da una idea de lo que hace la inversión de matriz por usted, y cómo encaja con el resto de su procesamiento, podríamos sugerir alternativas.

-Adam

Answer 7

La simple eliminación gaussiana funcionaría bien.

Depende de las bibliotecas / clases / estructuras que estés usando. Puede echar un vistazo a la GSL .

Answer 8

Si la matriz representa una transformación afín (muchas veces este es el caso con matrices 4x4 siempre que no introduzca un componente de escala), lo inverso es simplemente la transposición de la parte superior de rotación 3x3 con la última columna negada. Obviamente, si necesita una solución generalizada, buscar la eliminación de Gauss es probablemente la más fácil.