摘要：本篇博客將介紹該公式及其應(yīng)用，首先我們來看一下該公式的內(nèi)容及其證明。應(yīng)用循環(huán)三對(duì)角線性方程組的求解本篇博客將詳細(xì)講述公式在循環(huán)三對(duì)角線性方程組的求解中的應(yīng)用。

??Sherman-Morrison公式以 Jack Sherman 和 Winifred J. Morrison命名，在線性代數(shù)中，是求解逆矩陣的一種方法。本篇博客將介紹該公式及其應(yīng)用，首先我們來看一下該公式的內(nèi)容及其證明。
??(Sherman-Morrison公式)假設(shè)$Ainmathbb{R}^{n imes n}$為可逆矩陣，$u,vinmathbb{R}^{n}$為列向量，則$A+uv^{T}$可逆當(dāng)且僅當(dāng)$1+v^{T}A^{-1}u eq 0$, 且當(dāng)$A+uv^{T}$可逆時(shí)，該逆矩陣由以下公式給出：
$$(A+uv^{T})^{-1}=A^{-1}-{A^{-1}uv^{T}A^{-1} over 1+v^{T}A^{-1}u}.$$
證明：
$(Leftarrow)$當(dāng)$1+v^{T}A^{-1}u eq 0$時(shí)，令$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階單位矩陣。

$$ {displaystyle {egin{aligned} XY&=(A+uv^{T})left(A^{-1}-{A^{-1}uv^{T}A^{-1} over 1+v^{T}A^{-1}u} ight) &=AA^{-1}+uv^{T}A^{-1}-{AA^{-1}uv^{T}A^{-1}+uv^{T}A^{-1}uv^{T}A^{-1} over 1+v^{T}A^{-1}u} &=I+uv^{T}A^{-1}-{uv^{T}A^{-1}+uv^{T}A^{-1}uv^{T}A^{-1} over 1+v^{T}A^{-1}u} &=I+uv^{T}A^{-1}-{u(1+v^{T}A^{-1}u)v^{T}A^{-1} over 1+v^{T}A^{-1}u} &=I+uv^{T}A^{-1}-uv^{T}A^{-1} &=Iend{aligned}}} $$

同理，有$YX=I$.因此，當(dāng)$1+v^{T}A^{-1}u eq 0$時(shí)，$(A+uv^{T})^{-1}=A^{-1}-{A^{-1}uv^{T}A^{-1} over 1+v^{T}A^{-1}u}.$
$(Rightarrow)$當(dāng)$u=0$時(shí)，顯然有$1+v^{T}A^{-1}u=1 eq 0.$當(dāng)$u eq0$時(shí),用反正法證明該命題成立。假設(shè)$A+uv^{T}$可逆，但$1+v^{T}A^{-1}u = 0$，則有
$$(A+uv^{T})A^{-1}u=u+u(v^{T}A^{-1}u)=(1+v^{T}A^{-1}u)u=0.$$
因?yàn)?A+uv^{T}$可逆，故$A^{-1}$u=0,又因?yàn)?A^{-1}$可逆，故$u=0$,此與假設(shè)$u eq 0$矛盾。因此，當(dāng)$A+uv^{T}$可逆時(shí)，有$1+v^{T}A^{-1}u eq 0.$

??在Sherman-Morrison公式中，令$A=I$,則有：$I+uv^{T}$可逆當(dāng)且僅當(dāng)$1+v^{T}u eq 0$, 且當(dāng)$I+uv^{T}$可逆時(shí)，該逆矩陣由以下公式給出：
$$(I+uv^{T})^{-1}=I-{uv^{T} over 1+v^{T}u}.$$
??再令$v=u$,則$1+u^{T}u > 0$, 因此，$I+uu^{T}$可逆，且
$$(I+uu^{T})^{-1}=I-{uu^{T} over 1+u^{T}u}.$$

