元バイオ系

元バイオウェット系がデータサイエンスやらを勉強していくブログ。 基本自分用のまとめ。

EMアルゴリズムを理解したい【最尤法から】

よく忘れるので一旦ちゃんとEMアルゴリズム*1についてまとめておきます。
最尤法について振り返ってからEMアルゴリズムの説明に入っていきます。
あまり計算に自信はありませんので、間違いを見つけたら教えていただけると幸いです。

最尤法

観測事象の確率密度関数から尤度を構成し、その確率密度関数のパラメータの推測を尤度の最大化によって行う方法です。
さっそく説明をしていきたいところですがその前に、尤度について振り返りましょう。

尤度についてはこちらをどうぞ

hotoke-x.hatenablog.com

それでは最尤法についてですが、実は一度正規分布を例に一度記事を書いています。

hotoke-x.hatenablog.com

ここでは特定の分布を仮定せず、一般的な最尤法について書いていきます。


 \thetaを母数、 x を確率変数 \displaystyle Xの実現値、 p({\boldsymbol x}|{\boldsymbol \theta})を母集団が従う確率密度関数とします。
この時、 n個の観測値 x_1, x_2, \ldots , x_nの同時確率関数は  \displaystyle {\boldsymbol x} = \left(x_1, \ldots, x_n \right)^Tとして、

$$ \begin{align} p({\boldsymbol x} | {\boldsymbol \theta}) = \prod_{i=1}^{n} p(x_i|{\boldsymbol \theta}) \end{align} $$

で与えられます。これを、同時確率分布ではなく、母数 {\boldsymbol \theta}の関数として考えると

$$ \begin{align} L_n ({\boldsymbol \theta} |{\boldsymbol x}) = \prod_{i=1}^{n} p(x_i | {\boldsymbol \theta}) \end{align} $$

と書け、これを尤度関数と呼びます。

最尤法は、式(2)の尤度を最大化する \displaystyle {\boldsymbol \theta}を推定する方法です。この時の \displaystyle {\boldsymbol \theta}最尤推定値と呼び、 \displaystyle \hat {\boldsymbol \theta}と書くのが慣習的となっています。また、 \displaystyle {\boldsymbol x}の関数、 \displaystyle \hat {\boldsymbol \theta} ({\boldsymbol x})最尤推定量と呼びます。忘れがちなのでちゃんと区別しましょう。

現実問題では、最尤推定量を陽に得ることが難しいことが多いため、対数尤度

$$ \begin{align} l \left({\boldsymbol \theta}, {\boldsymbol x} \right) = \log L \left({\boldsymbol \theta}, {\boldsymbol x} \right) = \sum_{i=1}^{n} \log p \left( x_i | {\boldsymbol \theta} \right) \end{align} $$

を数値的に最適化して \displaystyle \hat {\boldsymbol \theta} ({\boldsymbol x})を求めることが多く、その最適化法の一つがEMアルゴリズムがあります。最適化手法はもちろん他にも色々あります。その他についてはこちら。

hotoke-x.hatenablog.com

EMアルゴリズム

実現値が得られない観測に対する確率変数を \displaystyle Yとすると、観測データの確率密度関数

$$ \begin{align} p ({\boldsymbol x}|{\boldsymbol \theta}) = \int p ({\boldsymbol x},{\boldsymbol y}|{\boldsymbol \theta}) \mathrm{d}{\boldsymbol y} \end{align} $$

と書けます。「観測できない点は考えられる点を全部調べる」という気持ちが積分に表れていますね。

なぜ観測できていない点をわざわざ式に埋め込むかというと、不完全データの対数尤度関数より完全データの対数尤度関数の方が扱いやすいからです。

式の上では完全データとはいえ、実際には観測できてない確率変数が混ざっているのでちょっと扱いを工夫します。それがEMアルゴリズムです。

「なぜそれで良いのか」という細かいことは置いといて、まずはEMアルゴリズムの枠組みをおさらいしましょう。

  1. パラメータの初期値 \displaystyle {\boldsymbol \theta}^{(0)}を適当な値に設定(スケールくらいは当たりをつけて置くと良い)。
  2. 観測値 \displaystyle {\boldsymbol x}とパラメータの推定値 \displaystyle {\boldsymbol \theta}^{(0)}から、観測できない \displaystyle Yについて条件付き期待値を算出。これを疑似的な \displaystyle Yの観測値とする。疑似観測 \displaystyle Yから疑似完全データを構成(Eステップ)
  3. 疑似完全データに基づいて対数尤度 \displaystyle l ({\boldsymbol \theta}, {\boldsymbol x})を最大化するパラメータ \displaystyle {\boldsymbol \theta}^{(1)}を推定(Mステップ)
  4.  \displaystyle {\boldsymbol \theta}^{(1)}を新たなパラメータ初期値としてEステップ、Mステップを繰り返す。

さて、いろいろ気になる点があります。

  1. 疑似的な観測値を構成するのに条件付き期待値を使ってよいのか
  2. 疑似的な完全データの尤度を最大化することが、観測データの尤度を最大にしたことになるのか
  3. 疑似的な完全データによる推測と、完全な観測による推測とで乖離はないのか(バイアスがあるのではないかということ)

1つ目は今回の記事では解説しきれそうにないので後日追記ということにして、2つ目の疑問から確認していきます。

疑似的な完全データの尤度を最大化することが、観測データの尤度を最大にしたことになるのか

観測 \displaystyle {\boldsymbol X}={\boldsymbol x}が与えられたときの尤度 \displaystyle l \left({\boldsymbol \theta}, {\boldsymbol x} \right)

$$ \begin{align} l \left({\boldsymbol \theta}, {\boldsymbol x} \right) &= \log p \left({\boldsymbol x}|{\boldsymbol \theta} \right) \\ &= \log p \left({\boldsymbol x, \boldsymbol y}|{\boldsymbol \theta} \right) - \log p \left({\boldsymbol y}|{\boldsymbol x, \boldsymbol \theta} \right) \end{align} $$

となります。しかし困りました。観測できない \displaystyle {\boldsymbol y}が式に含まれています。仕方がないので適当な \displaystyle {\boldsymbol \theta}={\boldsymbol \theta^{(k)}}として完全データの条件付き期待値を計算しましょう(つまり観測していない部分を周辺化するということ)。適当に \displaystyle {\boldsymbol \theta^{(k)}}を使っても問題ない理由もこの後わかります。

$$ \begin{align} l \left({\boldsymbol \theta}, {\boldsymbol x} \right) &= \int p \left({\boldsymbol y} \middle| {\boldsymbol x, \boldsymbol \theta^{(k)}} \right) \left\{ \log p \left({\boldsymbol x, \boldsymbol y}|{\boldsymbol \theta} \right) - \log p \left({\boldsymbol y}|{\boldsymbol x, \boldsymbol \theta} \right) \right\} \mathrm{d}{\boldsymbol y} \\ &= \mathrm{E}_{\boldsymbol \theta^{(k)} } \left[ \log p \left({\boldsymbol x, \boldsymbol y}|{\boldsymbol \theta} \right) \right] - \mathrm{E}_{\boldsymbol \theta^{(k)} } \left[ \log p \left({\boldsymbol y} \middle| {\boldsymbol x, \boldsymbol \theta} \right) \right] \end{align} $$

これは \displaystyle \boldsymbol \theta^{(k)} \displaystyle \boldsymbol \thetaの関数になっているので、新たに関数 \displaystyle Q, Hを導入して

\begin{align} l \left({\boldsymbol \theta}, {\boldsymbol x} \right) &= Q \left({\boldsymbol \theta}, {\boldsymbol \theta^{(k)}} \right) - H \left( {\boldsymbol \theta}, {\boldsymbol \theta^{(k)}} \right) \end{align}

と書くことができます。また、 \displaystyle Q, Hは、 \displaystyle \boldsymbol \theta^{(k)}によって定まる関数  p\left({\boldsymbol y} \middle| {\boldsymbol x, \boldsymbol \theta^{(k)}} \right) を引数にとる汎関数であるという見方もできます。

さて、ここで \displaystyle Qの式を改めて眺めてみます。 \displaystyle Q

$$ \begin{align} Q \left({\boldsymbol \theta}, {\boldsymbol \theta^{(k)}} \right) &= \int p \left({\boldsymbol y} \middle| {\boldsymbol x, \boldsymbol \theta^{(k)}} \right) \log p \left({\boldsymbol x, \boldsymbol y}|{\boldsymbol \theta} \right) \mathrm{d}{\boldsymbol y} \end{align} $$

です。これは、事前分布に基づいて完全データの対数尤度の期待値を計算しているだけです。完全データの条件付き期待値を用いた対数尤度の最適化ができれば、実際に完全データが得られた場合の対数尤度の値へ収束することがわかります。初期値に適当な \displaystyle {\boldsymbol \theta^{(k)}}を使って問題ないのは、「完全データの条件付き期待値の対数尤度が最適化可能である」という前提からなのでした(実際は局所解へ陥る可能性があることに注意)。

いやいや関数 \displaystyle Hは放っておいていいのかよと思われるかもしれませんが、この後説明します。

疑似的な完全データによる推測と、完全な観測による推測とで乖離はないのか(バイアスがあるのではないかということ)

当然バイアスはあります。実は、疑似データを用いたことによる現実との乖離度合いを表しているのが \displaystyle Hになります。式を眺めてみましょう。

$$ \begin{align} H \left( {\boldsymbol \theta}, {\boldsymbol \theta^{(k)}} \right) &= - \int p \left({\boldsymbol y} \middle| {\boldsymbol x, \boldsymbol \theta^{(k)}} \right) \log p \left({\boldsymbol y}|{\boldsymbol x, \boldsymbol \theta} \right) \mathrm{d}{\boldsymbol y} \\ \end{align} $$

ある程度勉強してきた方ならピンとくるのではないでしょうか。クロスエントロピーの式になっていますね。つまり、理想的な最適解が求まったとすればこの値は自己エントロピーに一致し最小となります。

対数尤度を最大にするのにこの項は最小になってよいのかと思われる方がいるかもしれませんが、その分、完全データの条件付き期待値の対数尤度が大きくなります。そもそも「疑似データを用いたことによる現実との乖離度合いを表している」だけなのです。

疑似的な観測値を構成するのに条件付き期待値を使ってよいのか(pending)

まず、これは個人の見解であり正確な解釈ではない可能性があることをご理解下さい。
実はEMアルゴリズムは暗に指数分布族と呼ばれる確率密度関数を使うことを仮定しているのだと思います。
指数分布で考えると、対数尤度をとるときにexponentialをきれいに打ち消せるのです。
すると全体が見えてきて、十分統計量をその条件付き期待値で置き換えて考える問題に帰着します。
今回は実用重視の記事にしたいので、数理統計の細かいところは次の機会もしくは追記という形にしたいと思います。

参考

小西 貞則、越智 義道、大森 裕浩(2008)「計算統計学の方法」朝倉書店

qiita.com

*1:Dempster, A. P., Laird, N. M. and Rubin, D. B. (1977), Maximum Likelihood from Incomplete Data Via the EM Algorithm. Journal of the Royal Statistical Society: Series B (Methodological), 39: 1-22. doi:10.1111/j.2517-6161.1977.tb01600.x

neovimでAnacondaの環境ごとにpython補完させる

環境:Ubuntu18.04 LTS

Anaconda (Miniconda) のインストール

neovim用環境の作成
conda create -n neovim python=3
pip install neovim

neovimのインストール

sudo add-apt-repository ppa:neovim-ppa/stable
sudo apt-get update
sudo apt-get install neovim

Deinの導入

mkdir -p $XDG_CONFIG_HOME/nvim/dein/toml
touch $XDG_CONFIG_HOME/nvim/dein/toml/init.viml
touch $XDG_CONFIG_HOME/nvim/dein/toml/dein.toml
touch $XDG_CONFIG_HOME/nvim/dein/toml/dein_lazy.toml

init.vimを編集。 典型的なinit.vimなのでコピペで動くと思います。

let g:python3_host_prog = $HOME.'/miniconda3/envs/neovim3/bin/python'

" initialize
if &compatible
  set nocompatible
endif
set runtimepath+=~/.cache/dein/repos/github.com/Shougo/dein.vim
if dein#load_state('~/.cache/dein')
  call dein#begin('~/.cache/dein')
  call dein#load_toml('$XDG_CONFIG_HOME/nvim/dein.toml', {'lazy': 0})
  call dein#load_toml('$XDG_CONFIG_HOME/nvim/dein_lazy.toml', {'lazy': 1})
  call dein#end()
  call dein#save_state()
endif
filetype plugin indent on
syntax enable

if dein#check_install()
  call dein#install()
endif

Deopleteの導入 dein.tomlへ以下を記入

[[plugins]] repo = 'Shougo/deoplete.nvim'
hook_add = 'let g:deoplete#enable_at_startup = 1'

Deoplete-jediの導入 activateしているAnacondaのpython_pathを自動取得してdeopleteへ反映

dein_lazy.tomlへ以下を記入(遅延読み込みさせる)

[[plugins]]
repo = 'deoplete-plugins/deoplete-jedi'
on_ft = 'python'
hook_add = """
let g:deoplete#sources#jedi#python_path = system("which python")[:-2]
"""

system("which python") でpython pathを取得し、[:-2]で改行文字を落としています。

source activateまたはconda activateで環境に入ってneovimを起動すれば、自動的にその環境での補完をしてくれます。

Julia tips #10: マクロを使って動的に構造体を定義する

Juliaでゼロから学ぶDeep Learningという記事を書きました。

hotoke-x.hatenablog.com

実はこちらの実装、ほぼ本家のPythonそっくりになるように書かれています。

というのも私がDeep Learning初心者なので、自己流で進めすぎるとどこが間違っているのかわからなくなってしまうためです。

「Juliaらしさ」も未だつかめておらず、見比べたらわかるんですがほぼPythonです。

実装に四苦八苦する中でマクロについてちょっと調べる機会があったので記事にしておこうと思った次第です。

ただし、Documentを全て読んだわけではないですし、今でも色々わかっていないということはご了承下さい。

(忘れないうちに記事にして消化しておきたかった)

構造体を動的に定義したい

とりあえず何も考えずに以下のコードを実行(ver. 1.1.0です)。

macro make_struct(struct_name, schema...)
    fields=[:($(entry.args[1])::$(entry.args[2])) for entry in schema]
    field_names=[:($(entry.args[1])) for entry in schema[1:2]]
    obj = :OBJ
    __init = [:($obj.$name=$name) for name in field_names] 
    esc(quote mutable struct $struct_name
        $(fields...)
            function (::Type{$(struct_name)})($(field_names...))
                $obj = new()
                $(__init...)
                $obj
            end
        end
    end)
end
println(@macroexpand @make_struct STRUCT_NAME (x,Integer) (y, Integer) (z, Integer))
@make_struct STRUCT_NAME (x, Integer) (y, Integer) (z, Integer)
STRUCT_NAME(1,2)

STRUCT_NAMEという構造体が定義されたはずです。

このマクロは前半でコードの部品生成、後半部分で部品の埋め込みを行っています。

部品生成部分にコメントをつけました。

