# About the problem of seeking matrix inverse

Recently I have solved the inverse matrix in three ways, found that their work efficiency has obvious differences.

This is the case:

typedef signed char int8;
typedef signed short int16;
typedef signed int int32;
typedef unsigned char uint8;
typedef unsigned short uint16;
typedef unsigned int uint32;
typedef float float32;
typedef double float64;

enum
{
MINORS,
COFACTORS,
};

template <int32 COLS, int32 ROWS, typename T>
class Matrix
{
puiblic:

T elements[COLS * ROWS];
//
Matrix() :
_columns(COLS),
_rows(ROWS),
_size(COLS * ROWS) {
identity();
}

//
void set(int32 col, int32 row, T value) {
elements[_columns * row + col] = value;
}

//
const T& get(int32 col, int32 row) const {
return elements[_columns * row + col];
}

// Minors, Cofactors and Adjugate solving the inverse matrix:

Matrix<COLS, COLS, T> getConversion(uint32 mode) const {

Matrix<COLS, COLS, T> m;
Matrix<COLS, COLS, T> temp;

for (int32 i = 0; i < COLS; i++) {

for (int32 j = 0; j < COLS; j++) {
// calculate the matrix of minors.
cofactor(temp, j, i, COLS);

T det = temp.determinantAid(COLS - 1);

// calculate the matrix of cofactors.
if (mode > MINORS)
det *= ((i + j) % 2 == 0) ? 1 : -1;

// calculate the matrix of adjugate.
m.set(i, j, det);
else
m.set(j, i, det);

}
}

return m;
}

// Calculates the inverse of matrix.
bool invert() {

assert(_columns == _rows);

T det = determinant();
if (det == 0) {
identity();
return false;
}

det = 1.0f / det;

// find inverse using formula "inverse(A) = adj(A)/det(A)".
for (int32 i = 0; i < COLS; i++) {
for (int32 j = 0; j < COLS; j++)
set(i, j, adj.get(i, j) * det);

}

return true;
}

// Gaussian elimination solves the inverse matrix：

bool invert2(){
assert(_columns == _rows);

int32 i, j, k, swap;
float32 t;
Matrix<COLS, COLS, T> temp;

for (i = 0; i < COLS; i++) {
for (j = 0; j < COLS; j++)
temp.set(j, i, get(j, i));

}

identity();

for (i = 0; i < COLS; i++) {
// look for largest element in column
swap = i;
for (j = i + 1; j < COLS; j++) {
if (abs(temp.get(i, j)) > abs(temp.get(i, i)))
swap = j;

}

if (swap != i) {
// swap rows.
for (k = 0; k < COLS; k++) {

t = temp.get(k, i);
temp.set(k, i, temp.get(k, swap));
temp.set(k, swap, t);

t = get(k, i);
set(k, i, get(k, swap));
set(k, swap, t);
}
}

if (temp.get(i, i) == 0)
// no non-zero pivot.
// the matrix is singular, which shouldn't
// happen.  This means the user gave us a bad matrix.
return false;

t = temp.get(i, i);
for (k = 0; k < COLS; k++) {
temp.set(k, i, temp.get(k, i) / t);
set(k, i, get(k, i) / t);
}

for (j = 0; j < COLS; j++) {
if (j != i) {
t = temp.get(i, j);
for (k = 0; k < COLS; k++) {
temp.set(k, j, temp.get(k, j) - temp.get(k, i) * t);
set(k, j, get(k, j) - get(k, i) * t);
}
}
}
}
return true;
}

private:
T determinantAid(int32 n) const {
switch (n)
{
case 1:
return get(0, 0);

case 2:
return get(0, 0) * get(1, 1) - get(1, 0) * get(0, 1);

default:
T det = 0;

Matrix<COLS, COLS, T> temp;

int32 sign = 1;

// iterate for each element of first row.
for (int32 i = 0; i < n; i++) {
cofactor(temp, 0, i, n);
det += sign * get(0, i) * temp.determinantAid(n - 1);

// terms are to be added with alternate sign.
sign = -sign;
}

return det;
}
}

void cofactor(Matrix<COLS, ROWS, T>& temp,
int32 p, int32 q, int32 n) const {

int32 i = 0, j = 0;

// looping for each element of the matrix.
for (int32 col = 0; col < n; col++) {

for (int32 row = 0; row < n; row++) {
// copying into temporary matrix only those element
// which are not in given row and column.
if (col != p && row != q) {

temp.set(j++, i, get(col, row));
// row is filled, so increase row index and
// reset col index.
if (j == n - 1) {
j = 0;
i++;
}
}
}
}
}
}

