Sherman-Morrison式とその応用
Sherman-Morrison式
Sherman‐Morrisson式はJack ShermanとWinifred J.Morrisonによって命名され,線形代数では逆行列を解く方法である.このブログでは、この式とその応用について説明します.まず、その式の内容と証明を見てみましょう. (Sherman-Morrison式)仮定A∈Rn×n A ∈ R n × nが可逆行列であり、u,v∈Rnu,v∈Rnが列ベクトルである場合、A+uvTA+uvTは可逆的であり、1+vTA−1 u≠0 1+vTA−1 u≠0のみであり、A+uvTA+uvTが可逆的である場合、この逆行列は以下の式によって与えられる.
(A+uvT)−1=A−1−A−1uvTA−11+vTA−1u. ( A + u v T ) − 1 = A − 1 − A − 1 u v T A − 1 1 + v T A − 1 u .
証明:
(⇐)当
1+vTA−1 u≠0 1+vT A−1 u≠0の場合
X=A+uvT,Y=A−1−A−1 uvTA−11+vTA−1 uX=A+uvT,Y=A−1−A−1 uvT−A−1 uvT−1+vT A−1 uであれば証明するだけである
XY=YX=IX Y=Y X=Iであればよく、そのうち
Iはn次単位行列である.
XY=(A+uvT)(A−1−A−1uvTA−11+vTA−1u)=AA−1+uvTA−1−AA−1uvTA−1+uvTA−1uvTA−11+vTA−1u=I+uvTA−1−uvTA−1+uvTA−1uvTA−11+vTA−1u=I+uvTA−1−u(1+vTA−1u)vTA−11+vTA−1u=I+uvTA−1−uvTA−1=I 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 Y X = I .したがって
1+vTA−1 u≠0 1+vT A−1 u≠0の場合、
(A+uvT)−1=A−1−A−1uvTA−11+vTA−1u. ( A + u v T ) − 1 = A − 1 − A − 1 u v T A − 1 1 + v T A − 1 u .
(⇒)(⇒)当
u=0 u=0の場合、明らかに
1+vTA−1u=1≠0. 1 + v T A − 1 u = 1 ≠ 0. 当
u≠0 u≠0の場合、この命題が成立したことをどうせ法で証明する.仮定
A+uvT A+uvTは可逆的であるが
1+vTA−1 u=01+vT A−1 u=0であると、
(A+uvT)A−1u=u+u(vTA−1u)=(1+vTA−1u)u=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+u v Tは可逆的であるため、
A−1 A−1 u=0である、また、
A−1 A−1は可逆的であるため
u=0 u=0と仮定
u≠0 u≠0矛盾.したがって
A+uvT A+uvT可逆の場合、
1+vTA−1u≠0. 1 + v T A − 1 u ≠ 0.
Sherman-Morrison式の応用
1:A=IA=Iを適用する場合のSherman-Morrison式
Sherman-Morrison式において、A=I A=Iとすると、I+uvTI+uvTが可逆的であり、1+vTu≠0 1+vTu≠0のみであり、I+uvTI+uvTが可逆的である場合、この逆行列は以下の式によって与えられる.
(I+uvT)−1=I−uvT1+vTu. ( I + u v T ) − 1 = I − u v T 1 + v T u .
再命令を下す.
v=u v=uの場合
1+uTu>0 1+uT u>0なので、
I+uuTI+uuTは可逆的であり、
(I+uuT)−1=I−uuT1+uTu. ( I + u u T ) − 1 = I − u u T 1 + u T u .
応用2：BFGSアルゴリズム
BFGSアルゴリズムにおける近似Hessian行列の逆を解くために用いられるBFGSアルゴリズムにおけるSherman‐Morrison式の応用.このブログはSherman-Morrison公式のBFGSアルゴリズムへの応用を与えるつもりはなく、BFGSアルゴリズムを紹介するブログを書き、その時にその公式の応用を与え、その後にブログのリンクを追加する(筆者はまだ書いていないため).
応用3:循環三対角線性方程式群の解
このブログでは,循環三対角線方程式群の解におけるSherman‐Morrison式の応用について詳細に述べる.まず理論知識の紹介部分を与える. A∈Rnに対して×n A ∈ R n × nは可逆行列、u,v∈Rnu,v∈Rnは列ベクトル、1+vTA-1 u≠0 1+vT A-1 u≠0であり、方程式(A+uvT)x=b(A+uvT)x=bを解く必要がある.これに対して、まず次の2つの方程式を解くことができます.
Ay=b,Az=u A y = b , A z = u
.
そして令
x=y−vTy 1+vTzz x=y−vTy 1+vTz、この解は元の方程式の解であり、以下のように検証される.
(A+uvT)x=(A+uvT)(y−vTy1+vTzz)=Ay+uvTy−vTy1+vTzAz−vTy1+vTzuvTz=b+uvTy−vTyu+vTyuvTz1+vTz=b+uvTy−(1+vTz)vTyu1+vTz=b+uvTy−uvTy=b ( 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は2つの方程式に分けられ,元の方程式をある程度簡略化できる.次に,サイクル三対角線方程式群の解を紹介する.
いわゆる循環三対角線性方程式群は、係数行列が以下の形式であることを指す.
A=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢b1a20⋮0cnc1b2⋱⋮⋯00c2⋱an−2⋯⋯⋯0⋱bn−2an−100⋮0cn−2bn−1ana10⋮0cn−1bn⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥ A = [ b 1 c 1 0 ⋯ 0 a 1 a 2 b 2 c 2 0 ⋮ 0 0 ⋱ ⋱ ⋱ 0 ⋮ ⋮ ⋮ a n − 2 b n − 2 c n − 2 0 0 ⋯ ⋯ a n − 1 b n − 1 c n − 1 c n 0 ⋯ 0 a n b n ]
サイクルの3対角線性方程式のグループは書くことができます
Ax=d A x=d、ここで
d=(d1,d2,...,dn)T. d = ( d 1 , d 2 , . . . , d n ) T .
この方程式の解について,我々は
u=(γ,0,0,...,cn)T,v=(1,0,0,...,a1γ)T u = ( γ , 0 , 0 , . . . , c n ) T , v = ( 1 , 0 , 0 , . . . , a 1 γ ) T,そして
A=A’+uvT A=A’+uvT
A’A’は以下の通りである.
A′=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢b1−γa20⋮00c1b2⋱⋮⋯00c2⋱an−2⋯⋯⋯0⋱bn−2an−100⋮0cn−2bn−1an00⋮0cn−1bn−a1cnγ⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥ A ′ = [ b 1 − γ c 1 0 ⋯ 0 0 a 2 b 2 c 2 0 ⋮ 0 0 ⋱ ⋱ ⋱ 0 ⋮ ⋮ ⋮ a n − 2 b n − 2 c n − 2 0 0 ⋯ ⋯ a n − 1 b n − 1 c n − 1 0 0 ⋯ 0 a n b n − a 1 c n γ ]
A’A’は三対角行列である.以上の理論知識に基づいて,我々は以下の2つの方程式を解くだけである.
A′y=d,A′z=u, A ′ y = d , A ′ z = u ,
そして、
y,z y,z求める
x x .以上の2つの方程式は3対角線形方程式群であり,追跡法(またはThomas法)で解くことができ,具体的なアルゴリズムはブログ:3対角線形方程式群(tridiagonal systems of equations)の解を参照することができる.
さらに,Sherman‐Morrison公式の考え方を用いて,循環三対角線性方程式群を三対角線形方程式群に変換して解くことができた.このアルゴリズムのPython言語実装を以下に示す.
Python実現
ジルコニウム我々が解くサイクル三対角線方程式のグループは以下の通りである.
⎡⎣⎢⎢⎢⎢⎢⎢4100314100014100014120014⎤⎦⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢x1x2x3x4x5⎤⎦⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢76668⎤⎦⎥⎥⎥⎥ [ 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 ] [ x 1 x 2 x 3 x 4 x 5 ] = [ 7 6 6 6 8 ]
この方程式を解く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]