macro make_struct(struct_name, schema...)
    # 部品(Expr)生成###################################
    fields=[:($(entry.args[1])::$(entry.args[2])) for entry in schema]   # [x::Integer, y::Integer, z::Integer]を生成
    field_names=[:($(entry.args[1])) for entry in schema[1:2]]   # [x::Integer, y::Integer]を生成(内部コンストラクタの引数に使う)
    obj = :OBJ   # このOBJという名前を内部コンストラクタで使う
    __init = [:($obj.$name=$name) for name in field_names]   # [OBJ.x=x, OBJ.y=y]を生成
    # 部品(Expr)生成###################################

    esc(quote mutable struct $struct_name
        $(fields...)
            function (::Type{$(struct_name)})($(field_names...))
                $obj = new()
                $(__init...)
                $obj
            end
        end
    end)
end
println(@macroexpand @make_struct STRUCT_NAME (x,Integer) (y, Integer) (z, Integer))
@make_struct STRUCT_NAME (x, Integer) (y, Integer) (z, Integer)
STRUCT_NAME(1,2)

埋め込み部分が何をしているかは

println(@macroexpand @make_struct STRUCT_NAME (x,Integer) (y, Integer) (z, Integer))

でわかります。

おわりに

ゼロから始めるDeep Learningを実装する中でJuliaの辞書が遅いということに気がついたので、「動的にstructを定義できたら速いのでは?」と思ったので今回のようなことを思いついたのでした。

ちなみに、「Juliaの辞書が遅い」というのは正確には間違っています。

Juliaの辞書は型指定無しもしくは抽象型で作った時にPythonより遅くなるようです。具象型であればPythonより高速です。

これは後ほど記事にする予定です。

Juliaでゼロから作るDeep Learning

もともとPython使いの自分には、ゼロから作るDeep LearningPythonでやっても面白くない!

ということでJuliaでやってみた。

速度が犠牲にならないように多少は考慮したつもりだが、本家Python版の方が速いかも。

なんにも考えなくても速度がでるNumpy凄い...

Juliaは数学的に自然な形でコードが書けて嬉しいが、行列演算をちょっと工夫して書かないとNumpyより遅くなる。

また、実装していて気がついたのだがどうやらdictionaryへのアクセスがJuliaは遅いらしい?

discourse.julialang.org

dictionaryは便利だけど、できるだけstructを使って回すのが良さそうです。

最後までちゃんと読まなかったようで、Juliaの方が早いとの事。失礼しました。

twitter.com

ただ、何にも工夫しないと結局遅いようなので後ほど試してまとめます。

以下にコードがあります。

商用・非商用問わず自己責任でご利用下さい。

github.com

Julia 小技集

完全に自分用
便利な小技をメモしていく

随時更新。

現在の環境

Julia Version 1.1.0
Commit 80516ca202 (2019-01-21 21:24 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: AMD Ryzen 7 1700 Eight-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, znver1)
Environment:
  JULIA_DEPOT_PATH = /opt/julia
  JULIA_PKGDIR = /opt/julia
  JULIA_VERSION = 1.1.0
  • 型を一段階アップグレードする関数
    例:Int64 -> Int128
wilden(x)
  • Array{Array{hoge, 1}, 1}を2次元配列に直したい
vcat(a...)
  • 配列の要素ごとに関数を適用
map([1,2,3]) do x
    何らかの処理
end

xに配列の要素を代入してdo-end間の処理をしてくれる。

  • 多次元配列の成分ごとに関数を適用
mapslices(関数, 対象配列,dims=次元の指定)

次元の指定では[1,2]のように複数の次元を指定できる。

Meta programming

動的に構造体を定義

discourse.julialang.org

例)

macro make_struct(struct_name, schema...)
    fields=[:($(entry.args[1])::$(entry.args[2])) for entry in schema]
    esc(quote mutable struct $struct_name
        $(fields...)
        end
    end)
end
println(@macroexpand @make_struct STRUCT_NAME (x,Int32) (y, Float32))
@make_struct STRUCT_NAME (x, Int32) (y, Float32)

その他

デバッグ

マルチGPU

  • 関数のベクトル、多次元配列への拡張
function f(x)
    xに対する処理
end


function g(x::Vector{T})  where {T<:AbstractFloat}
    xに対する処理
end

# 行列の成分ごとに上記の関数gを適用
function g(x::Matrix{T}) where {T<:AbstractFloat}
    mapslices(g, x, dims=次元の指定)
end

構造に対しては明示的に、数値型に関してはパラメトリック型で定義すると使いやすい(個人の感想) f(x)を要素ごとに適用したいならf.(x)とすれば良い。
@viewなどで参照を渡す場合は具体型ではなく抽象型で定義したほうが良い(SubArrayがわたされるので)

【初心者】Deep Learning用語まとめ

「ゼロから作るDeep Learning」を読んでいる。 一気読みできればいいが、そういう訳ににもいかないので簡単な用語を忘れることがある。 そんな時に本を開くのは億劫なので、ここに簡単な用語をまとめておく。


ニューラルネットワークに入る前に...

  • パーセプトロン
    多入力単出力の構造 信号は0, 1 (ステップ関数) 単純パーセプトロンでAND, NAND, ORゲートを作れる。(多層パーセプトロンにすればXORも可能)

  • バイアス
    上流からの重み付き和とは別に常に加算される値

  • 白色化(whitening)
    データ全体の分布の形状を均一にする方法

  • one-hot-label
    ラベルの表現方式のひとつ。
    ラベルが0~9まであり、正解ラベルが7の時、[0,0,0,0,0,0,1,0,0]のように表示する。

ニューラルネットワーク

  • 活性化関数
    入力信号の総和を出力信号へ変換する。単純パーセプトロンではステップ関数を使っていた。ステップ関数以外を使えばニューラルネットワークと呼ぶようになるらしい(ゼロから作るDeep Learning, p44より)。活性化関数には線形関数を用いてはならないことに注意(多層にする意味がなくなる。単純な合成写像を考えてみればよい。)。

    • ステップ関数
      閾値を超えたら1、それ以外は0を出力する関数。
      all-or-nothingな信号という点で一番ニューロンっぽい気がする(個人の感想)

    • シグモイド関数
       {} $$ \begin{align} \frac{1}{1+\exp(-x)} \end{align} $$ 値を[0, 1]に押し込めるのが嬉しい関数。
      滑らかな関数(微分可能)であるところもポイント。

    • ReLU (Rectified Linear Unit) 関数
      閾値を超えたら入力をそのまま出力、それ以外は0を出力する関数。

  • 出力層の活性化関数
    出力層だけ活性化関数を特別扱いする。

    • 恒等関数
      入力をそのまま出力する。回帰問題に用いる。

    • ソフトマックス関数
       \displaystyle y_kを入力信号、 \displaystyle y_kを出力信号として  {}
      $$ \begin{align} y_k=\frac{\exp(a_k)}{\sum_{i=1}^{n} \exp(a_i)} \end{align} $$
      これも出力を[0, 1]に押し込めるのがうれしい。
      オーバーフロー対策に定数  \displaystyle C を導入して  {} $$ \begin{align} y_k = \frac{\exp(a_k-C)}{\sum_{i=1}^{n}{\exp(a_i-C)}} \end{align} $$
      とするのが良いらしい。 \displaystyle Cには \displaystyle \{a_1, a_2, \ldots, a_n\}の最大値を持ってくることが一般的。

損失関数

