Monday, February 17, 2014

The Sherman-Morrison Formula

In a recent post, I mentioned how the upfront cost of LU decomposition pays for itself when solving a system of linear equations \[Ax=b_i, \quad \quad i = 1, 2, ... N.\] Now suppose the (presumably large) matrix \(A\) was modified slightly to form a new matrix \(A'\). Let us assume that we had a LU factorization of the original unmodified matrix, which would have allowed us to solve \(Ax=b\) very efficiently.

If we now wanted to solve \(A'x=b\), can we somehow avoid having to perform a brand new LU factorization on the modified matrix \(A'\)?

Recall that the factorization itself doesn't come very cheap.

The Sherman-Morrison formula (SMF) says "sure can!" as long as the modification can be expressed as an outer-product of two vectors. That is, \[A' = A + uv^T.\] As an example, suppose that the element in the \((j,k)\) position \(a_{jk}\) was changed to \(a'_{kj}\). This modification can be expressed as,\[A' = A + (a'_{jk} - a_{jk}) \mathbf{e}_{j} \mathbf{e}_{k}^T,\] where \(\mathbf{e}_{j}\) is a vector with "1" in the row \(j\) and zeroes everywhere else.

The Sherman-Morrison formula says \[(A')^{-1} = (A+uv^T)^{-1} = A^{-1} - {A^{-1}uv^T A^{-1} \over 1 + v^T A^{-1}u}.\] The formula can be readily verified, as shown on wikipedia.

Note: The scalar analog to question the SMF seeks to address is: given \(1/a\), how can you find \(1/(a+b)\)?

The SMF answer would be: \[1/a \times \left(1- {b/a \over 1 + b/a} \right).\] Returning to the original problem \(A'x = b \implies x = (A')^{-1} b\). Thus \[x = A^{-1}b - {A^{-1}uv^T A^{-1}b \over 1 + v^T A^{-1}u}.\] We can set \(y = A^{-1}b\), and \(z=A^{-1}u\). Since a LU factorization is available for \(A\), we can solve for \(y\) and \(z\) in \(\mathcal{O}(m^2)\) operations, where \(m\) is the size of the square matrix \(A\).

This gives us \[x = y - \left({v^T y \over 1 + v^T z}\right) z.\] In general, this can be performed in \(\mathcal{O}(3 m^2)\) operations, as opposed to to \(\mathcal{O}(m^3)\) for a fresh factorization.

No comments: