# 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. Commented Jun 15, 2012 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
Commented Jun 15, 2012 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.
– rm -rf
Commented Jun 15, 2012 at 21:33

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.

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. Commented Nov 22, 2012 at 11:29