Matplotlib でお絵かきアプリっぽいものを作る

Matplotlib でお絵かきアプリっぽいものを作ります。

1. matplotlib のバージョン

GUIAPIの仕様は変わる可能性が高いため、matplotlib のバージョンを記載しておきます。
バージョン 2.2.2 で動作確認済みです。以下のようにするとバージョンの確認ができます。

import matplotlib
print(matplotlib.__version__)
# Out: '2.2.2'

2. 仕様

アプリの仕様は以下とします。

No 状態 動作
1 マウス左ボタン プレス ペンを下ろす
2 マウス左ボタン リリース ペンを上げて、軌跡を消去
3 マウス左ボタンをプレス中にマウスを動かす ペンで軌跡を描画
4 マウス左ボタンをリリース中にマウスを動かす なにもしない

簡単にいってしまうと、マウス左ボタンを押している間だけマウスの軌跡を描画するアプリです。ボタンを離すと軌跡は消えます。一筆書きをしている間は軌跡をみることができます。

3. コード

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, Arrow

fig, ax = plt.subplots()

being_pressed = False
pre_x = None
pre_y = None

def on_motion(event):
    global pre_x, pre_y
    if None in (event.xdata, event.ydata): # マウスが画面の外にいるときは何もしない
        #print('x=%d, y=%d, xdata=%f, ydata=%f' %(event.x, event.y, event.xdata, event.ydata))
        return
    if being_pressed:
        ax.add_patch(Circle((event.xdata, event.ydata), 0.01, color='r'))
        if None not in (pre_x, pre_y): # pre_x と pre_y が存在するとき矢印を描画
            ax.add_patch(
                Arrow(pre_x, pre_y, event.xdata-pre_x, event.ydata-pre_y,
                      linestyle='solid', color='blue', width=0.02))
        pre_x = event.xdata
        pre_y = event.ydata
        fig.canvas.draw() # 必須!

def on_button_press(event):
    print('button pressed')
    global being_pressed
    being_pressed = True

def on_button_release(event):
    print('button released')
    global being_pressed, pre_x, pre_y
    being_pressed = False
    pre_x = None
    pre_y = None
    # reset patches
    for patch in reversed(ax.patches):
        patch.remove()
    fig.canvas.draw() # 必須!


cid_1 = fig.canvas.mpl_connect('motion_notify_event', on_motion)
cid_2 = fig.canvas.mpl_connect('button_press_event', on_button_press)
cid_3 = fig.canvas.mpl_connect('button_release_event', on_button_release)

plt.show()

4. コードの説明

マウスの左ボタンを押しているかどうかを being_pressed 変数で管理することとしています。pre_x, pre_y はボタンが押されているときの、マウスを動かす前のマウスの座標を保存するために定義しています。(pre_x, pre_y はボタンが押されてないときは存在しないので None で初期化します)

ボタンが押されたときに独自関数(on_button_press関数)を呼んでほしいため

fig.canvas.mpl_connect('button_press_event', on_button_press)

としています。これにより canvas で button_press_event が発生すると on_button_press 関数が呼ばれるようになります。同様に、 ボタンが離されたとき(button_release_event イベントが発生したとき)には on_button_release、 マウスが動いたとき(motion_notify_event イベントが発生したとき)には on_motion 関数が呼ばれるように設定します。

話はそれますが、matplotlib で他にどのようなイベントが提供されているか興味のある方は、以下を参照してください。
https://matplotlib.org/users/event_handling.html より抜粋

Event name Class and description
'button_press_event' MouseEvent - mouse button is pressed
'button_release_event' MouseEvent - mouse button is released
'draw_event' DrawEvent - canvas draw (but before screen update)
'key_press_event' KeyEvent - key is pressed
'key_release_event' KeyEvent - key is released
'motion_notify_event' MouseEvent - mouse motion
'pick_event' PickEvent - an object in the canvas is selected
'resize_event' ResizeEvent - figure canvas is resized
'scroll_event' MouseEvent - mouse scroll wheel is rolled
'figure_enter_event' LocationEvent - mouse enters a new figure
'figure_leave_event' LocationEvent - mouse leaves a figure
'axes_enter_event' LocationEvent - mouse enters a new axes
'axes_leave_event' LocationEvent - mouse leaves an axes

話を戻して、以下のコードで現在のマウスのx, y座標に丸を描画します。0.01 は 丸の大きさです。何度か試して 0.01 が良さそうだったのでそうしています。

ax.add_patch(Circle((event.xdata, event.ydata), 0.01, color='r'))

以下のコードで、前のマウスのx, y座標と現在の座標に矢印を描画します。0.02 は線の太さです。これも試した結果これくらいがよさそうだったのでそうしています。

ax.add_patch(
    Arrow(pre_x, pre_y, event.xdata-pre_x, event.ydata-pre_y,
    linestyle='solid', color='blue', width=0.02))

これら 0.01 や 0.02 を計算でだせるようにするのは私の宿題としておきます。

5. 実行結果

f:id:keimina:20190705002816g:plain
matplotlibによるお絵かきアプリっぽいもの

6. まとめ・考察

matploblib でお絵かきアプリっぽいものを作りました。
マウスがクリックされたら〜するとか、マウスが動いたら〜するなどいろいろ応用できるのではないかと思います。筆跡を時系列データとして機械学習して分類するのもおもしろそうです、夏休みの自由研究にもできそうではないでしょうか。

7. 今後の予定

個人的にはこれを使って軌跡に囲まれた領域のデータの点だけを取得して、それらのデータの可視化したいです(実はもうしたけどそれは次回以降、書きます)。

8. 宿題

matplotlib の Arrow の太さとか Circle の大きさを計算できるように仕様を把握する

それでは、おやすみなさい。

勾配降下法により PCA(主成分分析)の第2主成分を数値解析的に求め、PCAについて理解を深める

1. はじめに

前回記事で、第2主成分は「第1主成分と直交する&分散が最大となる結合係数を求めるということだと理解しています」と書きました。本記事は、その理解で正しいか検証したときの記録です。検証した結果、その理解で正しいということがわかりました。

2. 前回記事

keimina.hatenablog.jp

3. サンプルデータについて

前回同様、使用するサンプルデータは 3次元の特徴量をもつ3つのデータとします。
現実にこのサイズのデータで PCA を行うことは稀かと思われますが、PCA の理解には十分だと思います。

4. 実装の方針

第2主成分が第1主成分と直交するというのは第1主成分をベクトル Wa 、 第2主成分をベクトル Wb とすると Wa と Wb の内積が 0 になるということです。 Wa と Wb が以下の時

Wa = [Wa1, Wa2, Wa3]
Wb = [Wb1, Wb2, Wb3]

内積は以下のように定義されます。

inner_product = Wa1*Wa2 + Wa2*Wb2 + Wa3*Wb3


Wa と Wb が直交しているとき、内積が 0 なので

Wa1*Wb1 + Wa2*Wb2 + Wa3*Wb3 = 0

を Wa, Wb は満たします。ちなみに前回記事で第1主成分 Wa は求めることができましたので、 Wa1, Wa2, Wa3 には数値が入ることがわかります。上記の式を変形して Wb1 について解くと

Wb1 = -(Wa2*Wb2 + Wa3*Wb3) / Wa1

となります。 これを Wb の Wb1 に代入すると

Wb = [-(Wa2*Wb2 + Wa3*Wb3) / Wa1, Wb2, Wb3]

と書くことができます(細かいですが Wa1が 0 でないときに限ります)。Wa と Wb が直交するという条件があるので、それを反映させたときの Wb がこの式になります。この Wb をサンプルデータに適用し適用後のデータ間の分散が最大になる Wb が第2主成分であると理解していますので、その Wb と sklearn.decomposition.PCA を使用して求めた第2主成分を比較し、同じ方向を向いているか検証します。

※当たり前ですが Wb を求めるときに Wa2, Wa3 すなわち Wa の情報が必要です。このように依存関係がありますのでコーディング時に注意します。

5. ソースコード

ソースコードについてですが、前回は行列のシンボルとして sympy.MatrixSymbol を使用していたのですが、今回は sympy.Symbol と np.array を使うことにしました。理由は、前回記事のソースコードでは、Matrixによるキャストが大量発生してしまい、期待していたほど可読性が上がるわけではなかったからです。(MatrixSymbol を使用すると行列の各要素の変数を自動的に作成してくれて便利だったのですが残念です)そのため本記事では MatrixSymbol の代わりに、 Symbol と np.array を使うことにしました。

# -*- coding:utf-8 -*-
# MatrixSymbol を使わないで np.array で統一した & Second component 追加したバージョン

import numpy as np
from sympy import symbols, diff, Eq, solve
from sklearn.decomposition import PCA
from itertools import chain

# 1行目が観測データ1、2行目が観測データ2、3行目が観測データ3
# 1列目が特徴量1、2列目が特徴量2、3列目が特徴量3 とする
data = np.array([[3,5,7],[11,13,17],[19,23,29]], dtype=np.float64)

################################################################
# 独自手法で
# first component (第一主成分)を求めるための処理
################################################################

# 観測データ入力用の 3x3 の行列のシンボルを作成
# 1列目が特徴量1、2列目が特徴量2、3列目が特徴量3
Data = np.array([[symbols('Data%d%d'%(i, j)) for j in range(1, 4)] for i in range(1, 4)])

# first component 用のシンボルを作成(サイズ=3x1)
Wa = np.array([[symbols('Wa%d1'%i)] for i in range(1, 4)])

# second component 用のシンボルを作成(サイズ=3x1)
Wb = np.array([[symbols('Wb%d1'%i)] for i in range(1, 4)])