??Sherman-Morrison公式在BFGS算法中的應(yīng)用，可用來求解BFGS算法中近似Hessian矩陣的逆。本篇博客并不打算給出Sherman-Morrison公式在BFGS算法中的應(yīng)用，將會(huì)再寫篇博客介紹BFGS算法，到時(shí)再給出該公式的應(yīng)用，并會(huì)在之后補(bǔ)上該博客的鏈接（因?yàn)楣P者還沒寫）。

??本篇博客將詳細(xì)講述Sherman-Morrison公式在循環(huán)三對(duì)角線性方程組的求解中的應(yīng)用。
??首先給給出理論知識(shí)介紹部分。
??對(duì)于$Ainmathbb{R}^{n imes n}$為可逆矩陣，$u,vinmathbb{R}^{n}$為列向量，$1+v^{T}A^{-1}u eq 0$，需要求解方程$(A+uv^{T})x=b.$對(duì)此，我們可以先求解以下兩個(gè)方程：
$$Ay=b,qquad Az=u$$.
然后令$x=y-frac{v^{T}y}{1+v^{T}z}z$,該解即為原方程的解，驗(yàn)證如下：

$$ {displaystyle {egin{aligned} (A+uv^{T})x&=(A+uv^{T})(y-frac{v^{T}y}{1+v^{T}z}z) &=Ay+uv^{T}y-frac{v^{T}y}{1+v^{T}z}Az-frac{v^{T}y}{1+v^{T}z}uv^{T}z &=b+uv^{T}y-frac{v^{T}yu+v^{T}yuv^{T}z}{1+v^{T}z} &=b+uv^{T}y-frac{(1+v^{T}z)v^{T}yu}{1+v^{T}z} &=b+uv^{T}y-uv^{T}y &=bend{aligned}}} $$

??這樣將原方程$ (A+uv^{T})x=b$分成兩個(gè)方程，可以在一定程度上簡(jiǎn)化原方程。接下來，我們將介紹循環(huán)三對(duì)角線性方程組的求解。
??所謂循環(huán)三對(duì)角線性方程組，指的是系數(shù)矩陣為如下形式：

$$ A=egin{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_nend{bmatrix} $$

循環(huán)三對(duì)角線性方程組可寫成$Ax=d$,其中$d=(d_{1},d_2,...,d_{n})^{T}.$
??對(duì)于此方程的求解，我們令$u=(gamma, 0,0,...,c_{n})^{T}, v=(1,0,0,...,frac{a_1}{gamma})^{T}$, 且$A=A^{"}+uv^{T}$,其中$A^{"}$如下：

$$ A^{"}=egin{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} $$

$A^{"}$為三對(duì)角矩陣。根據(jù)以上的理論知識(shí)，我們只需要求解以下兩個(gè)方程
$$A^{"}y=d,qquad A^{"}z=u，$$
然后，就能根據(jù)$y,z$求出$x$.而以上兩個(gè)方程為三對(duì)角線性方程組，可以用追趕法（或Thomas法）求解，具體算法可以參考博客：三對(duì)角線性方程組(tridiagonal systems of equations)的求解 。
??綜上，我們利用Sherman-Morrison公式的思想，可以將循環(huán)三對(duì)角線
性方程組轉(zhuǎn)化為三對(duì)角線性方程組求解。我們將會(huì)在下面給出該算法的Python語言實(shí)現(xiàn)。

??我們要解的循環(huán)三對(duì)角線性方程組如下：

$$ {egin{bmatrix} 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} end{bmatrix}} {egin{bmatrix}{x_{1}}{x_{2}}{x_{3}}{x_{4}} {x_{5}}end{bmatrix}}={egin{bmatrix}{76 668}end{bmatrix}} $$

??用Python實(shí)現(xiàn)解該方程的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()

輸出結(jié)果如下：

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

https://en.wikipedia.org/wiki...

http://wwwmayr.in.tum.de/konf...

https://scicomp.stackexchange...

https://blog.csdn.net/jclian9...