ニューラルネットワークの学習に用いる損失関数いろいろ。 以下ではyを出力(ラベルごとの確率)、tを正解ラベルのone-hot-labelとする。kはデータの次元数とする。

  • 二乗和誤差
    説明不要。お馴染み。
     {} $$ \begin{align} E = \frac{1}{2} \sum_{k} (y_k-t_k)^2 \end{align} $$

  • 交差エントロピー誤差
     \displaystyle t_kに関する \displaystyle -\log y_kの期待値を計算していることになる。情報理論の観点で言えば、 \displaystyle t_kに関する \displaystyle y_kの平均符号長を計算していることになる。情報理論ででてくるエントロピーと異なり自然対数を用いているが、本質的には同じものなので注意。今回は細かい説明はしないので、ちゃんと理解したい場合は符号化理論を勉強しましょう。
     {} $$ \begin{align} E = - \sum_{k} t_k \log y_k \end{align} $$

ニューラルネットワークの学習

  • ミニバッチ学習
    ニューラルネットワークの学習において、全データについて損失関数を計算するのは現実的でない。
    そこで一部のデータをランダムサンプリングして損失関数を求め、それを全データの代表値とみなして学習を行う方法のこと。
    ランダムサンプリングしているといっても、そもそもどうやってデータが集められたのかには注意する必要があるなと思った(個人の感想)。

  • 勾配法
    勾配法については以下を参照

hotoke-x.hatenablog.com

  • 確率的勾配降下法SGD; Stochastic Gradient Descent)
    ミニバッチとしてランダムサンプリングされたデータを使用した勾配降下法のこと。
    以下の式に従ってパラメータ \boldsymbol{W}の更新を行う。  {} $$ \begin{align} \boldsymbol{W} \leftarrow \boldsymbol{W}- η\frac{\partial L}{\partial \boldsymbol{W}} \label{s_g_d} \end{align} $$  Lは損失関数、 \etaは学習率である。

  • エポック(epoch)
    学習の単位のこと。
    訓練データをすべて使い切った時の回数に対応。
    1000個のデータに対して50個のミニバッチで学習する場合、20回のミニバッチ学習=1 epochとなる。
    エポックごとに訓練データの認識精度とテストデータの認識精度を比較する。この時、認識精度に差があると過学習が起きていることになる。

誤差逆伝播

数値微分を用いたニューラルネットワークの学習は実装が容易。一方で計算に時間がかかる。誤差逆伝播法で勾配の計算を効率よく行う。

  • 計算グラフで考える
    ノードごとの局所的な計算で済む。 局所的な計算結果を保持できるので、逆伝播によって微分が計算できる。
    ただし、バッチごとに計算を行っている場合行列演算になるので、行列積が可換でないことに注意。
    縮約記法を使えば簡単(以下の記事を参照)。

hotoke-x.hatenablog.com

hotoke-x.hatenablog.com

hotoke-x.hatenablog.com

  • 勾配確認
    数値微分による逆誤差伝播の実装と、解析解による逆誤差伝播の実装の値を比べる事。
    解析解による実装を行った時、解が間違っていないことを確認できる。
Momentum

SGDの\eqref{s_g_d}に慣性項を導入した。前の状態を引き継ぐことでSGDの欠点の一つであるパラメータの過剰更新を抑制する。パラメータ更新の式は以下の通り。

 {} $$ \begin{align} \boldsymbol{v}_{t+1} &\leftarrow \alpha \boldsymbol{v}_{t} - η \frac{\partial L}{\partial \boldsymbol{W}} \\ \boldsymbol{W} &\leftarrow \boldsymbol{W} + \boldsymbol{v}_{t+1} \end{align} $$

AdaGrad

Adaptive Gradientの略?
学習率を減衰させていく方法。初めは大雑把に、時間が経つにつれ小さく探索する。パラメータの更新式は以下の通り。

 {} $$ \begin{align} \boldsymbol{h} &\leftarrow \boldsymbol{h} + \frac{\partial L}{\partial \boldsymbol{W}} \odot \frac{\partial L}{\partial \boldsymbol{W}} \\ \boldsymbol{W} &\leftarrow \boldsymbol{W} -η \frac{1}{\sqrt{\boldsymbol{h}}} \frac{\partial L}{\partial \boldsymbol{W}} \end{align} $$

ここで、 \displaystyle \odotアダマール積(要素ごとの積)である。 \displaystyle \frac{\partial L}{\partial \boldsymbol{W}}が大きい要素は、次の学習率がより小さくなるように調整される。ループによって学習率が下がる一方なので、最終的には学習が進まなくなるのが欠点。これを改善したものにRMSPropという方法があるらしい。

Adam

MomentumとAdaGradの融合(AdaptiveなMomentumということだと思う)。
詳細は原著*1を読まないとわからない。

重みの初期値問題

小さい値をランダムに与える。
過学習を防ぐために小さい値で重みを学習したいという願いが込められている。
PRMLでみた正則化に似ていると思った(制約を与えたわけでは無いが、事前分布に願いを込める意味で)。

勾配消失

例えばシグモイド関数の値が0, 1付近だと、微分値は小さくなる。アクティベーションが0, 1に偏っていると、逆誤差伝播の時に勾配が消失し、学習が進まなくなる。特定の値にアクティベーションの値が集中するのも良くない。ノード全体が同じ値を出力するということなので、ノードが複数存在している意味が無くなってしまう。

  • Xavierの初期値
    前層ノードの個数 \displaystyle nを使って \displaystyle \frac{1}{\sqrt{n}}標準偏差を持つ分布を用いる(なんだか中心極限定理と関係がありそう)。 活性化関数が線形であることが前提らしい。sigmoid関数、tanh関数は対称で中央付近が線形とみなせるとのこと(それでいいの??)。詳細は原著参照*2

  • Heの初期値
    ReLUを用いる場合はHeの初期値が良いとの事。 \displaystyle \frac{2}{\sqrt{n}}標準偏差を持つ分布を用いる。ReLUは負の値が全部0に押し付けれられるので、Xavierの初期値よりも広がりを持たせるということらしい(ReLUが青天井の値を持つことを考えると3でも5でも良い気がする...)。 こちらも詳細は原著参照*3

  • Batch Normalization
    各層ごとにアクティベーション分布に広がりを持たせたいなら、各層ごとに調整したらいいじゃん。と思ったら書いてあった。データと各層ごとの出力を \mathcal{N}(0,1)に標準化する。 活性化関数の前と後どちらでBatch Normalizationするかは議論になっているらしい。
    個人的には活性化関数の後でBatch Normalizationの方が自然な気がしている。理由は、各層の入力に広がりを持たせたいのが目的なら、その直前で分布の調整をした方が自然に感じるから。活性化関数の前で標準化してしまうと、活性化関数によって再び分布が偏る可能性が排除できない。だからこそXavierの初期値では線形関数である前提が必要だったのだと思う。活性化関数の後なら活性化関数が非線形でも使えるかもしれない(という期待を込めた感想)。
    histogram equalizationとか使ったらどうなるんだろう..。アルゴリズム的に破綻するかもしれないけど。

正則化

過学習を抑制する方法。

過学習

訓練データのみに最適化されてしまって、テストデータに全然対応できない状態。

Weight decay

損失関数に、重み \displaystyle \boldsymbol{W}のL2ノルム正則化 {} $$ \begin{align} \frac{1}{2}\lambda\boldsymbol{W}^2 \end{align} $$ を加算して学習を行う。

 \displaystyle \lambdaはハイパーパラメータ。

(ABICのノリで \displaystyle \lambdaの自動決定とかできないだろうか)

Dropout

ノードをランダムに除去しながら学習する方法。ノードを除去することで毎回異なるモデルを生成するので疑似的なアンサンブル学習として見れるとのこと。

検証データ(validation data)

ハイパーパラメータを調整するためのデータ。 テストデータでハイパーパラメータを調整すると、ハイパーパラメータがテストデータに関して過学習を起こすので別に用意する必要がある。