# second component は first component と直行する(内積=0となる)条件を使用し、1変数を置き換える
inner_prod = Wa.T.dot(Wb)[0][0]
# solve 関数を使用して Wa と 直行する Wb の Wb[0][0] を求める
# solve を使用するとリストが戻るため、要素を取り出すため インデックスを指定して取り出す
new_Wb11 = solve(Eq(inner_prod), Wb[0, 0])[0]
# 全要素の数式に含まれる Wb11 を new_Wb11 と置き換える関数を定義
replace_Wb11 = np.vectorize(lambda x: x.subs({Wb[0, 0]: new_Wb11}))

# Data を以下のように変換し、次元を1次元に削減する
Y = Data.dot(Wa) # shape=(3, 1)
Y2 = Data.dot(Wb)

# Y の特徴量の分散を求める
Ymean = np.mean(Y, axis=0, keepdims=True)
Yhat = (Y - Ymean)   # 重心を (0, 0, 0) にする、shape=(3, 1)
V = Yhat.T.dot(Yhat) # 分散を求める、 shape=(1, 1)

# Y2 特徴量の分散を求める
Y2mean = np.mean(Y2, axis=0, keepdims=True)
Y2hat = (Y2 - Y2mean)   # 重心を (0, 0, 0) にする、shape=(3, 1)
V2 = Y2hat.T.dot(Y2hat) # 分散を求める、 shape=(1, 1)

# 分散 V が最大になるようにする
# ここでは、数値解析的に勾配降下法を使用しVの最大値を求めるため微分する
diff_V = np.vectorize(lambda x: diff(V[0, 0], x))
diff_V2 = np.vectorize(lambda x: diff(V2[0, 0], x))

dVdWa = diff_V(Wa) #shape=(3, 1)
dV2dWb = diff_V2(Wb) #shape=(3, 1)

# first components の値を保持する変数を用意しランダムな数値で初期化
np.random.seed(51) # 他の人がやっても同じ結果になるようにする
Wa_val = np.random.random((3,1))
Wb_val = np.random.random((3,1))

# ハイパーパラメータ(学習係数)を定義
b = 0.05

# 勾配を求める関数を定義
dVdWa_func = lambda Wa_val: np.array([w[0].subs(dict(zip(chain(Data.flatten(), Wa.flatten()),
                                       chain(data.flatten(), Wa_val.flatten()))))
                                      for w in dVdWa], dtype=np.float64).reshape((3, 1))
dV2dWb_func = lambda Wb_val, Wa_val: np.array([w[0].subs(dict(zip(chain(Data.flatten(), Wb.flatten(), Wa.flatten()),
                                       chain(data.flatten(), Wb_val.flatten(), Wa_val.flatten()))))
                                               for w in replace_Wb11(dV2dWb)], dtype=np.float64).reshape((3, 1))

print('===========================================')
print('================ 第1主成分 ================')
print('===========================================')
# 第1主成分を学習する
for i in range(101):
    # 勾配を求める
    dVdWa_val = dVdWa_func(Wa_val)
    # ハイパーパラメータ b をかけて Wa_val を更新する
    Wa_val += b*dVdWa_val
    # 正規化する
    Wa_val /= np.linalg.norm(Wa_val)
    if i%10==0:
        print('================ 第1主成分 学習 i = ', i, "================")
        v = V[0, 0].subs(dict(zip(chain(Data.flatten(), Wa.flatten()),
                                       chain(data.flatten(), Wa_val.flatten()))))
        print("分散V = ", v)
        print("Wa = ", Wa_val.T)



print('===========================================')
print('================ 第2主成分 ================')
print('===========================================')

# 第2主成分を学習する
calc_Wb = np.vectorize(lambda x: x.subs(dict(zip(chain(Wb.flatten(), Wa.flatten()),
                                                 chain(Wb_val.flatten(), Wa_val.flatten())))))
for i in range(101):
    # Wb を Wa と直交するようにする
    Wb_val = calc_Wb(replace_Wb11(Wb)).astype(np.float64)
    # 勾配を求める
    dV2dW_val = dV2dWb_func(Wb_val, Wa_val)
    # ハイパーパラメータ b をかけて Wb_val を更新する
    Wb_val += b*dV2dW_val
    # 正規化する
    Wb_val /= np.linalg.norm(Wb_val)
    if i%10==0:
        print('================ 第2主成分 学習 i = ', i, "================")
        v2 = V2[0, 0].subs(dict(zip(chain(Data.flatten(), Wb.flatten(), Wa.flatten()),
                                       chain(data.flatten(), Wb_val.flatten(), Wa_val.flatten()))))
        print("分散V2 = ", v2)
        print("Wb = ", Wb_val.T)

################################################################
# sklearn.decomposition.PCA で
# first component (第一主成分)を求めるための処理
################################################################
# solving by sklearn.decomposition.PCA
pca = PCA(n_components=3) # 1-3 components を求める
pca.fit(data) # 学習する

################################################################
# 独自手法 と sklearn.decomposition.PCA で
# 求めた first component (第一主成分)を表示するための処理
################################################################
print('\n================ 検証 ================')
# 独自手法 で求めた first component を表示(比較しやすいようにflattendを使用)
print('--- 独自手法 ---')
print('第1主成分 = ', Wa_val.flatten())
print('第2主成分 = ', Wb_val.flatten())

# sklearn.decomposition.PCA で求めた frist component を表示
w1_pca, w2_pca, w3_pca = pca.components_
print('--- sklearn.decomposition.PCA ---')
print('第1主成分 = ', w1_pca)
print('第2主成分 = ', w2_pca)

# sklearn.decomposition.PCA で求めた frist component と
# コサイン類似度を求める、w1_pca と Wa_val はノルムが1なので内積を
# とるだけで cosθ が求まる
print('--- コサイン類似度 ---')
cos_theta1 = w1_pca.dot(Wa_val)
cos_theta2 = w2_pca.dot(Wb_val)
print(' 第1主成分 cosθ = {} ({}°)'.format(*cos_theta1, *np.rad2deg(np.arccos(cos_theta1))))
print(' 第2主成分 cosθ = {} ({}°)'.format(*cos_theta2, *np.rad2deg(np.arccos(cos_theta2))))

# 第1主成分と第2主成分における分散の比率
print('--- 第1主成分と第2主成分における分散の比率(sklearn.decomposition.PCA) ---')
print('pca.explained_variance_ratio_[0] / pca.explained_variance_ratio_[1] \n =', pca.explained_variance_ratio_[0] / pca.explained_variance_ratio_[1])
print('--- 第1主成分と第2主成分における分散の比率(独自手法) ---')
print('V / V2 \n = ', v / v2)

6. ソースコードの説明

ソースコードの説明を以下に記載します。

diff_V = np.vectorize(lambda x: diff(V[0, 0], x))
diff_V2 = np.vectorize(lambda x: diff(V2[0, 0], x))

dVdWa = diff_V(Wa) #shape=(3, 1)
dV2dWb = diff_V2(Wb) #shape=(3, 1)

上記は、V の Wa に関する微分と V2 の Wb に関する微分を求めています。微分を行う関数である diff 関数をブロードキャストできるようにするため np.vectorize 関数を使用しています。

new_Wb11 = solve(Eq(inner_prod), Wb[0, 0])[0]

上記は、 Eq(inner_prod) で inner_prod = 0 となる数式を構築し、 solve(数式, Wb[0, 0]) で 数式をWb[0, 0] について解いています。solve 関数は解が複数ある場合に備えてリストで解を戻します。今回は、数式をみれば解が一つであることは自明なので、インデックスに 0 を指定してリストから要素をを取り出しています。ここで求めた new_Wb11 を使用して Wb の要素 Wb_11 を置き換えます。

replace_Wb11 = np.vectorize(lambda x: x.subs({Wb[0, 0]: new_Wb11}))

上記は、全要素の数式に含まれる Wb11 を new_Wb11 と置き換える関数を定義しています。

dVdWa_func = lambda Wa_val: np.array([w[0].subs(dict(zip(chain(Data.flatten(), Wa.flatten()),
                                       chain(data.flatten(), Wa_val.flatten()))))
                                      for w in dVdWa], dtype=np.float64).reshape((3, 1))

上記は Wa の勾配を求める関数を定義しています。シンボルに値を代入するために、key にシンボルを value に シンボルに対応する値を入れた辞書を渡しています。(TensorFlow がわかる人には placeholder に値を代入しているというと伝わるかもしれません。)また、前回記事では、微分をするのにjacobian を使用していましたが今回は for文で回して、各要素を微分しています。

dV2dWb_func = lambda Wb_val, Wa_val: np.array([w[0].subs(dict(zip(chain(Data.flatten(), Wb.flatten(), Wa.flatten()),
                                       chain(data.flatten(), Wb_val.flatten(), Wa_val.flatten()))))
                                               for w in replace_Wb11(dV2dWb)], dtype=np.float64).reshape((3, 1))

上記は Wb の勾配を求める関数を定義しています。 Wb を求めるときに Wa の情報が必要でコーディング時に注意と冒頭で説明しましたが、 上記がそのコードになります。 dict の中身をみると Wb だけでなく Wa の情報も使用していることがわかります。 また、 dV2dWb に直交条件を反映させるため replace_Wb11 関数を適用しています。これにより、 Wa に直交する Wb における分散の勾配を求めることができます。

7. 出力結果

出力結果は以下になります。

