Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have a very large sparse matrix and I need to obtain its approximate inverse and save it as an sparse matrix too. Any of my efforts as could be seen in what follows fail for large n

ClearAll[n, s, f, aInv]
n = 100000;
s = SparseArray[{Band[{1, 120}] -> -2., Band[{950, 1}] -> -1., 
    Band[{1, 1}] -> 20., Band[{1, 100}] -> 2., 
    Band[{6, 800}] -> 1.1}, {n, n}, 0.];
f = LinearSolve[s];
aInv = f[SparseArray[{Band[{1, 1}] -> 1.}, {n, n}, 
     0.]]; // AbsoluteTiming

The accuracy of the approximate inverse is not that important to me for the time being, if it could be obtained and saved in the sparse form then I can improve it by the existing iteration methods. Maybe the built-in pre-conditioners for large sparse matrices in Mathematica 8 could be good. I would be grateful if someone give me some hints to obtain an approximate inverse for such large sparse matrices.

share|improve this question
    
Please consider registering your account. This way you'll have access to all your question at the same place and you'll be able to edit them and clarify things. I just rejected an edit to this post, which may or may not have been done by you, because it would have changed the question slightly and invalidate the existing answer. –  Szabolcs Jun 15 '12 at 16:33
2  
@Fazlollah Next time you edit a question please try to use code formatting features. I would also suggest first posting a comment to ask if the original poster agrees with the change in focus of the example. –  Jens Jun 15 '12 at 21:26
1  
@Jens They don't have enough rep to comment, which is why they made the edit (see edit comment, where they ask how to comment). I disagree with them changing focus without consent and I'm rolling this one back. They had a similar change previously suggested and was rejected by J. M. Typically users w/o rep are advised to wait till they get 50 rep and Fazlollah should have to wait just like everyone else. –  rm -rf Jun 15 '12 at 21:33

2 Answers 2

A possible way for obtaining an approximate inverse is based on Newton-Schulz iterative method, which is given by

$$V_{k+1}=V_{k}(2I-AV_{k}),$$

wherein $I$ is the identity matrix and it converges, when the eigenvalues of $AV_{0}$ are less than one. Following this formula one may obtain an approximate inverse. A starting value could be constructed as Daniel Lichtblau gave above.

In general, the only difficulty is to construct an inital sparse matrix, which preserves convergence in all cases.

share|improve this answer
2  
To make your answer more interesting you should provide a Mathematica code fulfilling your ideas. –  Artes Nov 22 '12 at 11:29

If your matrix is diagonally dominant (in the example it is) then you can do as follows. Start with a diagonal matrix comprised of the reciprocals of the diagonal of your original matrix. Find the residual and use that to form a correction. iterate as long as needed.

Here is your example, scaled down.

ClearAll[n, s, f, aInv]
n = 1000;
s = SparseArray[{Band[{1, 120}] -> -2., Band[{950, 1}] -> -1., 
    Band[{1, 1}] -> 20., Band[{1, 100}] -> 2., 
    Band[{6, 800}] -> 1.1}, {n, n}, 0.];
f = LinearSolve[s];
aInv = f[SparseArray[{Band[{1, 1}] -> 1.}, {n, n}, 
     0.]];

inv1 = SparseArray[Band[{1, 1}] -> 1/Normal[Diagonal[s]], {n, n}, 0.];
sparseIden = SparseArray[Band[{1, 1}] -> 1., {n, n}, 0.];

In[1559]:= residual1 = inv1.s - sparseIden;
Max[Abs[residual1]]

Out[1560]= 0.1

Now we get our first correction term, find the new residual, etc.

delta1 = -residual1.inv1;
inv2 = inv1 + delta1;
residual2 = inv2.s - sparseIden;
Max[Abs[residual2]]

Out[1570]= 0.02

Repeat one more time.

delta2 = -residual2.inv2;
inv3 = inv2 + delta2;
residual3 = inv3.s - sparseIden;
Max[Abs[residual3]]

Out[1566]= 0.0006

In[1573]:= Max[Abs[inv3 - aInv]]

Out[1573]= 0.00003

At this point the residual and the difference between approx and actual inverses, component wise, are bounded by 10^-(3) and 10^(-4) respectively.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.