Inverse of 6X6 Matrix or NXN Matrix

Question

• Hi all,

I'm using Visual studio 2005 and C# as back end for my project.I have a 36 double values i need to store them as a 6X6 matrix and get inverse of that matrix for output of my project.I just came across visual studio don't support matrices beyond 3X3 size.

I'd really appreciate if someone can help me out by sending me some sample code to inverse a 6X6 or NXN matrix.I tried searching online all the assembly(.dll files) i found online giving me inaccurate values.Thanks in advance,awaiting your reply.

Thanks,
Sri.

Monday, April 14, 2008 1:41 PM

• Here is a rather basic method using Gauss-Jordan elimination (the same way you learn to calculate inverses by hand in an entry level Linear Algebra course).  Although it's not super optimal, for small enough matricies this should work just fine.  This method will perform rather poorly on large matricies in which case high performance numerical libraries should be used (they'll tend to use LU factorization which is much more efficient but also much more complicated)  but for smaller matricies this method should work quite nicely:

(Note that the helper methods take in a 'cols' parameter.  While this is not necessary as it can be retrieved from the array.GetLength method, by taking it as a parameter the call to GetLength only needs to happen once.  Also since these are private methods and not callable by outside code, there is no input validation)

Code Snippet

/// <summary>
/// Calculate the inverse of a matrix using Gauss-Jordan elimination.
/// </summary>
/// <param name="data">The matrix to invert.</param>
/// <param name="inverse">The inverse of the matrix.</param>
/// <returns>false of the matrix is singular, true otherwise.</returns>
public static bool Invert(double[,] data, double[,] inverse)
{
if (data == null)
throw new ArgumentNullException("data");
if (inverse == null)
throw new ArgumentNullException("null");

int drows = data.GetLength(0),
dcols = data.GetLength(1),
irows = inverse.GetLength(0),
icols = inverse.GetLength(1),
n, r;
double scale;

//Validate the matrix sizes
if (drows != dcols)
throw new ArgumentException(

"data is not a square matrix", "data");
if (irows != icols)
throw new ArgumentException(

"inverse is not a square matrix", "inverse");
if (drows != irows)
throw new ArgumentException(

"data and inverse must have identical sizes.");

n = drows;

//Initialize the inverse to the identity
for (r = 0; r < n; ++r)
for (int c = 0; c < n; ++c)
inverse[r, c] = (r == c) ? 1.0 : 0.0;

//Process the matrix one column at a time
for (int c = 0; c < n; ++c)
{

//Swap rows if the current value is too close to 0.0

if (Math.Abs(data[c, c]) <= 2.0 * double.Epsilon)
{
for (r = c + 1; r < n; ++r)
if (Math.Abs(data[r, c]) <= 2.0 * double.Epsilon)
{
RowSwap(data, n, c, r);
RowSwap(inverse, n, c, r);
break;
}
if (r >= n)
return false;
}
scale = 1.0 / data[c, c];
RowScale(data, n, scale, c);
RowScale(inverse, n, scale, c);

//Zero out the rest of the column

for (r = 0; r < n; ++r)
{
if (r != c)
{
scale = -data[r, c];
}
}
}

return true;
}

/// <summary>
/// Swap 2 rows in a matrix.
/// </summary>
/// <param name="data">The matrix to operate on.</param>
/// <param name="cols">
/// The number of columns in <paramref name="data"/>.
/// </param>
/// <param name="r0">One of the 2 rows to swap.</param>
/// <param name="r1">One of the 2 rows to swap.</param>
private static void RowSwap(double[,] data, int cols,
int r0, int r1)
{
double tmp;

for (int i = 0; i < cols; ++i)
{
tmp = data[r0, i];
data[r0, i] = data[r1, i];
data[r1, i] = tmp;
}
}

/// <summary>
/// Perform scale and add a row in a matrix to another
/// row:  data[r1,] = a*data[r0,] + data[r1,].
/// </summary>
/// <param name="data">The matrix to operate on.</param>
/// <param name="cols">
/// The number of columns in <paramref name="data"/>.
/// </param>
/// <param name="a">
/// The scale factor to apply to row <paramref name="r0"/>.
/// </param>
/// <param name="r0">The row to scale.</param>
/// <param name="r1">The row to add and store to.</param>
private static void RowScaleAdd(double[,] data, int cols,
double a, int r0, int r1)
{
for (int i = 0; i < cols; ++i)
data[r1, i] += a * data[r0, i];
}

/// <summary>
/// Scale a row in a matrix by a constant factor.
/// </summary>
/// <param name="data">The matrix to operate on.</param>
/// <param name="cols">The number of columns in the matrix.</param>
/// <param name="a">
/// The factor to scale row <paramref name="r"/> by.
/// </param>
/// <param name="r">The row to scale.</param>
private static void RowScale(double[,] data, int cols,
double a, int r)
{
for (int i = 0; i < cols; ++i)
data[r, i] *= a;
}

On a side note, is calculating the inverse actually necessary?  If it's simply a middle step used for solving matrix equations, there are far more efficient and numerically stable ways of doing this that don't involve calculating the inverse.

Monday, April 14, 2008 7:29 PM

All replies

•  DesiRocks wrote:
 Hi all,I'm using Visual studio 2005 and C# as back end for my project.I have a 36 double values i need to store them as a 6X6 matrix and get inverse of that matrix for output of my project.I just came across visual studio don't support matrices beyond 3X3 size.I'd really appreciate if someone can help me out by sending me some sample code to inverse a 6X6 or NXN matrix.I tried searching online all the assembly(.dll files) i found online giving me inaccurate values.Thanks in advance,awaiting your reply.Thanks,Sri.