===========================================
================ 第1主成分 ================
===========================================
================ 第1主成分 学習 i =  0 ================
分散V =  532.873871148815
Wa =  [[0.50212335 0.54326168 0.67285874]]
================ 第1主成分 学習 i =  10 ================
分散V =  533.003128892760
Wa =  [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i =  20 ================
分散V =  533.003128892760
Wa =  [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i =  30 ================
分散V =  533.003128892760
Wa =  [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i =  40 ================
分散V =  533.003128892760
Wa =  [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i =  50 ================
分散V =  533.003128892760
Wa =  [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i =  60 ================
分散V =  533.003128892760
Wa =  [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i =  70 ================
分散V =  533.003128892760
Wa =  [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i =  80 ================
分散V =  533.003128892760
Wa =  [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i =  90 ================
分散V =  533.003128892760
Wa =  [[0.48958803 0.55232128 0.67471829]]
================ 第1主成分 学習 i =  100 ================
分散V =  533.003128892760
Wa =  [[0.48958803 0.55232128 0.67471829]]
===========================================
================ 第2主成分 ================
===========================================
================ 第2主成分 学習 i =  0 ================
分散V2 =  0.282286132046230
Wb =  [[-0.85629139  0.15854369  0.49155768]]
================ 第2主成分 学習 i =  10 ================
分散V2 =  0.303318086521561
Wb =  [[-0.86855083  0.24061133  0.43327317]]
================ 第2主成分 学習 i =  20 ================
分散V2 =  0.315596549516175
Wb =  [[-0.87190187  0.30165001  0.38573876]]
================ 第2主成分 学習 i =  30 ================
分散V2 =  0.322411743461676
Wb =  [[-0.87100083  0.34594949  0.3488216 ]]
================ 第2主成分 学習 i =  40 ================
分散V2 =  0.326088827550761
Wb =  [[-0.86849373  0.37775612  0.32096564]]
================ 第2主成分 学習 i =  50 ================
分散V2 =  0.328042467143946
Wb =  [[-0.8656847   0.40051315  0.30029855]]
================ 第2主成分 学習 i =  60 ================
分散V2 =  0.329071954623274
Wb =  [[-0.86312883  0.41679385  0.28511667]]
================ 第2主成分 学習 i =  70 ================
分散V2 =  0.329612106418076
Wb =  [[-0.86100686  0.42845606  0.27403029]]
================ 第2主成分 学習 i =  80 ================
分散V2 =  0.329894869159100
Wb =  [[-0.85933015  0.43682358  0.26596403]]
================ 第2主成分 学習 i =  90 ================
分散V2 =  0.330042715591538
Wb =  [[-0.85804391  0.44283655  0.26010852]]
================ 第2主成分 学習 i =  100 ================
分散V2 =  0.330119970975593
Wb =  [[-0.85707561  0.44716326  0.25586408]]

================ 検証 ================
--- 独自手法 ---
第1主成分 =  [0.48958803 0.55232128 0.67471829]
第2主成分 =  [-0.85707561  0.44716326  0.25586408]
--- sklearn.decomposition.PCA ---
第1主成分 =  [0.48958803 0.55232128 0.67471829]
第2主成分 =  [ 0.85440069 -0.45835872 -0.24475855]
--- コサイン類似度 ---
 第1主成分 cosθ = 1.0 (0.0°)
 第2主成分 cosθ = -0.9998720868190456 (179.0835687095191°)
--- 第1主成分と第2主成分における分散の比率(sklearn.decomposition.PCA) ---
pca.explained_variance_ratio_[0] / pca.explained_variance_ratio_[1] 
 = 1614.1609966447247
--- 第1主成分と第2主成分における分散の比率(独自手法) ---
V / V2 
 =  1614.57402082519

8. 考察

出力結果を見ると、第1主成分は前回と同じで、問題なく動作していることが確認できます。第2主成分については i = 0 のとき分散が 0.28 で i = 100 のとき 0.33 と増加していることがわかりますので、学習はうまく進んだといえます。検証の箇所で sklearn.decomposition.PCA で求めた第2主成分とのコサイン類似度を出していますがおおよそ cosθ = -1 となりました。これは、 ベクトルが反対方向を向いているということで、 分散の計算でこのベクトルの符号は打ち消されますので、この方向でも反対方向でも分散が最大なのは変わりません。ちなみにどちらの方向になるかは Wb_val の初期状態できまります。最後に、 第2主成分の分散の値が 0.33 と少なく思えるため、分散の計算があっているか確かめるのに、 第1主成分の分散と第2主成分の分散の比率を sklearn.decomposition.PCA で導出したものと比較しています。これにより、分散の比率は 1614 とほとんど同じですので、分散の計算はあっていることがわかります。このサンプルデータにおいては第2主成分は分散が少ないですので、ほとんど意味がない情報といえそうです。以上にて、第2主成分について私の理解について問題ないことが確認できました。最後に、 3x3 のサンプルデータにおける最後の主成分の第3主成分ですが、これは第1主成分と第2主成分に直行するノルム1のベクトルですので連立方程式をとくと解が求まります。連立方程式を解いて第3主成分( Wc とします)を求めるコードを以下に記載し終わりとします。

9. 第3主成分を求めるコード

連立方程式を解いて第3主成分( Wc とします)を求めるコードを以下に記載し終わりとします。

# 連立方程式を解いて 最後の主成分である第3主成分を求める
Wc = np.array([[symbols('Wc%d1'%i)] for i in range(1, 4)])
wc_wa_inner = Wc.T.dot(Wa)[0, 0]
wc_wb_inner = Wc.T.dot(Wb)[0, 0]
wc_norm = Wc.T.dot(Wc)[0, 0]

eq1 = Eq(wc_wa_inner)
eq2 = Eq(wc_wb_inner)
eq3 = Eq(wc_norm - 1)

# eq1 かつ eq2 かつ eq3 の条件を満たす Wc の各要素の変数についてシンボリックに解く
kais = np.array(solve([eq1, eq2, eq3], *Wc.flatten()))

# 求めたシンボリックな解に Wa と Wb の値を代入して実際の値を求める関数を定義
subs = np.vectorize(lambda x: x.subs(dict(zip(chain(Wa.flatten(), Wb.flatten()),chain(Wa_val.flatten(), Wb_val.flatten())))))

# 解を表示
print('第3主成分 = ', subs(kais)[0])

# 連立方程式で求めた第3主成分と sklearn.decomposition.PCA で求めたものとの内積を求めて値が 1 に近いこと(方向が同じであること)を確認する
cos_theta3 = float(np.inner(subs(kais)[0], w3_pca))
print('--- コサイン類似度 ---')
print(' 第3主成分 cosθ = {} ({}°)'.format(cos_theta3, np.rad2deg(np.arccos(cos_theta3))))

10. 出力結果(第3主成分を求めるコード)

3主成分 =  [0.160390052085319 0.703552574067647 -0.692306873225203]
--- コサイン類似度 ---
 第3主成分 cosθ = 0.9998720868190454 (0.9164312904816843°)

第3主成分についても、コサイン類似度をみると、少し誤差はあるものの、sklearn.decomposition.PCA のものとほとんど同じであることが確認できました。誤差を少なくするには、勾配降下法で第1主成分と第2主成分を求めるときに学習率を大きくするなどの工夫が必要です。

11. まとめ

3x3 のサンプルデータについて第1主成分と第2主成分を数値解析的に求め、第3主成分を代数的に求めました。それぞれ sklearn.decomposition.PCA で求めた値ととほぼ同じ値になることを確認しました。PCAの理解を深めました。

以上です。おやすみなさい。

点を移動させる行列と、その固有ベクトルについて可視化して理解を深める

1. はじめに

2x2の行列を点(x, y)に適応すると、新しい点に移動します。本記事は、この行列による点の変換を可視化し行列の理解を深める記録です。本題ではありませんが、行列の固有ベクトルも表示し、固有ベクトルについても理解を深めました。

2. 固有ベクトルについて

固有ベクトルは本記事の本題ではありませんが、固有ベクトルとは何なのかというと、行列を適応しても定数倍の効果しか及ぼさないベクトルをその行列の固有ベクトルと呼ぶようです。詳細は以下を参照してください。

ja.wikipedia.org

3. 仕様

行列による点の変換を可視化しするために考えた仕様は以下になります。

  1. 行列適用前の格子状の点を表示する(グレー)
  2. 行列適用後の点を表示する(赤)
  3. 固有ベクトルを表示する(青)
  4. 固有ベクトルが存在しないとき(固有ベクトルの要素に複素数が含まれるとき)は固有ベクトルを表示しない
  5. 点の移動量を矢印で表示する(グレー)
  6. 変換行列をスライダで調整できるようにする

4. ソースコード

仕様に従い実装しました。ソースコードは以下の通りです。

# xy座標上の点 A(x1, y1) を点B(x2, y2)に移動(写像)する行列 M を考えます。
# 点B は行列 M と 点A を用いて以下のように表せます。
# (a, b は点A と点B を表す縦ベクトルとします)
# b = M・a
import numpy as np
from sympy import Matrix, MatrixSymbol

# 描画関連
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider, Button
from matplotlib.patches import Arrow

# デバッグ用
# plt.ion()

# 行列適用前の点
xs, ys = np.mgrid[-10:11:2, -10:11:2] / 10

# 行列演算の定義
N = 2
M = MatrixSymbol('M', N, N)
a = MatrixSymbol('a', N, 1)
b = M*a

# 行列適用後の点
def calc_point(mat, point):
    return b.subs({M:mat, a:point}).doit()

# スライダーの axis を作成
fig, (ax, ax_button, *ax_sldrs) = plt.subplots(6, 1, gridspec_kw = {'height_ratios':[28, 1, 1, 1, 1, 1]})

# スライダーの名前を定義
snames = ['p11', 'p12', 'p21', 'p22']

# axis にボタンを作成
button = Button(ax_button, 'show/hide vector')

# axis にスライダーを作成
sldrs = [Slider(ax_sldr, sname, -1.0, 1.0, valinit=0.0) for ax_sldr, sname in zip(ax_sldrs, snames)]

# 行列 M は
# M = [[sldrs[0].val, sldrs[1].val],
#      [sldrs[2].val, sldrs[3].val]]
# とする

# 矢印の表示の有無
patches_show_state = True

def init_sldrs():
    # スライダの初期値を設定する
    # 単位行列
    sldrs[0].set_val(1.0)
    sldrs[1].set_val(0.0)
    sldrs[2].set_val(0.0)
    sldrs[3].set_val(1.0)

def calc_points():
    points = []
    # スライダから値を取得
    p11, p12 = sldrs[0].val, sldrs[1].val
    p21, p22 = sldrs[2].val, sldrs[3].val
    for x, y in zip(xs.flatten(), ys.flatten()):
        mat = Matrix([[p11, p12], [p21, p22]])
        # 行列適用前の点の座標を取得
        point1 = Matrix([[x], [y]])
        # 行列を適応
        point2 = calc_point(mat, point1)
        # 行列適用後の座標をリストに保存
        points.append(point2)
    # 保存したリストを x座標のリスト、 y座標のリストに変換する
    x, y = list(zip(*points))
    return x, y

def init_ax():
    # 行列適用後の点を描画する
    x, y = calc_points()
    ax_line, = ax.plot(x, y, color='r', marker='.', linestyle='None')
    # オリジナルの点も描画する
    ax.plot(xs, ys, color='k', alpha=0.1, marker='.', linestyle='None')
    # 矢印も描画する
    for x1, y1, x2, y2 in zip(xs.flatten(), ys.flatten(), x, y):
        # Arrow(x, y, dx, dy) で dx と dy を float にしないとエラーになる
        ax.add_patch(Arrow(x1, y1, float(x2-x1), float(y2-y1)))
    # 虚数が含まれない場合、固有ベクトルを表示する
    p11, p12 = sldrs[0].val, sldrs[1].val
    p21, p22 = sldrs[2].val, sldrs[3].val
    w, v = np.linalg.eig([[p11, p12],[p21, p22]])
    if np.all(np.isreal(v)):
        ax.add_patch(Arrow(0, 0, *(2*v[:,0]),  linestyle='solid', color='blue', width=0.02))
        ax.add_patch(Arrow(0, 0, *(2*v[:,1]),  linestyle='solid', color='blue', width=0.02))

    return ax_line

def update_ax(ax_line):
    # 行列適用後の点を描画する
    x, y = calc_points()
    ax_line.set_xdata(x)
    ax_line.set_ydata(y)
    update_patches()

def update_patches():
    x, y = calc_points()
    # 矢印を更新
    # https://stackoverflow.com/questions/49791848/matplotlib-remove-all-patches-from-figure?noredirect=1&lq=1
    for patch in reversed(ax.patches):
        patch.remove()
    if patches_show_state:
        for x1, y1, x2, y2 in zip(xs.flatten(), ys.flatten(), x, y):
            # Arrow(x, y, dx, dy) で dx と dy を float にしないとエラーになる
            ax.add_patch(Arrow(x1, y1, float(x2-x1), float(y2-y1), linestyle='solid', color='grey', width=0.02))
        # 固有ベクトルを求める
        p11, p12 = sldrs[0].val, sldrs[1].val
        p21, p22 = sldrs[2].val, sldrs[3].val
        w, v = np.linalg.eig([[p11, p12],[p21, p22]])
        # 虚数が含まれない場合、固有ベクトルを表示する
        if np.all(np.isreal(v)):
            ax.add_patch(Arrow(0, 0, *(2*v[:,0]),  linestyle='solid', color='blue', width=0.02))
            ax.add_patch(Arrow(0, 0, *(2*v[:,1]),  linestyle='solid', color='blue', width=0.02))
    fig.canvas.draw()

def button_clicked(x):
    global patches_show_state
    patches_show_state = not patches_show_state
    update_patches()

# 描画範囲を固定する
ax.set_xlim(-2.0, 2.0)
ax.set_ylim(-2.0, 2.0)

# x軸のラベルを上に持ってくる
# https://stackoverflow.com/questions/14406214/moving-x-axis-to-the-top-of-a-plot-in-matplotlib
ax.xaxis.tick_top()

# アスペクト比を 1:1 にする
ax.set_aspect("equal")

# スライダを初期化
init_sldrs()

# グラフを初期化
ax_line = init_ax()

# スライダが変わった時に呼ばれる関数を登録
for sldr in sldrs:
    sldr.on_changed(lambda t: update_ax(ax_line))

# ボタンがクリックした時に呼ばれる関数を登録
button.on_clicked(button_clicked)

plt.show()

5. ソースコードの説明

calc_point(mat, point) 関数により、与えられた行列を点(ベクトル)の計算を行っています。行列とベクトルの計算の定義には sympy.MatrixSymbol を使用しています。calc_point 関数にすべての格子状の点に適応し、格子状の点を移動させています。スライダを動かすと update_ax(ax_line) 関数がよばれ、これらの計算処理と矢印の描画などが行われ、表示が更新されます。


固有ベクトルを計算するのに np.linalg.eigh を使用しました。
docs.scipy.org

矢印の表示には matplotlib.patches.Arrow を使用しました。矢印の再描画ではまったのですが、矢印を描画する前に矢印をすべて消去(remove)して再描画しています。矢印の消去には以下を参考にしました。

stackoverflow.com

6. 出力結果

出力されるグラフは以下になります。スライダで2x2行列の各要素の値を変更し変換行列の作用を見ることができます。

f:id:keimina:20190429125528g:plain

点や矢印の意味は以下の通りです。

No. 点や矢印 意味
1 赤の点 移動後の点
2 グレーの点 変換前の点
3 グレーの矢印 変換前と変換後の点を結んだ矢印
4 青の矢印 行列の固有ベクトル

7. 考察

行列の初期状態は単位行列ですので、この状態では変換前と変換後で点の位置はかわりません。行列は[[p11, p12], [p21, p22]] で p11, p12, p21, p22 のパラメータからなります。 p11 を動かすと x 方向に変形します。 p22 を動かすと y 方向に変形します。 p12 を動かすと平行四辺形のように x 平方向に変形します。 p21 を動かすと平行四辺形のように y 方向に変形します。p12 と p21 の符号が違うとき 点は回転するように変形します。そしてこのとき固有ベクトルの要素が複素数になります(そのため固有ベクトルの矢印は表示していません)。点の移動量をみると、あたかも固有ベクトルの方向を目指すように点が移動しているように見えます。このように固有ベクトルは、行列の作用を説明してくれているようにも見えます。

8. まとめ

2x2 行列による点の変換と、固有ベクトルについて可視化しました。スライダーを動かして変換行列と、それに対応する固有ベクトルの性質を見ることができました。

以上です。

勾配降下法により PCA(主成分分析)の第1主成分を数値解析的に求め、PCAについて理解を深める

2. はじめに

本記事は、主成分分析(以下PCA)について自分の理解があっているか確かめた時の記録です。PCAの第1主成分(First component)について、自分の理解があっているか検証を行いました。Sympy で独自実装した数値解析のアルゴリズムによって求めた第1主成分と sklearn.decomposition.PCA によって求めたものを比較し、誤差が少ないことを確認しました。

2019/07/04 更新:
続編は以下です。
(第2主成分、第3主成分についても検証しました)
keimina.hatenablog.jp

3. 参考文献

主成分分析の数学的な話については、以下を参考にしました。

www.hellocybernetics.tech

4. PCAの第一主成分について自分が理解していること

私の理解している、PCAの第一主成分の求め方を書きます。例えば3次元の特徴量をもつ複数のデータがあったとして、それらのデータは 点(x, y, z) のように3次元で表現できます。(説明を簡略化するため、これらの点の集合の重心は (0, 0, 0) とします。)点の集合の要素のi番目の特徴量を、x[i], y[i], z[i] としその線形結合を Y[i](次元削減後の特徴量) とし、線形結合の係数を w1, w2, w3 とすると

Y[i] = w1*x[i] + w2*y[i] + w3*z[i]

とかけます。
ここで、例えば、 データが2つあったとして、

data[0] = (1, 2, 3)
data[1] = (4, 5, 6)

とします。すると Y は、

Y[0] = w1*1 + w2*2 + w3*3
Y[1] = w1*4 + w2*5 + w3*6

となります。当たり前ですが、以下のような変換が行われることがわかります。

data[0] (3次元) → Y[0] (スカラ)
data[1] (3次元) → Y[1] (スカラ)

この変換後の Y を使用し以下の V を定義します。

V = Y[0]*Y[0] + Y[1]*Y[1]

つまり、 V は Y[i] の分散です。

そして V が最大になるような w1, w2, w3 が第1主成分です。

5. 検証の方針

上記の理解で本当に正しいのか検証するため、以下の方針に従い実装しました。

  1. V を定式化する
  2. 勾配降下法により V が最大になるような第1主成分(線形結合の係数)を求める
  3. sklearn.decomposition.PCA により第1主成分を求める
  4. 2. と 3. で求めた第1主成分を比較し差分が少ないことを確認する

6. ソースコード

方針に従い実装したソースコードは以下です。
※実装には主に Sympy を使用しています。(TensorFlowでもできますが、Sympy だとprint文で数式の表示ができデバッグが楽なため Sympy を使用しました。)

import numpy as np
from sympy import MatrixSymbol, Matrix
from sklearn.decomposition import PCA

# 1行目が観測データ1、2行目が観測データ2、3行目が観測データ3
# 1列目が特徴量1、2列目が特徴量2、3列目が特徴量3 とする
data = np.array([[3,5,7],[11,13,17],[19,23,29]], dtype=np.float64)

################################################################
# 独自手法で
# first component (第一主成分)を求めるための処理
################################################################

# 特徴量を定義する
# 1行目が観測データ1、2行目が観測データ2、3行目が観測データ3
# 1列目が特徴量1、2列目が特徴量2、3列目が特徴量3
Data = MatrixSymbol('Data', 3, 3) # 観測データ入力用の 3x3 の行列のシンボルを作成
Data_a = np.array(Data) # numpy 表現にする

# 1行3列の行列(縦ベクトル)を作成(first components の変数)
W = MatrixSymbol('W', 3, 1) # first components 入力用の 3x1 の行列のシンボルを作成
W_a = np.array(Matrix(W))  # numpy 表現にする

# Data_a を以下のように変換し、次元を1次元に削減する
Y = Data_a.dot(W_a) # shape=(3, 1)

# 変換後の特徴量の分散を求める
Ymean = np.mean(Y, axis=0, keepdims=True)
Yhat = (Y - Ymean)   # 重心を (0, 0, 0) にする、shape=(3, 1)
V = Yhat.T.dot(Yhat) # 分散を求める、 shape=(1, 1)

# 分散 V が最大になるようにする
# ここでは、数値解析的に勾配降下法を使用しVの最大値を求めるため微分する
# Matrix(V).diff(W[0, 0])
# Matrix(V).diff(W[1, 0])
# Matrix(V).diff(W[2, 0])
# 上記を求めるのに jacobian() というものが便利だったため使用しました。
dVdW = Matrix(V).jacobian(W).T # shape=(3, 1)
dVdW_a = np.array(dVdW)  # numpy 表現にする

# first components の値を保持する変数を用意しランダムな数値で初期化
np.random.seed(51) # 他の人がやっても同じ結果になるようにする
W_a_val = np.random.random((3,1))

# ハイパーパラメータ(学習係数)を定義
b = 0.00003

# 勾配降下法により first component の最大値を求める。
for i in range(101):
    # 勾配を求める
    dVdW_a_val = np.array(dVdW.subs({W:Matrix(W_a_val), Data:Matrix(data)})).astype(np.float64)
    # ハイパーパラメータ b をかけて W_a_val を更新する
    W_a_val += b*dVdW_a_val
    # 正規化する
    W_a_val /= np.linalg.norm(W_a_val) #np.sum(W_a_val.T.dot(W_a_val)) でも同じ
    if i%10==0:
        print('学習 i = ', i)
        print("分散=", Matrix(V).subs({W:Matrix(W_a_val), Data:Matrix(data)}))
        print("W=", Matrix(W_a_val))


################################################################
# sklearn.decomposition.PCA で
# first component (第一主成分)を求めるための処理
################################################################
# solving by sklearn.decomposition.PCA
pca = PCA(n_components=3) # 1-3 components を求める
pca.fit(data) # 学習する

################################################################
# 独自手法 と sklearn.decomposition.PCA で
# 求めた first component (第一主成分)を表示するための処理
################################################################
print('\n================ first component ================')
# 独自手法 で求めた first component を表示(比較しやすいようにflattendを使用)
print('--- 独自手法 ---')
print(W_a_val.flatten())

# sklearn.decomposition.PCA で求めた frist component を表示
w1_pca = pca.components_[0]
print('--- sklearn.decomposition.PCA ---')
print(w1_pca)

# 差分 を表示
print('--- 差分 ---')
print(w1_pca - W_a_val.flatten())

# sklearn.decomposition.PCA で求めた frist component と
# コサイン類似度を求める、w1_pca と W_a_val はノルムが1なので内積を
# とるだけで cosθ が求まる
print('--- コサイン類似度 ---')
cos_theta = w1_pca.dot(W_a_val)
print('cosθ = {} rad ({}°)'.format(*cos_theta, *np.rad2deg(np.arccos(cos_theta))))

ソースコードの解説についてはソースコードのコメントをみてください。最後に誤差が少ないか確認するため、コサイン類似度を求めています。コサイン類似度については以下を参考にしてください。

ja.wikipedia.org

7. 出力結果

学習 i =  0
分散= Matrix([[326.893552407606]])
W= Matrix([[0.885069441686711], [0.0711731365114243], [0.459985291103423]])
学習 i =  10
分散= Matrix([[398.915653114400]])
W= Matrix([[0.828195710513823], [0.186379419417903], [0.528540043046834]])
学習 i =  20
分散= Matrix([[452.012116765482]])
W= Matrix([[0.765356021528304], [0.282142685524890], [0.578468378837745]])
学習 i =  30
分散= Matrix([[486.534391735230]])
W= Matrix([[0.705988103517652], [0.356151889455055], [0.612157356673220]])
学習 i =  40
分散= Matrix([[507.182771621072]])
W= Matrix([[0.655167743479454], [0.410912766518736], [0.633960508404080]])
学習 i =  50
分散= Matrix([[518.921766869242]])
W= Matrix([[0.614225060024774], [0.450555959412570], [0.647863336727261]])
学習 i =  60
分散= Matrix([[525.403750244901]])
W= Matrix([[0.582431990459451], [0.479009309986766], [0.656751899453244]])
学習 i =  70
分散= Matrix([[528.925385131233]])
W= Matrix([[0.558289211467764], [0.499398761609065], [0.662505874133985]])
学習 i =  80
分散= Matrix([[530.821835476805]])
W= Matrix([[0.540205689228539], [0.514032980713747], [0.666294160310340]])
学習 i =  90
分散= Matrix([[531.838240751787]])
W= Matrix([[0.526776505606956], [0.524565081220116], [0.668833304123724]])
学習 i =  100
分散= Matrix([[532.381592136768]])
W= Matrix([[0.516858344804962], [0.532166525145605], [0.670564121408182]])

================ first component ================
--- 独自手法 ---
[0.51685834 0.53216653 0.67056412]
--- sklearn.decomposition.PCA ---
[0.48958803 0.55232128 0.67471829]
--- 差分 ---
[-0.02727032  0.02015476  0.00415417]
--- コサイン類似度 ---
cosθ = 0.9994164292655139 (1.9575168306292177°)

出力結果を見ると 「学習 i = 0」の時、分散が326だったのが、「学習 i = 100」で 532 に増えていることがわかります。最大化の方向に学習がうまく進んでいるようです。最後に「--- 独自手法 ---」に独自手法で求めた第1主成分と「--- sklearn.decomposition.PCA ---」にsklearn.decomposition.PCA で求めた第1主成分を表示しています。若干の誤差はあるもののほぼ等しいことがわかります。等しさの指標としてベクトルの差分とコサイン類似度も表示しました。コサイン類似度を見るとベクトルの方向がほとんど一致していることがわかります。

※学習の経過を見るため、学習係数bはあえて小さくしています。学習係数bを大きくするとすぐに sklearn.decomposition.PCAで求めた値と同じ値へ収束します。

8. まとめ・考察

PCAの第1主成分(First component)について、自分の理解があっているか検証しました。PCA ではデータの特徴量を線形変換し次元削減した後それらの分散が最大となるような線形変換の結合係数を求めているということがわかりました。ただ、PCAは特徴量に線形変換を施すことが前提としてありますので、特徴量に非線形の関係性がある場合は、たとえ分散を最大化しても、元の特徴をうまく表せるとは限らないということも言えそうです。PCA の理解が深まりました。後は第2主成分以降の成分なのですが、これについては、第1主成分と直交する&分散が最大となる結合係数を求めるということだと理解しています。余裕があれば、こちらも記事にしたいと思います。

以上です。おやすみなさい。

SciPy Japan 2019 に参加しました

前に申し込んでいた SciPy Japan 2019 に参加してきました。

予定通り4/23, 4/24 に開催されました。無事に開催されて良かったです。開催の内容は、午前中はチュートリアル、午後は Enthoughtさん による公演や Science に関する問題を(一部) Python を使って解決したというようなトークのセッションでした。私は個人で参加したのですが、企業の一員として来られている方が多い印象でした。あっさりと書いていますが、トークの一部は私が求めていたような内容そのもので(トークの選別した人グッジョブです)感動しました、発表者の方が非常にわかりやすく説明してくれましたし、ワクワク感を味わえました。どれくらいワクワクしたかというと、私が小学生だったらこの発表を聞いたらサイエンスの分野に足を踏み入れたくなるくらいです。

SciPy のコミュニティを推進している人たち(主にEnthoughtさん)の発表を聞いて思ったのは、 SciPy の根底にあるものは Science だということです。Science に主眼を置いて、 Python やそれに付随する各種ライブラリやフレームワークは、あくまで目的達成のための手段に過ぎないということです。それなら SciPy ではなく SciCode でいいじゃないかと思うのですが、 コードを書くには Python の方が都合が良かったため SciPy と名付けたということです(と私は解釈した)。実際その決断は間違ってなさそうで、例えば大学の講師がデータ処理にPythonを主要言語に選択したなどのような話は、その決断の正しさに一票を追加するものでしょう。今はPython言語 の勢いがありますが、SciPy コミュニティは、この勢いに負けずに、 Science を根底とするというコンセプトは変えないで欲しいと思います。

f:id:keimina:20190424184727j:plain

今の所、SciPy は 2020年も日本に来てくれるようです、個人的に、それまでに以下を行いたいと思っています。

1. 英語力(聞く、話す)を身に付ける
2. SciPy っぽいコミュニティを作る or 参加する

以上

sympy でニューラルネットワークの重みの更新に使用する式を計算グラフにしてみる

1. はじめに

sympy(https://www.sympy.org/en/index.html) を使用し、ニューラルネットワークの重みの更新に使用する式を計算グラフにしてみます。具体的には、 sympy を使って、損失関数の勾配の計算を計算グラフとして表し、それを graphviz で画像に出力します。sympy と graphviz を使用すると、 数式を簡単に計算グラフに変換できます。ちなみに計算グラフとは、以下のようなものです。(実際にこれらを使用し作成したものです)

f:id:keimina:20190325203008p:plain
図1. c = a + b の計算グラフ

このように計算式をグラフとして表現したものを計算グラフと呼びます。

2. 使用するモジュールと、そのインストール

本記事では、以下の二つのモジュールを使用します。グラフの描画に graphviz というモジュールを使用します。sympy は計算グラフをテキストで出力できますが、単体では画像で出力することができないためグラフ描画アプリ(graphviz)を使用します。

  1. sympy
  2. graphviz

graphviz は以下のコマンドでインストールしました。

pip install graphviz

コマンドの出力結果

Collecting graphviz
  Downloading https://files.pythonhosted.org/packages/1f/e2/ef2581b5b86625657afd32030f90cf2717456c1d2b711ba074bf007c0f1a/graphviz-0.10.1-py2.py3-none-any.whl
distributed 1.21.8 requires msgpack, which is not installed.
Installing collected packages: graphviz
Successfully installed graphviz-0.10.1
You are using pip version 10.0.1, however version 19.0.3 is available.
You should consider upgrading via the 'pip install --upgrade pip' command.

コマンド実行後、"requires msgpack" と出力されていますが、 最終的に"Successfully installed" と出力されているため、問題なしとみなし気にせず先に進みました。

今回は、計算グラフを画像で出力したいため、 さらに以下の手順行う必要がありました。

  1. MacPorts のインストール(インストール手順は以下のサイトを参照)
  2. 以下のコマンドの実行しアプリをインストール
sudo port install graphviz

sympy のインストール手順について説明は省略します。

3. ニューラルネットワークのモデルと、順伝播の式

本記事で扱うニューラルネットワークのモデルと、順伝播の式は以下の通りです。

f:id:keimina:20190325203628p:plain
図2. ニューラルネットワークのモデルと、順伝播の式

4. ソースコード

3. のニューラルネットワークの損失関数の勾配を計算するコードにすると以下のようになります。
※本記事は、重みの更新を行う計算を、計算グラフに出力することを目的としているため、説明をわかりやすくするため損失関数 l は o1 としています。本来であれば、l = (o1 - y)**2 のような定義をすべきですが、パラメータが増えると、計算グラフが大きくなり見づらくなるため、簡略化しています。ご了承ください。
ソースコードで numpy を使用して行列計算していますが、sympy での行列計算の仕方がわからなかったためです。あまり真似しないように。。。

from sympy import symbols, Function, tanh
from sympy.printing import dot
import graphviz
import numpy as np

# sympy で使用する変数(シンボル)を定義する
# 入力値、中間層の値
i1, i2 = symbols("i1, i2")
h11, h12 = symbols("h11, h12")

# 重み
w111, w112, w121, w122 = symbols("w111, w112, w121, w122")
w211, w212 = symbols("w211, w212")

# バイアス
b11, b12 = symbols("b11, b12")
b21 = symbols("b21")

# numpy を利用して行列形式にする
I = np.array([[i1], [i2]])
W1 = np.array([[w111, w112], [w121, w122]])
W2 = np.array([[w211, w212]])
B1 = np.array([[b11], [b12]])
B2 = np.array([[b21]])

# 活性化関数
TANH = np.vectorize(tanh)

# 順伝播の数式を作成する
H1 = TANH(W1.dot(I) + B1)
O = TANH(W2.dot(H1) + B2)

# 行列からスカラを取り出す
o1 = O[0, 0]

# 最急降下法 によって重み w1 を更新する式(学習係数は 1 とする)
# 損失関数を定義する、重み更新を行う計算グラフを出したいだけなので l = o1  としてしまう
l = o1
dl_dw111 = l.diff(w111)

expr_forward = o1
src_forward = graphviz.Source(dot.dotprint(expr_forward))

expr_backward = dl_dw111
src_backward = graphviz.Source(dot.dotprint(expr_backward))

# 私の環境ではExecutableNotFound エラーが発生したので、以下のコードで回避しました
app_path = '/opt/local/bin/dot'
graphviz.ENGINES.add(app_path)
src_forward._engine = app_path
src_backward._engine = app_path

# 画像ファイルを作成する
filename_forward = "calc_forward_graph"
src_forward.render(filename_forward)
filename_backward = "calc_backward_graph"
src_backward.render(filename_backward)

順伝播の式が o1 です。勾配降下法によって損失関数 l を最小にするため重みパラメータ w111 を更新するときに行われる計算が dl_dw111 です。コードの最後の方で、sympy と graphviz を使用して、o1 と dl_dw111 の計算を行う計算グラフを出力しています。sympy.printing.dot で計算グラフを graphviz で使用可能な dot フォーマット形式のテキストに変換しています。そして、 dot フォーマット形式のテキストを graphviz.Source で graphviz に取り込み、 render メソッドを呼びだし pdf ファイルを作成しています。

5. 出力結果

出力結果は以下になります。図は上から順に、順伝播の計算を行うための計算グラフ(o1 の計算を行うグラフ)、損失関数lを最小化すべく w111 を勾配降下法により更新するための計算グラフ(dl_dw111 の計算を行うグラフ)、となっています

f:id:keimina:20190325204132p:plain
図3. 順伝播の計算グラフ("calc_forward_graph")

f:id:keimina:20190325204331p:plain
図4. w111 の重みを更新する数式の計算グラフ("calc_backward_graph")

楕円の中に書かれている Mul 、 Add 、 Pow 、 tanh はどのような計算を行うかを表していて、それぞれ、 乗算、加算、累乗、ハイパボリックタンジェントを表しています。ニューラルネットの学習時には、図中の w111、w112 、w121、 b11 、、、 などの変数のように見えるものは全て値が割り振られているため、 この計算グラフの計算を行うと、何かしらの値(スカラ値)が求まります。そしてこの値を使用して重み w111 を更新します。

6. まとめ

sympy と graphviz を使用して、ニューラルネットワークの重みの更新に使用する式を計算グラフにしました。このように、ニューラルネットワークのモデルが決まると、重みを更新するための式、すなわち計算グラフも決まります。

7. 参考文献

  1. ゼロから作るDeep Learning――Pythonで学ぶディープラーニングの理論と実装

Numpy 入門 (機械学習勉強会 in 新潟 2019/03/23)

Numpy 入門

Numpy 入門
(機械学習勉強会 in 新潟 2019/03/23)

1 はじめに

本記事は、機械学習勉強会 in 新潟 2019/03/23 の発表資料です。 numpy の入門者向けの記事になります。 numpy とは何か、何ができるかなどを説明します。numpy で簡単な数値計算を行えるようになることを目的とします。

※時間の都合上、 numpy の機能のほんの一部を紹介しています。 本記事は以下の参考文献を参考に書いています。 numpy についての詳細は以下の参考文献を参照してください。

2 numpy とは何か

Python数値計算のためのライブラリです。 numpy によってできることとメリットは以下になります

No. ライブラリ 説明 できること、メリット
1 numpy 数値計算ライブラリ 1. 行列の計算を行うコードを短く書ける
      2. 行列の計算、ベクトルの計算の高速化( 速度比較 参照)

3 numpy モジュールのインポート方法について

numpy を使用する場合、numpy モジュールのインポートは以下のように行います。

コード1. 「モジュールのインポート」:

import numpy as np

これにより、'np' の2文字で numpy を使用することができます。

コード2. 「モジュールの使用方法」:

import numpy as np
arr = np.array([1,2,3])
arr
Out: array([1, 2, 3])

使用する名前は numpy など np 以外の名前にもできますが、 np と書くのが一般的になっているためこちらを使用しましょう。

4 数値を作成する

はじめに、numpy で数値を作成しましょう。数値は、スカラ(Scalar)と呼ばれることがあります。 これらは、後述する配列の要素( 5.1.1 参照)となりますので、しっかり覚えましょう。 np.int32 を使用すると numpy の整数を作成できます。

コード3. 「numpy で数値の作成」:

import numpy as np
answer = np.int32(100)
answer
Out: 100

4.1 np.int32 とは

np.int32 はデータタイプ(dtype)として定義されています。 np.int32(100) を実行すると np.int32 データタイプのオブジェクトが作成されます。 np.int32 以外にも数値を表わすデータタイプが用意されています。

データタイプ 数値の種類 用途
np.uint8 8bit 符号無し整数 グレースケール画像
np.int32 32bit 整数 画像処理
np.float32 32bit 浮動少数 画像処理、平均値、確率
np.float64 64bit 浮動少数 数値計算(精度が求められるもの)
   

真偽値(True/False)を現す場合は np.bool データタイプが用意されています。

5 「1次元配列」の数値計算をする

配列の基本となる1次元配について学習しましょう。

5.1 「1次元配列」を作成する

1次元の配列を作成しましょう。1次元配列とは以下のようなものです。配列は'要素'と呼ばれるもので構成されています。

値1 値2 値3 ,,,

5.1.1 配列の'要素'とは

一般的に、上記 5.1 のような配列の中にある値1や値2などのことを配列の”要素”(Element)と呼びます。 1次元配列だけでなく、2次元配列や3次元配列も同様です。配列は要素によって構成されています。

※ numpy では配列を作成する関数の引数に dtype を指定することで、その配列の要素のデータタイプを指定できます。 引数 dtype を指定しない場合は、自動的に配列のデータタイプが決定されます。 numpy では配列は同一のデータタイプの要素から構成されます。

5.1.2 np.array 関数

np.array 関数を使用すると、標準 Python のリストを配列に変換することができます.

コード4. 「標準 Python のリストを配列に変換する」:

import numpy as np
arr = np.array([1, 2, 3])
arr
Out: array([1, 2, 3])

5.1.3 np.arange 関数

np.arange 関数を使用すると、数値の1次元配列を作成することができます。

コード5. 「np.arange 関数を使用し、数値の1次元配列を作成」:

import numpy as np
arr = np.arange(3)
arr
Out: array([0, 1, 2])

5.2 「1次元配列」の形状 shape について

配列には shape と呼ばれる、配列の形状を表す情報があります。 配列の形状 shape を調べるには以下のようにします。 この例の場合、配列は長さ3の1次元配列であることがわかります。

コード6. 「配列の形状 shape を調べる」:

import numpy as np
arr = np.arange(3)
answer = arr.shape
answer
Out: (3,)

5.3 「1次元配列」 + 「1次元配列」

1次元配列と1次元配列の四則演算の基本を学習します。 ここでは例として、1次元配列と1次元配列の足し算を行います。 以下のように、1次元配列を二つ作成し、足し算記号(+) を使用して求めます。 計算は、要素ごと(element wise)に行われます。

コード7. 「1次元配列と1次元配列の足し算」:

import numpy as np
arr_1 = np.array([1, 2, 3])
arr_2 = np.array([0, 0, 1])
answer = arr_1 + arr_2
answer
Out: array([1, 2, 4])

その他の演算子(-*/など)についても、同じように、配列の要素ごとに計算されます。 形状の異なる配列同士の計算は、エラーになるか、ならない場合があります。 形状の異なる配列同士の計算については 11 を参照してください。

5.4 配列の要素に数学関数を適用する

1次元配列の要素ごと(element wise)に数学関数を適用してみましょう。 ここでは例として平方根を適用してみましょう。

コード8. 「配列の各要素に平方根を適用」:

import numpy as np
arr = np.array([1, 2, 3])
answer = np.sqrt(arr)
answer
Out: array([1.        , 1.41421356, 1.73205081])

np.sqrt 関数以外にも様々な数学関数が用意されています。 それらの関数も、上記のように使用することで、要素ごとに適用することができます。 以下に数学関数の一部を記載します。

No. 関数 説明
1 np.exp 指数関数
2 np.log 対数関数
3 np.sin sin関数
4 np.arcsin arcsin関数

5.5 1次元配列の要素の合計を求める

合計を求める場合は、 np.sum 関数を使用します。

コード9. 「np.sum 関数を使用して1次元配列の要素の合計を求める」:

import numpy as np
arr = np.arange(11)
answer = np.sum(arr)
answer
Out: 55

6 「2次元配列」の数値計算をする

2次元配列について学習しましょう。

6.1 「2次元配列」を作成する

2次元の配列を作成しましょう。2次元配列とは以下のようなものです。

値11 値12 値13 ,,,  
値21 値22 値23 ,,,  
値31 値32 値33 ,,,  
,,,        

6.1.1 np.array 関数

np.array 関数を使用すると、2次元配列を作成することができます。(1次元配列の時と同様)

コード10. 「np.array 関数を使用して2次元配列を作成する」:

import numpy as np
arr = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
arr
Out: 
array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9]])

6.1.2 reshape メソッド

配列の reshape メソッドを使用すると配列の形状を変えることができます。 ここでは、 1次元配列の形状を変えて、2次元配列を作成してみましょう。

コード11. 「reshape メソッドを使用して1次元配列の形状を変えて、2次元配列を作成する」:

import numpy as np
arr = np.arange(9)
answer = arr.reshape((3, 3))
answer
Out: 
array([[0, 1, 2],
       [3, 4, 5],
       [6, 7, 8]])

6.2 「2次元配列」 + 「2次元配列」

2次元配列同士の四則演算、関数を適用した時の動作は、1次元配列の時と同様、要素ごとに計算が行われます。 (そのため2次元配列の四則演算については説明を省略します。)2次元配列を使用すると、行列の掛け算を行うことができます。

行列の掛け算とは、具体的には以下のような計算のことです。 行列 A と行列 B の掛け算(積)は A・B と書き、A と B と A・B は以下のような関係になります。 ※説明の簡略化のため行列の大きさは2x2としています。

  • 行列 A

    a b
    c d
  • 行列 B

    e f
    g h
  • 行列 A・B

    a*e+b*g a*f+b*h
    c*e+d*g c*f+d*h

行列の掛け算を行うには np.dot 関数を使用します。

コード12. 「np.dot 関数を使用して行列の掛け算を行う」:

import numpy as np
arr_1 = np.arange(9).reshape((3, 3))
arr_2 = np.arange(9).reshape((3, 3))
answer = np.dot(arr_1, arr_2)
answer
Out: 
array([[ 15,  18,  21],
       [ 42,  54,  66],
       [ 69,  90, 111]])

配列オブジェクトの dot メソッドを使用して以下のように書くこともできます。

コード13. 「dot メソッドを使用して行列の掛け算を行う」:

import numpy as np
arr_1 = np.arange(9).reshape((3, 3))
arr_2 = np.arange(9).reshape((3, 3))
answer = arr_1.dot(arr_2)
answer
Out: 
array([[ 15,  18,  21],
       [ 42,  54,  66],
       [ 69,  90, 111]])

注意: np.dot は '∗'による掛け算と同じではありません。 '∗' による掛け算は要素ごと(element wise)に掛け算が行わます。 一方 np.dot は行列の掛け算ですので、積和演算(掛け算と足し算)が行われます。

6.3 「2次元配列」の要素の合計を求める

'すべての要素の合計'を求める場合は、1次元配列の時と変わりませんので省略しますが、 2次元配列では、さらに、'列の合計'と'行の合計'を求めることができます。

6.3.1 列の合計

np.sum 関数の引数で axis=0 とすることで列の合計を求めることができます。

コード14. 「列の合計を求める」:

import numpy as np
arr = np.arange(9).reshape((3, 3))
answer = np.sum(arr, axis=0)
answer
Out: array([ 9, 12, 15])

6.3.2 行の合計

np.sum 関数の引数で axis=1 とすることで行の合計を求めることができます。

コード15. 「行の合計を求める」:

import numpy as np
arr = np.arange(9).reshape((3, 3))
answer = np.sum(arr, axis=1)
answer
Out: array([ 3, 12, 21])

6.3.3 配列の次元数を保つには

6.3.16.3.2 で np.sum を使用して求めた配列を見ると、2次元配列の入力に対して出力が1次元配列となっており、次元数が変わってしまうことに気付くかと思います。 np.sum 関数の引数で keepdims=True とすることで、配列の次元数を保つことができます。

コード16. 「np.sum 関数の引数で keepdims=True とし配列の構造をなるべく保つ」:

import numpy as np
arr = np.arange(9).reshape((3, 3))
answer = np.sum(arr, axis=1, keepdims=True)
answer
Out: 
array([[ 3],
       [12],
       [21]])

6.4 行列計算用の関数

2次元配列には行列計算用の関数を使用することができます。 行列計算用の関数は np.linalg にあります。 その一部を紹介します。

関数 説明
np.linalg.inv 逆行列を求めるための関数
np.linalg.norm 行列のノルムを求めるための関数

7 配列を作成する関数

配列を作成する方法は、5.16.1 で紹介した方法が全てではありません。 以下に、配列を作成するための関数の一部を紹介します。

7.1 固定された値を要素にもつ配列を作成する関数

No. 関数 説明
1 np.zeros 全ての要素が 0 の配列を作成します
2 np.ones 全ての要素が 1 の配列を作成します
3 np.full 全ての要素が x の配列を作成します(x は引数で指定)

コード17. 「np.zeros 関数を使用し全ての要素が 0 の配列を作成」:

import numpy as np
arr = np.zeros((3, 3))
arr
Out: 
array([[0., 0., 0.],
       [0., 0., 0.],
       [0., 0., 0.]])

7.2 単位行列を作成する関数

No. 関数 説明
1 np.eye 単位行列を表す配列を作成します

コード18. 「np.eye 関数を使用し単位行列を表す配列を作成」:

import numpy as np
arr = np.eye(3)
arr
Out:
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

7.3 ランダムな値を要素にもつ配列を作成する関数

No. 関数 説明
1 np.random.randint ランダムな整数の配列を作成します
2 np.random.random ランダムな少数の配列を作成します

コード19. 「np.random.randint 関数を使用しランダムな整数の配列を作成」:

import numpy as np
arr = np.random.randint(0, 10, (3, 3))
arr
Out:
array([[4, 0, 7],
       [2, 3, 4],
       [3, 0, 5]])

コード20. 「np.random.random 関数を使用しランダムな少数の配列を作成」:

import numpy as np
arr = np.random.random((3, 3))
arr
Out:
array([[0.23963676, 0.14041946, 0.67624721],
       [0.75209081, 0.3182057 , 0.18579186],
       [0.0888395 , 0.87384273, 0.24740674]])

7.4 似たような配列を作成する関数

配列の計算をプログラミングしているときに、ある配列と同じ形状をもつ配列だが、値は違うという配列を作成したい時があります。そのような場合は、以下の関数を使用します。

No. 関数 説明
1 np.zeros_like 引数で指定した配列と同じ形状
   
    全ての要素が 0 の配列を作成します
     
2 np.ones_like 引数で指定した配列と同じ形状
   
    全ての要素が 1 の配列を作成します

コード21. 「np.zeros_like 関数を使用し形状が同じ配列を作成」:

import numpy as np
arr = np.random.randint(0, 10, (3, 3))
arr_2 = np.zeros_like(arr)
arr_2
Out[4]: 
array([[0, 0, 0],
       [0, 0, 0],
       [0, 0, 0]])

7.5 ある範囲の数値を要素にもつ配列を作成する関数

No. 関数 説明
1 np.linspace 数直線の範囲を等幅で分割した1次元配列を作成します
    1変数関数を調べる時などに使用できます
2 np.meshgrid 2次元座標などの座標を配列を使い表現するときに使用します。
    2変数関数を調べる時などに使用できます

コード22. 「np.linspace 関数の使用例」:

import numpy as np
arr = np.linspace(0.0, 0.9, 10)
arr
Out: array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9])

コード23. 「np.meshgrid 関数の使用例」:

import numpy as np
x = np.linspace(0, 0.4, 5)
y = np.linspace(0, 0.4, 5)
x_mesh, y_mesh = np.meshgrid(x, y)
answer = x_mesh, y_mesh
answer
Out:
(array([[0. , 0.1, 0.2, 0.3, 0.4],
        [0. , 0.1, 0.2, 0.3, 0.4],
        [0. , 0.1, 0.2, 0.3, 0.4],
        [0. , 0.1, 0.2, 0.3, 0.4],
        [0. , 0.1, 0.2, 0.3, 0.4]]), array([[0. , 0. , 0. , 0. , 0. ],
        [0.1, 0.1, 0.1, 0.1, 0.1],
        [0.2, 0.2, 0.2, 0.2, 0.2],
        [0.3, 0.3, 0.3, 0.3, 0.3],
        [0.4, 0.4, 0.4, 0.4, 0.4]]))

8 配列の要素へのアクセス

8.1 配列の要素の取得

添字(インデックス)を使用することで、配列の要素を取得することができます。

8.1.1 要素の取得の基本

以下のようにすると、配列の要素の値を取得することができます。

コード24. 「1次元配列の配列の要素の取得の基本」:

import numpy as np
arr = np.array([1, 3, 5])
arr[2]
Out: 5

コード25. 「2次元配列の配列の要素の取得の基本(1)」:

import numpy as np
arr = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]])
arr[2][2]
Out: 23

上記コードは以下のようにも書くことができます。

コード26. 「2次元配列の配列の要素の取得の基本(2)」:

import numpy as np
arr = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]])
arr[2, 2]
Out: 23

