数値計算の中心課題の一つである，行列に関する演算について見ていきます． 多次元，大規模な行列に対する効率のよい計算法が多数開発されており， 多くの既存のライブラリが用意されています． 本章ではそれらの中心をなす， 逆行列(matrix inverse)と固有値(Eigen values)に関して具体的な計算方法を示します． 現実的な問題には既存のライブラリを使うのが上策ですが， それでも基礎となる原理の理解や，ちょっとした計算， ライブラリの結果の検証に使えるルーチンを示します．
逆行列は連立一次方程式を解くことと等価です． ルーチン的なやり方にガウスの消去法があります． これは上三角行列になれば代入を適宜おこなうことで 解が容易に求まることを利用します． さらに，初期値から始めて次々に解に近づけていく反復法があります． この代表例であるJacobi(ヤコビ)法と， 収束性を高めたGauss-Seidel(ガウス-ザイデル)法を紹介します．
上記の手法をより高速にした修正コレスキー分解と共役傾斜(共役勾配)法がありますが， 複雑になるのでパスります． 必要ならばNumRecipeを読んでください．
逆行列は連立一次方程式を解くことと等価です． すなわち，$A$を行列，$x$を未知数ベクトル，$b$を数値ベクトルとすると，
$$ \begin{array}{rl} Ax &= b \\ A^{-1}Ax &= A^{-1}b \\ x &= A^{-1}b \end{array} $$です． 未知数の少ない連立一次方程式では，適当に組み合わせて未知数を消していけばいいですが， 未知数が多くなってしまうと破綻します． 未知数の多い多元連立一次方程式で，ルーチン的に解を求めていく方法がガウス消去法で， 前進消去と後退代入という2つの操作からなります．
後退代入(Backward substitution)による解の求め方を先ず見ましょう． たとえば，
$$ \begin{array}{rl} x+y-2z & = -4 \\ -3y+3z & = 9\\ -2z & = -3 \end{array} $$では，下から順番に$z\rightarrow y\rightarrow x$と適宜代入することによって， 簡単に解を求めることが出来ます． 係数で作る行列でこのような形をした上三角行列にする操作を 前進消去あるいはガウスの消去法(Gaussian elimination)といいます． 下三角行列L(lower triangular matrix)と上三角行列U(upper triangular matrix)の積に分解する操作
$$ A = L.U $$をLU分解(LU decomposition)といいます． 例えば先に示した上三角行列を係数とする連立方程式は，
$$ \begin{array}{rl} x+y-2z&=-4 \\ x-2y+z&=5 \\ x-2y-z&=2 \end{array} $$を変形することで得られます．連立方程式からの変形を示してください．
係数行列(coefficient matrix)Aと定数項(b)を定義します．
from pprint import pprint
import scipy.linalg as linalg # SciPy Linear Algebra Library
import numpy as np
# A = np.array([[1,1,-2,-4],[1,-2,1,5],[2,-2,-1,2]])
A = np.array([[1,1,-2],[1,-2,1],[1,-2,-1]])
b = np.array([[-4],[5],[2]])
pprint(A)
pprint(b)
array([[ 1, 1, -2],
[ 1, -2, 1],
[ 1, -2, -1]])
array([[-4],
[ 5],
[ 2]])
単に逆行列を求める際は
inv_A = linalg.inv(A)
print(inv_A)
np.dot(inv_A,b)
[[ 0.66666667 0.83333333 -0.5 ] [ 0.33333333 0.16666667 -0.5 ] [ 0. 0.5 -0.5 ]]
array([[ 0.5],
[-1.5],
[ 1.5]])
拡大係数行列はhstack(horizontal stack)かcolumn_stackを使います．
e_A = np.hstack((A, b))
pprint(e_A)
pprint(np.column_stack((A,b)))
array([[ 1, 1, -2, -4],
[ 1, -2, 1, 5],
[ 1, -2, -1, 2]])
array([[ 1, 1, -2, -4],
[ 1, -2, 1, 5],
[ 1, -2, -1, 2]])
P, L, U = linalg.lu(e_A)
print("P:")
pprint(P)
print("L:")
pprint(L)
print("U:")
pprint(U)
P:
array([[ 1., 0., 0.],
[ 0., 1., 0.],
[ 0., 0., 1.]])
L:
array([[ 1., 0., 0.],
[ 1., 1., 0.],
[ 1., 1., 1.]])
U:
array([[ 1., 1., -2., -4.],
[ 0., -3., 3., 9.],
[ 0., 0., -2., -3.]])
元の拡大行列はこれらの掛け算で求められます．
print("e_A=P.L.U:")
pprint(np.dot(np.dot(P,L),U))
e_A=P.L.U:
array([[ 1., 1., -2., -4.],
[ 1., -2., 1., 5.],
[ 1., -2., -1., 2.]])
LU分解すれば線形方程式の解が容易に求まることは理解できると思う．具体的に$A$をLU分解する行列(消去行列と称す)T1,T2の係数は次のようにして求められる．
from pprint import pprint
import scipy.linalg as linalg # SciPy Linear Algebra Library
import numpy as np
A = np.array([[1.0,1,-2],[1,-2,1],[1,-2,-1]])
A0 = np.array([[1.0,1,-2],[1,-2,1],[1,-2,-1]])
b = np.array([[-4.0],[5],[2]])
n = 3
L = np.identity(n)
T = []
for i in range(n): #i行目
T.append(np.identity(n))
for j in range(i+1, n):
am = A[j,i]/A[i,i] #i行の要素を使って，i+1行目の先頭を消す係数を求める
T[i][j,i]=-am #i番目の消去行列に要素を入れる
L[j,i]=am #LTMの要素
for k in range(n):
A[j,k] -= am*A[i,k] #もとの行列をUTMにしていく
b[j] -= b[i]*am #数値ベクトルも操作
pprint(A)
pprint(b)
array([[ 1., 1., -2.],
[ 0., -3., 3.],
[ 0., 0., -2.]])
array([[-4.],
[ 9.],
[-3.]])
上のコードによって得られた消去行列．
pprint(T)
[array([[ 1., 0., 0.],
[-1., 1., 0.],
[-1., 0., 1.]]),
array([[ 1., 0., 0.],
[ 0., 1., 0.],
[ 0., -1., 1.]]),
array([[ 1., 0., 0.],
[ 0., 1., 0.],
[ 0., 0., 1.]])]
これを実際に元の行列$A0$に作用させると，UTMが求められる．
np.dot(T[1],np.dot(T[0],A0))
array([[ 1., 1., -2.],
[ 0., -3., 3.],
[ 0., 0., -2.]])
求められたLTM, UTMを掛けると
np.dot(L,A)
array([[ 1., 1., -2.],
[ 1., -2., 1.],
[ 1., -2., -1.]])
元の行列を得られる．L,Aに求めたい行列が入っていることを確認．
pprint(L)
pprint(A)
array([[ 1., 0., 0.],
[ 1., 1., 0.],
[ 1., 1., 1.]])
array([[ 1., 1., -2.],
[ 0., -3., 3.],
[ 0., 0., -2.]])
数値ベクトルも期待通り変換されている．
pprint(b)
array([[-4.],
[ 9.],
[-3.]])
luの出力で得られる拡張係数行列のUTMは，次の通り．
np.hstack((A, b))
array([[ 1., 1., -2., -4.],
[ 0., -3., 3., 9.],
[ 0., 0., -2., -3.]])
ガウス消去法で困るのは，割ろうとした対角要素が0の場合である．しかし，この場合にも，方程式の順序を，行列の行と右辺の値をペアにして入れ替えれば解決する．この割る ほうの要素をピボット要素あるいはピボット(pivot,バスケの軸足を動かさずにくるくる回すやつ)と呼ぶ．この操作は，変数の並びを変えたわけではなく，単に方程式の 順番を変更する操作に相当する．
さらに対角要素の数値が厳密に0でなくとも，極端に0に近づいた場合にも，その数で割った数値が大きくなり他の数との差を取ると以前に示した情報落ちの可能性が出てくる． この現象を防ぐためには，絶対値が最大のピボットを選んで行の入れ替えを毎回おこなうといい結果が得られることが知られている．
MapleのLUDecompositionコマンドをこのような行列に適用すると，置換行列(permutation matrix)Pが単位行列ではなく，ピボット操作に対応した行列となる．P.A=L.Uとなることに注意．
下の例では，$a_{11}$と$a_{31}$の数値を少し変えている． その結果，置換行列Pはidentity行列ではなく， 2，3行目でpivotが起こっていることを表している．
A = np.array([[2,1,-2],[1,-2,1],[2,-2,-1]])
b = np.array([[-4.0],[5],[2]])
e_A = np.hstack((A, b))
pprint(e_A)
P, L, U = linalg.lu(e_A)
print("P:")
pprint(P)
print("L:")
pprint(L)
print("U:")
pprint(U)
array([[ 2., 1., -2., -4.],
[ 1., -2., 1., 5.],
[ 2., -2., -1., 2.]])
P:
array([[ 1., 0., 0.],
[ 0., 0., 1.],
[ 0., 1., 0.]])
L:
array([[ 1. , 0. , 0. ],
[ 1. , 1. , 0. ],
[ 0.5 , 0.83333333, 1. ]])
U:
array([[ 2. , 1. , -2. , -4. ],
[ 0. , -3. , 1. , 6. ],
[ 0. , 0. , 1.16666667, 2. ]])
以下のような連立方程式を
$$ \left[ \begin {array}{c} 5\,x+y+z+u\\ x+3\,y+z+u\\ x-2\,y-9\,z+u\\ x+3\,y-2\,z+5\,u\end {array} \right] = \left[ \begin {array}{c} -6\\ 2\\ -7\\ 3\end {array} \right] $$形式的に解くと
$$ x=\frac{-6-(y+z+u)}{5} $$となる．他の未知数も，
| y= | |
|---|---|
| z= | |
| u= |
となる．適当に初期値($x_0,y_0,z_0,u_0$)をとり，下側の方程式に代入すると，得られた出力($x_1,y_1,z_1,u_1$)はより精解に近い値となる．これを繰り返すことによって精解が得られる．これをヤコビ(Jacobi)法と呼び，係数行列の対角要素が非対角要素にくらべて大きいときに適用できる．多くの現実の問題ではこの状況が成り立っている．
Gauss-Seidel法はJacobi法の高速版である．$n$番目の解の組が得られた後に一度に次の解の組に入れ替えるのではなく，得られた解を順次改良した解として使っていく．これにより，収束が早まる．以下にはヤコビ法のコードを示した．x1[i]の配列を変数に換えるだけで，Gauss-Seidel法となる．
# Jacobi iterative method for solving Ax=b by bob
import numpy as np
np.set_printoptions(precision=6, suppress=True)
A=np.array([[5,1,1,1],[1,3,1,1],[1,-2,-9,1],[1,3,-2,5]])
b=np.array([-6,2,-7,3])
n=4
x0=np.zeros(n)
x1=np.zeros(n)
for iter in range(0, 20):
for i in range(0, n):
x1[i]=b[i]
for j in range(0, n):
x1[i]=x1[i]-A[i][j]*x0[j]
x1[i]=x1[i]+A[i][i]*x0[i]
x1[i]=x1[i]/A[i][i]
for j in range(0, n):
x0[j]=x1[j]
print(x0)
[-1.2 0.666667 0.777778 0.6 ] [-1.608889 0.607407 0.562963 0.751111] [-1.584296 0.764938 0.54749 0.782519] [-1.618989 0.751429 0.518705 0.676892] [-1.589405 0.807797 0.506116 0.680422] [-1.598867 0.800956 0.497269 0.635649] [-1.586775 0.821983 0.492763 0.638108] [-1.590571 0.818635 0.489707 0.621271] [-1.585923 0.826531 0.488159 0.622816] [-1.587501 0.824982 0.487092 0.61653 ] [-1.585721 0.82796 0.486563 0.617348] [-1.586374 0.82727 0.48619 0.614993] [-1.585691 0.828397 0.486009 0.615389] [-1.585959 0.828098 0.485878 0.614503] [-1.585696 0.828526 0.485817 0.614684] [-1.585805 0.828398 0.485771 0.61435 ] [-1.585704 0.828561 0.48575 0.61443 ] [-1.585748 0.828508 0.485734 0.614304] [-1.585709 0.82857 0.485727 0.614338] [-1.585727 0.828548 0.485721 0.61429 ]
もう少しpython(numpy)の機能を使うとおしゃれな書き方ができる． 以下は，ネットから拾ってきたcode. 上(by bob)とやってることの本質的なところは同じですが， numpyのindexを自動的に判断する機能を使って簡潔に書いています． こういうのを読める，書けるようになれば，あなたもpython使い． 私はまだpythonに馴染めてなくて，いちいち頭の中で治しています．．．
# https://www.quantstart.com/articles/Jacobi-Method-in-Python-and-NumPy
from pprint import pprint
from numpy import array, zeros, diag, diagflat, dot
np.set_printoptions(precision=6, suppress=True)
def jacobi(A,b,N=25,x=None):
if x is None:
x = zeros(len(A[0]))
D = diag(A)
R = A - diagflat(D)
# Iterate for N times
for i in range(N):
x = (b - dot(R,x)) / D
print(x)
return x
A = array([[5,1,1,1],[1,3,1,1],[1,-2,-9,1],[1,3,-2,5]])
b = array([-6,2,-7,3])
sol = jacobi(A,b,N=10)
[-1.2 0.666667 0.777778 0.6 ] [-1.608889 0.607407 0.562963 0.751111] [-1.584296 0.764938 0.54749 0.782519] [-1.618989 0.751429 0.518705 0.676892] [-1.589405 0.807797 0.506116 0.680422] [-1.598867 0.800956 0.497269 0.635649] [-1.586775 0.821983 0.492763 0.638108] [-1.590571 0.818635 0.489707 0.621271] [-1.585923 0.826531 0.488159 0.622816] [-1.587501 0.824982 0.487092 0.61653 ]
次の行列AをLU分解し，上・下三角行列を求めよ．さらに連立方程式の解を求めよ. $$ \left[ \begin {array}{c} x+4\,y+3\,z\\ x-2\,y+z \\ 2\,x-2\,y-z \end {array} \right] = \left[ \begin {array}{c} 11\\ 11\\ 11 \end {array} \right] $$
次の連立方程式の拡大行列をLU分解せよ．pivot操作が必要となる． P.A = L.Uを確かめよ． $$ \left[ \begin {array}{c} 3\,w+2\,x+2\,y+z\\ 3\,w+2\,x+3\,y+z\\ w-2\,x-3\,y+z\\ 5\,w+3\,x-2\,y+5\,z\end {array} \right] = \left[ \begin {array}{c} -6\\ 2\\ -9\\ 2\end {array} \right] $$
Jacobi法のプログラムを参照してGauss-Seidel法のプログラムを作れ．
A=np.array([[5,1,1,1],[1,3,1,1],[1,-2,-9,1],[1,3,-2,5]])
b=np.array([-6,2,-7,3])
n=4
に対して計算を行い，Jacobi法と収束性を比べよ．
最初の数回はこんな感じ．
[-1.2 1.066667 0.407407 0.362963]
[-1.567407 0.932346 0.436763 0.528779]
[-1.579578 0.871345 0.46739 0.580064]
[-1.58376 0.845435 0.478382 0.600844]
[-1.584932 0.835236 0.482827 0.608976]