三次元グラフのプロットで行列計算を用いてfor文を回避するやり方（ガウス分布）

行列計算を用いたfor文回避

今回は授業資料で配られてなんだこれと思ったコードの数学的な意味と実装について考えてみました。

Python,Octave,Matlabなどfor文を回すより行列演算をした方が早い言語において、いかに行列計算を工夫してfor文を回避するかが重要になってきます。

　今回は、二次元ガウス分布のプロットをするときにいかに行列演算を使うかを考えます。
初めに一般の二変数関数のプロットを高速化させる方法を述べ、次に二次形式で場合さらに簡単なコードがかけることを述べ、最後に実際に実装します。

二次元ガウス分布

　二次元ガウス分布は
${ p(\boldsymbol{x|\boldsymbol{\mu}, \Sigma}) = \frac{1}{2\pi \sqrt{|\Sigma|}} exp \left\{ - \frac{1}{2} ( \boldsymbol{x} -\boldsymbol{\mu} )^T {\Sigma}^{-1} ( \boldsymbol{x} -\boldsymbol{\mu} ) \right\} }$
で表されます。
(ただし、 ${ \boldsymbol{x},\boldsymbol{\mu}は2次元のベクトル\Sigmaは2\times2の行列で\boldsymbol{x}= \left( \begin{array}{c} x\\ y\\ \end{array} \right) 。 }$ )

例えば、x=-3～3,y=-4～4の範囲でこの関数をプロットすることを考えると、
何も考えずに書くならばx=-3～3,y=-4～4で作られる平面を0.1の幅で区切って各点でのz座標をfor文を二回回して計算するということが考えられます。
そのfor文を避ける方法をこれから説明します。

for文を使う場合

今回の話の中心になるのが、指数部分の ${ \left\{ - \frac{1}{2} ( \boldsymbol{x} -\boldsymbol{\mu} ) {\Sigma}^{-1} ( \boldsymbol{x} -\boldsymbol{\mu} ) \right\} = \left\{ - \frac{1}{2} \left( \alpha_{11} x^2 + \alpha_{12} xy + \alpha_{21} xy + \alpha_{22} y^2 \right\} \right)…① }$ の部分です。
(ここで、 ${ {\Sigma}^{-1}= \begin{pmatrix} \alpha_{11} & \alpha_{12}\\ \alpha_{21} & \alpha_{22} \end{pmatrix} }$ と置いた)
ここを行列計算を用いてfor文なしで計算します。
for文を使うときは
${ \begin{array}{ccccc} (-3.0,-4.0) && (-2.9,-4.0) && \cdots && (2.9,4.0) && (3.0,4.0) \\ (-3.0,-3.9) && (-2.9,-3.9) && \cdots && (2.9,3.9) && (3.0,3.9) \\ \vdots && \vdots && \ddots && \vdots && \vdots \\ (3.0,3.9) && (2.9,3.9) && \cdots && (2.9,3.9) && (3.0,3.9) \\ (3.0,4.0) && (2.9,4.0) && \cdots && (2.9,4.0) && (3.0,4.0) \\　…② \end{array} }$
という61×81個の全部の(x,y)の組み合わせを①の式に順に代入していき対応する61×81個のz座標を求めてプロットすることになります。

行列計算を用いる場合

結論から先に書きますと、
${ \boldsymbol{X}(:)=(-3.0,-2.9, \cdots,2.9,3.0,-3.0,-2.9, \cdots,2.9,3.0, \cdots \cdots \cdots, -3.0,-2.9, \cdots,2.9,3.0) …③ \\ \boldsymbol{Y}(:)=(-4.0, \cdots,-4.0,-3.9, \cdots,-3.9, \cdots \cdots \cdots, 3.9, \cdots,3.9,4.0, \cdots,4.0)…④ }$
という二つの長さ61×81のベクトルを作り、X{:}を①の式のxにY{:}を①の式のyに代入すると、長さ61×81の各点でのz座標を表すベクトルを返してくれます。
ここで
X(:)は $(-3.0,-2.9, \cdots,2.9,3.0)$ というx座標を表す長さ61のベクトルを81個（y座標の数分）横に並べたベクトルで、長さは61×81です。
Y(:)は $(-4.0,-3.9, \cdots,3.9,4.0)$ というy座標を表す長さ81のベクトルの各成分を61個に増殖させたベクトルで、これも長さは61×81です。

