p1_ex12

三角行列

モジュールのインポート

import numpy as np
import scipy as sp
from scipy.linalg import solve_triangular # 上三角あるいは下三角行列に対する連立方程式の計算
from scipy.linalg import lu # LU分解
from scipy.linalg import cholesky # コレスキー分解

関数 solve_triangular

　次のような行列表示の連立方程式を考える． $$ \left( \begin{array}{cccc} 1 & 0 & 0 & 0\\ 2 & 3 & 0 & 0\\ 4 & 5 & 6 & 0\\ 7 & 8 & 9 & 1\\ \end{array} \right) \boldsymbol{x} = \left( \begin{array}{c} 1 \\ 8 \\ 32 \\ 54 \\ \end{array} \right) $$

a = np.array([[1,0,0,0],[2,3,0,0],[4,5,6,0],[7,8,9,1]])
x = np.array([1,2,3,4])
a,x

(array([[1, 0, 0, 0],
        [2, 3, 0, 0],
        [4, 5, 6, 0],
        [7, 8, 9, 1]]),
 array([1, 2, 3, 4]))

b = a@x
b

array([ 1,  8, 32, 54])

三角行列（Triangular matrices）の場合，
https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve_triangular.html

x_sol = solve_triangular(a, b, lower=True) # lower=Trueのとき，下三角行列を用いて計算する
x_sol

array([1., 2., 3., 4.])

a.dot(x) # Check the result

array([ 1,  8, 32, 54])

a.dot(x) == b

array([ True,  True,  True,  True])

LU分解

　行列$A$を $$A=PLU$$ の形に分解する方法として，LU分解（LU decomposition）がある．ただし．$P$は順序置換行列（permutation matrix），$L$は下三角行列（lower triangular matrix），$U$は上三角行列（ppper triangular matrix）を示す．
https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lu.html

a = np.array([[2,5,7,8],[5,2,2,8],[7,5,6,6],[5,4,4,8]])
p, l, u = lu(a)

順序置換行列は，

p # 順序置換行列

array([[0., 1., 0., 0.],
       [0., 0., 0., 1.],
       [1., 0., 0., 0.],
       [0., 0., 1., 0.]])

下三角行列は，

l # 下三角行列

array([[ 1.        ,  0.        ,  0.        ,  0.        ],
       [ 0.28571429,  1.        ,  0.        ,  0.        ],
       [ 0.71428571,  0.12      ,  1.        ,  0.        ],
       [ 0.71428571, -0.44      , -0.04347826,  1.        ]])

上三角行列は，

u # 上三角行列

array([[ 7.        ,  5.        ,  6.        ,  6.        ],
       [ 0.        ,  3.57142857,  5.28571429,  6.28571429],
       [ 0.        ,  0.        , -0.92      ,  2.96      ],
       [ 0.        ,  0.        ,  0.        ,  6.60869565]])

行列$A$は，

p@l@u

array([[2., 5., 7., 8.],
       [5., 2., 2., 8.],
       [7., 5., 6., 6.],
       [5., 4., 4., 8.]])

array([[2, 5, 7, 8],
       [5, 2, 2, 8],
       [7, 5, 6, 6],
       [5, 4, 4, 8]])

解の確認をして，

np.allclose(a - p@l@u, np.zeros((4, 4)))

True

x = np.array([1,2,3,4])
b = a@x
b

array([65, 47, 59, 57])

lu = sp.linalg.lu_factor(a)
sp.linalg.lu_solve(lu, b)

array([1., 2., 3., 4.])

コレスキー分解

　行列$A$が正定値対称行列（固有値が全て正であるような対称行列）ならば，下三角行列とその複素・転置より $$A = L L^{T*}$$ あるいは上三角行列とその複素・転置より $$A = U^{T*} U$$ に分解できる．これをコレスキー分解（Choleski decomposition）という．
https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.cholesky.html

　ここでは，実数について簡単な計算例を示そう．

a = np.array([[1,0.2],[0.2,1]])
a

array([[1. , 0.2],
       [0.2, 1. ]])

l = cholesky(a,lower=True) # Choleski decomposition
l

array([[1.       , 0.       ],
       [0.2      , 0.9797959]])

l@l.T

array([[1. , 0.2],
       [0.2, 1. ]])

l@l.T - a

array([[ 0.00000000e+00,  0.00000000e+00],
       [ 0.00000000e+00, -1.11022302e-16]])

array([[1. , 0.2],
       [0.2, 1. ]])

x = np.array([1,2])
b = a@x
b

array([1.4, 2.2])

さて， $$a x = b = L L^T x = L t$$ とおいて， $$L t = b$$ をまず解く．

t =np.linalg.solve(l, b)
x =np.linalg.solve(l.T, t)
x

array([1., 2.])

t2 = solve_triangular(l, b, lower=True) # lower=Trueのとき，下三角行列を用いて計算する
x2 = solve_triangular(l.T, t2) # lower=Trueのとき，上三角行列を用いて計算する
x2

array([1., 2.])

array([[1. , 0.2],
       [0.2, 1. ]])

u = cholesky(a) # Choleski decomposition
u

array([[1.       , 0.2      ],
       [0.       , 0.9797959]])

u.T@u

array([[1. , 0.2],
       [0.2, 1. ]])

u.T@u - a

array([[ 0.00000000e+00,  0.00000000e+00],
       [ 0.00000000e+00, -1.11022302e-16]])