8.1.2 要素を範囲で指定して取得

以下のようにすると、指定した範囲にある配列の要素を取得することができます。

コード27. 「1次元配列の要素を範囲で指定して取得」:

import numpy as np
arr = np.array([1, 3, 5])
arr[1:]
Out: array([3, 5])

コード28. 「2次元配列の要素を範囲で指定して取得」:

import numpy as np
arr = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]])
arr[1:, 1:]
Out: 
array([[11, 13],
       [21, 23]])

8.1.3 逆順の配列の取得

以下のようにすると、逆順の配列を取得できます。

コード29. 「逆順の配列の取得」:

import numpy as np
arr = np.arange(10)
arr[::-1]
Out:
array([9, 8, 7, 6, 5, 4, 3, 2, 1, 0])

8.1.28.1.3 で取得した arr[1:, 1:] や arr[::-1] は、あたかも配列のコピーのように思われるかもしれませんが、実際は配列のコピーではありません。これらはビューと呼ばれるものであり、実態は arr のデータを参照しているものです。扱いには注意しましょう。 8.2 でも説明します。

8.2 配列の要素の書き換え

8.2.1 要素の書き換えの基本

配列の要素を添字(インデックス)で指定して別の値に書き換えることができます。 以下のようにすると、配列の要素の値を書き換えることができます。

