Python で行列演算 = numpy

python で行列演算をしたい.

ちょうど待ち行列MAP/G/1の系内客数分布を数値計算しようとしていて, どうせなので汎用性高い物にしようと思ったの動機.
まぁありきたりな話, 新しいライブラリを使うと大抵最初はまる.

・numpy のインストール
と思ったらすでにnumpy は入っていた. バージョンは1.2.
古いようなので, easy-install でupgrade する.

% sudo easy_install --upgrade numpy

とここで, 鬼のようにwarning と error っぽいメッセージが. 萎える.

・numpy を使う

matrixってのが行列の定義のようだ.

>>> C_matrix = [
...  [-0.2,0.2],
...  [0.0, -0.2]
...  ]
>>> C_matrix = matrix(C_matrix)
>>> D_matrix = [
...     [0,0],
...     [0.2,0]
...     ]
>>> 
>>> C_matrix = matrix(C_matrix)
>>> D_matrix = matrix(D_matrix)
>>> C_matrix 
matrix([[-0.2,  0.2],
        [ 0. , -0.2]])
>>> D_matrix
matrix([[ 0. ,  0. ],
        [ 0.2,  0. ]])

のように定義できた. array ってのもあるが, これは完全に転置っぽいので直感的でないため, やめる.
(注) matrix で[[1,2],[3,4]] は array では [[1,3],[2,4]] と列で定義するとうこと.

演算もらくらくで便利そうだ.

>>> E_VEC = matrix([[1],]*2)
>>> E_VEC
matrix([[1],
        [1]])

>>> C_matrix + D_matrix
matrix([[-0.2,  0.2],
        [ 0.2, -0.2]])
>>> (C_matrix + D_matrix) * E_VEC
matrix([[ 0.],
        [ 0.]])
>>> dot((C_matrix + D_matrix), E_VEC)
matrix([[ 0.],
        [ 0.]])

* は各成分の掛け算になるらしい. dot()で行列のかけざん.


で, solve ってのが連立方程式の解のよう.
とりあえず, 行ベクトルPI * (C + D) = 0, PI*e = 0 を解いてCおよびDで構成される, (C+Dをジェネレータとして持つ)マルコフ連鎖の定常状態確率を解くことにする.



ここではまった.


solve は正方行列でないと扱えないらしく, おそらくフルランクでないと無理だろう.
とりあえず, C+Dは明らかにフルランクではなく, PI*e = 0と組み合わせることでようやく解ける.
どうしようかと思ったけどとりあえず暫定的に, C+Dの一列を抜いてeへ替える. そして右辺は[0,...,0,1]とおく.

以下のようになる.

from numpy.linalg import solve

def calculate_pi():
    sumC_D = C_matrix + D_matrix 
    pi_solver = []
    for column in sumC_D.tolist():
        column.append(1)
        pi_solver.append(column[1:])

    zeros_plus1 = [0] * (M_SIZE -1) + [1]    
    pi_solver = mat(pi_solver)
    zeros = mat(zeros_plus1)

あまりに煩雑すぎるが, 特にメソッドもなさそうなので暫定的に...
pi_solverは最左行を省く. (まぁこれがのちのち問題だと思う)
zeros はさっきかいた[0,...,0, 1]の行ベクトル.

で, 早速solve.

>>> pipi = solve(t,zeros)
Traceback (most recent call last):
  File "", line 1, in 
  File "/Library/Python/2.6/site-packages/numpy-1.5.1-py2.6-macosx-10.6-universal.egg/numpy/linalg/linalg.py", line 316, in solve
    raise LinAlgError, 'Incompatible dimensions'
numpy.linalg.linalg.LinAlgError: Incompatible dimensions

はーなんでやーと思って, 色々やってみると, 同じことを matrix のかわりに arrayなら通る.
で, もしかたしらとちょっと不安だったが, どうやら

行* 行列 = 行

はだめで

行列 * 列 = 列

しかやってくれないぽい.

ふざけんなよーと思って, またもや暫定的に

def calculate_pi():
    sumC_D = C_matrix + D_matrix 
    pi_solver = []
    for column in sumC_D.tolist():
        column.append(1)
        pi_solver.append(column[1:])
    zeros_plus1 = [0] * (M_SIZE -1) + [1]
        
    pi_solver = mat(pi_solver)
    zeros_plus1 = mat(zeros_plus1)
    pi_vec = solve(pi_solver.T,zeros_plus1.T)

    return  pi_vec.T

とした. Tで転置らしい.


これは先が思いやられるな...