The Data Science Lab
Computing the Determinant of a Matrix Using Gaussian Elimination to Row Echelon Form with C#
Dr. James McCaffrey presents a complete end-to-end demonstration of computing the determinant of a matrix using the C# language. In machine learning scenarios, computing the determinant of a matrix is typically used during model training to determine if a matrix has an inverse or not.
One of the fundamental operations in machine learning is computing the inverse of a square matrix. But not all matrices have an inverse. The most common way to check if a matrix has an inverse or not is to compute a value called the determinant of the matrix. If the determinant of a matrix is 0, the matrix does not have an inverse. Any other value of the determinant means the source matrix does have an inverse.
There are several algorithms that can be used to compute the determinant of a matrix. One of the simplest techniques is called Gaussian elimination to row echelon form. This is quite a mouthful so the name of the technique is often shortened in various ways, such as to just "row echelon." This article presents a complete code implementation of computing the determinant of a matrix, using the row echelon form technique (aka Gaussian elimination), with the C# language.
Note: A closely related technique that can be used to compute the determinant of a matrix is called Gauss-Jordan elimination. The Gaussian elimination technique presented in this article converts a matrix to "row echelon" form. Gauss-Jordan elimination converts a matrix to "reduced row echelon" form.
[Click on image for larger view.] Figure 1: Computing the Determinant of a Matrix in Action
Gaussian elimination converts the source matrix to a version called row echelon form. The determinant is calculated as a side effect during the conversion process.
A good way to see where this article is headed is to take a look at the screenshot of a demo program in Figure 1. The demo begins by setting up a 2-by-2 matrix and computing its determinant:
Source matrix:
2.0 9.0
1.0 8.0
Determinant = 7.0000
Because the determinant is not 0, the source matrix has an inverse. As you'll see shortly, computing the determinant of a 2-by-2 matrix is very easy. Next, the demo sets up and computes the determinant of a 3-by-3 matrix:
Source matrix:
5.0 9.0 7.0
3.0 8.0 -1.0
-2.0 1.0 6.0
Determinant = 234.0000
Computing the determinant of a 3-by-3 matrix is also an easy problem. Next, the demo tackles a 4-by-4 matrix:
Source matrix:
1.0 -4.0 3.0 0.0
2.0 8.0 6.0 -1.0
9.0 7.0 5.0 -1.0
-1.0 0.0 2.0 3.0
Determinant = -1103.0000
Computing the determinant of a 4-by-4 or larger matrix is surprisingly difficult. The demo concludes by showing an example where the source matrix has a determinant of 0, and therefore the matrix does not have an inverse:
Source matrix:
1.0 4.0 3.0 0.0
2.0 8.0 6.0 0.0
9.0 7.0 5.0 1.0
-1.0 0.0 -2.0 -3.0
Determinant = 0.0000
This article assumes you have intermediate or better programming skill but doesn't assume you know anything about the Gaussian elimination to row echelon form technique. The demo is implemented using C#, but you should be able to refactor the demo code to another C-family language if you wish. Most normal error checking has been removed to keep the main ideas as clear as possible.
The source code for the demo program is a bit too long to be presented in its entirety in this article. The complete code is available in the accompanying file download, and is also available online.
Understanding Gaussian Elimination to Row Echelon Form
Maybe the best way to understand the Gaussian elimination technique is to look at a concrete example. Here is a 7-by-7 matrix that's in row echelon form:
3.8 3.7 8.2 7.6 1.5 2.2 8.3
0.0 0.0 5.4 3.8 1.8 0.0 4.9
0.0 0.0 0.0 8.4 3.3 6.3 2.6
0.0 0.0 0.0 0.0 0.0 5.7 0.5
0.0 0.0 0.0 0.0 0.0 2.3 8.1
0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
A square matrix is in row echelon form when a.) if there are any rows of all zeros, they are all on the bottom, b.) if you start at the upper left corner, and draw a border so that the 0.0 entries are below, then you get a staircase shape (that isn't necessarily symmetric). The word "echelon" can be traced back to the Latin word for "stairs."
The source matrix is converted to row echelon form using three operations:
- Any two rows can be swapped.
- Any row (or a multiple of it) can be added to another row.
- Any row can be multiplied by a constant.
If two rows are swapped, the determinant is multiplied by -1 (that is, changes sign). If you add a multiple of one row to another row, the determinant does not change. If you multiply a row by a constant k, the determinant is also multiplied by k.
Suppose one row of a 4-by-4 matrix is r1 = [3, 5, 6, 2] and another row is r2 = [1, 4, 7, 8]. If you add -3 times r2 to r1, r1 becomes [0, -7, -15, -22], which is one step closer to row echelon form.
The Demo Program
I used Visual Studio 2022 (Community Free Edition) for the demo program. I created a new C# console application and checked the "Place solution and project in the same directory" option. I specified .NET version 8.0. I named the project MatrixDeterminantEchelon. I checked the "Do not use top-level statements" option to avoid the weird program entry point shortcut syntax.
The demo has no significant .NET dependencies and any relatively recent version of Visual Studio with .NET (Core) or the older .NET Framework will work fine. You can also use the Visual Studio Code program if you like.
After the template code loaded into the editor, I right-clicked on file Program.cs in the Solution Explorer window and renamed the file to the more descriptive MatrixDeterminantEchelonProgram.cs. I allowed Visual Studio to automatically rename class Program.
The overall program structure is presented in Listing 1. All the control logic is in the Main() method in the Program class. The matrix determinant functionality is contained in a single static method named MatDeterminantEchelon(). The program also has a MatShow() method to display a matrix, and a MatLoad() method to load a matrix into memory from data stored in a text file.
Listing 1: Overall Program Structure
using System;
using System.IO;
using System.Collections.Generic;
namespace MatrixDeterminantEchelon
{
internal class MatrixDeterminantEchelonProgram
{
static void Main(string[] args)
{
Console.WriteLine("Begin matrix determinant using" +
" row echelon reduction technique ");
// first example (2-by-2)
// second example (3-by-3)
// third example (4-by-4) with non-zero determinant
// fourth example (4-by-4) with zero determinant
Console.WriteLine("End demo ");
Console.ReadLine();
} // Main()
static double MatDeterminantEchelon(double[][] A) { . . }
static void MatShow(double[][] A, int dec, int wid) { . . }
static double[][] MatLoad(string fn, int[] usecols,
char sep, string comment) { . . }
} // class Program
} // ns
The demo program begins with a simple 2-by-2 example:
double[][] A = new double[2][];
A[0] = new double[] { 2, 9 };
A[1] = new double[] { 1, 8 };
Console.WriteLine("Source matrix: ");
MatShow(A, 1, 5);
The demo creates an array-of-arrays style matrix. It is possible to define a Matrix class of some sort, but this creates unnecessary complexity in my opinion. The arguments to MatShow() mean display the matrix elements with 1 decimal, using 5 spaces.
The demo hard-codes the source matrix. In a non-demo scenario, you will often load values into a matrix from a text file. For example, suppose you have a data file like:
# the_data.txt
#
2.0, 9.0
1.0, 8.0
Then the statement double[][] A = MatLoad("the_data.txt", new int[] { 0, 1 }, ',', "#") means load from file "the_data.txt" the values in columns 0 and 1, where values are separated by the comma character, and lines beginning with "#" are comments to be ignored.
The determinant is computed and displayed using these statements:
double detA = MatDeterminantEchelon(A);
Console.WriteLine("Determinant = " +
detA.ToString("F4"));
The demo sets up three additional matrices and computes their determinants using the same pattern:
double[][] B = new double[3][];
B[0] = new double[] { 5, 9, 7 };
B[1] = new double[] { 3, 8, -1 };
B[2] = new double[] { -2, 1, 6 };
. . .
double detD = MatDeterminantEchelon(D);
Console.WriteLine("Determinant = " +
detD.ToString("F4"));
Console.WriteLine("End demo ");
The demo program presents a standalone function to compute a matrix determinant. You can either check the determinant value before calling a matrix inverse function, or you can integrate the determinant code inside a matrix inverse function. Design decisions like this can be surprisingly tricky and you should consider them carefully.
Computing the Determinant
The implementation of matrix determinant function is presented in Listing 2. The function begins by checking for simple cases of a 1-by-1 matrix with a single value, or a 2-by-2 matrix:
static double MatDeterminantEchelon(double[][] A)
{
int n = A.Length; // assume A is square
if (n == 1)
return A[0][0]; // usually
if (n == 2)
return (A[0][0] * A[1][1]) - (A[0][1] * A[1][0]);
. . .
If a 2-by-2 matrix has values
a b
c d
the determinant is (a * d) - (b * c). Very easy.
A simultaneous advantage and disadvantage of implementing matrix operations from scratch rather than using a code library is that you can omit error checking when you can be sure it won't be needed. Here, the demo code assumes the source matrix is square.
Another advantage of from-scratch code is that you can easily modify your code for special cases. The demo code assumes that a 1-by-1 matrix does have a determinant, but for some scenarios you might want to throw an exception.
Listing 2: Matrix Determinant Function Implementation
static double MatDeterminantEchelon(double[][] A)
{
// using Gaussian elimination to row echelon form
int n = A.Length; // assume A is square
// small n shortcuts
if (n == 1)
return A[0][0]; // usually
if (n == 2)
return (A[0][0] * A[1][1]) - (A[0][1] * A[1][0]);
if (n == 3)
{
double a = (A[1][1] * A[2][2]) - (A[1][2] * A[2][1]);
double b = (A[1][0] * A[2][2]) - (A[1][2] * A[2][0]);
double c = (A[1][0] * A[2][1]) - (A[1][1] * A[2][0]);
double result =
(A[0][0] * a) - (A[0][1] * b) + (A[0][2] * c);
return result;
}
// row echelon for n = 4 or larger
double det = 1.0;
// make a copy of A
double[][] B = new double[n][];
for (int i = 0; i < n; ++i)
B[i] = new double[n];
for (int i = 0; i < n; ++i)
for (int j = 0; j < n; ++j)
B[i][j] = A[i][j];
for (int i = 0; i < n; ++i) // main loop
{
int pr = i; // find pivot row
for (int r = i + 1; r < n; ++r)
{
if (Math.Abs(B[r][i]) > Math.Abs(B[pr][i]))
pr = r;
}
if (pr != i) // swap rows pr and i
{
for (int k = 0; k < n; ++k)
{
double tmp = B[i][k];
B[i][k] = B[pr][k];
B[pr][k] = tmp;
}
det *= -1.0; // swap flips sign
}
if (Math.Abs(B[i][i]) < 1.0e-8)
return 0.0;
det *= B[i][i]; // theory
for (int r = i + 1; r < n; ++r) // elim below pivot
{
double scale = B[r][i] / B[i][i]; // not 0
for (int k = i; k < n; ++k)
B[r][k] -= scale * B[i][k];
}
} // main loop
// MatShow(B, 4, 9); // B is in row echelon form
return det;
}
The implementation continues:
if (n == 3)
{
double a = (A[1][1] * A[2][2]) - (A[1][2] * A[2][1]);
double b = (A[1][0] * A[2][2]) - (A[1][2] * A[2][0]);
double c = (A[1][0] * A[2][1]) - (A[1][1] * A[2][0]);
double result =
(A[0][0] * a) - (A[0][1] * b) + (A[0][2] * c);
return result;
}
. . .
This code computes the determinant directly, using cofactors along the top row, instead of Gaussian elimination. If a 3-by-3 matrix is:
a b c
d e f
g h i
then the determinant is
1 * a * ((e * i) - (f * h)) +
-1 * b * ((d - i) * (f - g)) +
1 * c * ((d - h) * (e - g))
So, the determinant of a 3-by-3 matrix can be computed easily and directly.
At this point, the implementation is dealing with a matrix that is 4-by-4 or larger, and begins Gaussian elimination to row echelon form:
double det = 1.0;
double[][] B = new double[n][];
for (int i = 0; i < n; ++i)
B[i] = new double[n];
for (int i = 0; i < n; ++i)
for (int j = 0; j < n; ++j)
B[i][j] = A[i][j];
. . .
The result determinant is initially set to 1.0 and will be modified as the source matrix is converted to row echelon form using the Gaussian elimination technique.
The code creates a copy of the source matrix. Because an array-of-arrays style matrix is a reference object, working directly on it when passed as a parameter will change the matrix. This might be OK in some situations, but it's usually a better idea to make a copy of the source matrix. Next:
for (int i = 0; i < n; ++i) // main loop
{
int pr = i; // find pivot row
for (int r = i + 1; r < n; ++r) {
if (Math.Abs(B[r][i]) > Math.Abs(B[pr][i]))
pr = r;
}
. . .
Gaussian elimination uses a main loop that iterates over each row. The pivot row in each iteration is the row that starts with the largest magnitude value. The next few statements perform a row swap if necessary:
if (pr != i) { // swap rows pr and i
for (int k = 0; k < n; ++k) {
double tmp = B[i][k];
B[i][k] = B[pr][k];
B[pr][k] = tmp;
}
det *= -1.0; // swap flips sign
}
. . .
Recall that if a row swap occurs, the sign of the determinant flips. Now the code checks to see if the current diagonal element of the row echelon form is 0:
if (Math.Abs(B[i][i]) < 1.0e-8)
return 0.0;
det *= B[i][i];
. . .
When implementing any matrix code, an important issue is defining what value is zero. The demo hard-codes a value of 0.00000001 so that anything smaller in magnitude is defined to be zero. In many matrix code libraries, there is a global value for zero. In some code libraries, the zero value is passed as a parameter to all functions, often with a default value.
The demo code is ready to multiply rows below the current row i by the constant in B[i][i]. Recall that this will multiply the determinant by that constant. The code performs the multiplication at this point, but it could be done after the row multiplications.
The determinant function concludes with:
. . .
for (int r = i + 1; r < n; ++r) { // elim below pivot
double scale = B[r][i] / B[i][i];
for (int k = i; k < n; ++k)
B[r][k] -= scale * B[i][k];
}
} // main loop
return det;
}
Notice there is a division by the B[i][i] value. This cannot be zero because the short circuit return earlier would have caught it.
Wrapping Up
In machine learning, there are several algorithms that require computing the inverse of a matrix. One common example is kernel ridge regression using kernel matrix inverse training. In such situations, you usually want to check if the matrix involved has an inverse by computing the determinant of the matrix. If the determinant is 0, then the matrix does not have an inverse and you must use an alternative training technique, typically stochastic gradient descent.