コード30. 「1次元配列の配列の要素の書き換えの基本」:

import numpy as np
arr = np.array([1, 3, 5])
arr[2] = 100
arr
Out: array([  1,   3, 100])

コード31. 「2次元配列の配列の要素の書き換えの基本(1)」:

import numpy as np
arr = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]])
arr[2][2] = 100
arr
Out:
array([[  1,   3,   5],
    [  7,  11,  13],
    [ 17,  21, 100]])

上記コードは以下のようにも書くことができます。

コード32. 「2次元配列の配列の要素の書き換えの基本(2)」:

import numpy as np
arr = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]])
arr[2, 2] = 100
arr
Out:
array([[  1,   3,   5],
       [  7,  11,  13],
       [ 17,  21, 100]])

8.2.2 要素を範囲で指定して書き換え

以下のようにすると、指定した範囲にある配列の要素を書き換えることができます。

コード33. 「1次元配列の要素を範囲で指定して書き換え」:

import numpy as np
arr = np.array([1, 3, 5])
arr[1:] = 100
arr
Out: array([  1, 100, 100])

コード34. 「2次元配列の要素を範囲で指定して書き換え」:

import numpy as np
arr = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]])
arr[1:, 1:] = 0
arr
Out: 
array([[ 1,  3,  5],
       [ 7,  0,  0],
       [17,  0,  0]])

