Home Forums HCL Reviews Tutorials Articles Register Search Today's Posts Mark Forums Read
 LinuxQuestions.org How to find inverse of a sparse matrix?
 Linux - Newbie This Linux forum is for members that are new to Linux. Just starting out and have a question? If it is not in the man pages or the how-to's this is the place!

Notices

 06-24-2010, 07:01 AM #1 smp LQ Newbie   Registered: Jun 2009 Posts: 21 Rep: How to find inverse of a sparse matrix? Hi, I have a matrix of dimensions 100 x 100. It is a sparse matrix. It has 460 non-zero values (9540 are zero values). I want to find inverse of this matrix. Any help will be appreciated. Regards, smp
 06-24-2010, 07:20 AM #2 alli_yas Member   Registered: Apr 2010 Location: Johannesburg Distribution: Fedora 14, RHEL 5.5, CentOS 5.5, Ubuntu 10.04 Posts: 559 Rep: This sounds suspiciously like homework...
 06-24-2010, 07:25 AM #3 johnsfine LQ Guru   Registered: Dec 2007 Distribution: Centos Posts: 5,286 Rep: What are your other constraints? 100x100 is pretty tiny for a sparse matrix. Do you need an algorithm to invert some enormous number of matrices like that efficiently? Or do you just need to get one (or a few) inverted? Have you installed Octave yet? You should. Info here: http://www.gnu.org/software/octave/about.html But try finding a package from your distribution (Red Hat?) or from Fedora, before falling back to the generic Linux install methods offered by gnu.org If you have problems like inverting that matrix that are in support of some other activity (studying or practicing some branch of engineering, science or math that needs operations like that), learn how to use a tool such as Octave to do the work for you. But if alli_yas guessed right and this is just an example problem in a class were the point is for you to learn programming, don't ask us to do your work for you. Octave is still a nice tool to have to check your results, but obviously your instructor won't accept using Octave as your answer if the point was for you to learn how to do the programming. Last edited by johnsfine; 06-24-2010 at 07:30 AM.
 06-24-2010, 05:11 PM #4 pixellany LQ Veteran   Registered: Nov 2005 Location: Annapolis, MD Distribution: Arch/XFCE Posts: 17,802 Rep: What is a "sparse matrix" and how would that affect the algorithm for inverting? I have to assume that the web is awash in subroutines for inverting matrices.
 06-24-2010, 05:14 PM #5 pixellany LQ Veteran   Registered: Nov 2005 Location: Annapolis, MD Distribution: Arch/XFCE Posts: 17,802 Rep: PS; I am always amused / curious when someone posts a homework question like this and then logs out 1 minute later. Here we are--10 hours later--and the OP has not been back. Sometimes, I have an image of a struggling student sitting in a classroom plugging test questions into their Smart Phone. In fairness, the OP here does not fit the typical homework pattern... Last edited by pixellany; 06-24-2010 at 05:17 PM.
06-24-2010, 11:35 PM   #6
smp
LQ Newbie

Registered: Jun 2009
Posts: 21

Original Poster
Rep:
Hi,

(I am not student. I am not asking you to teach me C programming.
And this is not a homework question.)

There are Numerical recipes functions namely ludcmp() and lubksb() to invert the matrix. I used them for a simple 3x3 matrix. I got the inverse. Then I mulitiplied that inverse with the original matrix. As expected I got the identity matrix.

Now the same thing I tried for my 100x100 sparse matrix. I got the inverse. So to check, I multiplied this inverse by original matrix. But I did not get Identity matrix.

So I posted my problem.

Quote:
 What are your other constraints? 100x100 is pretty tiny for a sparse matrix. Do you need an algorithm to invert some enormous number of matrices like that efficiently? Or do you just need to get one (or a few) inverted?
At present, I need algorithm to invert only one sparse matrix.
We know that not every matrix is invertible. So do I need to check first that whether my sparse matrix is invertible? If yes, what I have to check actually ? One theorem from Linear Algebra says that if the matrix is invertible, then the columns of that matrix form a linearly independent set. Is it helpful for such check? Is there easy way to check it?

By the way, I am sorry if I had asked wrong question at wrong place.

Regards,
smp

06-25-2010, 05:40 AM   #7
johnsfine
LQ Guru

Registered: Dec 2007
Distribution: Centos
Posts: 5,286

Rep:
Quote:
 Originally Posted by pixellany What is a "sparse matrix" and how would that affect the algorithm for inverting?
A "sparse matrix" is stored differently than a "dense matrix". In a dense matrix you store the value of each position. In a sparse matrix, a more complicated data structure is used to identify which positions hold nonzero values and to store only the values and positions of the nonzero positions. Typically you need to store one index and value for each nonzero plus N extra indexes for an NxN sparse matrix. You don't need to store two indexes for each nonzero.

