# True divide and conquer inversion of large matrices

In http://math.stackexchange.com/questions/2735/solving-very-large-matrices-in-pieces there is a way shown to solve matrix inversion in pieces. Is it possible to turn it into a true divide and conquer algorithm?

I am refering to the answer in the referenced question that starts with:

$\textbf{A}=\begin{pmatrix}\textbf{E}&\textbf{F}\\\\ \textbf{G}&\textbf{H}\end{pmatrix}$

And computes the inverse via:

$\textbf{A}^{-1}=\begin{pmatrix}\textbf{E}^{-1}+\textbf{E}^{-1}\textbf{F}\textbf{S}^{-1}\textbf{G}\textbf{E}^{-1}&-\textbf{E}^{-1}\textbf{F}\textbf{S}^{-1}\\\\ -\textbf{S}^{-1}\textbf{G}\textbf{E}^{-1}&\textbf{S}^{-1}\end{pmatrix}$, where $\textbf{S} = \textbf{H}-\textbf{G}\textbf{E}^{-1}\textbf{F}$

In this algorithm I will have some steps:

...
sequential
1. Compute inverse $\textbf{E}^{-1}$
2. Compute inverse $(\textbf{H}-\textbf{G}\textbf{E}^{-1}\textbf{F})^{-1}$
...

Would be possible to device an algorithm that has steps:

...
parallel
1. Compute inverse $\textbf{X}^{-1}$
1. Compute inverse $\textbf{Y}^{-1}$
...

i.e. to identify two smaller matrices $\textbf{X}$ and $\textbf{Y}$ which are independent from their inverse?

Best Regards

-
There are several parallel algorithms ($\mathrm{NC}^2$, or more precisely, $\mathrm{TC}^1$) for the computation of determinants (and therefore matrix inversion), such as the one by Berkowitz dl.acm.org/citation.cfm?id=495 . – Emil Jeřábek Oct 17 '11 at 17:06
Not yet sure whether this will get me, where I want to go. Does TC^1 imply true divide and conquer? – Countably Infinite Oct 17 '11 at 18:26
Probably not, but that’s hard to tell without a formal definition of what “true divide and conquer” means. I think you’d have hard time inventing a sensible definition which would disallow algorithms that first ask for $X^{-1}$ and $Y^{-1}$ for some bogus $X,Y$, and then compute $A^{-1}$ without actually using this information. Where the Berkowitz algorithm gets you is an efficiently parallelizable algorithm, that’s all. – Emil Jeřábek Oct 17 '11 at 19:01
The above original algorithm is recursive. You pick E, F, G, H. And then for example if you need E^-1, then you apply the algorithm again. But problem is algorithm must work sequentially, it can only invoke itself again for S^-1, after it has computed E^-1. So if you can find some bogus X^-1 and Y^-1 in parallel and then construct A^-1 from it, without invoking inversion again, then we would be done. – Countably Infinite Oct 17 '11 at 22:30
Just to mention that you find this algorithm in my book on Matrices (Springer GTM 216). I use it to show that the complexity of matrix inversion is the same as the complexity of matrix multiplication. I don't discuss the parallelization. I don't know the origin of this algorithm but it is not very new. – Denis Serre Oct 18 '11 at 5:26

This type of divide and conquer approach is often used within H-matrix arithmetic, but it is known not very amenable to parallelism (see Ronald Kriemann's dissertation on parallel H-matrix arithmetic). However, for sparse matrices, it is possible to use the basic framework of the multifrontal method in order to simultaneously factor submatrices [start on slide 26 of http://www.hlnum.org/publications/winterschool2009-2.pdf]. Notice that your above scheme is trivially parallelizable for triangular matrices, so being able to form an LU decomposition with a parallel divide and conquer approach would effectively solve your problem.

If your matrix has structure, you may want to look into Newton-Schulz iteration, which is essentially a Newton iteration for the matrix inverse function, and it is known to be stable (see, for example, Higham's "Accuracy and Stability of Numerical Algorithms").

-
Thanks for the many information. Have first to look at each piece. – Countably Infinite Oct 17 '11 at 21:49