# Is there any way to obtain an approximate inverse for very large sparse matrices?

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.

-
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
@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
@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. – R. M. Jun 15 '12 at 21:33

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.

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

-