The complexity of identifying which positions are non zero typically makes each non zero almost twice as expensive to store as in a dense matrix and much more than twice as expensive to process. But then the zeroes are entirely free, neither stored nor processed. So if the whole matrix has almost all positions zero, the sparse matrix will be overwhelmingly less expensive to store and process.

Any algorithm for operating on a sparse matrix is more complex than an analogous algorithm on a dense matrix simply because of the extra complexity of the storage. In addition, many algorithms on sparse matrices are done with major extra complexity beyond the above in order to preserve the low ratio of non zero to zero positions through various intermediate steps.

Matlab includes a very powerful collection of functions for operating on sparse matrices. These are far more efficient, reliable and numerically stable than anything an amateur might code. I don't have experience with Octave, but I trust various claims I have read online that Octave is an acceptable substitute for Matlab.

Quote:
 Originally Posted by smp At present, I need algorithm to invert only one sparse matrix. We know that not every matrix is invertible. So do I need to check first that whether my sparse matrix is invertible? If yes, what I have to check actually ? One theorem from Linear Algebra says that if the matrix is invertible, then the columns of that matrix form a linearly independent set. Is it helpful for such check? Is there easy way to check it?
If you want to invert a sparse matrix, or you want to know whether a sparse matrix can be correctly inverted and if not, why not or any similar operations or questions, the obvious approach is to use a tool such as Octave.