上記コードの arr[1:] や arr[1:, 1:] は配列のコピーではないことに注意しましょう。 このようなものを numpy ではビュー(View)と呼びます。実態は arr のデータを参照しているものです。 このような ビュー arr[1:] や arr[1:, 1:] に対する値の代入は、ビューの参照先である arr に対する値の代入であることを意味します。

よくあるのが、値を保持しておきたい配列のビューに代入を行なってしまいその配列を破壊してしまうことです。 値を代入したい場合は、以下のように、配列のコピーを作成して、そのコピーに対して代入しましょう。

コード35. 「配列のコピーを作成して、そのコピーに対して代入する」:

import numpy as np
arr_orig = np.array([[1, 3, 5], [7, 11, 13], [17, 21, 23]])
arr = arr_orig.copy()
arr[1:, 1:] = 0
arr_orig
Out: 
array([[ 1,  3,  5],
       [ 7, 11, 13],
       [17, 21, 23]])

9 真偽値(True/False)について

'==' や '<' などの比較演算子を使用すると配列の要素ごとに比較が行われます。

コード36. 「比較演算子による配列の比較」:

import numpy as np
arr_1 = np.arange(9).reshape((3, 3))**2
answer = arr_1 >= 10
answer
Out: 
array([[False, False, False],
       [False,  True,  True],
       [ True,  True,  True]])

