The Sherman-Morrison-Woodbury Formula

The Sherman-Morrison-Woodbury Formula#

The Sherman-Morrison-Woodbury formula provides an efficient and powerful method for calculating the inverse of a matrix that has been updated by a low-rank matrix. Its primary advantage is that it allows you to find the new inverse without having to perform a full, computationally expensive matrix inversion from scratch, provided you already know the inverse of the original matrix. 💡

The Sherman-Morrison Formula (Rank-One Update)#

The simplest case of this identity is the Sherman-Morrison formula, which applies to a rank-one update. A rank-one update means we are perturbing an invertible matrix \(A\) by adding an outer product of two vectors, \(uv^T\).

If \(A\) is an invertible matrix, \(u, v\) are vectors, and \(1 + v^T A^{-1} u \neq 0\), then the inverse of the updated matrix \(A + uv^T\) is given by:

\[(A + uv^T)^{-1} = A^{-1} - \frac{A^{-1}uv^T A^{-1}}{1 + v^T A^{-1} u}\]

The condition on the scalar is not a technicality: \(A + uv^T\) is invertible precisely when it holds. This follows from the matrix determinant lemma,

\[\det(A + uv^T) = (1 + v^T A^{-1} u)\,\det(A),\]

so the determinant vanishes exactly when the denominator does.

The remarkable efficiency of this formula comes from the denominator, \(1 + v^T A^{-1} u\). Since \(v^T\), \(A^{-1}\), and \(u\) are a row vector, a matrix, and a column vector, respectively, their product is a scalar. This means we avoid a complex matrix inversion and instead only need to perform a simple scalar division.

Application: Solving a Linear System#

In practice, we often use the formula to efficiently solve a linear system of the form \((A + uv^T)x = b\), especially when \(A^{-1}\) is known or systems with \(A\) are easy to solve.

To compute \(x = (A + uv^T)^{-1} b\), we can apply the formula and group the operations intelligently:

\[x = \left( A^{-1} - \frac{A^{-1}u(v^T A^{-1})}{1 + v^T A^{-1} u} \right) b = A^{-1}b - A^{-1}u \left( \frac{v^T A^{-1} b}{1 + v^T A^{-1} u} \right)\]

This leads to the following efficient computational steps:

  1. Solve for an intermediate vector \(y = A^{-1}b\).

  2. Solve for another intermediate vector \(z = A^{-1}u\).

  3. Compute the scalar value in the denominator: \(\beta = 1 + v^T z\).

  4. Compute the scalar value in the numerator’s right part: \(\alpha = v^T y\).

  5. Combine these results to find the final solution: \(x = y - z (\frac{\alpha}{\beta})\).

This procedure replaces a full \(O(n^3)\) matrix inversion with two solves against \(A\) and a handful of vector operations. The saving depends on already having a factorization of \(A\): with an LU factorization in hand, each of steps 1 and 2 is a pair of \(O(n^2)\) triangular solves. Computing that factorization from scratch would itself cost \(O(n^3)\).

Numerical Accuracy#

The formula is exact, but computers round numbers during the calculation, and these small errors can sometimes lead to an inaccurate answer. A large correction \(z(\alpha/\beta)\) is not a problem by itself. The concern is when it nearly cancels \(y\): errors already present in the two vectors can then be large compared with their difference. Similarly, when \(\beta = 1 + v^T z\) is close to zero, a small rounding error in \(v^T z\) can be large compared with \(\beta\), affecting the division by \(\beta\).

We will discuss these accuracy questions in the next chapter, Solving Linear Systems.

The Woodbury Matrix Identity (Generalization)#

The Sherman-Morrison formula is a special case of the more general Woodbury matrix identity, which handles updates of a higher rank.

Theorem 1 (The Woodbury Matrix Identity)

Let \(A\) be an invertible \(n \times n\) matrix, \(U\) be an \(n \times k\) matrix, \(C\) be an invertible \(k \times k\) matrix, and \(V\) be a \(k \times n\) matrix. The Woodbury matrix identity states:

\[(A + UCV)^{-1} = A^{-1} - A^{-1}U(C^{-1} + VA^{-1}U)^{-1}VA^{-1}\]