代入後の式は
${ \left\{ - \frac{1}{2} \left( \alpha_{11} \boldsymbol{X}(:)*\boldsymbol{X}(:) + \alpha_{12} \boldsymbol{X}(:)*\boldsymbol{Y}(:) + \alpha_{21} \boldsymbol{X}(:)*\boldsymbol{Y}(:) + \alpha_{22} \boldsymbol{Y}(:)*\boldsymbol{Y}(:) \right\} \right)…⑤ }$
となります。ただし、*は内積ではなく要素ごとの積を表します。

このことからわかるように、これはガウス分布のプロットに限った話ではなくz=f(x,y)で表される一般の二変数関数の時にも、f(X{:},Y{:})というように代入すればプロットができることがわかります。

（なぜこれで上手くいくのかを注釈で説明して置きます。*1）

二次形式の時はさらに簡単に

ガウス分布の指数部分は二次形式のため、さらにコードを簡単にすることができます。
今からやることは、先ほどまでは①の右辺（二次形式をバラバラにしたもの）を使って計算していたところを①の左辺を使って計算することで二次形式のまま計算するということです。これによりコードを書くのが楽になるだけでなく、①の右辺をわざわざ求めるという計算の手間も省けます。

まず初めに、
${ \boldsymbol{A}= \left( \begin{array}{c} \boldsymbol{X}(:)\\ \boldsymbol{Y}(:)\\ \end{array} \right) }$
という2行61×81列の行列を用意します。
これと、 ${{\Sigma}^{-1}}$ との内積を取ると、
${ {\Sigma}^{-1} \cdot \boldsymbol{A}= \left( \begin{array}{c} \alpha_{11}\boldsymbol{X}(:)+\alpha_{12}\boldsymbol{Y}(:)\\ \alpha_{21}\boldsymbol{X}(:)+\alpha_{22}\boldsymbol{Y}(:)\\ \end{array} \right) }$

これと ${\boldsymbol{X}(:)}$ とで要素同士の積を計算すると、
${ \boldsymbol{A}*({\Sigma}^{-1} \cdot \boldsymbol{A}　)= \left( \begin{array}{c} \boldsymbol{X}(:)\\ \boldsymbol{Y}(:) \end{array} \right)* \left( \begin{array}{c} \alpha_{11}\boldsymbol{X}(:)+\alpha_{12}\boldsymbol{Y}(:)\\ \alpha_{21}\boldsymbol{X}(:)+\alpha_{22}\boldsymbol{Y}(:)\\ \end{array} \right) = \left( \begin{array}{c} \alpha_{11}\boldsymbol{X}(:)\boldsymbol{X}(:)+\alpha_{12}\boldsymbol{X}(:)\boldsymbol{Y}(:)\\ \alpha_{21}\boldsymbol{X}(:)\boldsymbol{Y}(:)+\alpha_{22}\boldsymbol{Y}(:)\boldsymbol{Y}(:)\\ \end{array} \right) }$
この行列に１行目と２行目を足し算するという関数を使えば、⑤式と一致する。
よって、 ${\boldsymbol{A}}$ を初めに作り、 ${\boldsymbol{A}*({\Sigma}^{-1} \cdot \boldsymbol{A})}$ を計算して列方向に足し算すれば二次元ガウス分布の指数部分が完成し、あとは定数倍したりexpにかませたりすればプロットができます。

実装

Python

まずpythonで書いてみました。（ただpythonならnumpyのmeshgridとか使えば今までのようなことをやらずに済むようです…）

# -*- coding: utf-8 -*-
"""
Created on Wed May 18 19:25:50 2016

"""

