Numpy を用いた行列演算の基礎(対角化まで)

筆者が初めてnumpyで行列についていくつか処理をする際に学んだことを，こちらにまとめておきたいと思います．対角化を目標にはしていますが対角化は後半にあり，内容のほとんどは行列に関する基本的な処理ですので，初めてnumpyで行列処理をする方にも理解できるよう心がけたつもりです．．．

行列を作成(定義)

numpy.ndarray型

一番よくつかわれて扱いやすく，困ったときの解説もネットに多い型です．おそらく，python で行列演算する際に最もよく使われる型ではないのでしょうか．なお，型とオブジェクトは厳密には違いますが，区別なく話したいと思いますので，型という言葉がたくさん出てきます．

要素を指定すると同時に，行列定義

import numpy as np
A = np.array([[0,1,-2],[-3,7,-3],[3,-5,5]])

最初の行でnumpyを使う宣言をしています．2行目で行列の要素を決めるとともに定義していて，
[math]
A = \left(
\begin{array}{ccc}
0 & 1 & -2 \\
-3 & 7 & -3 \\
3 & -5 & 5
\end{array}
\right)
[/math]
という行列を定義しています．この行列の型は「numpy.ndarray」です．よく見ると，np.array の丸括弧()ので中で多次元配列(リスト)を作成し，それをnp.arrayに変換するというのが2行目の仕組みであることがわかります．numpy.ndarray型はリスト型とは異なることに注意しよう．

行列の動的確保，その後要素を代入

上と同じ行列を，行列の大きさだけ先に定義しあとで要素を代入する方法で作成してみます．

import numpy as np
A = np.zeros((3, 3))

A[0][1] = 1
A[0][2] = -2
A[1][0] = -3
A[1][1] = 7
A[1][2] = -3
A[2][0] = 3
A[2][1] = -5
A[2][2] = 5

「np.zeros( (n,m) )」で「すべての行列要素が0である，n行m列の行列(numpy.array型)を作成する」という指示になります．これで希望の大きさを持った行列を定義することができます．それより後の行でゼロでない部分の要素を代入しています．「A[0][2]」は「行列Aの1行3列目」を指定しています．

リスト型(list)の多次元配列として行列を動的に作成

numpy.array型でなく，リストとして扱った方が都合がいいのであればこちらを用いましょう．次のサイトを参考にしました

Pythonで2次元配列の静的確保と動的確保 – sonickun.log

A = [[0 for x in range(3)] for y in range(3)]

作成時点で要素はすべて0になっています．

numpy.matrixを用いて作成

A = np.matrix([[0,1,-2],[-3,7,-3],[3,-5,5]])

この後に print( type(A) )と書いてみるととなります．この部分は筆者も勉強不足でまだ確実には言えないのですが，この型はあまり使われていないようです．詳しくはこちらを読んで勉強したいと思っています．

私訳「暫定的 NumPy チュートリアル」- naoya_t@hatenablog

複素行列の場合

虚数の設定(complex型)：正しい例

import numpy as np
A = np.array([[1.j,1,-2],[-3,7,-3],[3,-5,5]])

これによって作成される行列は
[math]
A = \left(
\begin{array}{ccc}
i & 1 & -2 \\
-3 & 7 & -3 \\
3 & -5 & 5
\end{array}
\right)
[/math]
です．単位虚数をpythonで表す際にははiでなくjで表すことに注意しましょう．この虚数はpythonで「complex型」として扱われます．

ちなみに，次のようなやり方でも同様の虚数を設定できます．

import numpy as np
A = np.array([[complex(0,1),1,-2],[-3,7,-3],[3,-5,5]])

実数a,bを用いて「complex(a,b)」とすることで[math]a+ib[/math]という虚数を作成できます．こちらの方がエラーなく安全に使えそうですね．このように作成しても先ほどの場合と同じくcomplex(a,b)はcomplex型になります．complex型はnumpyの機能でなくpythonの機能であることに注意しましょう．

虚数の設定：間違った例

虚数を要素として持つ行列を作る際の注意点です．まずは次のコードを見てみましょう

import numpy as np
A = np.array([[1*j,1,-2],[-3,7,-3],[3,-5,5]])

虚数の部分を見ていただくと一見もっともらしいですが，これを実行すると「NameError: name ‘j’ is not defined」というエラーが吐かれます．jを単位虚数と認識しないようです．

また，次の例も同じエラーを吐きます．

import numpy as np
A = np.array([[j,1,-2],[-3,7,-3],[3,-5,5]])

これらの虚数の設定については，numpyの問題ではなくpythonそもそもの問題です．虚数の扱いについては今後数値計算の投稿にも書きたいと思っていますが，それまではこちらを参照していただければよいでしょう．

Complex numbers usage in python – stackoverflow

numpy.arrayを使うことが安全

リストやmatrixなどの型があると紹介しましたが，行列演算に特に条件がない限りはnumpy.arrayを使うことが安全でしょう．例えばmatrix型をいったん作成をすると，後になって A[2][1]=3.0 のように行列要素にアクセスして変更することができませんでした．作成された時点で固定されるようです．また，この投稿で後に説明する「行列の演算，性質の確認」に書いてある機能もnumpy.arrayについての機能になり，ほかの型ではうまくいかないことがあります(正確にどれがうまくいかないかは未確認ですが，書いている途中になんかエラーが出るなと思ったらmatrix型でダメだったということがありました)．numpyの機能を安全かつ最大に活用するため，numpy.ndarray型を積極的に用いましょう．

おまけ：要素をすべて0に

あまり使う場面はないかと思いますが，次のような面白いコマンドがあります．

import numpy as np
A = np.array([[1.j,1,-2],[-3,7,-3],[3,-5,5]])
A.fill(0)