The Matrix class is for 3-D rotations only. If you want to operate with higher ranks you have to either code your own routine or purchase a package. Go to wikipedia.com or www.mathworld.wolfram.com and get an expression for inverse matrices of arbitrary dimensions. You will have to compute determinants first.

There are math packages online, quite a few, perhaps three or four, you will have to search. they are for numerical solutions of math problems. Obviously they contain much more than you need here. It is my recollection that they run around 400 bucks or so, maybe more.

Inverting a matrix is a simple undertaking code wise.

Monday, April 14, 2008 2:01 PM
• Here is a rather basic method using Gauss-Jordan elimination (the same way you learn to calculate inverses by hand in an entry level Linear Algebra course).  Although it's not super optimal, for small enough matricies this should work just fine.  This method will perform rather poorly on large matricies in which case high performance numerical libraries should be used (they'll tend to use LU factorization which is much more efficient but also much more complicated)  but for smaller matricies this method should work quite nicely:

(Note that the helper methods take in a 'cols' parameter.  While this is not necessary as it can be retrieved from the array.GetLength method, by taking it as a parameter the call to GetLength only needs to happen once.  Also since these are private methods and not callable by outside code, there is no input validation)

Code Snippet

/// <summary>
/// Calculate the inverse of a matrix using Gauss-Jordan elimination.
/// </summary>
/// <param name="data">The matrix to invert.</param>
/// <param name="inverse">The inverse of the matrix.</param>
/// <returns>false of the matrix is singular, true otherwise.</returns>
public static bool Invert(double[,] data, double[,] inverse)
{
if (data == null)
throw new ArgumentNullException("data");
if (inverse == null)
throw new ArgumentNullException("null");

int drows = data.GetLength(0),
dcols = data.GetLength(1),
irows = inverse.GetLength(0),
icols = inverse.GetLength(1),
n, r;
double scale;

//Validate the matrix sizes
if (drows != dcols)
throw new ArgumentException(

"data is not a square matrix", "data");
if (irows != icols)
throw new ArgumentException(

"inverse is not a square matrix", "inverse");
if (drows != irows)
throw new ArgumentException(

"data and inverse must have identical sizes.");

n = drows;

//Initialize the inverse to the identity
for (r = 0; r < n; ++r)
for (int c = 0; c < n; ++c)
inverse[r, c] = (r == c) ? 1.0 : 0.0;

//Process the matrix one column at a time
for (int c = 0; c < n; ++c)
{

//Swap rows if the current value is too close to 0.0

if (Math.Abs(data[c, c]) <= 2.0 * double.Epsilon)
{
for (r = c + 1; r < n; ++r)
if (Math.Abs(data[r, c]) <= 2.0 * double.Epsilon)
{
RowSwap(data, n, c, r);
RowSwap(inverse, n, c, r);
break;
}
if (r >= n)
return false;
}
scale = 1.0 / data[c, c];
RowScale(data, n, scale, c);
RowScale(inverse, n, scale, c);

//Zero out the rest of the column

for (r = 0; r < n; ++r)
{
if (r != c)
{
scale = -data[r, c];
}
}
}

return true;
}

/// <summary>
/// Swap 2 rows in a matrix.
/// </summary>
/// <param name="data">The matrix to operate on.</param>
/// <param name="cols">
/// The number of columns in <paramref name="data"/>.
/// </param>
/// <param name="r0">One of the 2 rows to swap.</param>
/// <param name="r1">One of the 2 rows to swap.</param>
private static void RowSwap(double[,] data, int cols,
int r0, int r1)
{
double tmp;

for (int i = 0; i < cols; ++i)
{
tmp = data[r0, i];
data[r0, i] = data[r1, i];
data[r1, i] = tmp;
}
}

/// <summary>
/// Perform scale and add a row in a matrix to another
/// row:  data[r1,] = a*data[r0,] + data[r1,].
/// </summary>
/// <param name="data">The matrix to operate on.</param>
/// <param name="cols">
/// The number of columns in <paramref name="data"/>.
/// </param>
/// <param name="a">
/// The scale factor to apply to row <paramref name="r0"/>.
/// </param>
/// <param name="r0">The row to scale.</param>
/// <param name="r1">The row to add and store to.</param>
private static void RowScaleAdd(double[,] data, int cols,
double a, int r0, int r1)
{
for (int i = 0; i < cols; ++i)
data[r1, i] += a * data[r0, i];
}

/// <summary>
/// Scale a row in a matrix by a constant factor.
/// </summary>
/// <param name="data">The matrix to operate on.</param>
/// <param name="cols">The number of columns in the matrix.</param>
/// <param name="a">
/// The factor to scale row <paramref name="r"/> by.
/// </param>
/// <param name="r">The row to scale.</param>
private static void RowScale(double[,] data, int cols,
double a, int r)
{
for (int i = 0; i < cols; ++i)
data[r, i] *= a;
}

On a side note, is calculating the inverse actually necessary?  If it's simply a middle step used for solving matrix equations, there are far more efficient and numerically stable ways of doing this that don't involve calculating the inverse.

Monday, April 14, 2008 7:29 PM
• Hey Thanks a ton for  sending me the sample code.
Monday, April 14, 2008 8:17 PM