この真偽値を配列のインデックスに指定すると、 真偽値のTrue の箇所の要素だけ集めた1次元配列を求めることができます。

コード37. 「真偽値の配列による要素の収集」:

import numpy as np
arr_1 = np.arange(9).reshape((3, 3))**2
answer = arr_1[arr_1 >= 10]
answer
Out: array([16, 25, 36, 49, 64])

以下のように、 np.where を使用すると値が True のインデックスを求めることができます。 np.where の戻り値は、「インデックスを要素にもつ配列」を要素にもつタプル(tuple)です。

コード38. 「np.where による条件を満たす配列のインデックスの取得」:

import numpy as np
arr_1 = np.arange(9).reshape((3, 3))**2
answer = np.where(arr_1 >= 10)
answer
Out: (array([1, 1, 2, 2, 2]), array([1, 2, 0, 1, 2]))

10 配列の制約

ここまで1次元配列、2次元配列の計算方法について説明しました。 numpy の配列は、標準Python のリストと違い、一度作成した配列に要素を追加したり、削除したりすることはできません。 一度作成した配列のデータタイプ(dtype)を変更することもできません。これらの操作を行いたい場合は、新たに配列を作成してください。 これらの制約を設けることにより高速(C言語並)な計算が可能になっています。

11 配列の長さや、形状が異なる場合に行われる計算について