The algorithms are difficult. I wouldn't have the patience to explain even the simplest of them here even if I didn't have an extra constraint (it could get close to topics I can't post about without violating intellectual property rights).

If you don't explain a good reason to attempt these algorithms yourself, rather than using Octave, I assume you are just ignoring the best advice.

If you do have a good reason, and especially if you expect to do any aspect of it better than Octave does, I think you have a very difficult project in front of you.

Last edited by johnsfine; 06-25-2010 at 07:22 AM. Reason: fixed minor typo

 06-25-2010, 06:30 AM #8 syg00 LQ Veteran   Registered: Aug 2003 Location: Australia Distribution: Lots ... Posts: 15,620 Rep: Matric manipulation was another of those subjects on the (seemingly endless) list of useless things regarded as absolutely essential by my University mentors. Maybe I just fell into the wrong branch of IT. From a position of of monumental ignorance I'd be prepared to second johnsfine.
06-25-2010, 07:32 AM   #9
johnsfine
LQ Guru

Registered: Dec 2007
Distribution: Centos
Posts: 5,286

Rep:
Quote:
 Originally Posted by syg00 Matric manipulation was another of those subjects on the (seemingly endless) list of useless things regarded as absolutely essential by my University mentors. Maybe I just fell into the wrong branch of IT. From a position of of monumental ignorance I'd be prepared to second johnsfine.
But I'll disagree with what you're apparently saying. Specifically that "useless" part.

If the thing you are learning or working with requires the use of matrices, then I think it is appropriate to skip learning the underlying detailed mechanics of sparse matrices and instead use an existing tool.

But if you are learning software engineering, one of the many things you ought to learn is the underlying mechanics of sparse matrices, because that is a very well researched topic containing concepts and issues that will help you be more effective in other areas of software engineering even if you will never do any real work with sparse matrices.

In either case, anyone in the field ought to learn basic concepts of matrices, including what inverting a matrix means, how you would do it, why you would do it, and why/what you would do something else instead if the size and/or condition number makes inverting the matrix an ineffective intermediate step toward your actual goal.

 06-25-2010, 07:43 AM #10 syg00 LQ Veteran   Registered: Aug 2003 Location: Australia Distribution: Lots ... Posts: 15,620 Rep: O.K., maybe I phrased it wrong, or maybe I just got it wrong. Most of the (pure) mathematics I learned were bugger all use to me earning a living in the real world. And I think I've made a pretty good living out of it all things considered. Perhaps the "ivory towers" meets "real world" conundrum. All the pressure to do post-graduate and masters seems bloody silly viewed from this distance.
06-25-2010, 07:53 AM   #11
johnsfine
LQ Guru

Registered: Dec 2007
Distribution: Centos
Posts: 5,286

Rep:
Quote:
 Originally Posted by syg00 Most of the (pure) mathematics I learned were bugger all use to me earning a living in the real world. And I think I've made a pretty good living out of it all things considered. Perhaps the "ivory towers" meets "real world" conundrum. All the pressure to do post-graduate and masters seems bloody silly viewed from this distance.
Now I've managed to maneuver myself into the opposite side of the big argument from my usual position. In general, I'm on your side. I almost always grab for the common sense approach over the approach supported by the ivory towers theorists and I'm usually right in that choice.

But there are exceptions. Some of the BS that is researched and taught provides important insights or even directly usable algorithms for the practical software engineer.

 06-25-2010, 08:07 AM #12 saikee Senior Member   Registered: Sep 2005 Location: Newcastle upon Tyne UK Distribution: Any free distro. Posts: 3,398 Blog Entries: 1 Rep: Rather talk about it I post a subroutine I wrote Code: Public Sub Gauss01(aa(), No_row) Rem local i,j,k,L,s Rem Standard Gaussian elimination of a full size matrix [a]{x}=y{} Rem aa(i,j) is the matrix of size No_row, vector {x} is the unknown Rem y(i) is appended to left hand side of [a] and stored as aa(i,n+1) Rem on exit the solution {x} of the matrix is store as aa(i,n+1) Rem Developed by DR. Sam Watt of Merz and Mclellan 1998 Rem Two progress bars 1 & 2 may be activated to show progress if needed For i = 1 To No_row Rem ProgressBar1 = 100 * i / No_row For j = i + 1 To No_row If aa(j, i) <> 0 Then S = aa(j, i) / aa(i, i) For K = i To No_row + 1 If aa(i, K) <> 0 Then aa(j, K) = aa(j, K) - aa(i, K) * S Next K End If Next j Next i Rem ----------------------------------back substitution----- For i = 1 To No_row Rem ProgressBar2 = 100 * i / No_row L = No_row + 1 - i S = 0 If L <> No_row Then For j = L + 1 To No_row S = S + aa(L, j) * aa(j, No_row + 1) Next j End If aa(L, No_row + 1) = (aa(L, No_row + 1) - S) / aa(L, L) Next i End Sub This may not answer the OP as I use a Gaussian elimination which is mathemtically exact but it doesn't deal with the sparseness. The size of the matrix can be as big as your operating system permitted. The power of the routine is all loadcases can be solved in one sweep. The code was originally written in Fortran and now in VB6. I dealt with only symmetric sparse matrices with which I could store only one half of it. Also the matrix can be aligned vertically with maximum band width thus reducing the storage requirement. The above rountine store the full matrix. AFAIK in the normal sense the sparseness of a matrix we could assemble part of it, carry out the necessary mathemical operations, write part of the matrix onto a backing storage unit, shift the lower half upward, add additional section and then repeat the whole process until the whole matrix is done. Since only a portion of the matrix is assembled and solved so the bandwidth with the known zero elements has no needed to be assembled, stored or operated on, making the process more efficient and faster. In the extreme case the entire matrix can be stored in a vector of one element wide but the solution is based on the enclosed steps. Unless the properties of the sparse matrix is accurately defined it could be difficult to write a computer code to make use of the festures in the inversion. Last edited by saikee; 06-25-2010 at 08:34 AM.
06-25-2010, 08:19 AM   #13
johnsfine
LQ Guru

Registered: Dec 2007
Distribution: Centos
Posts: 5,286

Rep:
Quote:
 Originally Posted by saikee Code: S = aa(j, i) / aa(i, i)
Why is it safe for this code to assume that aa(i,i) is nonzero?

Last edited by johnsfine; 06-25-2010 at 09:42 AM.

 06-25-2010, 08:28 AM #14 saikee Senior Member   Registered: Sep 2005 Location: Newcastle upon Tyne UK Distribution: Any free distro. Posts: 3,398 Blog Entries: 1 Rep: johnsfine , Because the line before it will do a bypass if it isn't. The matrice I used to do were based on a physical structural system which should never have zero in the diagonal otherwise the system is a mechanism and unstable. Last edited by saikee; 06-25-2010 at 08:31 AM.
06-25-2010, 09:41 AM   #15
johnsfine
LQ Guru

Registered: Dec 2007
Distribution: Centos
Posts: 5,286

Rep:
Quote:
 Originally Posted by saikee Because the line before it will do a bypass if it isn't.
False. The line before tested j,i not i,i.

Quote:
 The matrice I used to do were based on a physical structural system which should never have zero in the diagonal otherwise the system is a mechanism and unstable.
No zeroes on the diagonal is less of a constraint than "diagonally dominant". No zeroes on the diagonal still allows an intermediate value on the diagonal to be zero due to cancellations, making your code unsound.

If you wrote and posted an algorithm that is only valid for diagonally dominant input (as I think you have). You should have warned anyone who might use your code and might be unaware of that severe restriction.

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is Off HTML code is Off Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post frenchn00b Linux - Desktop 2 08-20-2009 10:00 AM brazilnut Programming 4 10-17-2008 06:30 PM Asuralm Programming 1 05-22-2008 10:52 AM johnpaulodonnell Programming 4 04-30-2008 01:45 PM ztdep Programming 2 04-01-2007 11:10 PM

All times are GMT -5. The time now is 08:26 PM.

 Contact Us - Advertising Info - Rules - LQ Merchandise - Donations - Contributing Member - LQ Sitemap -