Sherman-Morrison公式及其應用

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≠0時，令X=A+uvT,Y=A−1−A−1uvTA−11+vTA−1u,則只需證明XY=YX=I即可，其中I為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=I.因此，當1+vTA−1u≠0時，(A+uvT)−1=A−1−A−1uvTA−11+vTA−1u.
(⇒)當u=0時，顯然有1+vTA−1u=1≠0.當u≠0時,用反正法證明該命題成立。假設A+uvT可逆，但1+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+uvT可逆，故A−1u=0,又因為A−1可逆，故u=0,此與假設u≠0矛盾。因此，當A+uvT可逆時，有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=u,則1+uTu>0, 因此，I+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+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 = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 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=d,其中d=(d1,d2,...,dn)T.
對於此方程的求解，我們令u=(γ,0,0,...,cn)T,v=(1,0,0,...,a1γ)T, 且A=A′+uvT,其中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' y = d, A' z = u ，

然後，就能根據y,z求出x.而以上兩個方程為三對角線性方程組，可以用追趕法（或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()

輸出結果如下：

The solution is [1.0, 1.0, 1.0, 1.0, 1.0]

Sherman-Morrison公式及其應用