Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Matrix inverse accuracy

I have a large world, about 5,000,000 x 1,000,000 units. The camera can be near some object or far enough as to see the whole world.
I get the mouse position in world coordinates by unprojecting (Z comes from depth buffer). The problem is that it involves a matrix inverse. When using big and small numbers (e.g. translating away from origin and scaling to see more world) at the same time, the calculations become unstable.

Trying to see the accuracy of this inverse matrix I look at the determinant. Ideally it never will be zero, due to the nature of transformation matrices. I know that being 'det' a small value means nothing on its own, it can be due to small values in the matrix. But it can also be a sign of numbers becoming wrong.

I also know I can calculate the inverse by inverting each transformation and multiplicaying them. Does it provide more accuracy?

How can I tell if my matrix is getting degenerated, suffer numerical issues?

like image 301
Ripi2 Avatar asked Feb 16 '17 22:02

Ripi2


People also ask

How can you tell if the inverse of a matrix is correct?

Multiply your test-inverse with your original matrix. This should give the identity matrix. Check it for more and more of the rows for increasing level of certainty.

Is matrix inversion expensive?

Inverting a matrix But the basic idea is that the inversion itself is not faster than the LU decomposition, which effectively solves for x. As we can see inverting a matrix to solve Ax=b is roughly three times as expensive as directly solving for x.

Why is LU decomposition better than inverse?

Normally you don't need/want to obtain the inverse of an matrix, because it is costly and many times unnecessary. The LU decomposition can be used with its necessary to solve a problem with many right hand sides. This situation often arise when solving PDE's using the finite differences method.

When can matrix not be inverted?

If the determinant of the matrix is zero, then it will not have an inverse; the matrix is then said to be singular. Only non-singular matrices have inverses.


1 Answers

