4
$\begingroup$

Suppose I have a covariance matrix $A_{n\times n}$ and $A^{-1}_{n\times n}$ is its inverse. Then I randomly exclude the $i$-th row and column, where $1\le i\le n \in \mathbb{N}$, obtaining the new symmetric matrix $A_{n-1\times n-1}$.

Is there any way to calculate the inverse of $A_{n-1\times n-1}$ using the inverse already calculated $A^{-1}_{n\times n}$ to simplify the calculus?

$\endgroup$
5
$\begingroup$

When $i=n,$ write $\mathbb{A}$ in block matrix form

$$\mathbb A = \pmatrix{A & B \\ C & D}$$

where $A$ is the $n-1 \times n-1$ matrix obtained by omitting the last row and column of $\mathbb{A},$ $B = C^\prime$ is the first $n-1$ entries in the last column, and $D = \mathbb{A}_{nn}$ is a nonzero number because $\mathbb{A}$ is an invertible definite symmetric matrix.

Similarly write $$\mathbb{A}^{-1} = \pmatrix{a & b \\ c & d}$$ in block matrix form. We are looking for an efficient formula for $A^{-1}$ in terms of $a,b,c,d.$

By definition, the product of a matrix and its inverse is the $n\times n$ identity. Let's compute it using block matrix operations:

$$\pmatrix{\mathbb{I}_{n-1} & 0 \\ 0 & 1} = \mathbb{I}_n = \mathbb{A}\,\mathbb{A}^{-1} = \pmatrix{Aa + Bc & Ab + Bd \\ Ca + Dc & Cb+Dd}.$$

In the upper left block we almost have the result we would like: it says $Aa + Bc$ is the identity. The trick will be to adjust $a$ to compensate for the $Bc$ term in that block.

In the upper right block we find $0 = Ab + Bd.$ Use this to re-express

$$Bc = (Bd)(d^{-1}c) = ((Ab + Bd) - Ab)(d^{-1}c) = (0 - Ab)(d^{-1}c) = -Abd^{-1}c.$$

Consequently (returning to the upper left block),

$$\mathbb{I}_{n-1} = Aa + Bc = Aa - Abd^{-1}c = A(a - bd^{-1}c),$$

demonstrating that

$$A^{-1} = a - b\,d^{-1}\,c.\tag{*}$$

The computational effort is minimal: the matrix product requires $2(n-1)^2$ multiplications and then $a$ is updated with $(n-1)^2$ subtractions. This $O(n^2)$ performance is the best possible because potentially all entries of $a$ will change and there are $O(n^2)$ of them.


The case for general $i$ is now readily solved by permuting the $i^\text{th}$ row and column into the last positions. The (self-inverse) permutation matrix $\mathbb{P}^{\,i;n}$ given by

$$\mathbb{P}^{\,i;n}_{jk} = \left\{\eqalign{1, & j=k\text{ and } j\notin \{i,n\} \\ 1, & \{j,k\}=\{i,n\} \\ 0 & \text{otherwise}}\right.$$

does the trick via conjugation: the matrix $$\mathbb{P}^{\,i;n}\, \mathbb{A}\, \mathbb{P}^{\,i;n}$$

has the $i^\text{th}$ row and column moved into the last row and column. Thus, apply formula $(*)$ to this permuted version of $\mathbb{A}.$

Note that this solution works for arbitrary invertible square matrices, whether or not they are covariance matrices, provided only that $\mathbb{A}_{i,i}\ne 0.$


As an example, take

$$\mathbb{A} = \pmatrix{2 & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & 2}$$

and $i=2.$ That is, we are given

$$\mathbb{A}^{-1} = \frac{1}{4}\pmatrix{3 & 2 & 1 \\ 2 & 4 & 2 \\ 1 & 2 & 3}$$

and we wish to update it to find the inverse of $\mathbb{A}$ with its second row and column removed; that is, we wish to compute

$$\pmatrix{2 & 0 \\ 0 & 2}^{-1} = \pmatrix{\frac{1}{2} & 0 \\ 0 & \frac{1}{2}}.$$

The permutation matrix is

$$\mathbb{P}^{2;3} = \pmatrix{1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0}.$$

Conjugating gives

$$\mathbb{P}^{2;3}\,\mathbb{A}^{-1}\,\mathbb{P}^{2;3} = \frac{1}{4}\pmatrix{3 & 1 & 2 \\ 1 & 3 & 2 \\ 2 & 2 & 4}.$$

From this we read off

$$a = \frac{1}{4}\pmatrix{3 & 1 \\ 1 & 3},\quad c = b^\prime = \frac{1}{4}\pmatrix{2 & 2},\quad d = \frac{1}{4}\pmatrix{4}=1.$$

Thus formula $(*)$ gives

$$A^{-1} = a - b d^{-1} c = \frac{1}{4}\pmatrix{3 & 1 \\ 1 & 3} - \frac{1}{2}\pmatrix{1 \\ 1}\,(1)^{-1}\,\frac{1}{2}\pmatrix{1 & 1} = \frac{1}{4}\pmatrix{2 & 0 \\ 0 & 2} = \pmatrix{\frac{1}{2} & 0 \\ 0 & \frac{1}{2}},$$

which is correct.


This R code implements the algorithm. It is followed by an example of its use. The permutation is implemented by R's native subscripting function [, which therefore is efficient.

inverse.update <- function(x, i) {
  a <- x[-i,-i, drop=FALSE]
  b <- x[-i,i, drop=FALSE]
  c <- x[i,-i, drop=FALSE]
  d <- x[i,i]
  a - b %*% c / d # For production code, should throw an error when d is 0.
}
#
# Example.
#
A <- matrix(c(2,-1,0, -1,2,-1, 0,-1,2), 3)
A.inv <- solve(A)
i <- 2
(x.1 <- solve(A[-i,-i]))           # The desired result, directly obtained
(x.0 <- inverse.update(A.inv, i))  # The result via an update
$\endgroup$
  • 1
    $\begingroup$ @Sycorax (re your question about developing an intuition for this kind of calculation): I spent a year or so as a math graduate student doing these kinds of calculations for my PhD research :-). I wouldn't recommend that to anybody, but I would say that some practice with manual matrix calculations--although somewhat painful--can be handy. I rarely need to do such calculations for statistical work, maybe once every several years, but the few times I have used them they led to huge optimizations in the code I was writing. $\endgroup$ – whuber Feb 18 at 18:01
  • 1
    $\begingroup$ Oh, you mean we have to put in hours of work to get better at our art? That sounds rather inconvenient. ;-) In all seriousness, thanks for the feedback. I haven't ever worked on a project that requires strenuous optimization of that kind because I only work on "productized" services, but I can certainly see why that work will pay dividends down the road. $\endgroup$ – Sycorax says Reinstate Monica Feb 18 at 18:27
  • 1
    $\begingroup$ @whuber what a great explanation! Your code worked perfectly for me. $\endgroup$ – Ga13 Feb 19 at 12:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.