# Updating the covariance matrix after deleting the i-th column and row

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?

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}.$$

$$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

• @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. – whuber Feb 18 at 18:01
• 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. – Sycorax says Reinstate Monica Feb 18 at 18:27
• @whuber what a great explanation! Your code worked perfectly for me. – Ga13 Feb 19 at 12:27