ハイパーパラメータの最適化

書籍を読んだ感じでは、グリッドサーチっぽかった。ただ、ハイパーパラメータ探索ではグリッドサーチよりもランダムに探索したほうが良い結果になるとの事。ハイパーパラメータごとに結果に与える影響度が異なるからだとか。ハイパーパラメータを探索するスケールが適当でないことが多いってことかな?

ベイズ的に最適化する方法もあるらしい*4MCMCマルコフ連鎖モンテカルロ法)かABC(近似ベイズ計算)だろうか。パラメトリックに最適化するなら事前分布の設定がまた問題になりそう。何はともあれ、これも原著を読まないことにはわからない。

畳み込みニューラルネットワーク(Convolutional neural network: CNN)

全結合層だと、データの位置情報を無視して学習してしまう。 CNNでは位置情報も学習できる(ということになっている)。

Convolution レイヤ

フィルタを通すことで、出力サイズを落とす(次元という言い方は正しくないかもしれない)。 例えばMNISTなら、28x28だった画像をフィルタを通すことで26x26にするような操作。 3次元データ(例えばRGB画像)に対してフィルタを適用する場合は、フィルタも3次元にし、出力される3層のz軸方向に和をとる(なので出力は2次元)。

フィルタを複数個用意すれば出力は3次元になる。バッチの個数分同じフィルタを用意すれば一括処理が可能になる。

パディング

Convolution レイヤばかり通していたら出力サイズが小さくなる一方なので、入力行列の外側を何らかの値(例えば0)で埋めて、畳み込みによる出力サイズの現象を防ぐ方法。

ストライド

フィルターを適用する位置の間隔

Pooling レイヤ

空間を小さくする演算処理。
Poolingのウィンドウサイズとストライドは同じ値にするのが通常。
学習パラメータがない。
チャネル数も変化しない。
微小な位置変化に対してロバスト。ちょっとぐらい入力データが動いても、フィルタの特性(フィルタ適用範囲の最大値をとってくるなど)によって出力結果が変化しにくくなる。

画像認識ではMax poolingが一般的らしい。Average poolingもあるが、確かに外れ値に引っ張られるのはちょっと使いにくそう(素人的感想)。 入力をFFTした画像とかにしたらどうなるのだろう。ノイズはノイズで学習しやすくなる気がする(これも素人的感想)

代表的なCNN

LeNet

  • 畳み込み層とプーリング層を重ねて、最後に全結合層を使う。
  • 活性化関数はReLU
    オリジナルではプーリング層はただサブサンプリングするだけで、活性化関数もsigmoidだった。現在はMax poolingが主流らしい。

AlexNet

  • 活性化関数にReLUを使用
  • 局所正規化を行う(上で書いたBatch Normalizationのようなものだと思う)。
  • Dropoutを使用する

Deep Learning

VGG*5

  • ReLUを使用
  • 重みの初期値はHeの初期値を使用(ReLUを使用するから?)
  • 小さなフィルタ(3x3)による畳み込み(些細な特徴でも拾いたいから?)
  • Adamでパラメータ更新
  • 全結合層の後にDropout(上層でメチャクチャ特徴を学習させて、最後にアンサンブル学習することで変な学習を回避??)

データを回転させたり動かしたりして疑似的にデータを増やす「データ拡張」も認識精度向上に貢献するとのこと。

層を増やす利点

  1. 単純にパラメータが減る。
    5x5のフィルターによる出力と、3x3のフィルタを2回通すのとでは、出力サイズは同じだが、後者の方がフィルタの要素数が小さくて済む。
  2. より少ないデータで学習ができる 層が深いと高度な情報が抽出できるため

他のモデルの例

  • GoogLeNet*6
    VGGとは違って、層のつながり自体が2次元的(インセプション構造というらしい)。1x1フィルターの畳み込みでサイズを削減することによって計算速度を向上させているとの事。
    ただ、層のつながりを2次元的にする恩恵はなんなのだろうか...。これも要原著参照ということで。

  • ResNet*7
    入力層を畳み込み層をすっ飛ばして、出力と合算させる構造が入っている。
    層を深くしすぎて抽出できる特徴がスッカスカになって来たころに、昔を思い出させようとするお気持ちなのだろうか。もしそうなら、層の深ささえ適切なら同じような結果が得られそうな気がした。また、合算処理によって新たな特徴マップが得られるということなら、データ拡張でも似たような結果が得られそうだとも思った。知らんけど。

物体検出の例

  1. R-CNN*8
    高速版もあるらしい*9

セグメンテーションの例

  1. FCN*10
    全結合層が畳み込み層からなる。

画像キャプション生成

代表的なものとしては、NIC*11(Recurrent Neural Netwotk (RNN)の一種)がある。 再帰的な構造があり、時系列や言語などの連続データに使われる。

画像スタイル変換

絵のテイストを変えたりする*12

画像生成

スタイル変換とは異なり、完全に新しい画像を生成する。
DCGAN*13

Deep Q-Network(DQN, 強化学習の一つ)

  1. エージェントが行動する
  2. 行動に対する報酬を得る これを繰り返してずーっと学習させる。
    DQNでは最適行動価値関数とかいう関数をCNNで近似するらしい。これも論文読まないとわからないかも*14

宿題

参考書籍

斎藤 康毅(2016)「ゼロから作るDeep Learning-Pythonで学ぶディープラーニングの理論と実装」株式会社オライリー・ジャパン

*1:D. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014.

*2:X. Glorot and Y. Bengio, “Understanding the difficulty of training deep feedforward neural networks,” Proceedings of the thirteenth international conference on artificial intelligence and statistics, 2010.

*3:K. He, X. Zhang, S. Ren, and J. Sun, “Delving deep into rectifiers: Surpassing human-level performance on imagenet classification,” Proceedings of the IEEE international conference on computer vision, 2015.

*4:J. Snoek, H. Larochelle, and R. Adams, “Practical bayesian optimization of machine learning algorithms,” Advances in neural information processing systems, 2012.

*5:K. Simonyan and A. Zisserman, “Very Deep Convolutional Networks for Large-Scale Image Recognition,” arXiv preprint arXiv:1409.1556, 2014.

*6:C. Szegedy, W. Liu, Y. Jia, P. Sermanet, S. Reed, D. Anguelov, D. Erhan, V. Vanhoucke, and A. Rabinovich, “Going deeper with convolutions,” Proceedings of the IEEE conference on computer vision and pattern recognition, 2015.

*7:K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” Proceedings of the IEEE conference on computer vision and pattern recognition, 2016

*8:R. Girshick, J. Donahue, T. Darrell, and J. Malik, “Rich feature hierarchies for accurate object detection and semantic segmentation,” Proceedings of the IEEE conference on computer vision and pattern recognition, 2014.

*9:S. Ren, K. He, R. Girshick, and J. Sun, “Faster r-cnn: Towards real-time object detection with region proposal networks,” Advances in neural information processing systems, 2015

*10:E. Shelhamer, J. Long, and T. Darrell, “Fully Convolutional Networks for Semantic Segmentation,” Ieee T Pattern Anal, vol. 39, no. 4, pp. 640–651, 2016.

*11:O. Vinyals, A. Toshev, S. Bengio, and D. Erhan, “Show and tell: A neural image caption generator,” Proceedings of the IEEE conference on computer vision and pattern recognition, 2015.

*12:L. Gatys, A. Ecker, and M. Bethge, “A neural algorithm of artistic style,” arXiv preprint arXiv:1508.06576, 2015.

*13:A. Radford, L. Metz, and S. Chintala, “Unsupervised representation learning with deep convolutional generative adversarial networks,” arXiv preprint arXiv:1511.06434, 2015.

*14:Mnih, K. Kavukcuoglu, D. Silver, A. Rusu, J. Veness, M. Bellemare, A. Graves, M. Riedmiller, A. Fidjeland, G. Ostrovski, S. Petersen, C. Beattie, A. Sadik, I. Antonoglou, H. King, D. Kumaran, D. Wierstra, S. Legg, and D. Hassabis, “Human-level control through deep reinforcement learning,” Nature, vol. 518, no. 7540, p. 529, 2015.

*15:D. Silver, A. Huang, C. Maddison, A. Guez, L. Sifre, G. van den Driessche, J. Schrittwieser, I. Antonoglou, V. Panneershelvam, M. Lanctot, S. Dieleman, D. Grewe, J. Nham, N. Kalchbrenner, I. Sutskever, T. Lillicrap, M. Leach, K. Kavukcuoglu, T. Graepel, and D. Hassabis, “Mastering the game of Go with deep neural networks and tree search,” Nature, vol. 529, no. 7587, p. 484, 2016.

Windows10+Anaconda+PyTorch

Chainerを使いたかったけどWindowsは公式サポートでない上に、CUDA、cuDNNのインストールがやたら面倒くさい...。 他のライブラリでもっと簡単に使えるものはないかと探していたら、PyTorchならAnaconda経由で環境が整うらしいのでメモ。
linuxでもmacOSでも同様の手順で行けると思います(試してない)。

ちなみにtensorflowもAnacondaだけで環境が整います。

hotoke-x.hatenablog.com


以下の記事によると「PyTorchはCUDAをマニュアルインストールする必要はないよ!」とのこと。

discuss.pytorch.org

インストール方法

前提としてAnacondaのインストールは終わっているものとします。

PyTorchの公式

で、以下のようにポチポチと自分の環境を選択します。

f:id:hotoke-X:20190213161505p:plain
pytorch

「CUDA」のところは、自分が使いたいバージョンを選べば良いようです。すると、「Run this Command:」にcondaコマンドが表示されるのでコピペして実行しましょう。これだけです。

正しくインストールできたかは以下のpythonコードで確認できます。

import torch 
print(torch.cuda.is_available())

Trueが返ってくれば正常にインストールできています。

CUDAのバージョンとかドライバのバージョンとか、ubuntuカーネルとか気にしなくてよいのが最高にうれしいですね。

テンソルを使いこなしたい3(行列の行列微分)

hotoke-x.hatenablog.com

hotoke-x.hatenablog.com

の続き。
本記事ではテンソル微分演算についてまとめる。行列の行列微分がややこしかったので、縮約記法で計算すればスッキリできるのではと思ったのが事の発端。ようやくここまで来た。

準備

スカラーのベクトル微分、ベクトルのスカラー微分、ベクトルのベクトル微分スカラーの行列微分を以下のように定義する。また、それぞれの最右辺は縮約記法での表現を表す。 {} $$ \begin{align} \frac{\mathrm{d}f}{\mathrm{d}\boldsymbol{x}} &= \left( \begin{array}{c} \frac{\mathrm{d}f}{\mathrm{d}x_{1}} \\ \frac{\mathrm{d}f}{\mathrm{d}x_{2}} \\ \vdots \\ \frac{\mathrm{d}f}{\mathrm{d}x_{n}} \end{array} \right) = \partial{}_{i} f \\ \frac{\mathrm{d}\boldsymbol{f}}{\mathrm{d}x} &=\left(\begin{array}{c} \frac{\mathrm{d}f_{1}}{\mathrm{d}x} \\ \frac{\mathrm{d}f_{2}}{\mathrm{d}x} \\ \vdots \\ \frac{\mathrm{d}f_{n}}{\mathrm{d}x} \end{array} \right) = \partial{} f_{i} \\ \frac{\mathrm{d}\boldsymbol{f}}{\mathrm{d}\boldsymbol{x}} &= \left(\begin{array}{cccc} \frac{\mathrm{d}f_{1}}{\mathrm{d}x_{1}} & \frac{\mathrm{d}f_{1}}{\mathrm{d}x_{2}} & \ldots & \frac{\mathrm{d}f_{1}}{\mathrm{d}x_{n}} \\ \frac{\mathrm{d}f_{2}}{\mathrm{d}x_{1}} & \frac{\mathrm{d}f_{2}}{\mathrm{d}x_{2}} & \ldots & \frac{\mathrm{d}f_{2}}{\mathrm{d}x_{n}}\\ \vdots & \vdots & \ddots & \vdots \\ \frac{\mathrm{d}f_{m}}{\mathrm{d}x_{1}} & \frac{\mathrm{d}f_{m}}{\mathrm{d}x_{2}} & \ldots & \frac{\mathrm{d}f_{m}}{\mathrm{d}x_{n}} \end{array} \right) = \partial{}_{j} f_{i} \label{vec_by_vec} \\ \frac{\mathrm{d}f}{\mathrm{d}\boldsymbol{X}} &= \left(\begin{array}{cccc} \frac{\mathrm{d}f}{\mathrm{d}x_{11}} & \frac{\mathrm{d}f}{\mathrm{d}x_{12}} & \ldots & \frac{\mathrm{d}f}{\mathrm{d}x_{1n}} \\ \frac{\mathrm{d}f}{\mathrm{d}x_{21}} & \frac{\mathrm{d}f}{\mathrm{d}x_{22}} & \ldots & \frac{\mathrm{d}f}{\mathrm{d}x_{2n}}\\ \vdots & \vdots & \ddots & \vdots \\ \frac{\mathrm{d}f}{\mathrm{d}x_{m1}} & \frac{\mathrm{d}f}{\mathrm{d}x_{m2}} & \ldots & \frac{\mathrm{d}f}{\mathrm{d}x_{mn}} \end{array} \right) = \partial{}_{ij} f \label{scalar_by_matrix} \end{align} $$

また、\eqref{vec_by_vec}より  {} $$ \begin{align} \partial{}_{j} x_{i}= \delta_{ij} \end{align} $$ である(この後使う)。

なお、\eqref{vec_by_vec}は転置した形として定義されていることもある。行列表現のまま計算を進めようとすると、行列積が可換でないことに注意する必要がある。縮約記法で計算すれば気にしなくて良くなる。

さて、ベクトルの行列微分、行列のベクトル微分、行列の行列微分はどうなるのだろう。行列の行列微分がわかれば他もノリでわかると思うので、上記と同様に縮約記法を使って行列の行列微分の定義から、ニューラルネットワークの逆誤差伝播の式を導出するところまでやってみる。

行列の行列微分

行列 \displaystyle Xの自身による微分は縮約記法を使って以下のように定義される。

 {} $$ \begin{align} \frac{\partial X_{kl}}{\partial X_{ij}} = \delta_{ik}\delta_{li} \end{align} $$

要は、 \displaystyle X_{kl} \displaystyle X_{kl}微分したときは \displaystyle 1それ以外は \displaystyle 0と言っているだけである。「じゃあ全成分 \displaystyle 1になるじゃん」とか一瞬思うかもしれない(自分は思った)がそれは大間違い。 \displaystyle k,lそれぞれについて \displaystyle i,jを走査して計算するので4階テンソルになっている(添字も4つあるし)。

逆誤差伝播

以下のような単純なニューラルネットワークを考える。 Lは損失関数でスカラーである。形式的に書いてるだけなので実態は何でも良い。

f:id:hotoke-X:20190210064732p:plain

これを式で書けば

$$ \begin{align} Y &= WX + B \\ X &= \left(\begin{array}{c} x_{1} \\ \vdots \\ x_{n} \end{array} \right), W = \left(\begin{array}{ccc} w_{11} & \ldots & w_{1n} \\ \vdots & \ddots & \vdots \\ w_{m1} & \ldots & w_{mn} \end{array} \right), B = \left(\begin{array}{c} b_{1} \\ \vdots \\ b_{m} \end{array} \right) \end{align} $$

である。計算グラフにして逆誤差伝播を赤で書き入れると以下のようになる。

f:id:hotoke-X:20190210082429p:plain

ナニコレ。特に \displaystyle \left( \frac{\partial Y}{\partial X}\right)^{\mathrm{T}} \frac{\partial L}{\partial Y}。そりゃ \displaystyle Xの次元に一致しなきゃいけないから逆算すればわかるけど気持ち悪い。一方で、 \displaystyle \frac{\partial L}{\partial Y} \frac{\partial Y}{\partial W}も連鎖律としてはスッキリしているが、ベクトルの行列微分がある。しかも \displaystyle X^{\mathrm{T}}についてもなんで転置されてるんだかよくわからない。

逆誤差伝播(縮約記法でスッキリ ver.)

では縮約記法で \displaystyle \frac{\partial L}{\partial X}, \frac{\partial L}{\partial W}を計算してみよう。

 {} $$ \begin{align} Y_{i} = W_{ij} X_{j} + B_{i} \end{align} $$ として、  {} $$ \begin{align} \frac{\partial L}{\partial X_{j}} &= \frac{\partial L}{\partial Y_{k}} \frac{\partial Y_{k}}{\partial X_{j}} = \frac{\partial L}{\partial Y_{k}} \frac{\partial \left(W_{kl} X_{l} + B_{k} \right)}{\partial X_{j}} = \frac{\partial L}{\partial Y_{k}} W_{kl} \delta_{lj} \\ &= \frac{\partial L}{\partial Y_{k}} W_{kj} = W_{kj} \frac{\partial L}{\partial Y_{k}} = W^{\mathrm{T}} \frac{\partial L}{\partial Y} \\ \frac{\partial L}{\partial W_{ij}} &= \frac{\partial L}{\partial Y_{k}} \frac{\partial Y_{k}}{\partial W_{ij}} = \frac{\partial L}{\partial Y_{k}} \frac{\partial \left(W_{kl} X_{l} + B_{k} \right)}{\partial W_{ij}} = \frac{\partial L}{\partial Y_{k}} X_{l} \delta_{ik}\delta_{lj} \\ &= \frac{\partial L}{\partial Y_{i}} X_{j} = \frac{\partial L}{\partial Y} X^{\mathrm{T}} \end{align} $$

最右辺は行列表現で書いた。今まで悩んでたのが馬鹿らしくなるくらい簡単に計算できた。

まとめ

もう全部これで良くない?

とはいえ、実装と実行速度を考えると行列表現が得られた方が実用的であることは明白。なので、縮約記法で計算してから行列表現に戻してやるのが実務的なやり方なんだと思う。

追記:
前回記事でまとめたエディントンのイプシロンは全く必要なかった(笑)。まぁでも知っていることは良い事だ。

参考

http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3274/pdf/imm3274.pdf

テンソルを使いこなしたい2(レヴィ=チヴィタ記号)

hotoke-x.hatenablog.com

の続き。一気に書くとまとまりがないので分けた。

本記事ではエディントンのイプシロン(レヴィ=チヴィタ記号)についてまとめる。これにより、ベクトル積、行列式、三重積などがスッキリ表現できる。

エディントンのイプシロン(レヴィ=チヴィタ記号)

添字が (1, 2, 3)の偶置換なら 1、奇置換なら -1、それ以外は0とするもの。 Wikipediaではテンソルであると言い切っているが、 レヴィ・チヴィタの記号 [物理のかぎしっぽ]では、疑テンソルであると説明している(ただ、あまり実用上は気にしなくて良さそう)。以下のように定義されている。

$$ \begin{align} \varepsilon_{i j k}=\left\{\begin{array}{ll} 1 & \left( (i, j, k)\in \left\{(1,2,3), (2,3,1), (3,1,2)\right\} \right) \\ -1 & \left( (i, j, k)\in \left\{(1,3,2), (3,2,1), (2,1,3)\right\} \right) \\ 0 & \left( \text {otherwise} \right) \end{array}\right. \end{align} $$

すなわち

$$ \begin {align} \varepsilon_{i j k}=\varepsilon_{j k i}=\varepsilon_{k i j}=-\varepsilon_{i k j}=-\varepsilon_{j i k}=-\varepsilon_{k j i} \label{prop_epsilon} \end {align} $$

3つの直交した単位ベクトル( \boldsymbol{e}_{i}, \boldsymbol{e}_{j}, \boldsymbol{e}_{k} )を使えば、スカラー三重積でも表現できる。

$$ \begin{align} \left| \begin{array}{ccc} \boldsymbol{e}_{i} & \boldsymbol{e}_{j} & \boldsymbol{e}_{k} \end{array} \right| :=\boldsymbol{e}_{i} \cdot\left(\boldsymbol{e}_{j} \times \boldsymbol{e}_{k}\right)=\varepsilon_{i j k} \label{triple} \end{align} $$

スカラー三重積の性質からも、\eqref{prop_epsilon}が成り立っていることがわかる。さらに、\eqref{triple}より また、クロネッカーのデルタを用いれば以下のような性質があることが示せる。

$$ \begin{align} \varepsilon_{i j k} \varepsilon_{l m n} &=\operatorname{det} \left[ \begin{array}{ccc} {\delta_{i l}} & {\delta_{i m}} & {\delta_{i n}} \\ {\delta_{j l}} & {\delta_{j m}} & {\delta_{j n}} \\ {\delta_{k l}} & {\delta_{k m}} & {\delta_{k n}} \end{array} \right] \\ &=\delta_{i l} \left(\delta_{j m} \delta_{k n}-\delta_{j n} \delta_{k m} \right)+\delta_{i m} \left(\delta_{j n} \delta_{k l}-\delta_{j l} \delta_{k n}\right)+\delta_{i n}\left(\delta_{j l} \delta_{k m}-\delta_{j m} \delta_{k l} \right) \label{mul_epsilon} \end{align} $$

証明は具体的に i, j, k, l, m, nに数字を入れて計算してみれば良い。例えば i=l=2の時、考えられる (j,k,m,n)のパターンは

$$ \begin{align} \left(j, k, m, n \right) \in \left\{\left(1, 3, 1, 3\right), \left(1, 3, 3, 1\right), \left(3,1, 1, 3\right), \left(3,1, 3,1\right)\right\} \end{align} $$

となる。なお、エディントンのイプシロンの定義より値が0となる組み合わせは無視した。 この時、\eqref{mul_epsilon}は 0となる項を省略すると

$$ \begin{align} \delta_{il} \left(\delta_{jm}\delta_{kn} - \delta_{jn}\delta_{km} \right) &= \delta_{22} \left(\delta_{jm}\delta_{kn} - \delta_{jn}\delta_{km} \right) \end{align} $$

となるので、\eqref{mul_epsilon}が成り立つことを簡単に確認できる。他の条件でも同様である。 また、\eqref{mul_epsilon}より以下のような性質も導ける*1

$$ \begin{align} \sum_{k} \varepsilon_{i j k} \varepsilon_{l m k} &=\operatorname{det} \left[ \begin{array}{ll} {\delta_{i l}} & {\delta_{i m}} \\ {\delta_{j l}} & {\delta_{j m}}\end{array}\right] \\ &=\delta_{i l} \delta_{j m}-\delta_{i m} \delta_{j l} \\ \sum_{j, k} \varepsilon_{i j k} \varepsilon_{l j k} &=2 \delta_{i l} \\ \sum_{i, j, k} \varepsilon_{i j k} \varepsilon_{i j k}&=6 \end{align} $$

ちょっと実践

ベクトルの外積がシンプルになる。

$$ \begin{align} \boldsymbol{a} \times \boldsymbol{b} = \sum_{i, j} \varepsilon_{i j k} a_j b_k \end{align} $$

3次元の行列式スカラー三重積)もコンパクトに表現できる。