import numpy as np
import scipy as sc
import matplotlib.pyplot as plt
import math
from mpl_toolkits.mplot3d import Axes3D

mu=np.array([0,0],dtype=np.float)           #average　vector
sigma=np.array([[2,1],[1,2]],dtype=np.float)#covariance matrix
det=np.linalg.det(sigma)                    #determinant of sigma
sigmainv=np.linalg.inv(sigma)               #inverse of sigma
constant=1.0/(2*math.pi*math.sqrt(det))
x=np.array(np.arange(-3,3.1,0.1),dtype=np.float,ndmin=2) #ndmin=2にする事でベクトルでなく１行６１列の行列にする
y=np.array(np.arange(-4,4.1,0.1),dtype=np.float,ndmin=2) #そうすれば、転置したときに縦ベクトルとして扱える

X=np.tile(x,len(y[0])) 
Y=np.repeat(y.T,len(x[0]),axis=0) 
Y=Y.T
A=np.r_[X,Y]                                #列方向に連結
exponential=-0.5*np.sum((A*(sigmainv.dot(A))),axis=0)#np.sumで列方向に足し算
Z=constant*np.exp(exponential)

X=X.reshape((len(y[0]),len(x[0])))
Y=Y.reshape((len(y[0]),len(x[0])))
Z=Z.reshape((len(y[0]),len(x[0])))


fig = plt.figure()


ax=Axes3D(fig)
ax.plot_wireframe(X,Y,Z,rstride=1,cstride=1)
plt.show()

するとこんなグラフができます。
f:id:TakeToHazama:20160521115130p:plain

Octave

Octaveで書くと相当コードが短くて済みますね。非常に楽です（笑）

Mu=[0;0];
Sigma=[2 1;1 2];
x=[-3:0.1:3];#横ベクトル
y=[-4:0.1:4];#横ベクトル
X=repmat(x-Mu(1),[length(y) 1]);
Y=repmat(y-Mu(2),[1 length(x)]);
A=[X(:) Y(:)]';#ブログ内でのX(:)と違いコード内ではこれは縦ベクトルとして出てきます。
constant=1/(2*pi*sqrt(det(Sigma)));
Z=exp(-sum(A.*(inv(Sigma)*A))/2)*constant;#「.」が付くと要素動詞の演算を表す
z=reshape(Z,[length(y) length(x)]);
figure(1); clf; surf(x,y,z); view(45,60)
figure(2);clf;contour(x,y,z);#真上から見た図

figure(1)
f:id:TakeToHazama:20160520092009p:plain
figure(2)
f:id:TakeToHazama:20160520091341p:plain

*1:これらのベクトルを①の式の右辺に代入しますと入れますと、「X(:)の1番目の要素とY(:)の1番目の要素（すなわち行列②の1行1列要素）を取ってきて計算したものを結果のベクトルの1番目の要素に格納」 …… 「X(:)の61番目の要素とY(:)の61番目の要素（すなわち行列②の1行61列要素）を取ってきて計算したものを結果のベクトルの61番目の要素に格納」「X(:)の62番目の要素とY(:)の62番目の要素（すなわち行列②の2行1列要素）を取ってきて計算したものを結果のベクトルの61番目の要素に格納」 …… 「X(:)の122番目の要素とY(:)の122番目の要素（すなわち行列②の2行61列要素）を取ってきて計算したものを結果のベクトルの61番目の要素に格納」 ………… 「X(:)の61×80+1番目の要素とY(:)の61×80+1番目の要素（すなわち行列②の81行1列要素）を取ってきて計算したものを結果のベクトルの61×80+1番目の要素に格納」 …… 「X(:)の61×81番目の要素とY(:)の61×81番目の要素（すなわち行列②の81行61列要素）を取ってきて計算したものを結果のベクトルの61×81番目の要素に格納」という作業が行われ、欲しかったzの値が求まります。

隙間ブログ

情報系M1のメモ。機械学習、画像処理、考えた事などなど。