Pythonプログラミング(ステップ8・行列の計算)

数学(線形代数)で学んだ行列の計算をPythonで行う。

1.行列式の計算

数学(線形代数)の授業で行列式の計算を学んだ際に、何とも手間のかかることか、と、感じた諸君も多いのではなかろうか。 そういった面倒な計算はコンピュータにやらせるに限る。 まずは、効率は度外視して、「教科書通り」に計算する方法をPythonに移植してみよう。

例えば、第1列目に着目して行列式を余因子展開すると、 $$ \begin{eqnarray} \begin{vmatrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1n} \\ a_{21} & a_{22} & a_{23} & \cdots & a_{2n} \\ a_{31} & a_{32} & a_{33} & \cdots & a_{3n} \\ \vdots & & & & \\ a_{n1} & a_{n2} & a_{n3} & \cdots & a_{nn} \end{vmatrix} & = & (-1)^{1+1} a_{11} \begin{vmatrix} a_{22} & a_{23} & \cdots & a_{2n} \\ a_{32} & a_{33} & \cdots & a_{3n} \\ \vdots & & & \\ a_{n2} & a_{n3} & \cdots & a_{nn} \end{vmatrix} \\ & + & (-1)^{2+1} a_{21} \begin{vmatrix} a_{12} & a_{13} & \cdots & a_{1n} \\ a_{32} & a_{33} & \cdots & a_{3n} \\ \vdots & & & \\ a_{n2} & a_{n3} & \cdots & a_{nn} \end{vmatrix} \\ & + & (-1)^{3+1} a_{31} \begin{vmatrix} a_{12} & a_{13} & \cdots & a_{1n} \\ a_{22} & a_{23} & \cdots & a_{2n} \\ \vdots & & & \\ a_{n2} & a_{n3} & \cdots & a_{nn} \end{vmatrix} + \cdots \end{eqnarray} $$ といった具合になり、サイズ$n$の行列の行列式が、サイズ$n-1$の行列式$n$個の計算に還元できる。

この様子を愚直にPythonコードに翻訳した例を以下に示す。

行列はPythonの二次元リストによって表現している。 リストの添字は0から始まるので、数学の教科書とは1だけ表記がずれているかもしれないが、符号の計算($(-1)^{i+j}$)のところの表式は($i+j$のズレが2なので)そのままで構わない。

プログラム中では2つの関数を定義している。 cofactor(i,j,A)は、行列$A$の$i$行$j$列を取り除いた部分行列を生成するための関数、そして、 determinant(A)が行列式を計算する関数で、公式通りに、小行列式を再帰呼び出しで求めながら「総和」しているだけのシンプルな構造である。

余因子展開による行列式計算
(determinant.py)

# coding: utf-8

def submatrix(i,j,A):
    n = len(A[0])
    B = [ ]
    for k in range(n):
        if k==i:
            continue
        x=[ ] # gen row vector
        for l in range(n):
            if l==j:
                continue
            x.append(A[k][l])
        B.append(x)
    return B

# cofactor expansion
def determinant(A):
    n=len(A)
    if n==1:
        return A[0][0]
    det=0
    for i in range(n):
        det = det + (-1)**(i+0) * A[i][0] * determinant(submatrix(i,0,A))
    return det

# メイン部

X=[[1,3,2,4,5],
   [2,4,3,1,3],
   [3,1,2,3,4],
   [4,5,3,2,1],
   [5,2,1,4,3]]

d=determinant(X)

print(d)

icon-pc 練習:行列式の計算量

上のコードでサイズ$n$の行列を計算する場合、関数determinant( )は何回呼び出されることになるか。

icon-pc 練習:数式処理

上のコードを元にして

X=[['a','b','c','d'],
   ['e','f','g','h'],
   ['i','j','k','l'],
   ['m','n','o','p']]

のように、文字の行列を入力としてして与えたとき、記号の式の形で行列式を出力するよう改造してみなさい。

icon-hint ヒント

元のコードで「計算」を行っている箇所はdet = det + (-1)**(i+0) * A[i][0] * determinant(cofactor(i,0,A)) のみである。 ここで、加算や乗算の代わりに、"+何々"、"*(何々)"といった文字列を連結していけば良さそうだ。

できるだけ不要な符号マークや括弧を出力しないよう工夫すると良い。[['a','b','c'],['d','e','f'],['g','h','i']]の場合の出力例:

a*(e*i-h*f)-d*(b*i-h*c)+g*(b*f-e*c)

2.NumPyライブラリの利用

Pythonには、行列の計算をはじめとして、一般的・汎用的な計算を行うための「外付け」機能がいろいろと用意されている。 行列計算でよく使われているのがNumPy(ナン・パイ)である。 この節では、NumPyの基本的な使い方を紹介する。

準備

もしまだなら、ウェブ等の情報を手がかりに、使っているシステムにまずNumPyをインストールする。そして、Pythonコードのはじめに

import numpy as np

と書く。こうしておくことで、np.何々 という名前でNumPyの機能が使えるようになる。

ndarrayオブジェクトの生成

NumPyでは行列を ndarray という形式(オブジェクト)で管理する。これはリストに似ているものの、別物である。 ベクトルや行列を表現・生成するには、以下のようにする:

