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?
 A: 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