Next I created a Maxtrix4x4 class using Minors, Cofactors and Adjugate solutions：

class Matrix4x4
{
public:
float32 m11, m12, m13, m14,
m21, m22, m23, m24,
m31, m32, m33, m34,
m41, m42, m43, m44;

Matrix4x4(float32 m11 = 1.0f, float32 m12 = 0.0f, float32 m13 = 0.0f, float32 m14 = 0.0f,
float32 m21 = 0.0f, float32 m22 = 1.0f, float32 m23 = 0.0f, float32 m24 = 0.0f,
float32 m31 = 0.0f, float32 m32 = 0.0f, float32 m33 = 1.0f, float32 m34 = 0.0f,
float32 m41 = 0.0f, float32 m42 = 0.0f, float32 m43 = 0.0f, float32 m44 = 1.0f) {
.........
}

bool invert() {

float32 a = m33 * m44 - m43 * m34;
float32 b = m32 * m44 - m42 * m34;
float32 c = m32 * m43 - m42 * m33;
float32 d = m31 * m44 - m41 * m34;
float32 e = m31 * m43 - m41 * m33;
float32 f = m31 * m42 - m41 * m32;

float32 g = m22 * a - m23 * b + m24 * c;
float32 h = m21 * a - m23 * d + m24 * e;
float32 i = m21 * b - m22 * d + m24 * f;
float32 j = m21 * c - m22 * e + m23 * f;

// calculate the determinant.
float32 det = m11 * g - m12 * h + m13 * i - m14 * j;

// close to zero, can't invert.
if (det == 0) {
identity();
return false;
}

float32 m = m23 * m44 - m43 * m24;
float32 n = m22 * m44 - m42 * m24;
float32 o = m22 * m43 - m42 * m23;
float32 p = m21 * m44 - m41 * m24;
float32 q = m21 * m43 - m41 * m23;
float32 r = m21 * m42 - m41 * m22;

float32 s = m23 * m34 - m33 * m24;
float32 t = m22 * m34 - m32 * m24;
float32 u = m22 * m33 - m32 * m23;
float32 v = m21 * m34 - m31 * m24;
float32 w = m21 * m33 - m31 * m23;
float32 x = m21 * m32 - m31 * m22;

det = 1.0f / det;

float32 t11 = g * det;
float32 t12 = -h * det;
float32 t13 = i * det;
float32 t14 = -j * det;

float32 t21 = (-m12 * a + m13 * b - m14 * c) * det;
float32 t22 =  (m11 * a - m13 * d + m14 * e) * det;
float32 t23 = (-m11 * b + m12 * d - m14 * f) * det;
float32 t24 =  (m11 * c - m12 * e + m13 * f) * det;

float32 t31 =  (m12 * m - m13 * n + m14 * o) * det;
float32 t32 = (-m11 * m + m13 * p - m14 * q) * det;
float32 t33 =  (m11 * n - m12 * p + m14 * r) * det;
float32 t34 = (-m11 * o + m12 * q - m13 * r) * det;

float32 t41 = (-m12 * s + m13 * t - m14 * u) * det;
float32 t42 =  (m11 * s - m13 * v + m14 * w) * det;
float32 t43 = (-m11 * t + m12 * v - m14 * x) * det;
float32 t44 =  (m11 * u - m12 * w + m13 * x) * det;

m11 = t11;
m12 = t21;
m13 = t31;
m14 = t41;

m21 = t12;
m22 = t22;
m23 = t32;
m24 = t42;

m31 = t13;
m32 = t23;
m33 = t33;
m34 = t43;

m41 = t14;
m42 = t24;
m43 = t34;
m44 = t44;

return true;    }
}


Next, I tested the time it took to calculate：

#include <cstdio>
#include <iostream>
#include <windows.h>
#include <time.h>

float64 timeSpent = 0;
LARGE_INTEGER nFreq;
LARGE_INTEGER nBeginTime;
LARGE_INTEGER nEndTime;

Matrix4x4 m4;
m4.m11 = 2.1018f;   m4.m12 = -1.81754f; m4.m13 = 1.2541f;  m4.m14 = 2.442f;
m4.m21 = 0.54194f;  m4.m22 = 2.75391f;  m4.m23 = -0.1167f; m4.m24 = 0.0f;
m4.m31 = -5.81652f; m4.m32 = -7.9381f;  m4.m33 = 4.2816f;  m4.m34 = 23.33819f;
m4.m41 = 9.5076f;   m4.m42 = 10.9058f;  m4.m43 = 2.0f;     m4.m44 = 4.8239f;

Matrix<4, 4, float32> mat1;
Matrix<4, 4, float32> mat2;
float32* e1 = mat1.elements;
e1[0] = 2.1018f;   e1[1] = -1.81754f; e1[2] = 1.2541f;  e1[3] = 2.442f;
e1[4] = 0.54194f;  e1[5] = 2.75391f;  e1[6] = -0.1167f; e1[7] = 0.0f;
e1[8] = -5.81652f; e1[9] = -7.9381f;  e1[10] = 4.2816f; e1[11] = 23.33819f;
e1[12] = 9.5076f;  e1[13] = 10.9058f; e1[14] = 2;       e1[15] = 4.8239f;

mat2 = mat1;

QueryPerformanceFrequency(&nFreq);
QueryPerformanceCounter(&nBeginTime); // start timer

for (i = 0; i < 100000; i++)
mat1.invert();

QueryPerformanceCounter(&nEndTime); // end timer

printf("time1:%f\n", timeSpent);

QueryPerformanceCounter(&nBeginTime);

for (i = 0; i < 100000; i++)
mat2.invert2();

QueryPerformanceCounter(&nEndTime);

printf("time2:%f\n", timeSpent);

QueryPerformanceCounter(&nBeginTime);

for (i = 0; i < 100000; i++)
m4.invert();

QueryPerformanceCounter(&nEndTime);

printf("time3:%f\n", timeSpent);


Output：

time1:0.122209
time2:0.014264
time3:0.002685
........
........
........
time1:0.100992
time2:0.014209
time3:0.002736
time1:0.101950
time2:0.014248
time3:0.002731
.......
.......


Obvious difference!

What do you think?

Regarding possible performance differences: yes, they are possible. But I'm also pretty sure you are not measuring it correctly. You need to compile with optimizations turned on, that is e.g. for gcc / clang you'll need to at least add -O3 to your compiler flags, also note that figures you get will strongly depend on architecture, used compiler etc. You may want to inspect some benchmark frameworks for C++, e.g. Quick Bench.

That being said, there are some other things you can improve. Since this is tagged as C++, you should use standard library, e.g. std::array from header array as the storage for your matrix class, you should probably also get familiar with chrono if you insist on measuring time "by hand", particularly steady_clock might be of interest to you (this part is a guess, since you have not included actual implementation for the QueryPerformance... functions).
I'd also recommend making the Matrix class template more container like- define proper constructors, member access functions, iterators (possibly) inside the class. On the other hand, I'd remove the remaining stuff from the class and make those functions free and operating on the matrix class e.g. the invert declaration could look like this:

template<std::size_t Dim, typename T>
using SquareMatrix = Matrix<Dim, Dim, T>;

template<std::size_t Dim, typename T>
SquareMatrix<Dim, T> invert (SquareMatrix<Dim, T> const& mat);


Also note that there is a very small number of types on which the algorithms work, this could be manifested using some SFINAE (that is probably a bit more advanced topic, but there is plenty of resources online on how to use that technique).

The Matrix4x4 is not needed- this is a very bad implementation of Matrix constructor implementation, should be handled by proper Matrix constructor instead.

Use enum class instead of plain enum, it's type-safe and prohibits implicit conversions from underlying numeric type (this will also force you to write the algorithm choice part properly- the ifs around the underlying values are not pretty).

Lastly, force your compiler to do some work for you: add -Werror -Wall -Wpedantic to your compiler flags, so that your code won't compile if you produce warnings (-Wall is a must to avoid a ton of undefined / unspecified behavior). If you are using MVSC compiler then you'll need /W3 instead.