LinuxQuestions.org

LinuxQuestions.org (/questions/)
-   Linux - Newbie (https://www.linuxquestions.org/questions/linux-newbie-8/)
-   -   How to find inverse of a sparse matrix? (https://www.linuxquestions.org/questions/linux-newbie-8/how-to-find-inverse-of-a-sparse-matrix-816111/)

smp 06-24-2010 07:01 AM

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

alli_yas 06-24-2010 07:20 AM

This sounds suspiciously like homework...

johnsfine 06-24-2010 07:25 AM

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.

pixellany 06-24-2010 05:11 PM

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.

pixellany 06-24-2010 05:14 PM

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...

smp 06-24-2010 11:35 PM

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

johnsfine 06-25-2010 05:40 AM

Quote:

Originally Posted by pixellany (Post 4014058)
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 (Post 4014339)
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.

syg00 06-25-2010 06:30 AM

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.

johnsfine 06-25-2010 07:32 AM

Quote:

Originally Posted by syg00 (Post 4014615)
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.

syg00 06-25-2010 07:43 AM

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.

johnsfine 06-25-2010 07:53 AM

Quote:

Originally Posted by syg00 (Post 4014678)
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.

saikee 06-25-2010 08:07 AM

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.

johnsfine 06-25-2010 08:19 AM

Quote:

Originally Posted by saikee (Post 4014710)
Code:

        S = aa(j, i) / aa(i, i)

Why is it safe for this code to assume that aa(i,i) is nonzero?

saikee 06-25-2010 08:28 AM

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.

johnsfine 06-25-2010 09:41 AM

Quote:

Originally Posted by saikee (Post 4014727)
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.


All times are GMT -5. The time now is 04:40 PM.