Sherman-Morrison公式及其应用

Sherman-Morrison公式

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

证明：
(⇐)(⇐)当 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阶单位矩阵。

同理，有 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+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=IA=I时的Sherman-Morrison公式

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

  再令 v=uv=u,则 1+uTu>01+uTu>0, 因此， I+uuTI+uuT可逆，且

应用2：BFGS算法

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

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

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

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

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

循环三对角线性方程组可写成 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′A′为三对角矩阵。根据以上的理论知识，我们只需要求解以下两个方程

然后，就能根据 y,zy,z求出 xx.而以上两个方程为三对角线性方程组，可以用追赶法（或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]

参考文献

1. https://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula
2. http://wwwmayr.in.tum.de/konferenzen/Jass09/courses/2/Soldatenko_paper.pdf
3. https://scicomp.stackexchange.com/questions/10137/solving-system-of-linear-equations-with-cyclic-tridiagonal-matrix
4. https://blog.csdn.net/jclian91/article/details/80251244