Kei Minagawa's Blog

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

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

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

SciPy Japan 2019 に参加します

SciPy Japan 2019 に参加します。

www.scipyjapan2019.scipy.org

開催は4/23, 4/24の二日間、場所は東京。チュートリアルトーク、ポスター発表が行われるようです。チュートリアルも二日間に渡って行われるようです。

英語のリスニングができずに、ほとんど理解できないかもしれず、無謀なのですが、せっかくの機会ですので、参加してみます。果たしてどうなることか、、、

Wikipediaの文書を使って文の長さの平均を求める

Wikipediaの文書を使って文の長さの平均値と最頻値を求めました。機械学習勉強会 in 新潟で発表した内容です。Wikipedia の約60億文字以上のXMLファイルから、独自に定義した文を抜き出し、その文に対して統計量を求めました。使用言語はいつも通りPythonです。

1. 動作環境について

コードは以下の環境で動作確認を行っています。

  • CPU:6コア
  • メモリ:256GB
  • HDD:50GB

※環境はGCE(Google Compute Engine)で行っています。実行する場合はメモリ量に注意してください。

2. Wikipedia 文書のアーカイブについて

Wikipedia の文書は以下からダウンロードできます。 "jawiki-latest-pages-articles.xml.bz2"の中にあるXMLが今回解析に使用したファイルになります。

dumps.wikimedia.org

3. 発表スライド

発表スライドは、一番下に置いておきましたので、そちらをご覧ください。

4. やってみた感想

まず、このような分析は、自然言語処理の分野では頻度分析と呼ばれるらしいです。統計ですね。最近は統計のビジネス書をコツコツ読んで勉強していたたため、おかげでサンプルデータの分布を見た瞬間パラメタ推定できると気づくことができました。
パラメタ推定は本で読んでいたけど実際やってみると、ニューラルネットワークなど必要とせず非常にシンプルなのですが、強力でした。今回は、分布が対数正規分布にフィットしたのですがその理由を考察するときに、仮定した関数の性質が役に立ちました。結果、個人的に納得感のある考察を得ることができました(対数正規分布に従うのはなぜか?の考察です)。
最後ですが、Wikipediaを日常的に使用していたため、その中に法則性が存在するという事実に驚きました。おそらく解析しないと一生気付かなかったと思います。個人的に、その法則性の存在を客観的に示せたので非常に有意義だったと思います。また、求めた最頻値の値自体にも、価値があると思います。この値はWikipediaで遭遇する可能性が一番高い文の文字列の長さなので、それは、もしかすると我々が最も親しんでいる文字列長かもしれません。そのため24文字を目指して文を書くとわかりやすさの改善に繋がるつかもしれません。ビッグデータの時代ですが、手付かずのデータが山ほどあるとのことなので、まだ誰も知らない事実がデータの中に大量に眠っていると推測されます。今後それらのデータの中から、想像もしなかった事実が出てくることを期待しています。

5. 発表会で上がった質問や意見

発表会で上がった質問をここにメモしておきます。これらの質問や意見については、今後扱うかもしれませんし、扱わないかもしれません。個人的には、最頻値を求めて十分満足してしまったため、あとやったら面白そうなことは、コーパス教師なし学習クラスタリングなど)をさせてテキストマイニングするくらいしか思いつかないです。あと人工無脳

1. 他の分布で推定してみたらどうなるのか?例えばポアソン分布でやってみたらどうなるか?
→ ググったら一番最初に出てきたキーワードが対数正規分布だったのでそれでやってみた、ポアソン分布ではパラメタ推定していない。(ただし、対数正規分布でパラメタ推定を行ったきっかけはどうであれ、実際は、考察に記載したように対数正規分布を仮定するもっともらしい理由があるため、ポアソン分布や他の分布でパラメタ推定は行う必要性を感じなかった)
2. 集めた文(コーパスとも言える)を今後どういったことに使えるのか?
→ オートエンコーダによる教師なし学習(クラスタリング)を行いテキストマイニングを行うことが考えられる。
3. Wikipediaではなくニュースサイトなどの文の場合は違う分布になるのか?
→ 本などをみるとコーパスの分野が変わると違う分布になったり、あるコーパスで学習した機械学習のモデルが、別の分野のコーパスを入力して予測すると精度が落ちるなどの事例がある様子。そのため、ニュースサイトの分布とは異なるものになる可能性が高い。
4. XMLの読み込みにかかる時間はどれくらいか
→ 約15分くらい
5. 正規表現の適用にかかる時間はどれくらいか
→ 6コアCPUで分散処理しているため早い、約5分くらい。CPUを増やせばほぼリニアに早くなる(←発表では10分くらいと言ったかもしれませんが訂正します)

発表スライド

発表したスライドになります。正規表現が冗長になっており(スクエアブラケットの中のバックスラッシュのこと)、もっとコンパクトにできますが、そこは認識しています。時間がないため修正せずにあげています。ご指摘ありましたらコメントをお願いいたします。
以上。

f:id:keimina:20190225193006j:plain
スライド1
f:id:keimina:20190225193114j:plain
スライド2
f:id:keimina:20190225193134j:plain
スライド3
f:id:keimina:20190225193152j:plain
スライド4
f:id:keimina:20190225193210j:plain
スライド5
f:id:keimina:20190225193306j:plain
スライド6
f:id:keimina:20190225193403j:plain
スライド7
f:id:keimina:20190225193417j:plain
スライド8
f:id:keimina:20190225193434j:plain
スライド9
f:id:keimina:20190225193448j:plain
スライド10
f:id:keimina:20190225193501j:plain
スライド11
f:id:keimina:20190225193520j:plain
スライド12
f:id:keimina:20190225193609j:plain
スライド13