您现在的位置是：首页 > 文章详情

Sherman-Morrison公式及其应用

日期：2018-05-08点击：941收藏

Sherman-Morrison公式

Sherman-Morrison公式以 Jack Sherman 和 Winifred J. Morrison命名，在线性代数中，是求解逆矩阵的一种方法。本篇博客将介绍该公式及其应用，首先我们来看一下该公式的内容及其证明。
(Sherman-Morrison公式)假设A∈Rn×n为可逆矩阵，u,v∈Rn为列向量，则A+uvT可逆当且仅当1+vTA−1u≠0, 且当A+uvT可逆时，该逆矩阵由以下公式给出：

(A + u v T) - 1 = A - 1 - A - 1 u v T A - 1 1 + v T A - 1 u .

证明：

(⇐)(⇐)当

1+vTA−1u≠01+vTA−1u≠0时，令

X=A+uvT,Y=A−1−A−1uvTA−11+vTA−1uX=A+uvT,Y=A−1−A−1uvTA−11+vTA−1u,则只需证明

XY=YX=IXY=YX=I即可，其中

II为n阶单位矩阵。

X Y = (A + u v T) (A - 1 - A - 1 u v T A - 1 1 + v T A - 1 u) = A A - 1 + u v T A - 1 - A A - 1 u v T A - 1 + u v T A - 1 u v T A - 1 1 + v T A - 1 u = I + u v T A - 1 - u v T A - 1 + u v T A - 1 u v T A - 1 1 + v T A - 1 u = I + u v T A - 1 - u ( 1 + v T A - 1 u ) v T A - 1 1 + v T A - 1 u = I + u v T A - 1 - u v T A - 1 = I

同理，有

YX=IYX=I.因此，当

1+vTA−1u≠01+vTA−1u≠0时，

(A+uvT)−1=A−1−A−1uvTA−11+vTA−1u.(A+uvT)−1=A−1−A−1uvTA−11+vTA−1u.

(⇒)(⇒)当

u=0u=0时，显然有

1+vTA−1u=1≠0.1+vTA−1u=1≠0.当

u≠0u≠0时,用反正法证明该命题成立。假设

A+uvTA+uvT可逆，但

1+vTA−1u=01+vTA−1u=0，则有

(A + u v T) A - 1 u = u + u (v T A - 1 u) = (1 + v T A - 1 u) u = 0.

因为

A+uvTA+uvT可逆，故

A−1A−1u=0,又因为

A−1A−1可逆，故

u=0u=0,此与假设

u≠0u≠0矛盾。因此，当

A+uvTA+uvT可逆时，有

1+vTA−1u≠0.1+vTA−1u≠0.

Sherman-Morrison公式的应用

应用1：A=I时的Sherman-Morrison公式

在Sherman-Morrison公式中，令A=I,则有：I+uvT可逆当且仅当1+vTu≠0, 且当I+uvT可逆时，该逆矩阵由以下公式给出：

(I + u v T) - 1 = I - u v T 1 + v T u .

再令

v=uv=u,则

1+uTu>01+uTu>0, 因此，

I+uuTI+uuT可逆，且

(I + u u T) - 1 = I - u u T 1 + u T u .

应用2：BFGS算法

Sherman-Morrison公式在BFGS算法中的应用，可用来求解BFGS算法中近似Hessian矩阵的逆。本篇博客并不打算给出Sherman-Morrison公式在BFGS算法中的应用，将会再写篇博客介绍BFGS算法，到时再给出该公式的应用，并会在之后补上该博客的链接（因为笔者还没写）。

应用3：循环三对角线性方程组的求解

本篇博客将详细讲述Sherman-Morrison公式在循环三对角线性方程组的求解中的应用。
首先给给出理论知识介绍部分。
对于A∈Rn×n为可逆矩阵，u,v∈Rn为列向量，1+vTA−1u≠0，需要求解方程(A+uvT)x=b.对此，我们可以先求解以下两个方程：

A y = b, A z = u

.
然后令

x=y−vTy1+vTzzx=y−vTy1+vTzz,该解即为原方程的解，验证如下：

(A + u v T) x = (A + u v T) (y - v T y 1 + v T z z) = A y + u v T y - v T y 1 + v T z A z - v T y 1 + v T z u v T z = b + u v T y - v T y u + v T y u v T z 1 + v T z = b + u v T y - ( 1 + v T z ) v T y u 1 + v T z = b + u v T y - u v T y = b

这样将原方程

(A+uvT)x=b(A+uvT)x=b分成两个方程，可以在一定程度上简化原方程。接下来，我们将介绍循环三对角线性方程组的求解。
所谓循环三对角线性方程组，指的是系数矩阵为如下形式：

A = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ b 1 a 2 0 ⋮ 0 c n c 1 b 2 ⋱ ⋮ \dots 0 0 c 2 ⋱ a n - 2 \dots \dots \dots 0 ⋱ b n - 2 a n - 1 0 0 ⋮ 0 c n - 2 b n - 1 a n a 1 0 ⋮ 0 c n - 1 b n ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

