Matrix Inverse Algorithm
- Thread starter Tsotne
- Start date
I have defined a linear algebra and statistics library in C# and gonna do the same for the rest of mathematics. Then Ill extend those codes to C++. I have developed many inverse algorithms but one problem is that the difficulty increases exponentially when adding one dimension to a matrix. So if we are at 200x200 matrix and have some difficulty n. Once we add the 1 more dimension, (201x201) the difficulty rises to n^3 that is directly translated into the processing time. The fastest algorithm I developed once was only able to calculate the 100x100 matrix inverse in 30 seconds. MATLAB is pretty faster in that but still, the dimension is restricted there also.
A great explanation on why computing a matrix inverse is usually a bad idea: http://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/
Daniel Duffy
C++ author, trainer
- Joined
- 10/4/07
- Messages
- 10,469
- Points
- 648
Matrix inversion is something that fell out of fashion in the 19th century. I think what you want is to solve
Ax = b
or am I missing something?
Thanks for your replies. Just a hypothetical case when you want to analyze the whole market including all the stocks traded on NYSE. (Covariance/Correlation matrix) I'm not trying to solve the system at this moment. I can do it many different ways and have done so. For the library, I need to include the matrix inverse function in the class. (Not for a particular purpose) So I'm trying to get suggestions about the fastest one. I have defined it many ways in programming languages C++/C#. When I execute the program, the method takes 15 seconds with 100x100 matrix. And when I add the dimensions the processing time increases much. I have learned C++/C# from the mathematical side so I'm more concerned and aware of those languages in math point of view rather than the programming constructs themselves. So what I'm actually interested in is: Can I write the code in some manner that breaks into the programming convictions that could decrease the compiling time? Many thanks for replies
Daniel Duffy
C++ author, trainer
- Joined
- 10/4/07
- Messages
- 10,469
- Points
- 648
Actually, I thought about this a bit more,. In some cases it is useful to precompute the inverse to avoid having to solve the same system at each iteration.
But the danger is inverse(A) is a dense matrix even when A is sparse.
There are many methods to compute inverse(A) but I cannot comment as I have not examined them.
But the danger is inverse(A) is a dense matrix even when A is sparse.
There are many methods to compute inverse(A) but I cannot comment as I have not examined them.
Converting matrix into row-echelon form
Thanks Mr Daniel for Responses. I have defined the algorithm of converting the given matrix into a row-echelon form. I have been working on it since the last post so it is quite quick-made. It is working I checked for several times but has some difficulties when increasing the dimensions. (But don't pay attention to details for now, I did it in short time so I'll correct everything). You are not of course asked to evaluate this code into details since it is quite hard to pin down every detail there. Im just showing the difficulty related to my last question.
I know a lot in C# already and have such imagination on the compiler execution: When compiler starts working, it mobilizes all the resources to be prepared to the matching code during the execution phase. So my interest here is if it is possible to decline the compiler to mobilize and hold ready (loaded) all the resources while executing and force it to drop the non-useful resources in that particular code.
It's like drawing the graph and checking for objects (or dependence of objects) if they are rooted while the garbage collection is to be performed. In the similar manner, im interested if we can "root" the
resources and drop ones not needed that could increase the processing time. Thanks
Thanks Mr Daniel for Responses. I have defined the algorithm of converting the given matrix into a row-echelon form. I have been working on it since the last post so it is quite quick-made. It is working I checked for several times but has some difficulties when increasing the dimensions. (But don't pay attention to details for now, I did it in short time so I'll correct everything). You are not of course asked to evaluate this code into details since it is quite hard to pin down every detail there. Im just showing the difficulty related to my last question.
I know a lot in C# already and have such imagination on the compiler execution: When compiler starts working, it mobilizes all the resources to be prepared to the matching code during the execution phase. So my interest here is if it is possible to decline the compiler to mobilize and hold ready (loaded) all the resources while executing and force it to drop the non-useful resources in that particular code.
It's like drawing the graph and checking for objects (or dependence of objects) if they are rooted while the garbage collection is to be performed. In the similar manner, im interested if we can "root" the
resources and drop ones not needed that could increase the processing time. Thanks
Code:
static class Lin_Alg
{
public static double[,] Gauss_Jordan(double[,] A)
{
double[,] B = RemoveZeroRows(A);
int m = B.GetLength(0);
int n = B.GetLength(1);
int c = -1; int r = 0;
int t = 0; int k = 0;
double[] R = new double[n];
double[] K = new double[n];
double p = -1; double pp = -1;
while (k < m)
{
//Finding the first non-zero column index*******************************
for (int i = t; i < n; i++)
{
for (int j = k; j < m; j++)
{
if (B[j, i] != 0)
{
c = i;
break;
}
}
if (c == i)
break;
}
//***********************************************************************
//Finding the first non-zero entry in the column "c" to start iterating**
for (int i = k; i < m; i++)
{
if (B[i, c] != 0)
{
r = i;
break;
}
}
//*************************************************************************
//Interchanging rows(if neccessary)to get the first entry of the located...
//...column non-zero*******************************************************
if (r != k)
{
for (int i = 0; i < n; i++)
{
K[i] = B[k, i];
R[i] = B[r, i];
B[k, i] = R[i];
B[r, i] = K[i];
}
}
//****************************************************************************
//Checking for the first first entry from the previous iteration if it is 1...
//...if not divide the rows by the multiplicative inverse of that number******
p = B[k, c];
if (p != 1)
{
for (int i = 0; i < n; i++)
B[k, i] *= Math.Pow(p, -1);
}
//Then multiplying the first number times other non-zero entry rows to get all.
//...numbers equal to 0 below the selected number*********************
for (int i = 0; i < m; i++)
{
if (i == k)
continue;
else
{
if (B[i, c] != 0)
{
pp = B[i, c];
for (int j = 0; j < n; j++)
{
B[i, j] -= pp * B[k, j];
}
}
}
}
//********************************************************************
//Adjusting the indexes for the next iteration************************
t = c + 1; k++;
}
//************************************************************************
//Adding the removed zero rows if there were any**************************
if (GetZeroRows(A) != 0)
{
double[,] G = new double[m + GetZeroRows(A), n];
G = B;
for (int i = m + 1; i <= GetZeroRows(A); i++)
{
for (int j = 0; j < n; j++)
{
G[i, j] = 0;
}
}
return G;
}
else
return B;
}
private static double[,] RemoveZeroRows(double[,] A)
{
int m = A.GetLength(0);
int n = A.GetLength(1);
//int c = 0; int r = 0;
int q = 0;
int[] zeroIndexes = ZeroRowIndexes(A);
bool index = false;
//Remove the zero rows (if any) and free up the matrix with nonzero parameters*****
int w = GetZeroRows(A);
int l = m - w;
if (w != 0)
{
double[,] B = new double[l, n];
for (int i = 0; i < m; i++)
{
for (int j = 0; j < w; j++)
{
if (zeroIndexes[j] != i)
index = true;
else
{
index = false;
break;
}
}
if (index)
{
for (int k = 0; k < n; k++)
B[q, k] = A[i, k];
q++;
}
}
return B;
}
else
return A;
}
private static int[] ZeroRowIndexes(double[,] A)
{
int m = A.GetLength(0);
int n = A.GetLength(1);
int r = GetZeroRows(A);
int[] zeroRowIndexes = new int[r];
int q = 0; bool zero = false;
if (r != 0)
{
for (int i = 0; i < m; i++)
{
if (A[i, 0] == 0)
{
for (int j = 0; j < n; j++)
{
if (A[i, j] == 0)
{
zero = true;
zeroRowIndexes[q] = i;
}
else
{
zero = false;
break;
}
}
if (zero)
q++;
}
}
return zeroRowIndexes;
}
else
{
return null;
}
}
private static int GetZeroRows(double[,] A)
{
int m = A.GetLength(0);
int n = A.GetLength(1);
int r = m; int k = 0;
for (int i = 0; i < m; i++)
{
for (int j = 0; j < n; j++)
{
if (A[i, j] != 0)
{
k++;
break;
}
}
}
r -= k;
return r;
}The one thing the professor in numerical analysis nailed in my head is you never use matrix inverse to solve systems of equations b\c it is computationally expensive. Gauss-Jordan elimination is about 3 times slower than LU for solving systems of equations. I do not know other algorithms for matrix inverse.
Your 2nd and 3rd paragraph do not make sense. You want to speed up compilation or execution? Your compilation time should be alright, to speed up execution you can try to parallelize your code.
On a general note I would forget about C# for now and first implement the algorithms in C++, then if you won't run out of steam implement in them C#. Also, I would take a look at Boost library, its matrix library is primitive but its implemented as inline functions which is more convenient to use.
Your 2nd and 3rd paragraph do not make sense. You want to speed up compilation or execution? Your compilation time should be alright, to speed up execution you can try to parallelize your code.
On a general note I would forget about C# for now and first implement the algorithms in C++, then if you won't run out of steam implement in them C#. Also, I would take a look at Boost library, its matrix library is primitive but its implemented as inline functions which is more convenient to use.
Thank you vincegata for reply. I don't need to solve the system with the Gauss-Jordan algorithm. Ordinary Gaussian algorithm solves the system of linear equations if needed but now Im not focusing on the algorithm itself. I implemented it just for the sake of interest. The point of this thread was how to decrease the execution time. I said in the previous posts that increasing the number of dimension in a square matrix increases the execution time catastrophically. Matrix inverse is just the example. Generally many flow-heavy codes have similar problems. I can come up with many ways how to solve the problem mathematically. Algorithms to solve the system of equations(no matter whether they are linear or nonlinear) are nothing new to me. But any suggestion from the side of programming how to decrease the time of execution will be new. For example consider the simple world where Gauss has not yet been born and the only way to solve the system is matrix inverse. What now? Provided we have no mathematical tools, can we do anything with programming manoeuvre to compute such a matrix??? Thanks again
- Joined
- 8/4/06
- Messages
- 44
- Points
- 18
Try A = LU decomp first (if applicable, even go for Cholesky if possible). Then, (if L and U are invertable), inv(A) = inv(LU) = inv(U) * inv(L), and inverting L and U are quick. Well, this works for me, 0.13 sec for 100x100, 18 sec for 500x500, 260 sec for 1000x1000, and my laptop is a dinosaur.
Daniel Duffy
C++ author, trainer
- Joined
- 10/4/07
- Messages
- 10,469
- Points
- 648
Interesting approach,. Some issues
1. Choosing the 'stable' L and U; preprocessing might be needed.
2. Increased memory consumption (inv(L) is dense).
3. inv(L) and inv(U) could be done in parallel (e.g. parallel sections).
Is it possible to use QR to find inverse in the same way? QR is more efficient than LU.
Not sure if it is more efficient than LU. You might find this useful:
http://en.wikipedia.org/wiki/QR_decomposition
- Joined
- 3/21/08
- Messages
- 594
- Points
- 38
it is not only available RAM and/or HDD space, there is something else you should consider
- Joined
- 8/4/06
- Messages
- 44
- Points
- 18
I agree, especially with #3, thanks for pointing that out Daniel.
Matrix inversion is O(n^3), and Gauss-Jordan is the way to go if you want to do it by hand. Performance-wise, the only jumps you should see when increasing n should be when you reach points when your matrix is not able to fit cache (L1, then L2) anymore; otherwise the execution time should increase with n^3. If you want to do calculations by hand, use C or something alike - C# is not the language for doing numerical work. For further performance improvements, one would try to parallelize calculations - how to accomplish that depends on parallelization approach (specialized instructions like SSE, multi-threading, GPU computing, distributed computing etc.) you want to employ. If you want to check how fast good code could be, find some LAPACK library (I guess MKL would do), and check out DGETRI routine (it works by doing LU factorization, then inverting U, and finally solving the system A^(-1)*L=U^(-1)). Matlab is probably using the same routine under the hood, so Matlab execution times should be good indicator on how fast calculation could be on particular machine.
We can surely achieve the same speed and performance using C# as C++ with one of the tools you provided. Actually when we direct the compiler to use less resources (which are not needed) or multitask the current activity, we have the same effect(if done properly) as C/C++. C# compiler is loading LINQ, .NET4 new features different from C++ and executes the code such. That's like a man carrying 5 heavy bags. In C++, there are row, less features and it looks like another man carrying 3 heavy bags. If we tell that first man(C#) to drop unnecessary bags(LINQ,.NET4 features,new forms,just for that particular purpose of course) then they can surely go walk the same speed.
In case of arrays for example, we don't have such prepared arrays with any dimension as in C#. In VB C/C++ we have dynamic massive, which certainly is harder for use right? Also, does LINQ techniques help mathematicians perform their computations more effectively than in language lacking it??? So if you are talking about speed, I somehow agree, but in case of disadvantage in ability to perform numerical computations I strongly disagree. I have experience with both, C/C++ and now C# and see many differences made from math point of view. Something I could do in C++ in tens of lines are compacted and can do in 5 lines.
In case of arrays for example, we don't have such prepared arrays with any dimension as in C#. In VB C/C++ we have dynamic massive, which certainly is harder for use right? Also, does LINQ techniques help mathematicians perform their computations more effectively than in language lacking it??? So if you are talking about speed, I somehow agree, but in case of disadvantage in ability to perform numerical computations I strongly disagree. I have experience with both, C/C++ and now C# and see many differences made from math point of view. Something I could do in C++ in tens of lines are compacted and can do in 5 lines.
Similar threads
