transpuesta - matriz inversa math2me

Eficiente matriz inversa 4x4(transformada afín)

Question

Esperaba que alguien pudiera señalar una fórmula eficiente para la transformación de matriz afín 4x4. Actualmente mi código utiliza la expansión de cofactor y asigna una matriz temporal para cada cofactor. Es fácil de leer, pero es más lento de lo que debería ser.

Tenga en cuenta que esto no es tarea y sé cómo resolverlo manualmente utilizando la expansión de co-factor 4x4, es solo un dolor y no es realmente un problema interesante para mí. También busqué en Google y encontré algunos sitios que ya le dieron la fórmula ( http://www.euclideanspace.com/maths/algebra/matrix/functions/inverse/fourD/index.htm ). Sin embargo, este probablemente podría optimizarse aún más mediante el cálculo previo de algunos de los productos. Estoy seguro de que a alguien se le ocurrió la "mejor" fórmula para esto en un momento u otro.

Answer 1

Creo que la única forma de calcular un inverso es resolver n veces la ecuación: A x = y, donde y abarca los vectores unitarios, es decir, el primero es (1,0,0,0), el segundo es (0 , 1,0,0), etc.

(Usar los cofactores (regla de Cramer) es una mala idea, a menos que desee una fórmula simbólica para el inverso).

La mayoría de las bibliotecas de álgebra lineal le permitirán resolver esos sistemas lineales, e incluso calcular un inverso. Ejemplo en python (usando numpy):

from numpy.linalg import inv
inv(A) # here you go

Answer 2

Debería poder explotar el hecho de que la matriz es afín para acelerar las cosas en una inversión completa. A saber, si su matriz se ve así

A = [ M   b  ]
    [ 0   1  ]

donde A es 4x4, M es 3x3, b es 3x1, y la fila inferior es (0,0,0,1), luego

inv(A) = [ inv(M)   -inv(M) * b ]
         [   0            1     ]

Dependiendo de su situación, puede ser más rápido calcular el resultado de inv (A) * x en lugar de formar realmente inv (A). En ese caso, las cosas se simplifican a

inv(A) * [x] = [ inv(M) * (x - b) ]
         [1] = [        1         ] 

donde x es un vector 3x1 (generalmente un punto 3D).

Por último, si M representa una rotación (es decir, sus columnas son ortonormales), entonces puede usar el hecho de que inv (M) = transponer (M). Entonces, calcular el inverso de A es simplemente una cuestión de restar el componente de traducción y multiplicar por la transposición de la parte 3x3.

Tenga en cuenta que si la matriz es ortonormal es algo que debe saber del análisis del problema. Comprobarlo durante el tiempo de ejecución sería bastante caro; aunque es posible que desee hacerlo en compilaciones de depuración para verificar que sus suposiciones se cumplan.

Espero que todo eso esté claro ...

Answer 3

En caso de que alguien quiera guardar algo de mecanografía, aquí hay una versión de AS3 que escribí en la página 9 (versión más eficiente del Teorema de Expansión de Laplace) del enlace publicado arriba por phkahler:

public function invert() : Matrix4 {
    var m : Matrix4 = new Matrix4();

    var s0 : Number = i00 * i11 - i10 * i01;
    var s1 : Number = i00 * i12 - i10 * i02;
    var s2 : Number = i00 * i13 - i10 * i03;
    var s3 : Number = i01 * i12 - i11 * i02;
    var s4 : Number = i01 * i13 - i11 * i03;
    var s5 : Number = i02 * i13 - i12 * i03;

    var c5 : Number = i22 * i33 - i32 * i23;
    var c4 : Number = i21 * i33 - i31 * i23;
    var c3 : Number = i21 * i32 - i31 * i22;
    var c2 : Number = i20 * i33 - i30 * i23;
    var c1 : Number = i20 * i32 - i30 * i22;
    var c0 : Number = i20 * i31 - i30 * i21;

    // Should check for 0 determinant

    var invdet : Number = 1 / (s0 * c5 - s1 * c4 + s2 * c3 + s3 * c2 - s4 * c1 + s5 * c0);

    m.i00 = (i11 * c5 - i12 * c4 + i13 * c3) * invdet;
    m.i01 = (-i01 * c5 + i02 * c4 - i03 * c3) * invdet;
    m.i02 = (i31 * s5 - i32 * s4 + i33 * s3) * invdet;
    m.i03 = (-i21 * s5 + i22 * s4 - i23 * s3) * invdet;

    m.i10 = (-i10 * c5 + i12 * c2 - i13 * c1) * invdet;
    m.i11 = (i00 * c5 - i02 * c2 + i03 * c1) * invdet;
    m.i12 = (-i30 * s5 + i32 * s2 - i33 * s1) * invdet;
    m.i13 = (i20 * s5 - i22 * s2 + i23 * s1) * invdet;

    m.i20 = (i10 * c4 - i11 * c2 + i13 * c0) * invdet;
    m.i21 = (-i00 * c4 + i01 * c2 - i03 * c0) * invdet;
    m.i22 = (i30 * s4 - i31 * s2 + i33 * s0) * invdet;
    m.i23 = (-i20 * s4 + i21 * s2 - i23 * s0) * invdet;

    m.i30 = (-i10 * c3 + i11 * c1 - i12 * c0) * invdet;
    m.i31 = (i00 * c3 - i01 * c1 + i02 * c0) * invdet;
    m.i32 = (-i30 * s3 + i31 * s1 - i32 * s0) * invdet;
    m.i33 = (i20 * s3 - i21 * s1 + i22 * s0) * invdet;

    return m;
}

Esto produjo con éxito una matriz de identidad cuando multipliqué varias matrices de transformación 3D por el inverso devuelto por este método. Estoy seguro de que puede buscar / reemplazar para obtener esto en el idioma que desee.

Answer 4

IIRC puede reducir en gran medida el código y el tiempo precompilando un manojo (12?) 2x2 determinantes. Divida la matriz por la mitad verticalmente y calcule cada 2x2 tanto en la mitad superior como en la inferior. Uno de estos determinantes más pequeños se usa en cada término que necesitará para el mayor cálculo y cada uno se vuelve a utilizar.

Además, no use una función determinante separada: reutilice los subdeterminantes que calculó para el adjunto para obtener el determinante.

Oh, acabo de encontrar this.

Hay algunas mejoras que puedes hacer sabiendo que también es cierto tipo de transformación.

Answer 5

Para dar seguimiento a las excelentes respuestas de pkhaler y Robin Hilliard anteriores, aquí está el código ActionScript 3 de Robin convertido en un método C #. Con suerte, esto puede ahorrar algo de tipeo para otros desarrolladores de C #, así como para desarrolladores de C / C ++ y Java que necesitan una función de inversión de matriz 4x4:

public static double[,] GetInverse(double[,] a)
{
    var s0 = a[0, 0] * a[1, 1] - a[1, 0] * a[0, 1];
    var s1 = a[0, 0] * a[1, 2] - a[1, 0] * a[0, 2];
    var s2 = a[0, 0] * a[1, 3] - a[1, 0] * a[0, 3];
    var s3 = a[0, 1] * a[1, 2] - a[1, 1] * a[0, 2];
    var s4 = a[0, 1] * a[1, 3] - a[1, 1] * a[0, 3];
    var s5 = a[0, 2] * a[1, 3] - a[1, 2] * a[0, 3];

    var c5 = a[2, 2] * a[3, 3] - a[3, 2] * a[2, 3];
    var c4 = a[2, 1] * a[3, 3] - a[3, 1] * a[2, 3];
    var c3 = a[2, 1] * a[3, 2] - a[3, 1] * a[2, 2];
    var c2 = a[2, 0] * a[3, 3] - a[3, 0] * a[2, 3];
    var c1 = a[2, 0] * a[3, 2] - a[3, 0] * a[2, 2];
    var c0 = a[2, 0] * a[3, 1] - a[3, 0] * a[2, 1];

    // Should check for 0 determinant
    var invdet = 1.0 / (s0 * c5 - s1 * c4 + s2 * c3 + s3 * c2 - s4 * c1 + s5 * c0);

    var b = new double[4, 4];

    b[0, 0] = ( a[1, 1] * c5 - a[1, 2] * c4 + a[1, 3] * c3) * invdet;
    b[0, 1] = (-a[0, 1] * c5 + a[0, 2] * c4 - a[0, 3] * c3) * invdet;
    b[0, 2] = ( a[3, 1] * s5 - a[3, 2] * s4 + a[3, 3] * s3) * invdet;
    b[0, 3] = (-a[2, 1] * s5 + a[2, 2] * s4 - a[2, 3] * s3) * invdet;

    b[1, 0] = (-a[1, 0] * c5 + a[1, 2] * c2 - a[1, 3] * c1) * invdet;
    b[1, 1] = ( a[0, 0] * c5 - a[0, 2] * c2 + a[0, 3] * c1) * invdet;
    b[1, 2] = (-a[3, 0] * s5 + a[3, 2] * s2 - a[3, 3] * s1) * invdet;
    b[1, 3] = ( a[2, 0] * s5 - a[2, 2] * s2 + a[2, 3] * s1) * invdet;

    b[2, 0] = ( a[1, 0] * c4 - a[1, 1] * c2 + a[1, 3] * c0) * invdet;
    b[2, 1] = (-a[0, 0] * c4 + a[0, 1] * c2 - a[0, 3] * c0) * invdet;
    b[2, 2] = ( a[3, 0] * s4 - a[3, 1] * s2 + a[3, 3] * s0) * invdet;
    b[2, 3] = (-a[2, 0] * s4 + a[2, 1] * s2 - a[2, 3] * s0) * invdet;

    b[3, 0] = (-a[1, 0] * c3 + a[1, 1] * c1 - a[1, 2] * c0) * invdet;
    b[3, 1] = ( a[0, 0] * c3 - a[0, 1] * c1 + a[0, 2] * c0) * invdet;
    b[3, 2] = (-a[3, 0] * s3 + a[3, 1] * s1 - a[3, 2] * s0) * invdet;
    b[3, 3] = ( a[2, 0] * s3 - a[2, 1] * s1 + a[2, 2] * s0) * invdet;

    return b;
}