Copyright (C) 2020 Andreas Kloeckner
import numpy as np
import scipy.linalg as la
Let's set up some matrices and data for the rank-one modification:
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:
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:
def solveA(b):
return la.lu_solve((LU, piv), b)
la.norm(np.dot(A, solveA(b)) - b)
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:
xhat = la.solve(Ahat, b)
Next, apply Sherman-Morrison to find xhat2
:
xhat2 = solveA(b) - solveA(u)*np.dot(v, solveA(b))/(1+np.dot(v, solveA(u)))
la.norm(xhat - xhat2)
5.010784718882614e-15