for starters see Understanding 4x4 homogenous transform matrices

  1. Improving accuracy for cumulative matrices (Normalization)

    To avoid degeneration of transform matrix select one axis as main. I usually chose Z as it is usually view or forward direction in my apps. Then exploit cross product to recompute/normalize the rest of axises (which should be perpendicular to each other and unless scale is used then also unit size). This can be done only for orthonormal matrices so no skew or projections ... Orthogonal matrices must be scaled to orthonormal then inverted and then scaled back to make this usable.

    You do not need to do this after every operation just make a counter of operations done on each matrix and if some threshold crossed then normalize it and reset counter.

    To detect degeneration of such matrices you can test for orthogonality by dot product between any two axises (should be zero or very near it). For orthonormal matrices you can test also for unit size of axis direction vectors ...

    Here is how my transform matrix normalization looks like (for orthonormal matrices) in C++:

    double reper::rep[16]; // this is my transform matrix stored as member in `reper` class
    //---------------------------------------------------------------------------
    void reper::orto(int test) // test is for overiding operation counter
    {
        double   x[3],y[3],z[3]; // space for axis direction vectors
        if ((cnt>=_reper_max_cnt)||(test)) // if operations count reached or overide
        {
            axisx_get(x);      // obtain axis direction vectors from matrix
            axisy_get(y);
            axisz_get(z);
            vector_one(z,z);   // Z = Z / |z|
            vector_mul(x,y,z); // X = Y x Z  ... perpendicular to y,z
            vector_one(x,x);   // X = X / |X|
            vector_mul(y,z,x); // Y = Z x X  ... perpendicular to z,x
            vector_one(y,y);   // Y = Y / |Y|
            axisx_set(x);      // copy new axis vectors into matrix
            axisy_set(y);
            axisz_set(z);
            cnt=0;             // reset operation counter
        }
    }
    
    //---------------------------------------------------------------------------
    void reper::axisx_get(double *p)
    {
        p[0]=rep[0];
        p[1]=rep[1];
        p[2]=rep[2];
    }
    //---------------------------------------------------------------------------
    void reper::axisx_set(double *p)
    {
        rep[0]=p[0];
        rep[1]=p[1];
        rep[2]=p[2];
        cnt=_reper_max_cnt; // pend normalize in next operation that needs it
    }
    //---------------------------------------------------------------------------
    void reper::axisy_get(double *p)
    {
        p[0]=rep[4];
        p[1]=rep[5];
        p[2]=rep[6];
    }
    //---------------------------------------------------------------------------
    void reper::axisy_set(double *p)
    {
        rep[4]=p[0];
        rep[5]=p[1];
        rep[6]=p[2];
        cnt=_reper_max_cnt; // pend normalize in next operation that needs it
    }
    //---------------------------------------------------------------------------
    void reper::axisz_get(double *p)
    {
        p[0]=rep[ 8];
        p[1]=rep[ 9];
        p[2]=rep[10];
    }
    //---------------------------------------------------------------------------
    void reper::axisz_set(double *p)
    {
        rep[ 8]=p[0];
        rep[ 9]=p[1];
        rep[10]=p[2];
        cnt=_reper_max_cnt; // pend normalize in next operation that needs it
    }
    //---------------------------------------------------------------------------
    

    The vector operations looks like this:

    void  vector_one(double *c,double *a)
    {
        double l=divide(1.0,sqrt((a[0]*a[0])+(a[1]*a[1])+(a[2]*a[2])));
        c[0]=a[0]*l;
        c[1]=a[1]*l;
        c[2]=a[2]*l;
    }
    
    void  vector_mul(double *c,double *a,double *b)
    {
        double   q[3];
        q[0]=(a[1]*b[2])-(a[2]*b[1]);
        q[1]=(a[2]*b[0])-(a[0]*b[2]);
        q[2]=(a[0]*b[1])-(a[1]*b[0]);
        for(int i=0;i<3;i++) c[i]=q[i];
    }
    
  2. Improving accuracy for non cumulative matrices

    Your only choice is use at least double accuracy of your matrices. Safest is to use GLM or your own matrix math based at least on double data type (like my reper class).

    Cheap alternative is using double precision functions like

    glTranslated
    glRotated
    glScaled
    ...
    

    which in some cases helps but is not safe as OpenGL implementation can truncate it to float. Also there are no 64 bit HW interpolators yet so all iterated results between pipeline stages are truncated to floats.

    Sometimes relative reference frame helps (so keep operations on similar magnitude values) for example see:

    ray and ellipsoid intersection accuracy improvement

    Also In case you are using own matrix math functions you have to consider also the order of operations so you always lose smallest amount of accuracy possible.

  3. Pseudo inverse matrix

    In some cases you can avoid computing of inverse matrix by determinants or Horner scheme or Gauss elimination method because in some cases you can exploit the fact that Transpose of orthonormal rotational matrix is also its inverse. Here is how it is done:

    void matrix_inv(GLfloat *a,GLfloat *b) // a[16] = Inverse(b[16])
    {
        GLfloat x,y,z;
        // transpose of rotation matrix
        a[ 0]=b[ 0];
        a[ 5]=b[ 5];
        a[10]=b[10];
        x=b[1]; a[1]=b[4]; a[4]=x;
        x=b[2]; a[2]=b[8]; a[8]=x;
        x=b[6]; a[6]=b[9]; a[9]=x;
        // copy projection part
        a[ 3]=b[ 3];
        a[ 7]=b[ 7];
        a[11]=b[11];
        a[15]=b[15];
        // convert origin: new_pos = - new_rotation_matrix * old_pos
        x=(a[ 0]*b[12])+(a[ 4]*b[13])+(a[ 8]*b[14]);
        y=(a[ 1]*b[12])+(a[ 5]*b[13])+(a[ 9]*b[14]);
        z=(a[ 2]*b[12])+(a[ 6]*b[13])+(a[10]*b[14]);
        a[12]=-x;
        a[13]=-y;
        a[14]=-z;
    }
    

    So rotational part of the matrix is transposed, projection stays as was and origin position is recomputed so A*inverse(A)=unit_matrix This function is written so it can be used as in-place so calling

    GLfloat a[16]={values,...}
    matrix_inv(a,a);
    

    lead to valid results too. This way of computing Inverse is quicker and numerically safer as it pends much less operations (no recursion or reductions no divisions). Of coarse this works only for orthonormal homogenuous 4x4 matrices !!!*

  4. Detection of wrong inverse

    So if you got matrix A and its inverse B then:

    A*B = C = ~unit_matrix
    

    So multiply both matrices and check for unit matrix...

    • abs sum of all non diagonal elements of C should be close to 0.0
    • all diagonal elements of C should be close to +1.0
like image 60
Spektre Avatar answered Oct 03 '22 23:10

Spektre