Sherman-Morrison公式及其应用
Sherman-Morrison公式
Sherman-Morrison公式以 Jack Sherman 和 Winifred J. Morrison命名，在线性代数中，是求解逆矩阵的一种方法。本篇博客将介绍该公式及其应用，首先我们来看一下该公式的内容及其证明。
(Sherman-Morrison公式)假设\(A\in\mathbb{R}^{n\times n}\)为可逆矩阵，\(u,v\in\mathbb{R}^{n}\)为列向量，则\(A+uv^{T}\)可逆当且仅当\(1+v^{T}A^{-1}u\neq 0\), 且当\(A+uv^{T}\)可逆时，该逆矩阵由以下公式给出：
证明：
\((\Leftarrow)\)当\(1+v^{T}A^{-1}u\neq 0\)时，令\(X=A+uv^{T}, Y=A^{-1}-{A^{-1}uv^{T}A^{-1} \over 1+v^{T}A^{-1}u}\),则只需证明\(XY=YX=I\)即可，其中\(I\)为n阶单位矩阵。
同理，有\(YX=I\).因此，当\(1+v^{T}A^{-1}u\neq 0\)时，\((A+uv^{T})^{-1}=A^{-1}-{A^{-1}uv^{T}A^{-1} \over 1+v^{T}A^{-1}u}.\)
\((\Rightarrow)\)当\(u=0\)时，显然有\(1+v^{T}A^{-1}u=1\neq 0.\)当\(u\neq0\)时,用反正法证明该命题成立。假设\(A+uv^{T}\)可逆，但\(1+v^{T}A^{-1}u = 0\)，则有
因为\(A+uv^{T}\)可逆，故\(A^{-1}\)u=0,又因为\(A^{-1}\)可逆，故\(u=0\),此与假设\(u\neq 0\)矛盾。因此，当\(A+uv^{T}\)可逆时，有\(1+v^{T}A^{-1}u \neq 0.\)
Sherman-Morrison公式的应用
应用1：\(A=I\)时的Sherman-Morrison公式
在Sherman-Morrison公式中，令\(A=I\),则有：\(I+uv^{T}\)可逆当且仅当\(1+v^{T}u\neq 0\), 且当\(I+uv^{T}\)可逆时，该逆矩阵由以下公式给出：
再令\(v=u\),则\(1+u^{T}u > 0\), 因此，\(I+uu^{T}\)可逆，且
应用2：BFGS算法
Sherman-Morrison公式在BFGS算法中的应用，可用来求解BFGS算法中近似Hessian矩阵的逆。本篇博客并不打算给出Sherman-Morrison公式在BFGS算法中的应用，将会再写篇博客介绍BFGS算法，到时再给出该公式的应用，并会在之后补上该博客的链接（因为笔者还没写）。
应用3：循环三对角线性方程组的求解
本篇博客将详细讲述Sherman-Morrison公式在循环三对角线性方程组的求解中的应用。
首先给给出理论知识介绍部分。
对于\(A\in\mathbb{R}^{n\times n}\)为可逆矩阵，\(u,v\in\mathbb{R}^{n}\)为列向量，\(1+v^{T}A^{-1}u\neq 0\)，需要求解方程\((A+uv^{T})x=b.\)对此，我们可以先求解以下两个方程：
{\displaystyle
{\begin{aligned}
(A+uv{T})x&=(A+uv)(y-\frac{v{T}y}{1+vz}z)\
&=Ay+uv{T}y-\frac{vy}{1+v{T}z}Az-\frac{vy}{1+v{T}z}uvz\
&=b+uv{T}y-\frac{vyu+v{T}yuvz}{1+v^{T}z}\
&=b+uv{T}y-\frac{(1+vz)v{T}yu}{1+vz}\
&=b+uv{T}y-uvy\
&=b\end{aligned}}}
A=\begin{bmatrix}
b_1&c_1&0&\cdots&0&a_1\
a_2&b_2&c_2&0&\vdots&0\
0&\ddots&\ddots&\ddots&0&\vdots\
\vdots&\vdots&a_{n-2}&b_{n-2}&c_{n-2}&0\
0&\cdots&\cdots&a_{n-1}&b_{n-1}&c_{n-1}\
c_n&0&\cdots&0&a_n&b_n\end{bmatrix}
A^{'}=\begin{bmatrix}
b_1-\gamma&c_1&0&\cdots&0&0\
a_2&b_2&c_2&0&\vdots&0\
0&\ddots&\ddots&\ddots&0&\vdots\
\vdots&\vdots&a_{n-2}&b_{n-2}&c_{n-2}&0\
0&\cdots&\cdots&a_{n-1}&b_{n-1}&c_{n-1}\
0&0&\cdots&0&a_n&b_n-\frac{a_1c_n}{\gamma}\end{bmatrix}
然后，就能根据\(y,z\)求出\(x\).而以上两个方程为三对角线性方程组，可以用追赶法（或Thomas法）求解，具体算法可以参考博客：三对角线性方程组(tridiagonal systems of equations)的求解 。
综上，我们利用Sherman-Morrison公式的思想，可以将循环三对角线性方程组转化为三对角线性方程组求解。我们将会在下面给出该算法的Python语言实现。
Python实现
我们要解的循环三对角线性方程组如下：
用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()
输出结果如下：
The solution is [1.0, 1.0, 1.0, 1.0, 1.0]
参考文献
- https://en.wikipedia.org/wiki/Sherman–Morrison_formula
- http://wwwmayr.in.tum.de/konferenzen/Jass09/courses/2/Soldatenko_paper.pdf
- https://scicomp.stackexchange.com/questions/10137/solving-system-of-linear-equations-with-cyclic-tridiagonal-matrix
- https://blog.csdn.net/jclian91/article/details/80251244
注意：本人现已开通两个微信公众号： 用Python做数学（微信号为：python_math）以及轻松学会Python爬虫（微信号为：easy_web_scrape）， 欢迎大家关注哦~~
