Sherman-Morrison

In [1]:
import numpy as np
import scipy.linalg as la

Let's set up some matrices and data for the rank-one modification:

In [6]:
n = 5
A = np.random.randn(n, n)
u = np.random.randn(n)
v = np.random.randn(n)

b = np.random.randn(n)

Ahat = A + np.outer(u, v)

Let's start by computing the "base" factorization.

We'll use lu_factor from scipy, which stuffs both L and U into a single matrix (why can it do that?) and also returns pivoting information:

In [7]:
LU, piv = la.lu_factor(A)
print(LU)
print(piv)
[[ 1.08364401 -2.18363086 -1.07875426 -0.17119634  0.40477559]
 [ 0.94001453  2.11377723  0.49515679  0.04580586 -0.48416682]
 [ 0.11060844 -0.01756305 -1.72250622 -0.49623917 -1.42697035]
 [-0.2463909  -0.41778398  0.78182516  1.19130607  1.03773055]
 [-0.06177908 -0.79834934 -0.1559638   0.71449418 -0.04804558]]
[3 3 2 4 4]

Next, we set up a subroutine to solve using that factorization and check that it works:

In [9]:
def solveA(b):
    return la.lu_solve((LU, piv), b)

la.norm(np.dot(A, solveA(b)) - b)
Out[9]:
3.074360713233696e-15

As a last step, we try the Sherman-Morrison formula:

$$(A+uv^T)^{-1} = A^{-1} - {A^{-1}uv^T A^{-1} \over 1 + v^T A^{-1}u}$$

To see that we got the right answer, we first compute the right solution of the modified system:

In [11]:
xhat = la.solve(Ahat, b)

Next, apply Sherman-Morrison to find xhat2:

In [13]:
xhat2 = solveA(b) - solveA(u)*np.dot(v, solveA(b))/(1+np.dot(v, solveA(u)))
In [14]:
la.norm(xhat - xhat2)
Out[14]:
5.010784718882614e-15
  • What's the cost of the Sherman-Morrison procedure?