Many classical machine learning algorithms use matrix inversion, but inverting a matrix is a very difficult problem. I have a personal C# library of matrix functions, including a MatInverse() function.
My MatInverse() function uses Crout’s decomposition, which works for any size square matrix. However, decomposition is a bit of overkill if the source matrix is 2×2 or 3×3. For those special cases, it’s easier to implement a direct solution. If the source matrix is 4×4, a direct solution is possible but a bit ugly. See the demo code below.
The 2×2 case is:
public static double[][] MatInverse2x2(double[][] m)
{
double[][] result = MatCreate(2, 2);
double det = (m[0][0] * m[1][1]) -
(m[0][1] * m[1][0]); // determinant
result[0][0] = m[1][1] / det;
result[0][1] = -1 * m[0][1] / det;
result[1][0] = -1 * m[1][0] / det;
result[1][1] = m[0][0] / det;
return result;
}
where
public static double[][] MatCreate(int rows, int cols)
{
double[][] result = new double[rows][];
for (int i = 0; i "lt"; rows; ++i) // "less-than"
result[i] = new double[cols];
return result;
}
I rarely include the special case matrix inverse functions in my programs because for the problems I work with, 2×2 and 3×3 matrices are somewhat rare. The 4×4 case is more common but direct computation of the inverse has a lot of code.
In aircraft models, the most common small scales I built when I was a young man were 1/72, 1/48, 1/32. Here are three examples of beautiful WWII airplane models I found on the Internet. Left: A Bell P-39 Airacobra in 1/72 scale. I loved the mid-engine design and tricycle landing gear. Center: A North American P-51 Mustang in 1/32 scale. I loved the ventral air scoop and raised bubble canopy. Right: A Vought F4U Corsair in 1/32 scale. I loved the inverted gull wing and short rear fuselage.
Demo program. Replace “lt”, “gt”, “lte”, “gte” with Boolean operator symbols.
namespace SmallMatrixInverse
{
internal class SmallMatrixProgram
{
static void Main(string[] args)
{
Console.WriteLine("\nSmall matrix inverse C# demo ");
double[][] A = Utils.MatRandom(2, 2, -1.0, 1.0, 0);
Console.WriteLine("\nSource matrix: ");
Utils.MatShow(A, 4, 8);
double[][] Ainv1 = Utils.MatInverse(A);
Console.WriteLine("\nInverted using decomposition: ");
Utils.MatShow(Ainv1, 8, 12);
double[][] Ainv2 = Utils.MatInverse2x2(A);
Console.WriteLine("\nInverted using 2x2 specific: ");
Utils.MatShow(Ainv2, 8, 12);
Console.WriteLine("---------------------");
double[][] B = Utils.MatRandom(3, 3, -1.0, 1.0, 1);
Console.WriteLine("\nSource matrix: ");
Utils.MatShow(B, 4, 8);
double[][] Binv1 = Utils.MatInverse(B);
Console.WriteLine("\nInverted using decomposition: ");
Utils.MatShow(Binv1, 8, 12);
double[][] Binv2 = Utils.MatInverse3x3(B);
Console.WriteLine("\nInverted using 3x3 specific: ");
Utils.MatShow(Binv2, 8, 12);
Console.WriteLine("---------------------");
double[][] C = Utils.MatRandom(4, 4, -1.0, 1.0, 2);
Console.WriteLine("\nSource matrix: ");
Utils.MatShow(C, 4, 8);
double[][] Cinv1 = Utils.MatInverse(C);
Console.WriteLine("\nInverted using decomposition: ");
Utils.MatShow(Cinv1, 8, 12);
double[][] Cinv2 = Utils.MatInverse4x4(C);
Console.WriteLine("\nInverted using 4x4 specific: ");
Utils.MatShow(Cinv2, 8, 12);
Console.WriteLine("\nEnd demo ");
Console.ReadLine();
} // Main
} // Program
public class Utils
{
public static double[][] MatCreate(int rows, int cols)
{
double[][] result = new double[rows][];
for (int i = 0; i "lt" rows; ++i)
result[i] = new double[cols];
return result;
}
public static double[][] MatRandom(int rows, int cols,
double lo, double hi, int seed)
{
Random rnd = new Random(seed);
double[][] result = MatCreate(rows, cols);
for (int i = 0; i "lt" rows; ++i)
for (int j = 0; j "lt" cols; ++j)
result[i][j] = (hi - lo) * rnd.NextDouble() + lo;
return result;
}
public static double[][] MatInverse(double[][] m)
{
// assumes determinant is not 0
// that is, the matrix does have an inverse
int n = m.Length;
double[][] result = MatCreate(n, n); // make a copy
for (int i = 0; i "lt" n; ++i)
for (int j = 0; j "lt" n; ++j)
result[i][j] = m[i][j];
double[][] lum; // combined lower & upper
int[] perm; // out parameter
MatDecompose(m, out lum, out perm); // ignore return
double[] b = new double[n];
for (int i = 0; i "lt" n; ++i)
{
for (int j = 0; j "lt" n; ++j)
if (i == perm[j])
b[j] = 1.0;
else
b[j] = 0.0;
double[] x = Reduce(lum, b); //
for (int j = 0; j "lt" n; ++j)
result[j][i] = x[j];
}
return result;
}
static int MatDecompose(double[][] m,
out double[][] lum, out int[] perm)
{
// Crout's LU decomposition for matrix determinant
// and inverse.
// stores combined lower & upper in lum[][]
// stores row permuations into perm[]
// returns +1 or -1 according to even or odd number
// of row permutations.
// lower gets dummy 1.0s on diagonal (0.0s above)
// upper gets lum values on diagonal (0.0s below)
// even (+1) or odd (-1) row permutatuions
int toggle = +1;
int n = m.Length;
// make a copy of m[][] into result lu[][]
lum = MatCreate(n, n);
for (int i = 0; i "lt" n; ++i)
for (int j = 0; j "lt" n; ++j)
lum[i][j] = m[i][j];
// make perm[]
perm = new int[n];
for (int i = 0; i "lt" n; ++i)
perm[i] = i;
for (int j = 0; j "lt" n - 1; ++j) // note n-1
{
double max = Math.Abs(lum[j][j]);
int piv = j;
for (int i = j + 1; i "lt" n; ++i) // pivot index
{
double xij = Math.Abs(lum[i][j]);
if (xij "gt" max)
{
max = xij;
piv = i;
}
} // i
if (piv != j)
{
double[] tmp = lum[piv]; // swap rows j, piv
lum[piv] = lum[j];
lum[j] = tmp;
int t = perm[piv]; // swap perm elements
perm[piv] = perm[j];
perm[j] = t;
toggle = -toggle;
}
double xjj = lum[j][j];
if (xjj != 0.0)
{
for (int i = j + 1; i "lt" n; ++i)
{
double xij = lum[i][j] / xjj;
lum[i][j] = xij;
for (int k = j + 1; k "lt" n; ++k)
lum[i][k] -= xij * lum[j][k];
}
}
} // j
return toggle; // for determinant
} // MatDecompose
static double[] Reduce(double[][] luMatrix,
double[] b) // helper
{
int n = luMatrix.Length;
double[] x = new double[n];
//b.CopyTo(x, 0);
for (int i = 0; i "lt" n; ++i)
x[i] = b[i];
for (int i = 1; i "lt" n; ++i)
{
double sum = x[i];
for (int j = 0; j "lt" i; ++j)
sum -= luMatrix[i][j] * x[j];
x[i] = sum;
}
x[n - 1] /= luMatrix[n - 1][n - 1];
for (int i = n - 2; i "gte" 0; --i)
{
double sum = x[i];
for (int j = i + 1; j "lt" n; ++j)
sum -= luMatrix[i][j] * x[j];
x[i] = sum / luMatrix[i][i];
}
return x;
} // Reduce
//public static double MatDeterminant(double[][] m)
//{
// double[][] lum;
// int[] perm;
// double result = MatDecompose(m, out lum, out perm);
// for (int i = 0; i "lt" lum.Length; ++i)
// result *= lum[i][i];
// return result;
//}
// -------
public static double[][] MatInverse2x2(double[][] m)
{
double[][] result = MatCreate(2, 2);
double det = (m[0][0] * m[1][1]) -
(m[0][1] * m[1][0]); // determinant
result[0][0] = m[1][1] / det;
result[0][1] = -1 * m[0][1] / det;
result[1][0] = -1 * m[1][0] / det;
result[1][1] = m[0][0] / det;
return result;
}
public static double[][] MatInverse3x3(double[][] m)
{
// assumes determinant is not 0
// that is, the matrix does have an inverse
double[][] result = MatCreate(3, 3);
double det = (m[0][0] * m[1][1] * m[2][2]) +
(m[0][1] * m[1][2] * m[2][0]) +
(m[0][2] * m[1][0] * m[2][1]) -
(m[0][2] * m[1][1] * m[2][0]) -
(m[0][1] * m[1][0] * m[2][2]) -
(m[0][0] * m[1][2] * m[2][1]); // determinant
result[0][0] = 1 * (m[1][1] * m[2][2] - m[1][2] * m[2][1]);
result[0][1] = -1 * (m[0][1] * m[2][2] - m[0][2] * m[2][1]);
result[0][2] = 1 * (m[0][1] * m[1][2] - m[0][2] * m[1][1]);
result[1][0] = -1 * (m[1][0] * m[2][2] - m[1][2] * m[2][0]);
result[1][1] = 1 * (m[0][0] * m[2][2] - m[0][2] * m[2][0]);
result[1][2] = -1 * (m[0][0] * m[1][2] - m[0][2] * m[1][0]);
result[2][0] = 1 * (m[1][0] * m[2][1] - m[1][1] * m[2][0]);
result[2][1] = -1 * (m[0][0] * m[2][1] - m[0][1] * m[2][0]);
result[2][2] = 1 * (m[0][0] * m[1][1] - m[0][1] * m[1][0]);
for (int i = 0; i "lt" 3; ++i)
for (int j = 0; j "lt" 3; ++j)
result[i][j] /= det;
return result;
} // 3x3
public static double[][] MatInverse4x4(double[][] m)
{
// assumes determinant is not 0
// that is, the matrix does have an inverse
double[][] result = MatCreate(4, 4);
// double det =
//m[0][0] * m[1][1] * m[2][2] * m[3][3] +
//m[0][0] * m[1][2] * m[2][3] * m[3][1] +
//m[0][0] * m[1][3] * m[2][1] * m[3][2] -
//m[0][0] * m[1][3] * m[2][2] * m[3][1] -
//m[0][0] * m[1][2] * m[2][1] * m[3][3] -
//m[0][0] * m[1][1] * m[2][3] * m[3][2] -
//m[0][1] * m[1][0] * m[2][2] * m[3][3] -
//m[0][2] * m[1][0] * m[2][3] * m[3][1] -
//m[0][3] * m[1][0] * m[2][1] * m[3][2] +
//m[0][3] * m[1][0] * m[2][2] * m[3][1] +
//m[0][2] * m[1][0] * m[2][1] * m[3][3] +
//m[0][1] * m[1][0] * m[2][3] * m[3][2] +
//m[0][1] * m[1][2] * m[2][0] * m[3][3] +
//m[0][2] * m[1][3] * m[2][0] * m[3][1] +
//m[0][3] * m[1][1] * m[2][0] * m[3][2] -
//m[0][3] * m[1][2] * m[2][0] * m[3][1] -
//m[0][2] * m[1][1] * m[2][0] * m[3][3] -
//m[0][1] * m[1][3] * m[2][0] * m[3][2] -
//m[0][1] * m[1][2] * m[2][3] * m[3][0] -
//m[0][2] * m[1][3] * m[2][1] * m[3][0] -
//m[0][3] * m[1][1] * m[2][2] * m[3][0] +
//m[0][3] * m[1][2] * m[2][1] * m[3][0] +
//m[0][2] * m[1][1] * m[2][3] * m[3][0] +
//m[0][1] * m[1][3] * m[2][2] * m[3][0];
// slightly more efficient version
double det =
m[0][0] *
(m[1][1] * m[2][2] * m[3][3] +
m[1][2] * m[2][3] * m[3][1] +
m[1][3] * m[2][1] * m[3][2] -
m[1][3] * m[2][2] * m[3][1] -
m[1][2] * m[2][1] * m[3][3] -
m[1][1] * m[2][3] * m[3][2])
- m[1][0] *
(m[0][1] * m[2][2] * m[3][3] +
m[0][2] * m[2][3] * m[3][1] +
m[0][3] * m[2][1] * m[3][2] -
m[0][3] * m[2][2] * m[3][1] -
m[0][2] * m[2][1] * m[3][3] -
m[0][1] * m[2][3] * m[3][2])
+ m[2][0] *
(m[0][1] * m[1][2] * m[3][3] +
m[0][2] * m[1][3] * m[3][1] +
m[0][3] * m[1][1] * m[3][2] -
m[0][3] * m[1][2] * m[3][1] -
m[0][2] * m[1][1] * m[3][3] -
m[0][1] * m[1][3] * m[3][2])
- m[3][0] *
(m[0][1] * m[1][2] * m[2][3] +
m[0][2] * m[1][3] * m[2][1] +
m[0][3] * m[1][1] * m[2][2] -
m[0][3] * m[1][2] * m[2][1] -
m[0][2] * m[1][1] * m[2][3] -
m[0][1] * m[1][3] * m[2][2]);
result[0][0] = m[1][1] * m[2][2] * m[3][3] +
m[1][2] * m[2][3] * m[3][1] +
m[1][3] * m[2][1] * m[3][2] -
m[1][3] * m[2][2] * m[3][1] -
m[1][2] * m[2][1] * m[3][3] -
m[1][1] * m[2][3] * m[3][2];
result[0][1] = -m[0][1] * m[2][2] * m[3][3] -
m[0][2] * m[2][3] * m[3][1] -
m[0][3] * m[2][1] * m[3][2] +
m[0][3] * m[2][2] * m[3][1] +
m[0][2] * m[2][1] * m[3][3] +
m[0][1] * m[2][3] * m[3][2];
result[0][2] = m[0][1] * m[1][2] * m[3][3] +
m[0][2] * m[1][3] * m[3][1] +
m[0][3] * m[1][1] * m[3][2] -
m[0][3] * m[1][2] * m[3][1] -
m[0][2] * m[1][1] * m[3][3] -
m[0][1] * m[1][3] * m[3][2];
result[0][3] = -m[0][1] * m[1][2] * m[2][3] -
m[0][2] * m[1][3] * m[2][1] -
m[0][3] * m[1][1] * m[2][2] +
m[0][3] * m[1][2] * m[2][1] +
m[0][2] * m[1][1] * m[2][3] +
m[0][1] * m[1][3] * m[2][2];
result[1][0] = -m[1][0] * m[2][2] * m[3][3] -
m[1][2] * m[2][3] * m[3][0] -
m[1][3] * m[2][0] * m[3][2] +
m[1][3] * m[2][2] * m[3][0] +
m[1][2] * m[2][0] * m[3][3] +
m[1][0] * m[2][3] * m[3][2];
result[1][1] = m[0][0] * m[2][2] * m[3][3] +
m[0][2] * m[2][3] * m[3][0] +
m[0][3] * m[2][0] * m[3][2] -
m[0][3] * m[2][2] * m[3][0] -
m[0][2] * m[2][0] * m[3][3] -
m[0][0] * m[2][3] * m[3][2];
result[1][2] = -m[0][0] * m[1][2] * m[3][3] -
m[0][2] * m[1][3] * m[3][0] -
m[0][3] * m[1][0] * m[3][2] +
m[0][3] * m[1][2] * m[3][0] +
m[0][2] * m[1][0] * m[3][3] +
m[0][0] * m[1][3] * m[3][2];
result[1][3] = m[0][0] * m[1][2] * m[2][3] +
m[0][2] * m[1][3] * m[2][0] +
m[0][3] * m[1][0] * m[2][2] -
m[0][3] * m[1][2] * m[2][0] -
m[0][2] * m[1][0] * m[2][3] -
m[0][0] * m[1][3] * m[2][2];
result[2][0] = m[1][0] * m[2][1] * m[3][3] +
m[1][1] * m[2][3] * m[3][0] +
m[1][3] * m[2][0] * m[3][1] -
m[1][3] * m[2][1] * m[3][0] -
m[1][1] * m[2][0] * m[3][3] -
m[1][0] * m[2][3] * m[3][1];
result[2][1] = -m[0][0] * m[2][1] * m[3][3] -
m[0][1] * m[2][3] * m[3][0] -
m[0][3] * m[2][0] * m[3][1] +
m[0][3] * m[2][1] * m[3][0] +
m[0][1] * m[2][0] * m[3][3] +
m[0][0] * m[2][3] * m[3][1];
result[2][2] = m[0][0] * m[1][1] * m[3][3] +
m[0][1] * m[1][3] * m[3][0] +
m[0][3] * m[1][0] * m[3][1] -
m[0][3] * m[1][1] * m[3][0] -
m[0][1] * m[1][0] * m[3][3] -
m[0][0] * m[1][3] * m[3][1];
result[2][3] = -m[0][0] * m[1][1] * m[2][3] -
m[0][1] * m[1][3] * m[2][0] -
m[0][3] * m[1][0] * m[2][1] +
m[0][3] * m[1][1] * m[2][0] +
m[0][1] * m[1][0] * m[2][3] +
m[0][0] * m[1][3] * m[2][1];
result[3][0] = -m[1][0] * m[2][1] * m[3][2] -
m[1][1] * m[2][2] * m[3][0] -
m[1][2] * m[2][0] * m[3][1] +
m[1][2] * m[2][1] * m[3][0] +
m[1][1] * m[2][0] * m[3][2] +
m[1][0] * m[2][2] * m[3][1];
result[3][1] = m[0][0] * m[2][1] * m[3][2] +
m[0][1] * m[2][2] * m[3][0] +
m[0][2] * m[2][0] * m[3][1] -
m[0][2] * m[2][1] * m[3][0] -
m[0][1] * m[2][0] * m[3][2] -
m[0][0] * m[2][2] * m[3][1];
result[3][2] = -m[0][0] * m[1][1] * m[3][2] -
m[0][1] * m[1][2] * m[3][0] -
m[0][2] * m[1][0] * m[3][1] +
m[0][2] * m[1][1] * m[3][0] +
m[0][1] * m[1][0] * m[3][2] +
m[0][0] * m[1][2] * m[3][1];
result[3][3] = m[0][0] * m[1][1] * m[2][2] +
m[0][1] * m[1][2] * m[2][0] +
m[0][2] * m[1][0] * m[2][1] -
m[0][2] * m[1][1] * m[2][0] -
m[0][1] * m[1][0] * m[2][2] -
m[0][0] * m[1][2] * m[2][1];
for (int i = 0; i "lt" 4; ++i)
for (int j = 0; j "lt" 4; ++j)
result[i][j] /= det;
return result;
} // 4x4
public static void MatShow(double[][] m,
int dec, int wid)
{
for (int i = 0; i "lt" m.Length; ++i)
{
for (int j = 0; j "lt" m[0].Length; ++j)
{
double v = m[i][j];
if (Math.Abs(v) "lt" 1.0e-8) v = 0.0; // hack
Console.Write(v.ToString("F" +
dec).PadLeft(wid));
}
Console.WriteLine("");
}
}
} // class Utils
} // ns
.NET Test Automation Recipes
Software Testing
SciPy Programming Succinctly
Keras Succinctly
R Programming
2026 Visual Studio Live
2025 Summer MLADS Conference
2026 DevIntersection Conference
2025 Machine Learning Week
2025 Ai4 Conference
2026 G2E Conference
2026 iSC West Conference
You must be logged in to post a comment.