This identity holds provided that the matrix \((C^{-1} + VA^{-1}U)\) is invertible.

The key advantage is that instead of inverting the large \(n \times n\) matrix on the left, we only need to invert the much smaller \(k \times k\) matrix \((C^{-1} + VA^{-1}U)\) on the right. This is extremely beneficial when \(k\) is much smaller than \(n\).

Taking \(k = 1\) recovers the Sherman-Morrison formula. Set \(U = u\), \(C = [1]\) and \(V = v^T\). Then \(C^{-1} + VA^{-1}U\) is the \(1 \times 1\) matrix \(1 + v^T A^{-1} u\), inverting it is a scalar division, and the identity reads

\[(A + uv^T)^{-1} = A^{-1} - \frac{A^{-1}u v^T A^{-1}}{1 + v^T A^{-1} u}.\]

Proof. Woodbury Matrix Identity

The most straightforward way to prove the identity is to show that the product of the original matrix \((A + UCV)\) and its proposed inverse from the formula results in the identity matrix, \(I\).

Let’s multiply \((A + UCV)\) by the proposed inverse:

\[(A + UCV) \left( A^{-1} - A^{-1}U(C^{-1} + VA^{-1}U)^{-1}VA^{-1} \right)\]

We can expand this product using the distributive property:

\[\begin{split} \begin{aligned} = {} & A(A^{-1}) - A(A^{-1}U(C^{-1} + VA^{-1}U)^{-1}VA^{-1}) \\ & {} + UCV(A^{-1}) - UCV(A^{-1}U(C^{-1} + VA^{-1}U)^{-1}VA^{-1}). \end{aligned} \end{split}\]

Now, let’s simplify each term:

  • \(A A^{-1} = I\)

  • Again using \(A A^{-1} = I\),

    \[\begin{split} \begin{aligned} & A A^{-1} U (C^{-1} + VA^{-1}U)^{-1}VA^{-1} \\ & \qquad = U (C^{-1} + VA^{-1}U)^{-1}VA^{-1}. \end{aligned} \end{split}\]

This simplifies the expression to:

\[\begin{split} \begin{aligned} = {} & I - U(C^{-1} + VA^{-1}U)^{-1}VA^{-1} \\ & {} + UCVA^{-1} - UCVA^{-1}U(C^{-1} + VA^{-1}U)^{-1}VA^{-1}. \end{aligned} \end{split}\]

Now, let’s factor out the common term \(U\) on the left side of the last three terms:

\[\begin{split} = I + U \left( \begin{aligned} & - (C^{-1} + VA^{-1}U)^{-1}VA^{-1} + CVA^{-1} \\ & {} - CVA^{-1}U(C^{-1} + VA^{-1}U)^{-1}VA^{-1} \end{aligned} \right). \end{split}\]

Let’s focus on the expression inside the parentheses. We can factor out the term \((C^{-1} + VA^{-1}U)^{-1}VA^{-1}\) from the first and third parts:

\[\begin{split} = I + U \left( \begin{aligned} & CVA^{-1} \\ & {} - (I + CVA^{-1}U)(C^{-1} + VA^{-1}U)^{-1}VA^{-1} \end{aligned} \right). \end{split}\]

The key step is to simplify the term \((I + CVA^{-1}U)\). We can factor out \(C\) on the left:

\[I + CVA^{-1}U = C(C^{-1} + VA^{-1}U)\]

Now, we substitute this back into our main expression:

\[\begin{split} = I + U \left( \begin{aligned} & CVA^{-1} \\ & {} - C(C^{-1} + VA^{-1}U)(C^{-1} + VA^{-1}U)^{-1}VA^{-1} \end{aligned} \right). \end{split}\]

The matrix \((C^{-1} + VA^{-1}U)\) and its inverse \((C^{-1} + VA^{-1}U)^{-1}\) cancel out to become the identity matrix:

\[= I + U \left( CVA^{-1} - C(I)VA^{-1} \right)\]

This simplifies to:

\[= I + U \left( CVA^{-1} - CVA^{-1} \right) = I + U(0) = I\]

The product is the identity matrix. Both factors are square, so a one-sided inverse is automatically two-sided, and the proposed expression is therefore \((A + UCV)^{-1}\).