2019年8月17日土曜日

numpyを使った行列の計算

最近線形代数を復習している。学生の時の講義で単位は取ったはずだが、もう様々なことがうろ覚えになってしまっている。とりあえず、先月にごく簡単な本を一冊読んだ。色々とうろ覚えだったことがはっきりしてきた。ただ単に復習をしているだけのはずだが、結構楽しい(むしろ、学生の時よりも楽しいのはなぜだろう?)。これなら、まだ時間も記憶力もあった(はずの)20代の頃にもっと勉強しておけばよかったと悔やむ気持ちだ。

大学の講義では、章末問題などを解きながら計算を覚えていくが、今回は少し違ったアプローチを取って、傍らでpythonによる計算をしながら行列演算の概念を押さえようと思っている。pythonで行列の演算をするには、いくつか方法があるようだ。ごく簡単に行列っぽい扱いをするには、普通にリストが要素となっているリストを作れば良い。
List1 = [[1,2,3],[4,5,6],[7,8,9]]
print(List1)
# [[1, 2, 3], [4, 5, 6]]
これは今までにも使ってきたやり方だが、行列の演算をするには不便だ。pythonの場合、numpyのarrayやmatrixがあり、これを使えば計算ができるようだ。どちらも基本的に同じことができるようだが、名前の通り、numpy.matrixの方が二次元の行列に特化したクラスのようだ。

まずはnumpyをインポートする。定石の方法に従ってnpとしてインポートする。行列を作るときは以下のようにする。
import numpy as np
ExArr1 = np.array([[1,2,3],[4,5,6]])
print(ExArr1)
# [[1 2 3]
#  [4 5 6]]

ここでは、最初の[1,2,3]が一行目、[4,5,6]が二行目となる。または、matrixを使い、
ExMat1 = np.matrix([[1,2,3],[4,5,6]])
print(ExMat1)
# [[1 2 3]
#  [4 5 6]]
となる。
行列の積を求めるときは、numpy.dot()、numpy.matmul()、または@が使える。
Arr1 = np.array([[1,2,3],[4,5,6],[7,8,9]])
Arr2 = np.array([[1,0,0],[0,1,0],[0,0,1]])
print(np.dot(Arr1, Arr2))
# [[1 2 3]
#  [4 5 6]
#  [7 8 9]]
print(np.matmul(Arr1, Arr2))
# [[1 2 3]
#  [4 5 6]
#  [7 8 9]]
print(Arr1 @ Arr2)
# [[1 2 3]
#  [4 5 6]
#  [7 8 9]]
np.matrixについては、*でも行列の積が計算できる。注意するべきは、numpy.ndarrayでは、*を使うと行列の要素ごとの積になることだ。
Mat1 = np.matrix([[1,2,3],[4,5,6],[7,8,9]])
Mat2 = np.matrix([[1,0,0],[0,1,0],[0,0,1]])
print(Mat1*Mat2)
# [[1 2 3]
#  [4 5 6]
#  [7 8 9]]
print(Arr1*Arr2)
# [[1 0 0]
#  [0 5 0]
#  [0 0 9]]
転置行列は.Tで計算できる。
print(Arr1.T)
# [[1 4 7]
#  [2 5 8]
#  [3 6 9]]
逆行列も簡単に計算できる。**-1と.Iが使えるのはnumpy.matrixの方だけで、numpy.ndarrayでは.linalg.invだけが使える。
mat1 = np.matrix([[2,3],[1,4]])
print(mat1)
# [[2 3]
#  [1 4]]
print(mat1.I)
# [[ 0.8 -0.6]
#  [-0.2  0.4]]
print(mat1**-1)
# [[ 0.8 -0.6]
#  [-0.2  0.4]]
print(np.linalg.inv(mat1))
# [[ 0.8 -0.6]
#  [-0.2  0.4]]
こういうところが、pythonを横で書きながら勉強していった方が捗るような気がする。例えば、逆行列をかけると単位行列になるのは簡単に確かめることができる。
print(mat1 * mat1.I)
# [[1. 0.]
#  [0. 1.]]
print(mat1.I * mat1)
# [[1. 0.]
#  [0. 1.]]
こういう細かい確かめをしながら進めていくのは結構有効なのでは、と思ったのだ。教科書の例題をスクリプトに書き起こしながら読んでいくとか。そういう方法での授業や講義はされているのだろうか?

固有値問題の場合は、numpy.linalg.eig()を使う。
mat1 = np.matrix([[1,2],[-1,4]])
w, v = np.linalg.eig(mat1)
print(w)
# [2. 3.]
print(v)
# [[-0.89442719 -0.70710678]
#  [-0.4472136  -0.70710678]]
固有値がw、それぞれの固有値に対応する固有ベクトルがvの列になっている。これを使い、
print(v**-1 * mat1 * v)
# [[2. 0.]
#  [0. 3.]]
が簡単に計算してみることができる。
pythonの計算の仕方についてはこの方法で教科書の例をそのままスクリプトに書きおこす方式で一通りやってみようと思う。スクリプトを書いて計算する分には、証明を斜め読みにしながら進んでいくことも可能だろう。学び方として褒められたものではないかもしれないが。

参考:
https://note.nkmk.me
nkmkさんのサイトが充実していて、いつも参考にさせていただいている。