# Are there any decompositions of a symmetric matrix that allow for the inversion of any submatrix?

I am given a $J \times J$ symmetric matrix $\Sigma$. For a subset of $\{1, ..., J\}$, $v$, let $\Sigma_{v}$ denote the associated square submatrix. I am in need of an efficient scheme for inverting $\Sigma_v$ for potentially any subset $v$. Essentially, I will have to invert many different submatricies $\Sigma_v$, but I won't know which ones I need to invert in advance of running some program; I would rather invest in a good matrix decomposition at the outset, if one exists (or otherwise get whatever information necessary, if not a decomposition).

I've messed around with the eigen decomposition a little bit, but wasn't able to coerce the inverse of a submatrix out of it.

UPDATE: Apparently the term for the type of submatricies I want to invert is "principal submatrix". I'm wondering if I can't make some progress on this via the Cholesky decomposition. Suppose I'm content to calculate $\Sigma_{jj} ^ {-1}$ for any $j \in \{1, 2, ..., J\}$, where $\Sigma_{jj}$ denotes the submatrix obtained by deleting rows/columns greater than $j$. Write $\Sigma = LL^T$ and let $Q = L^{-1}$. Write $$L = \begin{pmatrix}L_1 & \mathbf 0 \\ B & L_2\end{pmatrix}, Q = \begin{pmatrix}Q_1 & \mathbf 0 \\ C & Q_2\end{pmatrix}$$ where $L_1$ and $Q_1$ and $j \times j$. It follows that $\Sigma_{jj} = L_1 L_1^T$ and $Q_1 = L_1 ^{-1}$ so that $\Sigma_{jj} ^{-1} = Q_1^T Q_1$. So, once I have the Cholesky decomposition I have the inverses of the leading principal submatricies. This doesn't solve the problem as stated since I may need to deal with other principal submatricies, but it should be a useful partial solution.

-
Would this be useful? It shows a way of updating the inverse when small changes are made to the matrix: alglib.net/matrixops/general/invupdate.php – Bitwise Oct 6 '12 at 0:21
Assuming your matrix is real matrix (all elements are real), you are guaranteed by spectral theorem that your matrix $A = QDQ^{T}$ where $D$ is diagonal and $Q$ is orthogonal. Getting this $Q$ is algorithmically Gram-Schimdt process (so rather simple). – s.b Oct 6 '12 at 0:22
@s.b yes, I'm aware; if that helps me invert the submatricies of $\Sigma$ (your $A$) then please tell me how, but as I said in my question I've already tried to figure out how to use the eigen decomposition to do it. – guy Oct 6 '12 at 0:54
> I can tell you that there is a rank one update formula for finding > the inverse of a submatrix of one less dimension, which is to say the > matrix attained after deletion of any desired row and any desired > column. If you think it would be useful I could recite it for you > after I look it up exactly. I'm also interested in inverting this kind of matrices. If you could post the answer i would be glad. Thanks in advance, P – user66612 Mar 13 at 18:48
@user66612: I would be very happy to see it. – Gil Mar 20 at 12:15
show 1 more comment

I imagine the best you can do is only one rank at a time building up or down to attain the desired sub-matrix portion. The simple formula is

$$\mathbf{A}^{-1} = E - \frac{\mathbf{f}\mathbf{g}^T}{h}$$

Where $$\pmatrix{\mathbf A & \mathbf b\\ \mathbf c^T & d}^{-1} = \pmatrix{\mathbf E & \mathbf f\\ \mathbf g^T & h}$$

To derive it we will use Bitwise's helpful comment; the formula $$(A + uv^T)^{-1}=A^{-1}-\frac{(A^{-1}u)(v^TA^{-1})}{1 + v^TA^{-1}u}$$ which shows the rank one aspect of the update relative to the already known inverse. If you use this with sub-matrices in mind, care would have to be given for how to reduce the dimension of $A$, this formula is not directly applicable for that; it does not show how to update a row and a column at once since that is not a rank one update. But it becomes necessary in the derivation.