最後にprint(A)を実行してみると，2行目で設定したAの要素がすべてゼロになっていることがわかります．このように「A.fill(x)」とすることで「行列の要素をすべてxにしろ」という指示になります．

ベクトルの作成

対角化をすると，ベクトルを扱うことがあるため，こちらについても触れておく．

ベクトルの動的作成

a = np.zeros(4)

これで4つの0を要素に持つゼロベクトルaが作成される．もちろん，リスト型でなくnumpy.array型である．

ベクトルの要素にアクセス，代入

a = np.zeros(4)
a[0] = 2

print(a) で出力すると[ 2. 0. 0. 0.]となります．つまり，numpy.array もリストの場合と同様に 0 から番号付けされていて，「a[0]=2」とすることで最初の要素に2を代入することができます．型はnumpy.arrayですが代入やアクセス方法はリストと同じです．

行列の演算，操作，性質の確認

行列の様々な演算や，行列式やランクの計算などの行列の性質を見る際には，numpy の linalg というところの機能をよく用います．

行列式(determinant)の計算

numpy.linalg.detという機能を使います．

import numpy as np
A = np.array([[0,1,-2],[-3,7,-3],[3,-5,5]])
print( np.linalg.det(A) )

この場合実行結果は「18.0」となりました．

ランク(rank)の計算

numpy.linalg.matrix_rankという機能を使います．

import numpy as np
A = np.array([[0,1,-2],[-3,7,-3],[3,-5,5]])
print( np.linalg.matrix_rank(A) )

実行結果は「3」となりました．

転置

行列のお尻に.Tという記号を付けるだけで実現できます．

import numpy as np
A = np.array([[0,1,-2],[-3,7,-3],[3,-5,5]])
B = A.T

BにAの転置行列が収められます．

複素共役

numpy.conjugateという機能を用います．

import numpy as np
A = np.array([[1.j,1,-2],[-3,7,-3],[3,-5,5]])
B = np.conjugate(A)

Bには，Aの要素についてすべて複素共役をとったものが収められます．確認のため最後の行に print(B-A)を追加して実行してみると，確かに複素数の値を持つ1行1列目だけ2.jという値を返しています．

ちなみに，このnumpy.conjugateは引数としてnumpy.ndarray型だけなく複素数単体を入れてもその複素共役を返してくれます．

エルミート共役

転置と複素共役の合わせ技で実現可能です．

import numpy as np
A = np.array([[1.j,1,-2],[-3,7,-3],[3,-5,5]])
B = np.conjugate(A.T)

このようにすることで，BにAのエルミート共役をとった行列が得られます．

行列の積

行列の積は numpy.dotという機能を使います．

import numpy as np
A = np.array([[1.j,1,-2],[-3,7,-3],[3,-5,5]])
B = np.array([[0,3,2],[-3,2,-5],[2,-1,9]])
C = np.dot(A,B)

「numpy.dot(A,B)」で「AとBの行列積を計算してください」という意味になります．なのでこのコードを実行するとCに[math]A\cdot B[/math]の結果が格納されます．行列の積は順番に依存しますが，numpy.dot の順番は掛け算の順番に等しいです．ちなみにこのnumpy.dotは行列とベクトルの積も計算できる代物になります．

逆行列

逆行列は「numpy.linalg.inv」という機能を使うことで求めてくれます．

import numpy as np
A = np.array([[1.j,1,-2],[-3,7,-3],[3,-5,5]])
B = np.linalg.inv(A)
print( np.dot(A,B) )

「numpy.linalg,inv(A)」で「Aという行列の逆行列を返してください」という意味になります．上のプログラムでは行列[math]A[/math]の逆行列を[math]B[/math]として，その二つをかけた際に単位行列になるかチェックしています．筆者の環境では丸め誤差で非対角項にも[math]10^{-14}[/math]程度の要素を持ってしまいました．

ベクトルの処理

ベクトルのノルム(norm)

import numpy as np
a = np.array([1, 3, -4])
print( np.linalg.norm(a) )

実行結果は5.09901951359でした．「numpy.linalg.norm(v)」で「ベクトルvの長さ(ノルム)を計算しろ」という指示になります．

ベクトルの要素をすべて足す

import numpy as np
a = np.array([1, 3, -4])
print( sum(a) )

実行結果は0でした．「sum(ベクトル)」で「ベクトルの要素をすべて足せ」という指示になります．

行列の対角化

実行列を例に挙げていますが，複素行列も可能です．

固有値を出す

import numpy as np
A = np.array([[0,1,-2],[-3,7,-3],[3,-5,5]])
l, P = np.linalg.eig( A )

lには固有値がnumpy.ndarray型で，Pには対角化された行列を[math]D[/math]としたとき，[math]D=P^{-1}AP[/math]となる[math]P[/math]がnumpy.ndarray型として格納されます．固有値の順番はランダムです(大きい順でも小さい順でもない)．また，lに格納されている固有値の順番とPは対応しています．つまり，Pの2列目はlの固有値の2番目，つまりl[1]の固有ベクトルとなっています．

固有ベクトルを出す

対角化で求まった[math]P[/math]から固有ベクトルをnumpy.ndarray型に収める方法です．[math]P[/math]が求まっている時点で少しひねるだけで実現可能です．

import numpy as np
A = np.array([[0,1,-2],[-3,7,-3],[3,-5,5]])
l, P = np.linalg.eig( A )
v = np.transpose(P)

上のようにすると，v[1]には固有値l[1]のときの(規格化された)固有ベクトルがnumpy.array型で格納されるようになります．ただ[math]P[/math]を転置をしてあげることでよかったんですね．