Need some facts about Newton-Schulz iterative method and its application to sparse matrices I am studying Newton-Schulz iterative method for obtaining an approximate inverse , 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. 
I want to know whether Newton's method work well for sparse matrices? Also, whether for the sparse matrix $A$, its inverse $V_k$ by Newton's method preserves sparsity? If not, what could be the way so that inverse matrix $V_k$ preserves sparsity for sparse $A$?
Any reference related with this discussion will also be helpful to me.
Thank you very much for your help. 
 A: It sometimes happens that someone re-invents the wheel. It's on page 1-3 of this paper:

Multigrid Calculus

Wrote it more than 10 years ago; cannot judge so well anymore if the text maybe helpful somehow. But anyway, here is a hopefully relevant summary of the first part of the paper.

Newton-Raphson algorithm
The Newton-Raphson method is a numerical algorithm for finding zeros $p$ of a
function $f(x)$. The gist of the method is to draw successive tangent lines and
determine where these lines intersect the x-axis (see figure below):
$$
   y - f(p_n) = f'(p_n).(x - p_n)
   \quad \mbox{where} \quad y = 0
   \quad \mbox{and} \quad x = p_{n+1}
   \quad \Longrightarrow
$$ $$
   p_{n+1} = p_n - \frac{f(p_n)}{f'(p_n)}
$$
Thereby assuming that the whole process will be convergent.
The method can be used for performing a division without actually performing a
division. The equation to be searched for its roots, in this case, is given by:
$$ \frac{1}{x} = a $$
Substitute $f(x) = 1/x - a$ in the above algorithm. Resulting in:
$$
   p_{n+1} = p_n - \frac{1/p_n-a}{-1.1/p_n^2} = p_n - (- p_n + a.p_n^2)
   \quad \Longrightarrow \quad p_{n+1} = 2.p_n - a.p_n^2
$$
It is well known that the Newton-Raphson method, if it converges, then it does
so quadratically, meaning that the inverse $1/a$ can be found rather quickly.
The algorithm has been used in the old days, on computers which had no floating
point division instruction available.



Thus by employing the Newton-Raphson algorithm, the inverse of a number can be
found with quadratic speed, by performing solely additions and multiplications.
Armed with this knowledge, let's make the transition from numbers to matrices.
Determining the inverse of a matrix, filled with many numbers, seems to be much
more like a challenge anyway.
Let the iterative process for matrices be defined by: $$
 P_0 = I  \quad \mbox{and} \quad  P_{n+1} = 2.P_n - A.P_n.P_n $$
Here $I$ is the unity matrix, $A$ is the matrix to be inverted and $P_n$ are
the successive iterands, which should converge to the inverse matrix $A^{-1}$.
Theorem.
$$
\mbox{Let}   \quad M = (I - A) \quad \mbox{or} \quad A = (I - M) \quad
\mbox{then:} \quad P_n = (I - M^{2^n}) . A^{-1}
$$
Proof by induction:
$$
\begin{array}{lll}
P_0 &=& I = A . A^{-1} = (I - M) . A^{-1} \\
P_{n+1} &=& \left(I - M^{2^{n+1}} \right) . A^{-1} \\
        &=& \left(I - M^{2^n.2} \right) . A^{-1}  \\
        &=& \left[I - (M^{2^n})^2 \right] . A^{-1} \\
        &=& \left(I - M^{2^n} \right) .
            \left(I + M^{2^n} \right) . A^{-1} \\
        &=& \mbox{because all matrices are mutually commutative:} \\
        &=& \left(I - M^{2^n} \right) . A^{-1} .
            \left[ 2.I - A.\left( I - M^{2^n} \right).A^{-1} \right] \\
        &=& P_n.(2.I - A.P_n) = 2.P_n - A.P_n.P_n
\end{array}
$$
Let $m = 2^n$, then:
$$
  P_n = (I - M^m) . (I-M)^{-1}
$$
For numbers, this would be the sum of a Geometric series :
$$
   (I - M^m).(I - M)^{-1} =
   (I + M + M^2 + M^3 + M^4 + M^5 + M^6 + ... + M^{m-1})
$$
For matrices, this Geometric Series turns out to be equivalent to an iterative
"incremental Jacobi" solution method, as will not be futher explained here.
But there also does exist a product of terms, called the Euler expansion
(don't remember where the naming comes from; a reference would be quite welcome):
$$ 
  (I - M^m).(I - M)^{-1} = \\
  (I + M^{m/2}).(I - M^{m/2}).(I - M)^{-1} = \\
  (I + M^{m/2}).(I + M^{m/4}).(I - M^{m/4}).(I - M)^{-1} = \\
  (I + M^{m/2}).(I + M^{m/4}). \, ... \, (I+M^2).(I+M).(I-M).(I-M)^{-1} = \\
  (I + M^{m/2}).(I + M^{m/4}). \, ... \, (I + M^8).(I + M^4).(I + M^2).(I + M)
$$
Let the system of equations to be solved be given by $A.w = b$. Using the above
sequences, we can write:
$$
    P_n = (I - M^{2^n}) . A^{-1} \quad \Longrightarrow \quad
    A^{-1} b = (I - M^{2^n})^{-1} . P_n b
$$ $$
 = (I - M^{2^n})^{-1} .
 (I + M^{m/2}).(I + M^{m/4}). \, ... \, (I + M^8).(I + M^4).(I + M^2).(I + M) b
$$
Another way of looking at this is the following:
$$
 (I-M)^{-1} = (I-M)^{-1}.(I+M)^{-1}.(I+M) = (I-M^2)^{-1}.(I+M) = \\
 (I-M^2)^{-1}.(I+M^2)^{-1}.(I+M^2).(I+M) = (I-M^4)^{-1}.(I+M^2).(I+M)
$$
With other words:
$$
 w = (I-M)^{-1}.b = (I-M)^{-1}.(I+M)^{-1}.(I+M).b = (I-M^2)^{-1}.(I+M).b
$$
Define $b := (I + M).b$ and $M := M^2$. Then again:
$$
 w = (I-M)^{-1}.b = (I-M)^{-1}.(I+M)^{-1}.(I+M).b = (I-M^2)^{-1}.(I+M).b
$$
And this process can be repeated until we find a way to determine $(I-M)^{-1}$,
preferrably without continuing the iterations ad infinitum.
Update. Quite in general, the inverse of a sparse matrix
is not sparse at all. And this, of course, will be
reflected in the iterands of the Newton-Schulz method for obtaining
an approximate inverse. The common remedy is simply not to calculate
the inverse of a large sparse matrix and use e.g. LU decomposition instead.
But, as argued in the abovementioned reference,
the Newton-Schulz method may be regarded as an Ansatz to MultiGrid
methods. And that is a quite different story, not to be elaborated
here.