You have as your knowns the conformably partitioned matrix and its inverse: $$\pmatrix{\mathbf A & \mathbf b\\ \mathbf c^T & d} \pmatrix{\mathbf E & \mathbf f\\ \mathbf g^T & h} = \pmatrix{ \mathbf I & \mathbf 0 \\ \mathbf 0^T & 1}$$ And you want to find the quantity $A^{-1}$. First see that $AE + \mathbf {bg}^T = \mathbf I$ Is the top left element of their multiple. Then:

\begin{align} AE + \mathbf {bg}^T &= \mathbf I \\ AE &= \mathbf I - \mathbf {bg}^T \\ AE\left(\mathbf I - \mathbf {bg}^T\right)^{-1} &= \mathbf{I} \\ \Rightarrow A^{-1} &= E\left(\mathbf I - \mathbf {bg}^T\right)^{-1} \end{align}

The inverse for $\mathbf{A}$ then is the matrix $E(\mathbf I-\mathbf {bg}^T)^{-1}$. The rank one formula may now be used to find the necessary elements. \begin{align} \mathbf{A}^{-1} & = E\left(\mathbf I - \mathbf {bg}^T\right)^{-1} \\ & =E\left(\mathbf I + \mathbf {(-b)g}^T\right)^{-1} \\ & =E\left(\mathbf I - \frac{\mathbf {(-b)g}^T}{1 + \mathbf{g}^T(-\mathbf{b})}\right) \\ & =E\left(\mathbf I + \frac{\mathbf {bg}^T}{1 - \mathbf{g}^T\mathbf{b}}\right) \\ & =E + \frac{(E\mathbf {b})\mathbf{g}^T}{1 - \mathbf{g}^T\mathbf{b}} \\ \end{align}

How about some further reducing: $$\pmatrix{\mathbf E & \mathbf f\\ \mathbf g^T & h}\pmatrix{\mathbf A & \mathbf b\\ \mathbf c^T & d} = \pmatrix{ \mathbf I & \mathbf 0 \\ \mathbf 0^T & 1}$$ \begin{align} E\mathbf{b} + \mathbf{f}d & = 0 \\ E\mathbf{b} &= -\mathbf{f}d \\ \text{and}\quad \mathbf{g}^Tb + hd &= 1 \\ \mathbf{g}^T\mathbf{b} &= 1 - hd \\ 1 - \mathbf{g}^T\mathbf{b} &= hd \\ \Rightarrow \mathbf{A}^{-1} &= E + \frac{-\mathbf{f}d\mathbf{g}^T}{hd} \\ &= E - \frac{\mathbf{f}\mathbf{g}^T}{h} \\ \end{align}

And you have your rank one update of the sub-matrix portion, repeated here for convenience:

$$\mathbf{A}^{-1} = E - \frac{\mathbf{f}\mathbf{g}^T}{h}$$

To go the other direction (adding a row and column) you can use what is called the bordering method as described in this answer

-
 Thanks (+1). The point of looking into the eigen decomposition was that I could perhaps do it up front if it gave made submatrix inversion as trivial as it makes the full matrix inversion, but it doesn't appear to. The logic was that, while it is harder, it would save time in the long run because all the work would be done up front at once. – guy Oct 6 '12 at 1:52 Yes, unfortunately (or fortunately else things might not be as interesting) eigenspace spectrums do not have easy update formulas. – adam W Oct 6 '12 at 2:02 Also, having the spectrum is the same as having a matrix with inverse, that is what the left and right eigenvectors are after all. So this method could be used on those vectors if desired. I think :) – adam W Oct 6 '12 at 2:10 Could you please recite the rank one update formula for finding the inverse? That would be extremely helpful for me. – Gil Mar 20 at 15:00 @Gil That would be the formula already given: $(A + uv^T)^{-1}=A^{-1}-\frac{(A^{-1}u)(v^TA^{-1})}{1 + v^TA^{-1}u}$. If you are referring to my "downdate" to lesser dimension as described in this answer, apply that formula to $E(\mathbf I-\mathbf {bg}^T)^{-1}$. So use $A=I$, $u=-\mathbf{b}$, and $v^\top = \mathbf{g}^\top$. When simplified, it becomes a rank one update of $\mathbf{E}$ – adam W Mar 20 at 17:48
show 1 more comment