I'm working on a program that requires calculating the inverse of an 8x8 matrix as fast as possible. Here's the code I wrote:

class matrix
{
public:
    int w, h;
    std::vector<std::vector<float>> cell;

    matrix(int width, int height)
    {
        w = width;
        h = height;
        cell.resize(width);
        for (int i = 0; i < cell.size(); i  )
        {
            cell[i].resize(height);
        }
    }
};

matrix transponseMatrix(matrix M)
{
    matrix A(M.h, M.w);

    for (int i = 0; i < M.h; i  )
    {
        for (int j = 0; j < M.w; j  )
        {
            A.cell[i][j] = M.cell[j][i];
        }
    }

    return A;
}

float getMatrixDeterminant(matrix M)
{
    if (M.w != M.h)
    {
        std::cout << "ERROR! Matrix isn't of nXn type.\n";
        return NULL;
    }
    float determinante = 0;
    if (M.w == 1)
    {
        determinante = M.cell[0][0];
    }
    if (M.w == 2)
    {
        determinante = M.cell[0][0] * M.cell[1][1] - M.cell[1][0] * M.cell[0][1];
    }
    else
    {
        for (int i = 0; i < M.w; i  )
        {
            matrix A(M.w - 1, M.h - 1);
            int cy = 0;
            for (int y = 1; y < M.h; y  )
            {
                int cx = 0;
                for (int x = 0; x < M.w; x  )
                {
                    if (x != i)
                    {
                        A.cell[cx][cy] = M.cell[x][y];
                        cx  ;
                    }
                }
                cy  ;
            }

            determinante  = M.cell[i][0] * pow(-1, i   0) * getMatrixDeterminant(A);
        }
    }

    return determinante;
}

float getComplementOf(matrix M, int X, int Y)
{
    float det;

    if (M.w != M.h)
    {
        std::cout << "ERROR! Matrix isn't of nXn type.\n";
        return NULL;
    }

    if (M.w == 2)
    {
        det = M.cell[1 - X][1 - Y];
    }
    else
    {
        matrix A(M.w - 1, M.h - 1);
        int cy = 0;
        for (int y = 0; y < M.h; y  )
        {
            if (y != Y)
            {
                int cx = 0;
                for (int x = 0; x < M.w; x  )
                {
                    if (x != X)
                    {
                        A.cell[cx][cy] = M.cell[x][y];
                        cx  ;
                    }
                }
                cy  ;
            }
        }
        det = getMatrixDeterminant(A);
    }

    return (pow(-1, X   Y) * det);
}

matrix invertMatrix(matrix M)
{
    matrix A(M.w, M.h);
    float det = getMatrixDeterminant(M);
    if (det == 0)
    {
        std::cout << "ERROR! Matrix inversion impossible (determinant is equal to 0).\n";
        return A;
    }

    for (int i = 0; i < M.h; i  )
    {
        for (int j = 0; j < M.w; j  )
        {
            A.cell[j][i] = getComplementOf(M, j, i) / det;
        }
    }

    A = transponseMatrix(A);

    return A;
}

While it does work, it does so way too slowly for my purposes, managing to calculate an 8x8 matrix's inverse about 6 times per second.

I've tried searching for more efficient ways to invert a matrix but was unsuccessfull in finding solutions for matrices of these dimensions. However I did find conversations in which people claimed that for matrices below 50x50 or even 1000x1000 time shouldn't be a problem, so I was wondering if I have missed something, either a faster method or some unnecessary calculations in my code.

Does anyone have experience regarding this and/or advice?

Sorry for broken english.

CodePudding user response：

Your implementation have problems as others commented on the question. The largest bottleneck is the algorithm itself, calculating tons of determinants.(It's O(n!)!)

If you want a simple implementation, just implement Gaussian elimination. See finding the inverse of a matrix and the pseudo code at Wikipedia. It'll perform fast enough for small sizes such as 8x8.

If you want a complex but more efficient implementation, use a library that is optimized for LU decomposition(Gaussian elimination), QR decomposition, etc.(Such as LAPACK or OpenCV.)