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.
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.)