形状の異なる配列同士の演算を行う場合、それ専用の規則(broadcasting rule)に従い処理が行われます。 その規則に従いエラーになることもあれば計算結果が求まることもあります。

11.1 エラーになるか計算結果が求まるかの判定方法

配列 arr1 と arr2 のそれぞれの配列の shape を一番右から左にかけて以下条件を満たす時、計算結果が求まります。 条件を満たさない場合はエラーとなります。また、エラーにならず、形状が異なる場合は、お互いに異なる形状を補うよう値が補完され、計算されます(詳細は割愛します)。

(以下は「Broadcasting」(https://docs.scipy.org/doc/numpy/user/basics.broadcasting.html) より引用)

  1. 両方の値が同じ、または、
  2. 片方の数値が1

11.2 Broadcasting rule の適用例

  1. OK

    配列 shape
    arr1 10
    arr2 1
    結果 10
  2. OK

    配列 shape  
    arr1 10 1
    arr2 1 10
    結果 10 10
  3. OK

    配列 shape    
    arr1 10 10 1
    arr2 1 1 3
    結果 10 10 3
  4. NG

    配列 shape    
    arr1 10 10 3
    arr2 1 1 2
    結果 エラー    

12 速度比較

numpy を使用すると、どれくらい速くなるのでしょうか。 100x100の要素からなる2次元配列の各要素を +1 する処理を numpy を使用する前と後で処理時間を比較してみました。 処理時間の計測には、マジックコマンド %timeit を使用しています。 その結果、numpy の使用により処理速度は 100倍 になりました。また可読性も向上したと思います。

コード39. 「速度比較 numpy 使用前」:

def func(x):
   y = x.copy()
   for i in range(100):
       for j in range(100):
           y[i][j] += 1
   return y

# リストの初期化
x = [[0 for i in range(100)] for j in range(100)]
# 処理時間計測
%timeit func(x)
785 µs ± 24.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

コード40. 「速度比較 numpy 使用後」:

import numpy as np

def func(x):
  y = x.copy()
  y += 1
  return y

# 配列の初期化
x = np.zeros((100, 100))
# 処理時間計測
%timeit func(x)
7.11 µs ± 284 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

13 まとめ

配列の作成方法、配列同士の計算や、配列への関数の適用方法について学習しました。

numpy を使って得られる恩恵は以下の2つです。

  • 行列計算のコード量が少なくなる
  • 処理速度が早くなる

14 最後に

既存の行列の計算のソースコードの可読性の向上や、処理速度向上などを目的に、 numpy を使用したいと思われるかもしれません。 まず計算の対象が 10 の制約の範囲内で処理が可能であることを確認しましょう。 次に、以下の表のように可読性と、処理速度を天秤にかけ、 numpy で実装するか標準 Python で実装するか考えると良いのではないかと思います。 (numpy を使ったからといって常に「可読性が高い&処理速度が早い」コードになる訳ではありません。)

可読性 処理速度 numpy を使用する?
高くなる 高くなる YES
高くなる 低くなる YES or NO
低くなる 高くなる YES or NO
低くなる 低くなる NO

numpy を使用してコーディングするとテクニカルで難解なコードになってしまうことがあります。 可読性重視であればテクニカルなコーディングはしないほうが良いでしょう。 numpy を使用してコーディングする際は、numpy を使用する当初の目的を見失わないようにしましょう。

それでは、よい numpy ライフを!

著者: Kei Minagawa (http://keimina.hatenablog.jp)

Created: 2019-03-23 Sat 11:25

Validate