$$ \begin{align} \operatorname{det} \left[ \begin{array}{lll} {a_{11}} & {a_{12}} & {a_{13}} \\ {a_{21}} & {a_{22}} & {a_{23}} \\ {a_{31}} & {a_{32}} & {a_{33}}\end{array}\right] &=\operatorname{det} \left[ \begin{array}{ccc} {\sum_{i} \delta_{1 i} a_{i 1}} & {\sum_{j} \delta_{1 j} a_{j 2}} & {\sum_{k} \delta_{1 k} a_{k 3}} \\ {\sum_{i} \delta_{2 i} a_{i 1}} & {\sum_{j} \delta_{2 j} a_{j 2}} & {\sum_{k} \delta_{2 k} a_{k 3}} \\ {\sum_{i} \delta_{3 i} a_{i 1}} & {\sum_{j} \delta_{3 j} a_{j 2}} & {\sum_{k} \delta_{3 k} a_{k 3}}\end{array}\right] \\ &= \sum_{i, j, k} \operatorname{det} \left[ \begin{array}{lll} {\delta_{1 i}} & {\delta_{1 j}} & {\delta_{1 k}} \\ {\delta_{2 i}} & {\delta_{2 j}} & {\delta_{2 k}} \\ {\delta_{3 i}} & {\delta_{3 j}} & {\delta_{3 k}}\end{array}\right] a_{i 1} a_{j 2} a_{k 3} \\ &=\sum_{i, j, k} \varepsilon_{i j k} a_{i 1} a_{j 2} a_{k 3} =\sum_{i, j, k} \varepsilon_{i j k} a_{1 i} a_{2 j} a_{3 k} \end{align} $$

最後の式変形は \operatorname{det} A = \operatorname{det} A^{\mathrm{T}}であることを利用した。

その他

実は添字の書き方には下付きと上付きで使い分けるルールがある。
本記事で扱う範囲を遥かに超えるので説明は省略した(というか説明できない)。
代わりに参考になるページ、資料へのリンクを貼っておく。

  1. 共変ベクトルと反変ベクトル [物理のかぎしっぽ]
    添字について解説してくれている
  2. http://takeno.iee.niit.ac.jp/~shige/math/lecture/basic2/data/detprf1.pdf
    導入は天下りだが、的p.9で添字の扱いに触れている。また、定理7で $$ \begin{align} |A B|=\left\{\begin{array}{ll} {0} & (n<m) \\ {|A||B|} & (n=m) \\ {\sum_{i_{1}<\cdots<i_{m}} \left|A_{i_{1}, \ldots, i_{m}}| \left|B^{i_{1}, \ldots, i_{m}}| \right.\right.} & (n>m) \end{array}\right. \end{align} $$ の証明をしている。 \varepsilon_{i j k} \varepsilon_{l m n} の計算の参考になる。
  3. ベクトル解析公式の証明 – 準備篇 | 高校物理の備忘録
    証明がたくさん書いてある。丁寧に途中式まである。

まとめ

使いこなすまでの道のりはまだまだ長そう...。

テンソルを使いこなしたい1(アインシュタインの縮約記法)

「ゼロから作るディープラーニング1」で、Affineレイヤの逆誤差伝播のために行列積の行列微分がシレっとでてきた。一旦書き下して納得はしたものの、行列の順序を意識しながら連鎖律を操るのはシンドイ。そこで、Einsteinの縮約記法でテンソルを扱えたらスッキリするのではと思ったので扱い方のメモ。以下がわかれば簡単な微分についても同様のノリで計算できるようになる。気が向いたら書く。

テンソルとは

ベクトルとか行列の多次元への拡張版。まとめてテンソルって呼ぶことにしましょうってノリで良いのでは(知らんけど)。

スカラー:0階テンソル
ベクトル:1階テンソル
行列  :2階テンソル

のような対応がある。

Einsteinの縮約記法

 \sum_iとか行列の次元とか、自明なものはどんどん省略しちゃおうぜってことらしい。

例)

 s    :0階テンソル
 \boldsymbol{u, v}  :1階テンソル
 \boldsymbol{T, A, B}:2階テンソル

とおくと、

$$ \begin{align} {\boldsymbol u}^T {\boldsymbol v} \longrightarrow \delta_{ij} u_{i} v_{j} = u_i v_i \label{tensor0} \end{align} $$

 \delta_{ij}クロネッカーのデルタ。直接対応しているわけではないが、この指標表示が意味するところは

$$ \begin{align} \mathrm{Tr}\hspace{1pt} \left({\boldsymbol u} {\boldsymbol v}^T \right) \end{align} $$

と同じようなものである。

$$ \begin{align} {\boldsymbol T} {\boldsymbol v} \longrightarrow T_{ij} v_j \label{tensor1} \end{align} $$

$$ \begin{align} {\boldsymbol u} {\boldsymbol v}^T &\longrightarrow u_i v_j \label{tensor2.1} \end{align} $$

$$ \begin{align} {\boldsymbol A} {\boldsymbol B} \longrightarrow A_{ij} B_{jk} \label{tensor2.2} \end{align} $$

ここで重要なのは、ベクトルや行列の成分を表す i, jなどの添字には特別な意味がない事で、次元の下限、上限を文脈から読み取ることになる。添字を表すのに無駄に文字を消費せずに済む。

また、計算結果に影響を与えないダミー指標がある。\eqref{tensor2.2}で言えば jがダミー指標で、計算の過程で総和をとるので計算結果に影響を与えない。つまり、計算のためだけに一時だけに導入された使い捨ての添字がダミー指標である。

テンソル積(Einsteinの縮約記法 ver.)

では、直積(クロス積*1、2つのベクトルのテンソル積)をEinsteinの縮約記法に書き直してみる。

$$ \begin{align} {\boldsymbol T} = {\boldsymbol u} \otimes {\boldsymbol v}^T \longrightarrow T_{ij} = {\boldsymbol u}_{i} {\boldsymbol v}_{j} = {\boldsymbol v}_{j} {\boldsymbol u}_{i} \label{outer} \end{align} $$

テンソルの成分を明示的に扱っているので行列積のように積の順序を気にしなくてよくなる。これだけでも大分スッキリする。

テンソル内積

そりゃそうだよねって形だけど念のため。

$$ \begin{align} {\boldsymbol A} \cdot {\boldsymbol B} = \delta_{ik} \delta_{jl} A_{ij} B_{kl} \end{align} $$

縮約記法から元々の記法へ

各成分へ基底をかけて和をとるだけ。
見たほうが早い。

$$ \begin{align} {\boldsymbol T} &= \left( \begin{array}{ccc} T_{11} & \ldots & T_{1n} \\ \vdots & \ddots & \vdots \\ T_{m1} & \ldots & T_{mn} \end{array} \right) \\ &= \sum_{i=1}^{m} \sum_{j=1}^{n} T_{ij} \left({\boldsymbol e}_i \otimes {\boldsymbol e}_j \right) \end{align} $$

なお、 {\boldsymbol e}_i,  {\boldsymbol e}_jは基底を表す。

まとめ

アインシュタイン先生に感謝

*1:内積に対して外積とも呼ばれることがある。ただし直積と外積は3次元では一致するが同義ではないことに注意が必要。また、outer productは直積の意味で使われるので日本語と英語の違いにも注意すること