循环三对角线性方程组可写成

Ax=dAx=d,其中

d=(d1,d2,...,dn)T.d=(d1,d2,...,dn)T.
对于此方程的求解，我们令

u=(γ,0,0,...,cn)T,v=(1,0,0,...,a1γ)Tu=(γ,0,0,...,cn)T,v=(1,0,0,...,a1γ)T, 且

A=A′+uvTA=A′+uvT,其中

A′A′如下：

A' = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ b 1 - γ a 2 0 ⋮ 00 c 1 b 2 ⋱ ⋮ \dots 0 0 c 2 ⋱ a n - 2 \dots \dots \dots 0 ⋱ b n - 2 a n - 1 0 0 ⋮ 0 c n - 2 b n - 1 a n 00 ⋮ 0 c n - 1 b n - a 1 c n γ ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

A′A′为三对角矩阵。根据以上的理论知识，我们只需要求解以下两个方程

A' y = d, A' z = u ，

然后，就能根据

y,zy,z求出

xx.而以上两个方程为三对角线性方程组，可以用追赶法（或Thomas法）求解，具体算法可以参考博客：三对角线性方程组(tridiagonal systems of equations)的求解。
综上，我们利用Sherman-Morrison公式的思想，可以将循环三对角线性方程组转化为三对角线性方程组求解。我们将会在下面给出该算法的Python语言实现。

Python实现

我们要解的循环三对角线性方程组如下：

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 4100314100014100014120014 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ x 1 x 2 x 3 x 4 x 5 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ 76668 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥

用Python实现解该方程的Python完整代码如下：

# use Sherman-Morrison Formula and Thomas Method to solve cyclic tridiagonal linear equation import numpy as np # Thomas Method for soling tridiagonal linear equation Ax=d # parameter: a,b,c,d are list-like of same length # b: main diagonal of matrix A # a: main diagonal below of matrix A # c: main diagonal upper of matrix A # d: Ax=d # return: x(type=list), the solution of Ax=d def TDMA(a,b,c,d): try: n = len(d) # order of tridiagonal square matrix # use a,b,c to create matrix A, which is not necessary in the algorithm A = np.array([[0]*n]*n, dtype='float64') for i in range(n): A[i,i] = b[i] if i > 0: A[i, i-1] = a[i] if i < n-1: A[i, i+1] = c[i] # new list of modified coefficients c_1 = [0]*n d_1 = [0]*n for i in range(n): if not i: c_1[i] = c[i]/b[i] d_1[i] = d[i] / b[i] else: c_1[i] = c[i]/(b[i]-c_1[i-1]*a[i]) d_1[i] = (d[i]-d_1[i-1]*a[i])/(b[i]-c_1[i-1] * a[i]) # x: solution of Ax=d x = [0]*n for i in range(n-1, -1, -1): if i == n-1: x[i] = d_1[i] else: x[i] = d_1[i]-c_1[i]*x[i+1] x = [round(_, 4) for _ in x] return x except Exception as e: return e # Sherman-Morrison Fomula for soling cyclic tridiagonal linear equation Ax=d # parameter: a,b,c,d are list-like of same length # b: main diagonal of matrix A # a: main diagonal below of matrix A # c: main diagonal upper of matrix A # d: Ax=d # return: x(type=list), the solution of Ax=d def Cyclic_Tridiagnoal_Linear_Equation(a,b,c,d): try: # use a,b,c to create cyclic tridiagonal matrix A n = len(d) A = np.array([[0] * n] * n, dtype='float64') for i in range(n): A[i, i] = b[i] if i > 0: A[i, i - 1] = a[i] if i < n - 1: A[i, i + 1] = c[i] A[0, n - 1] = a[0] A[n - 1, 0] = c[n - 1] gamma = 1 # gamma can be set freely u = [gamma] + [0] * (n - 2) + [c[n - 1]] v = [1] + [0] * (n - 2) + [a[0] / gamma] # modify the coefficient to form A' b[0] -= gamma b[n - 1] -= a[0] * c[n - 1] / gamma a[0] = 0 c[n - 1] = 0 # solve A'y=d, A'z=u by using Thomas Method y = np.array(TDMA(a, b, c, d)) z = np.array(TDMA(a, b, c, u)) # use y and z to calculate x # x = y-(v·y)/(1+v·z) *z # x is the solution of Ax=d x = y - (np.dot(np.array(v), y)) / (1 + np.dot(np.array(v), z)) * z x = [round(_, 3) for _ in x] return x except Exception as e: return e def main(): ''' equation: A = [[4,1,0,0,2], [1,4,1,0,0], [0,1,4,1,0], [0,0,1,4,1], [3,0,0,1,4]] d = [7,6,6,6,8] solution x should be [1,1,1,1,1] ''' a = [2, 1, 1, 1, 1] b = [4, 4, 4, 4, 4] c = [1, 1, 1, 1, 3] d = [7, 6, 6, 6, 8] x = Cyclic_Tridiagnoal_Linear_Equation(a,b,c,d) print('The solution is %s'%x) main()