Kei Minagawa's Blog

皆川圭(@keimina)のブログ、Pythonで試したことを書いていきます

(エクセルで作成した csv ファイルについて)一つのセルに書かれている複数の行を、テキストエディタ Emacs で一つの行で見れるようにする方法

1. はじめに

エクセルで作成した csv ファイルで一つのセルの中に複数の行の記載があるファイルは、よく見かけると思います。今回の記事では、この csv ファイルをテキストエディタ Emacs で開き、一つのセルに書かれている複数の行を一つの行で見れるように変換する方法を書きます。さらに、変換したあと逆変換し元に戻す方法も書きます。使用するエディタは Emacs です。

2. 注意事項、制限

以下の csv では正しく動作しません。
- CRLF(キャリッジリターン + ラインフィード)以外で、 CR(キャリッジリターン)を使用している csv ファイル

※変換、逆変換するときに CR を区切り文字として使用するため、 CR を単独で使用している csv ファイルの場合、正しく動作しません。

3. 入力ファイル

3.1. エクセル

以下の様なデータをエクセルで作成します。2行目のA列のセルと2行目のB列のセルが複数行になっています。これを csv で保存したものが入力ファイルです。

f:id:keimina:20190731224004p:plain

3.2. 入力ファイルの初期状態

入力ファイルを Emacs で開くと以下の様に表示されます。^M と表示されているのは CR(キャリッジリターン)を表す文字です。LF(ラインフィード)は改行として表示されています。

f:id:keimina:20190731224102p:plain

4. 方法

複数行から1行への変換、1行から複数行へ変換する具体的な方法は以下の通りです。Emacscsv ファイルを開き以下のような手順を実行すると複数行、CR(キャリッジリターン)区切りの1行にできます。理論については「5. 背景、理論」に記載しました。

4.1. 1行への変換

複数行から1行へ変換するには以下を行います。

4.1.1. 検索する正規表現

検索する正規表現を入力します。

f:id:keimina:20190731224137p:plain

4.1.2. 置換のための emacs lisp

置換するための文字列として emacs lisp 式を入力します。

f:id:keimina:20190731224208p:plain

4.1.3. 置換実行後

置換実行後の状態です。複数行になっていた箇所が、CR区切りでスッキリと1行にまとまりました。

f:id:keimina:20190731224234p:plain

4.2. 複数行への変換

1行から複数行へ変換(逆変換)するには以下を行います。

4.2.1. 検索する正規表現

検索する正規表現を入力します。

f:id:keimina:20190731224349p:plain

4.2.2. 置換のための emacs lisp

置換するための文字列として emacs lisp 式を入力します。

f:id:keimina:20190731224419p:plain

4.2.3. 置換実行後

置換実行後の状態です。入力ファイルの初期状態に戻りました。

f:id:keimina:20190731224445p:plain

5. 背景、理論

エクセルで作成した csv ファイルは、セルの中の改行は LF で、そうでない改行(個々のサンプルを隔てるための改行)は CRLF で表されるようです。今回は LF で記載された複数の行を一つにすることが目的なので、Emacs に表示されている CRLF を普通の改行に、普通の改行を CR に変換すればよさそうだと考えました。すなわち CRLF を LF に LF を CR に同時に変換するということです。そして、そうやってできたデータを元に戻す、逆の変換には、LF を CRLF に CR を LF に変換すればよさそうです。このような変換は、一般的にはプログラミング言語を使うと正規表現のグルーピングを用いた置換で実現できると思います。テキストエディタ Emacs では、正規表現で指定した文字列を置換するための文字列として emacs lisp 式を使用することができます(「バックスラッシュ+カンマ+式」の形式で記載します)ので、こちらを使用して置換しました。

6. まとめ

テキストエディタ Emacs を使用して、csv で一つのセルに書かれている複数の行を CR 区切りで一つの行として見れるように変換しました。さらに、それを逆変換し元に戻しました。

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

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で学ぶディープラーニングの理論と実装