# 横(行)ベクトル
x = np.array([0.4, 0.1, -1.2])	

# 行列
A = np.array([
   [1,3,2,4,5],
   [2,4,3,1,3],
   [3,1,2,3,4],
   [4,5,3,2,1],
   [5,2,1,4,3]]
   )

# 7行7列の単位行列を生成
E = np.identity(7)

# 3行4列の行列を生成(各要素は0)
Z = np.zeros((3,4))

このように生成したデータにアクセスするには、リストの場合と同様、A[i][j] のように要素を指定できるほか、 ndarrayでは A[i,j] という表現も可能である。

行列の計算

NumPyには基本的な行列計算機能が備わっている。使い方は以下の例を見れば一目瞭然であろう。 ここで、行列の積の計算には .dot() を用いる点に注意。A * B は、要素毎に積をとった行列、という意味になる。

A = np.array([[1,2],[3,4]])
B = np.array([[-2,1],[1,-1]])

# 行列の和
C = A + B

# 差
D = A - B

# 行列の積 
F = A.dot(B)
G = B.dot(A)

# 転置
H = A.transpose()
# または
H = A.T

† 一次元配列、例えば、np.array([1,2,3]) を転置しても [1,2,3] のままである。 横(行)ベクトル np.array([[1,2,3]]) を転置すると、縦(列)ベクトル[[1],[2],[3]] が得られる。

線形代数ライブラリ関数

さらに、NumPyの中の線形代数(linalg)ライブラリをさらに読み込むことで、以下の関数が使えるようになる。 関数を呼び出す前に、

from numpy.linalg import det,inv,solve

のように、使用する関数をあらかじめインポートしておく。

‡ 上記のようにimportする代わりに、以下の関数をnp.linalg.det( ) 等のように、.linalg.関数名( ) で呼び出すこともできる。

関数 機能
diag()
対角要素を成分とする1次元配列を返す
trace()
対角要素の和(トレース:跡)を返す
det()
行列式を返す
eig()
正方行列の固有値と固有ベクトルを返す
†固有ベクトルは縦(列)ベクトルとして返される点に注意
inv()
正方行列の逆行列を返す
pinv()
行列のMoore-Penrose疑似逆行列を返す
qr()
行列をQR分解した結果を返す
svd()
行列を特異値分解した結果を返す
solve()
連立1次方程式の解を求める
lstsq()
連立1次方程式の最小自乗解を求める

以下に、線形代数ライブラリ関数を使って、行列式の計算と、連立方程式の解を2通りの方法で求めるコードの例を示す。

# coding: utf-8

import numpy as np
from numpy.linalg import det,solve,inv

A = np.array([
   [1,3,2,4,5],
   [2,4,3,1,3],
   [3,1,2,3,4],
   [4,5,3,2,1],
   [5,2,1,4,3]]
   )
b = np.array([1,0,0,0,1])

# 行列式の計算
print(det(A))

# 連立方程式 A x = b を解く
x = solve(A,b)
print(x)

# 今度は、A^{-1} b で解を計算する
x2 = inv(A).dot(b)
print(x2)

# 両者の自乗差をチェックする
diff = (x-x2).dot((x-x2))
print(diff)

icon-pc 練習:行列式ふたたび

このページの最初の余因子展開による行列式の計算結果を、NumPyで検算してみなさい。

サイズの大きな行列(8程度以上)を用意し、両者の計算時間も比較してみなさい。

icon-hint ヒント

NumPyを使うと、要素の値がランダムな、任意のサイズの行列を簡単に生成できる。 以下は、要素が$(-1,1)$の一様乱数であるような $10 \times 10$の行列を生成する例である。

import numpy as np

A = np.random.uniform(-1, 1, (10, 10))

また、ndarrayの行列をリストに変換するには

X = A.tolist()

とすればよい。 X にAの内容と同じデータがリストとしてセットされる。

icon-teacher 解説: ランダム行列の固有値分布

上のヒントで、検算のためにランダムな行列を生成する方法を紹介したが、ランダムな行列の振る舞い方には興味深い「規則性」がいろいろと見つかっている。 例えば、各要素が、平均0で、分散 $1/n$ のガウス分布するようなランダムな値の $n \times n$ 正方行列の固有値を複素平面にプロットすると、 $n$ が大きい場合、半径1の円の内部に一様に分布することがわかっている (円則(circular law))。


サイズ100のランダムな行列を生成して、その100個の固有値を複素平面上に赤点でプロットした結果(右のコード)。

その様子を、NumPyとmatplotlibを使って確認するコードの例を以下に示す:

# coding: utf-8
import math
import numpy as np
from numpy.linalg import eig
import matplotlib.pyplot as plt

n = 100

A = np.random.normal(0, 1/math.sqrt(n), (n, n))

z,v = eig(A)

xdata = np.real(z)
ydata = np.imag(z)

plt.plot(xdata,ydata, '.', color=(1.0,0,0.0))
plt.xlim(-2,2)
plt.ylim(-2,2)
plt.xlabel('EIGENVALUE(REAL PART)')
plt.ylabel('EIGENVALUE(IMAG PART)')
plt.grid(True)
plt.show()

この例の中では、関数 np.random.normal(平均, 標準偏差, (行数, 列数))によって、正規(ガウス)分布に従うランダム行列を生成している。