読者です 読者をやめる 読者になる 読者になる

Gunosyデータ分析ブログ

Gunosyで働くデータエンジニアが知見を共有するブログです。

【Edward】MCMCの数学的基礎からStochastic Gradient Langevin Dynamicsの実装まで

機械学習 ベイズ統計 確率モデリング Edward

こんにちは。初めまして。
データ分析部新入りのmathetake(@mathetake)と申します。

先日個人ブログでこんなエントリを書いた人です:

mathetake.hatenablog.com

そんなこんなでTwitter就活芸人(?)として活動(?)してましたが、これからは真面目に頑張っていこうと思います。



今日はみんな大好きベイズモデリングおいて、事後分布推定に欠かせないアルゴリズム(群)の一つである*1

マルコフ連鎖モンテカルロ法(Markov chain Monte Carlo)

通称MCMCに関するエントリです。より具体的に、

MCMCの意義(§1.)から始め、マルコフ連鎖の数学的な基礎(§2.,3.,4.)、MCMCの代表的なアルゴリズムであるMetropolis-Hastings法(§5.)、その例の1つである*2Langevin Dynamics(§6.)、そして(僕の中で)絶賛大流行中のライブラリEdwardを使ってより発展的(?)なアルゴリズムであるStochastic Gradient Langevin Dynamicsの説明&実装(§7.,8.)していきたいと思います。


今までのデータ分析ブログのエントリと少しテイストが違うかもしれませんが、お楽しみいただけたら幸いです。

細心の注意を払ってはいますが、数学的正しさを保証する記事ではありませんので、詳細が気になる方は§Referenceにある資料を御覧ください。

はじめに

本題に入る前にまず、MCMCの意義について軽く触れておきます。

ベイズ統計におけるモデルパラメータの推定には事後分布からのサンプリングがかかせません。即ち、モデルのパラメータを{\theta \in \Theta},データの集合を{\mathcal{D}}として、ベイズの定理から得られる次のような確率分布

{ p \left( \theta \ |  \ \mathcal{D} \right) = \dfrac{p \left( \mathcal{D} \ | \ \theta  \right)p(\theta) }{ \displaystyle \int_\Theta p \left( \mathcal{D} \ | \ \theta  \right) p(\theta) d\theta } }

からサンプリングする必要がありますが、一般に分母( 正規化定数 )は解析的に求められず、サンプリングは困難です。

そこを、マルコフ連鎖の収束定理等の数学的基礎理論を使って上手くサンプリングする手法がMCMCで、広く応用されています。*3

ベイズ統計の数学的背景については、僕の個人ブログのエントリを御覧ください

mathetake.hatenablog.com

離散値マルコフ連鎖

{ \mathcal{X}:=\{X^0,X^1,X^2,\dots,X^t,\dots \} } を離散な集合{ \{1,2,3,4,..,k\}(k \in \mathbb{N}) } に値を取る、離散な確率変数の集合とします。確率的な値の時間発展だと考え、確率過程と呼ぶことにします。

(Definition 1.) {\mathcal{X} }マルコフ連鎖であるとは、全ての{ \{x^{s}\}_{s=0}^{t+1} \subset \{1,2,\dots,k \}}に対して

{ \mathbb{P} \left(X^{t+1} = x^{t+1} \ |  \ X^t=x^t ,\dots X^0=x^0 \right) = \mathbb{P}\left(X^{t+1} =x  \ | \ X^t=x^t  \right)\ \ \   \cdots (1) }

が成立する事。

感覚的には、次の時刻における値の分布は現在の値のみで決まり、それ以前の値には影響されない確率過程の事です。

数学的便宜上、次の定義も用意しておきます。

(Definition 2.) マルコフ連鎖{\mathcal{X} }斉時的であるとは、

{p_{i,j}^t:= \mathbb{P}\left(X^{t+1} =j \ | \ X^{t} = i \right) }

時間に依存しない事。またこの時 行列{ p_{i,j} := p_{i,j}^t }遷移行列 と呼ぶ。

つまり状態遷移の確率が時間に依らず一定なマルコフ連鎖の事で、状態遷移を表す行列を遷移行列 としています。


ずっと数学の話になってしまいましたが、最後に一つだけMCMCに関わる重要な命題とその系を述べておきます。

(Proposition 3.) {\mathcal{X} }を斉時的なマルコフ連鎖とする。この時 {X^t} の分布は 遷移行列{ p_{i,j} }{X^0}の分布(初期分布)により完全に決定される。

(Proof) {X^0,X^1,\dots} の分布を{\pi^0, \pi^1,\dots,}とし、{\pi^t_j}で確率変数{X^t}が値{j} を取る確率とする。マルコフ連鎖の性質(1)より

{ \begin{eqnarray} \pi^1_j &=& \sum_{i=1}^k \mathbb{P} \left( X^1=j \ | \ X^0=i \right) \mathbb{P} \left(  X^0 =i \right)  \\
&=& \sum_{i=1}^k \pi^0_i p_{i,j} \ \ \  \cdots (2) \end{eqnarray}}

また、

{ \begin{eqnarray} \pi^t_j &=& \sum_{i=1}^k \mathbb{P} \left( X^t=j \ | \ X^{t-1}=i \right) \mathbb{P} \left(  X^{t-1} =i \right)  \\
&=& \sum_{i=1}^k \pi^{t-1}_i p_{i,j}  \ \ \  \cdots (3) \end{eqnarray} }

が成立するので、帰納的に示される □

この命題の系として次が得られます

(Corollary 4.) 初期分布{\pi^{0}} と遷移行列 {T:= (p_{i,j})} を与えることで、式(2),(3)により斉時的マルコフ連鎖が得られる。

マルコフ連鎖収束定理とMCMC

あと少しMCMCの説明まで辿りつきます。もう少々数学にお付き合いください。

§1. で紹介したように、サンプリングしたい確率分布{ \pi =(\pi_i) }が手元(?)にあるとします。

(Definition 5.)遷移行列が{T:= (p_{i,j})} により与えられる、斉時的マルコフ連鎖{ \mathcal{X}}{\pi}不変分布に持つとは

{ \displaystyle \sum_{i=1}^k \pi_i p_{i,j} = \pi_j  \ \ \ ( \ \forall j \in \{ 1,\dots,k\})}

が成立する事。行列の式で書けば

{ \pi T = \pi  }

が成立すること。この時{\pi}{ \mathcal{X}}不変分布と言う。

さて、MCMCの肝となる定理は次のものです

(Theorem 6.(離散値マルコフ連鎖の収束定理))
{\pi } を不変分布に持つ斉時的マルコフ連鎖{ \mathcal{X}}(i) 非周期的 かつ (ii) 既約*4 である時、マルコフ連鎖は不変分布 {\pi } に収束する。即ち


{ \displaystyle \lim_{t \rightarrow \infty} \mathbb{P}\left( X^t = i \right)  =  \pi_i  }

が成立する。

ここまで来てやっと、MCMCの定義を与えることができます:

(Definition 7.) MCMC(Markov chain Monte Carlo)とは、サンプリングしたい確率分布{\pi} を不変分布とするような既約で非周期的なマルコフ連鎖を構築&サンプリング するアルゴリズムの事。

Theorem 6.により、MCMCにより生成されるサンプルの列{ \{x^t\}^\infty_{t=0} } は確率分布 {\pi}からのサンプルに収束し、目的を達成することができます。

既約性と非周期性を満たすようなマルコフ連鎖を構築するのはそんなに難しくはありません、が、サンプリングしたい確率分布{\pi} を不変分布とするようなマルコフ連鎖を構築するのは一般に困難です。

そこでよく用いられるのが*5詳細釣り合い条件と呼ばれる不変分布を持つための十分条件です:

(Definition 8.)遷移行列が{T:= (p_{i,j})} により与えられる、斉時的マルコフ連鎖{ \mathcal{X}}不変分布 {\pi} に対して詳細釣り合い条件を満たすとは


{ \displaystyle  \pi_i \  p_{i,j} = \pi_j \ p_{j,i}  \ \ \ \ (  \forall \ i, j \in \{ 1,\dots,k\})}


が成立する事。またこの時、マルコフ連鎖{ \mathcal{X}}{\pi} を不変分布に持つ。

連続値マルコフ連鎖の場合

今まで簡単のため、離散な値を持つマルコフ連鎖について話をしてきましたが、連続な確率変数の話に一般化する事ができます。

連続値( 便宜上{\mathbb{R}^d}値とする )なマルコフ連鎖{\mathcal{X} = \{ X^0,X^1,\dots,X^t,\dots \}}斉時的である時、離散値の場合の推移行列に対応する推移核{ T: \mathbb{R}^d \times \mathbb{R}^d \rightarrow \mathbb{R} }

{  \mathbb{P} \left( X^{t} \in A \  | \ X^{t-1} =x \right) = \displaystyle \int_A T(x,y) dy \ \ \ \ \ \ (A \subset \mathbb{R}^d ) }

を満たすものとして定めます。感覚的には現時刻で値 {x} を取る時、次の時刻の分布を表す密度関数です。

また、{ \mathcal{X}}が分布{\pi}を不変分布に持つとは

{ \displaystyle \int_{\mathbb{R}^d} \pi(x)T(x,y)dx = \pi(y) }

を満たすことであり、詳細釣り合い条件

{ \pi(x) \  T(x,y) = \pi(y) \ T(y,x)   \ \ \ \cdots (4)}

で与えられます。

(Theorem 9.(連続値マルコフ連鎖の収束定理))
{\pi } を不変分布に持つ斉時的マルコフ連鎖{ \mathcal{X}}(i) 非周期的 かつ (ii) 既約である時、{ \mathcal{X}}total variation distanceの意味で{\pi} に収束する。即ち、


{ \displaystyle \lim_{t \rightarrow \infty} \sup_{A \subset \mathbb{R}^d } \left| \pi^t(A) - \pi(A) \right| \ = \ 0 }


が成立する。

Metropolis-Hastings法

MCMCの代表的なアルゴリズム(群)である、Metropolis-Hastings法(以下M-H法)について説明します。

M-H法では、各 {x}に対して提案分布と呼ばれる確率分布 {q(x,y)dy}を用意し、採択確率と呼ばれる確率

{\alpha(x,y):= \min \left\{ 1, \dfrac{\pi(y)q(y,x) }{\pi(x)q(x,y)} \right\} }

を準備します。そして推移核 {T(x,y)}

{T(x,y) = \alpha(x,y)q(x,y) + A(x)\delta_{x,y} }

として定義&適当な初期分布{\pi^0 }を与えることで斉時的マルコフ連鎖 {\mathcal{X}_{MH} }を考えます。(ここで {\delta_{x,y} } はクロネッカーのデルタ関数で、{A(x)} は正規化定数を与える関数。)

(Theorem 10.) 斉時的マルコフ連鎖 {\mathcal{X}_{MH} }{\pi} を定常分布に持つ

(Proof) 定義から

{ \begin{eqnarray} 
\pi(x)q(x,y)\alpha(x,y) &=& \pi(x)q(x,y) \min \left\{ 1, \dfrac{\pi(y)q(y,x) }{\pi(x)q(x,y)} \right\}\\
&=& \min \left\{ \pi(x)q(x,y), \pi(y)q(y,x) \right\}\\
&=& \pi(y)q(y,x) \min \left\{ \dfrac{\pi(x)q(x,y)}{\pi(y)q(y,x)}, 1 \right\}\\
&=& \pi(y)q(y,x)\alpha(y,x)
\end{eqnarray}}

が従うので、クロネッカーのデルタの定義から(4)式が成立し詳細釣り合い条件が満たされる。□

Theorem 10.だけでは保証されない既約性や非周期性を満たすような推移核の具体的な設計は重要な課題ですが、ここではそのような性質が満たされ Theorem 6.(or 7.) が成立すると仮定しましょう。*6

この時マルコフ連鎖{\mathcal{X}_{MH}} からサンプリングする次のようなアルゴリズムを Metropolis-Hastings法と呼びます。

(Metropolis-Hastings法)
(1) 初期分布{ \pi^0} から {x^0} をサンプリングする:

{x^0 \sim \pi^0 }

(2) {t =0,1,2,3, \dots} に対して、以下を実行する
 (i) 標準一様分布から乱数 {u} を生成する:

{u \sim U(0,1) }
 (ii) {y} を提案分布 {q(x^t, y)} からサンプリングする:
{y \sim q(x^t, y) }
 (ii) 次の式により"次の点" {x^{t+1} }を決める:
{x^{t+1} \ := \ \displaystyle
  \begin{cases}
    y  \ \ \ \ \ if \ \  u < \alpha(x^t, y) \\
    x^t \ \ \ \ otherwise
  \end{cases}
}

このアルゴリズムが、実際に上で与えた{\mathcal{X}_{MH} } からサンプリングしている事は明らかでしょう。

§0. で述べたように採択確率の計算には目標の分布{\pi(x))}正規化定数が必要ない事に注目して下さい。

Langevin Dynamics

M-H法を実際に実行するためには、提案分布 {q(x,y)dy} を定義する必要があります。

ここではその例として、Langevin Dynamics法(Metropolis-adjusted Langevin Algorithm)(以下LD法) を紹介します。

LD法では提案分布を次のように定義します:

{q(x,y) := \mathcal{N} \left( x + \dfrac{\epsilon}{2} \dfrac{\partial \log \pi }{\partial x}(x),  \ \ \epsilon I \right)  }

ここで{\mathcal{N}\left( \mu, \Sigma \right)} は平均 {\mu}, 分散 {\Sigma}に従う正規分布で、{\epsilon}はstep size(またはlearning rate)と呼ばれるハイパーパラメータです。

ベイズモデルの事後分布に適応する場合において、決定論的に眺めると、LD法は事前分布で正規化した、ノイズ入り勾配降下法のようなものであると解釈することができます。

実際、事後分布推定において、{\pi} は観測データを {\mathcal{D = \{d_i \}}_{i=1}^N  }として

{ \pi(\theta) = p \left( \theta \ |  \ \mathcal{D} \right) = \dfrac{p(\theta) \prod_{i=1}^N p\left( d_i \ | \ \theta  \right) }{ \displaystyle \int_\Theta p \left( \mathcal{D} \ | \ \theta  \right) p(\theta) d\theta } }

で与えられ、その対数微分は

{\displaystyle \dfrac{\partial}{\partial \theta} \log \pi \left( \theta \right) =\dfrac{\partial}{\partial \theta} \log p \left( \theta \right)  + \dfrac{\partial}{\partial \theta} \sum_{i=1}^N \log p\left(d_i \ | \ \theta \right)  \ \ \ \cdots (5)}

のように計算できるので、上のような解釈ができます。

実はLG法は、Stanで有名になったHamiltonian Monte Carlo法の特別なケースと等価になっているので、気になる方は [1]や[2] を御覧ください。

Stochastic Gradient Langevin Dynamics(SGLD)

ここまでMCMCの数学的基礎からM-H法、そしてその具体例としてLG法を紹介しました。

LG法の問題点として、(5) 式の計算量がサンプル数が増えるほど膨大になっていく点があります。

近年はビッグデータと呼ばれるバズワードもあるように、サンプル数&パラメータ数が巨大なセッティングでモデリングする事が多いのでこのままLG法を適用する事はできません。*7


その問題点を克服するサンプリング手法として、ここで紹介するのが Stochastic Gradient Langevin Dynamics法 [8](以下SGLD法)です。

SGLD法では次のように初期分布からサンプルしたパラメータを更新していきます:

まず、次の2つの条件

{\displaystyle \sum_{t=0}^\infty \epsilon_t^2  < \infty \ , \ \ \ \ \ \ \sum_{t=0}^\infty \epsilon_t  = \infty}

を満たす数列{ \{\epsilon_t \}_{t=0}^\infty }を用意し、パラメータサンプル{ \{ \theta_t \}_{t=0}^\infty }を次の式によって取得して行きます

{\begin{eqnarray}
\displaystyle & \theta_{t+1} &\leftarrow  \ \ \ \theta_t  +  \dfrac{\epsilon_t}{2} \dfrac{\partial L}{\partial \theta} (\theta_t) + \eta_t  \ , \ \ \ \eta_t  \sim \mathcal{N}\left( 0 \ ,\epsilon_t \right) \\ \ \\
&\dfrac{\partial L}{\partial \theta} (\theta_t)& =  \dfrac{\partial \log p}{\partial \theta_t}  \left( \theta_t \right) +  \dfrac{N}{\left| \mathcal{S_t} \right|}\sum_{d \in \mathcal{S_t} } \dfrac{\partial \log p \left(d \ | \ \theta \right) }{\partial \theta}(\theta_t)
\end{eqnarray} }

ここで、{\mathcal{S_t}} はデータ{ \{d_i\}_{i=1}^N } からランダムに抽出された {N} より十分小さいミニバッチとします。

注意として、このアルゴリズムに対応するマルコフ連鎖は斉時的ではないので、上述の収束定理は適用できません。

ですが、例えば [9]で収束性に関する解析がされています。

EdwardでのSGLDの実装

最後に確率モデリング用ライブラリEdwardを用いて、SGLDをベイズ的線形回帰に適用してみようと思います。

Edwardの詳しい使い方は公式チュートリアルまたは次の論文

[1701.03757] Deep Probabilistic Programming

[1610.09787] Edward: A library for probabilistic modeling, inference, and criticism

をご覧ください。また質問等ありましたら@mathetakeまで気軽にリプライorDMください。

まず各種ライブラリをimportします。

import numpy as np
import tensorflow as tf
import edward as ed
from edward.models import Normal, Empirical
import time

次にデータセットを用意します。

N = 20000  # サンプル数
D = 50  # 特徴量の次元
N_ITER = 10000  # MCMCのiteration
MINI_BATCH_SIZE = 2500  #ミニバッチのサイズ 

# toy dataset. 切片=0はなし.
def build_toy_dataset(N, D, noise_std=0.1):
    w = np.random.randn(D).astype(np.float32)
    X = np.random.randn(N, D).astype(np.float32)
    Y = np.dot(X, w) + np.random.normal(0, noise_std, size=N)
    return w, X, Y

# データ生成。観測値のノイズの分散は既知とする。
w_true, X_data, Y_data = build_toy_dataset(N, D)

# ミニバッチを返す関数 
def next_batch(mini_batch_size=128):
    indexes = np.random.randint(N, size=mini_batch_size)
    return X_data[indexes], Y_data[indexes]

モデルを構築し, 推論のためのインスタンスを作ります。

# 観測データを挿入するためのデータを収めるplaceholder
x = tf.placeholder(tf.float32, [MINI_BATCH_SIZE, D])
y_ph = tf.placeholder(tf.float32, [MINI_BATCH_SIZE])

w = Normal(mu=tf.zeros(D), sigma=tf.ones(D))
b = Normal(mu=tf.zeros(1), sigma=tf.ones(1))
y = Normal(mu=ed.dot(x, w) + b, sigma=tf.ones(MINI_BATCH_SIZE)*0.1)

# 経験分布をposteriorの近似に使う
qw = Empirical(params=tf.Variable(tf.random_normal([N_ITER, D])))
qb = Empirical(params=tf.Variable(tf.random_normal([N_ITER, 1])))

# SGLD法用インスタンス
SGLD = ed.SGLD(latent_vars={w: qw, b: qb}, data={y: y_ph})


最後に推論を実行します。

# 推論GO 
# data辞書にはobservedな確率変数の観測データを送る。
# xの値は確率変数ではないので、updateの際feed_dictで送る。
SGLD = ed.SGLD(latent_vars={w: qw, b: qb}, data={y: y_ph})
SGLD.initialize(scale={y: float(N) / MINI_BATCH_SIZE}, step_size=0.00001, n_iter=N_ITER)

start = time.time()
init = tf.global_variables_initializer()
init.run()
for _ in tqdm(range(N_ITER)):
    X_batch, Y_batch = next_batch(MINI_BATCH_SIZE)
    _ = SGLD.update(feed_dict={x: X_batch, y_ph: Y_batch})
elapsed_time = time.time() - start
print("elapsed_time:{}".format(elapsed_time))


実行結果ですが、ミニバッチで勾配計算をしない通常のLD法と比較してみました:

アルゴリズム 実行時間(s) MSE of W Estimated b
SGLD法 26.4 {4.9 \times 10^{-5}} 0.0020
LD法 55.4 {5.3 \times 10^{-7}} 0.0001

といった感じです。精度はLD法には劣りますが、実行時間は約2倍短いと言った感じです。

より大きなデータセット&複雑なモデルの場合、この差はより顕著になるでしょう。

Appendix

(Definition) 斉時的マルコフ連鎖{ \mathcal{X}}非周期的であるとは、任意の{i \in \{1,\dots,k\} } に対して

{gcd \left\{ n \ | \ \mathbb{P} \left( X^n=i \  |  \ X^0 =i \right)  > 0 \right\} = 1 }

が成立する事。

(Definition) 斉時的マルコフ連鎖{ \mathcal{X}}既約であるとは, 任意の{i, j \in \{1,\dots,k \}} に対してある{t_{i,j} \in \mathbb{N} } が存在して

{ \mathbb{P}\left(X^{t_{i,j}}=j \ | \ X^0 =i  \right) >0   }

が成立する事。

*1:他には変分ベイズ法などがあります

*2:Hamiltonian Monte Carlo法の一種でもあります。

*3: 事後分布推定に限らず、”正規化定数が分からない分布からサンプリングする手法” として広く使われています。

*4: (i) 非周期的 かつ (ii) 既約 の定義はAppendixに付けておきます。

*5:詳細釣り合い満たさないようなMCMCの研究が最近流行っている、らしいです。例えば 詳細つりあいを満たさないマルコフ連鎖モンテカルロ法とその一般化を御覧ください。

*6:M-Hについて既約性や非周期性が満たされるための条件は論文[5]をご覧ください。

*7:LG法にかぎらず、M-H法は採択確率の計算が基本的にintractableです。ですが、例えばハミルトニアンモンテカルロ法をビッグデータに適用する手法については次の論文があります : Stochastic Gradient Hamiltonian Monte Carlo

世界を代表する8人の旬なトップ機械学習研究者たち (2017年上半期版)

論文 機械学習 学会

データ分析部の久保です。 最近行ったライブはAimerのAcoustic Live Tour 2017です。

早いもので2017年も3月になりましたが、機械学習分野は相変わらずとてもホットな分野です。 去年はAI、人工知能という言葉がディープラーニングとともにバズワードになり、その傾向は尚も続いています。

その流行の元となったのが機械学習なわけですが、今その最先端ではどういう人がどのような研究をしているのかをかなりざっくりと見ていきたいと思います。

調査方法は2013年に同様のことを行ったとき

qiita.com

と同じく、NIPSとICMLという機械学習の代表的国際会議の過去3年分を対象とし、1st authorの重要度をそれ以外の著者よりも重くしてスコアづけしました。具体的には複数人の著者がいる場合は1st authorを0.8として、残りの0.2を他の著者に分配、1人の場合は1としてあります。また1st authorと2nd authorの貢献度が同じと明記されている論文もありましたがそのあたりは考慮できていません。

結果5.0ポイント以上になった上位8人を今回は紹介したいと思います。 彼らの1st authorの論文も併記しておきます。

Prateek Jain (7.2ポイント)

https://www.microsoft.com/en-us/research/people/prajain/

  • マイクロソフトリサーチ(バンガロール, インド)の研究者
  • 専門は大規模非凸最適化および統計的学習理論
  • プライバシーやコンピュータビジョン、テキストマイニング、自然言語処理などへの応用も積極的
  • 主なNIPS, ICMLの論文
    • Alternating Minimization for Regression Problems with Vector-valued Outputs (NIPS 2015)
    • Predtron: A Family of Online Algorithms for General Prediction Problems (NIPS 2015)
    • Optimizing Non-decomposable Performance Measures: A Tale of Two Classes (ICML 2015)
    • Surrogate Functions for Maximizing Precision at the Top (ICML 2015) など

Ohad Shamir (6.8ポイント)

http://www.wisdom.weizmann.ac.il/~shamiro/

  • ワイツマン科学研究所(イスラエル)の研究者
  • 実用的な効率性と理論的な知見を組み合わせることを得意とする
  • 主なNIPS, ICMLの論文
    • Without-Replacement Sampling for Stochastic Gradient Methods: Convergence Results and Application to Distributed Optimization (NIPS 2016)
    • Convergence of Stochastic Gradient Descent for PCA (ICML 2016)
    • Fast Stochastic Algorithms for SVD and PCA: Convergence Properties and Convexity (ICML 2016) など

Cho-Jui Hsieh (6.0ポイント)

http://www.stat.ucdavis.edu/~chohsieh/rf/

  • カルフォルニア大学デービス校のアシスタント・プロフェッサー
  • 専門は大規模データに対する最適化
  • 主なNIPS, ICMLの論文
    • PASSCoDe: Parallel ASynchronous Stochastic dual Co-ordinate Descent (ICML 2015)
    • Fast Prediction for Large-Scale Kernel Machines (NIPS 2014)
    • QUIC & DIRTY: A Quadratic Approximation Approach for Dirty Statistical Models (NIPS 2014) など

Jacob Steinhardt (5.6ポイント)

http://cs.stanford.edu/~jsteinhardt/

  • スタンフォード大学の博士課程3年生(!)
  • 信頼できて人にとっても使いやすい機械学習とは何かを研究している
  • 安全なAIに関する長期的、また短期的なチャレンジについてのエッセイを書いていたりしている
  • 主なNIPS, ICMLの論文
    • Avoiding Imposters and Delinquents: Adversarial Crowdsourcing and Peer Prediction (NIPS 2016)
    • Unsupervised Risk Estimation Using Only Conditional Independence Structure (NIPS 2016)
    • Learning Fast-Mixing Models for Structured Prediction (ICML 2015) など

Jie Wang (5.0ポイント)

http://www-personal.umich.edu/~jwangumi/

  • ミシガン大学のアシスタント・プロフェッサー
  • 最適化全般、またバイオインフォマティクスや画像処理に対しても積極的
  • 主なNIPS, ICMLの論文
    • Multi-Layer Feature Reduction for Tree Structured Group Lasso via Hierarchical Projection (NIPS 2015)
    • Safe Screening for Multi-Task Feature Learning with Multiple Data Matrices (ICML 2015)
    • Two-Layer Feature Reduction for Sparse-Group Lasso via Decomposition of Convex Sets (NIPS 2014) など

Yining Wang (5.0ポイント)

http://yining-wang.com/

  • カーネギーメロン大学の博士課程3年生(!)
  • 行列補完や近似、部分空間クラスタリングに強み
  • 主なNIPS, ICMLの論文
    • Data Poisoning Attacks on Factorization-Based Collaborative Filtering (NIPS 2016)
    • Online and Differentially Private Tensor Decomposition (NIPS 2016)
    • A Theoretical Analysis of Noisy Sparse Subspace Clustering on Dimensionality-Reduced Data (ICML 2015) など

Elad Hazan (5.0ポイント)

https://www.cs.princeton.edu/~ehazan/

  • プリンストン大学の教授
  • 数理最適化やゲーム理論、計算複雑生など機械学習の中心的なところが研究分野
  • ちなみに2013年に同様の調査をしたときの唯一の連続ランクインで、当時はイスラエル工科大学の准教授であった
  • 主なNIPS, ICMLの論文
    • On Graduated Optimization for Stochastic Non-Convex Problems (ICML 2016)
    • Variance-Reduced and Projection-Free Stochastic Optimization (ICML 2016)
    • On Graduated Optimization for Stochastic Non-Convex Problems (ICML 2016) など

Kirthevasan Kandasamy (5.0ポイント)

http://www.cs.cmu.edu/~kkandasa/

  • カーネギーメロン大学の博士課程2年生(!!)
  • 統計とアルゴリズムの共通部分に興味があり、研究分野はバンディット問題、ベイズ最適化など
  • 主なNIPS, ICMLの論文
    • Gaussian Process Bandit Optimisation with Multi-fidelity Evaluations (NIPS 2016)
    • The Multi-fidelity Multi-armed Bandit (NIPS 2016)
    • Additive Approximations in High Dimensional Nonparametric Regression via the SALSA (ICML 2016) など

おまけ

今回のランキング方法でトップになる日本人(ぽい名前の)研究者は誰だろうと思いみてみると、東工大准教授の鈴木大慈先生(3.0ポイント)でした。

Taiji Suzuki's Homepage (鈴木大慈)

  • 主なNIPS, ICMLの論文
    • Minimax Optimal Alternating Minimization for Kernel Nonparametric Tensor Learning (NIPS 2016)
    • Convergence rate of Bayesian tensor estimator and its minimax optimality (ICML 2015)
    • Stochastic Dual Coordinate Ascent with Alternating Direction Method of Multipliers (ICML 2014)

また今回対象としたNIPSとICML以外にもAISTATS, WWW, WSDM, KDD, ICLR, IJCAI, COLTなど、さまざまな国際会議もあるので気になる人はチェックしてみてください。最近だとarxivにいきなり投稿されることも増えてきました。

所感

今回調査する前はやはり深層学習まわりの研究者が上位にくるのかなと思っていましたが、あまりそんなことはなく昔から着実に進んでいる研究分野の研究者が多かったという印象でした。

同じ年のNIPSやICMLに同時に1st authorとして複数論文を通している研究者もちらほらいるんですが、いったいどんな研究スタイルでそんなことできるのか個人的には気になるところです。。

機械学習はまだこれから発展していく研究分野であることは間違いないと思うので、これからも引き続き注目していきます。

Spark StreamingからAmazon Kinesis Analyticsへ移行する話

分析基盤 Kinesis AWS Spark

はじめに

こんにちは、データ分析部の森本です。主な業務は記事配信アルゴリズムの改善とログ基盤の整備です。
Gunosyでは、ユーザーへより良い記事を提供するためにアクセスログをストリーム処理し、集計結果を記事配信アルゴリズムに活用しています。 ストリームログ基盤にはSpark Streamingを利用していますが、現在Kinesis Analyticsへ移行中です。
この記事ではKinesis Analyticsへ移行する理由や運用上のTips等についてお話します。

Spark Streamingを利用したストリームログ基盤構成

f:id:moyomot:20170213185600p:plain

現在のストリームログ基盤はSpark Streamingで集計を行い、結果をRDSに保存しています。

なぜSpark StreamingからKinesis Analyticsへ移行するのか

サーバーコストと運用コストの削減を目的としています。

サーバーコストについて

Spark StreamingはEMR上で使用していますが、コストはそこそこします。 例えばサーバー(m4.xlarge)10台と考えると1ヶ月およそ約30万円になります。
($0.348 + $0.060) * ¥110 * 24hour * 30day * 10台
スポット価格で運用しサーバーコストを抑えてますが、スポット価格の変動に注意する必要があります。 これを複数クラスタ運用していたので、それなりのサーバーコストを要しました。 またクラスタの構築コスト、運用コストを考慮するとEMRが最適なため、自前でクラスタ環境の構築は行いませんでした。

運用コストについて

長期間運用していると、ノードがしばしば落ちます。
一時的に少数のノードが落ちてもクラスタ運用に問題ありませんが、長期的にはクラスタ状況をモニタリングし、問題が発生する前にノードの追加やクラスタの再構築が必要です。主観的に休日にクラスタ整備を行うことが多く、辛いです。(本当に土日よく落ちるのかノードダウンの統計情報を可視化してみたい)

Spark Streamingはこのアクセスログ集計でオーバースペックなところもありました。
そこで、新たな基盤となるKinesis Analyticsを利用したストリームログ基盤を紹介します。

Kinesis Analyticsを利用したストリームログ基盤構成

f:id:moyomot:20170213185624p:plain ログ集計はKinesis Analyticsが担当し、Kinesis FirefoseがElasticsearchに結果を保存しています。

  • Kinesis Analyticsとは
    ストリーミングデータをSQLクエリで集計できるサービスです。
    標準SQLを使用でき、容易に実装できました。
    集計はKinesis AnalyticsのSQLのみで完結するため、Sparkアプリケーションを実装するよりも開発コストを抑えることができました。 今回Kinesis AnalyticsのインプットにはKinesis Stream、アウトプットにはKinesis Firefoseを使用しています。

詳細はこちらの記事が詳しいです。

data.gunosy.io

運用してみて

現在Kinesis Analyticsを利用したログ基盤は運用テスト中で、大きな問題もなく順調に稼働しています。ここでは運用上のTipsについて紹介します。

マスターデータとJOINできる

Kinesis Analyticsは、SQLクエリ上で静的ファイル(リファレンスデータ)とJOINできます。
マスターデータ的なファイルをS3に保存しておき、必要なAPIを呼ぶことでKinesis Analytics上でJOINできます。 これは非常に便利な機能で重宝しています。
アクセスログ集計アプリケーションではマスターデータの更新を1日に1度Data Pipelineで行い、S3に再アップロードしています。 リファレンスデータの更新はLambdaで行っています。

AddApplicationReferenceDataSource - Amazon Kinesis Analytics

監視はDatadogとAWS Lambdaで実施

Gunosyの監視はCloudWatch、Datadog、Papertrail、PagerDutyなどの各種サービスを必要に応じて利用しています。 ログストリーム基盤の監視では最終的なアウトプット先であるElasticsearch Serviceに重きをおいて監視することにしました。 Elasticsearchの各種メトリックをDatadogでモニタリングし、アラート先をSlackに設定しています。 また直近N分間に登録されているドキュメント数もモニタリングしたいため、ドキュメント数をカウントし、Cloudwatchに登録するオリジナルスクリプトをLambdaで実装しています。(なんでも雑にLambdaで監視し、Slackに通知したくなってしまいます)

Kinesis Analyticsの対応しているリージョンが限定的

2017年2月現在、Kinesis Analyticsの対応しているリージョンは、アイルランド、オレゴン、バージニアのみとなっています。 調査したところ東京からレイテンシが最も小さいリージョンはオレゴンであったため、現在はオレゴンにKinesis Analytics環境を構築しています。 この場合、東京リージョンのログストリーム環境と比べて、オレゴンまでの転送分、遅延が多少発生します。(数秒レベル) 今回の用途では、そこまで遅延に対してシビアではないため利用に踏み切りましたが、利用用途に応じて検証する必要があります。

その他の所感

  • ElasticsearchのクエリはSQLと比較すると難しいが、中央値算出などいろいろな関数がサポートされているのが良い(なんでMySQLでは中央値を簡単に算出できないの…)
  • Kinesis streamのシャード数見積もり結構難しい、多めに見積もったほうが無難
  • Kinesis streamのシャード数分割も画面でできるようになって欲しい

Kinesis Analyticsのコスト

実際にそれなりのコストを抑えることに成功しましたが、具体的にどのくらい安くなったのかお伝えできないことが残念です。(すいません)
Kinesisの見積もり方法は下記リンクを参考にしてください。 また、オレゴンリージョンで環境を構築しているため、リージョンを越えるデータ転送料も発生しています。

料金 - Amazon Kinesis ストリーム | AWS

Amazon Kinesis Analytics 料金表 – アマゾン ウェブ サービス (AWS)

料金 - Amazon Kinesis Firehose | AWS

まとめ

  • コストの観点からアクセスログ集計はSparkStreamingからKinesis Analyticsに移行します
  • Kinesis Analyticsはフルマネージドで運用コストを抑えられます
  • Kinesis AnalyticsはSpark Streamingのような自由度はないため要件に応じた利用が求められます

ABテストの対象をいい感じに割り振る方法

Python 分析基盤 分析ノウハウ ABテスト

こんにちは、データ分析部の石塚 (@ij_spitz) です。 最近聴いている曲は久保田利伸さんのLA・LA・LA LOVE SONGです。 ロンバケ最高でした、月曜9時はOLが街から消えるというのも納得です。

Gunosyではプロダクト改善のためにABテストを用いて意思決定を行っています。 今回はタイトルにもある通り、ABテストを実現させる上で必要となる対象の割り振り方法を、Gunosyで以前使っていた従来の手法と半年ほど前に新しく導入した手法の2つをご紹介します。 いい感じってなんだよと思われるかもしれませんが、従来の手法の課題を解決するようにいい感じに割り振る方法と理解していただければと思います。 それぞれの運用上で気づいたメリット・デメリットなども合わせてご紹介します。

従来の手法

以前はユーザIDを100で割った余りを使用していました。 例えば、全ユーザの1%でテストしたいという時には、100で割った余り0を対象にしてテストを実施し、100で割った余り1と比較するという仕組みです。

LPなどのABテストでは1回切りの訪問のため、ユーザが複数回訪問してきたときに同じテストが割り当たることは考慮しなくてもいいですが、グノシーは毎日のようにユーザが継続して利用するニュースアプリなので、テスト期間中は同じユーザが同じテストに割り当たることが必須条件となります。 次にこの手法のメリットとデメリットを見ていきます。

メリット

  • 集計がラク
    • 標準SQLには MOD という剰余を計算する関数があるので、集計する際にはログが格納されているのがBigQueryでもRedshiftでも簡単に行うことができます。
    • そのためログにユーザが割り当たっているABテストの情報を追加せずとも集計することができます。
  • 計算コストが気にならない
  • 管理しやすい
    • ユーザIDさえわかれば、どのユーザがどのABテストに割り当たっているかをすぐ確認できるので、管理しやすいです。
    • Gunosyではスプレッドシートを使ってテストを管理していました。

デメリット

  • テスト対象を選ぶのが面倒
    • この手法ではテスト対象と比較対象の2つを自分で選ぶ必要があります。
    • テストの数が2個や3個だと対象を選ぶのも苦ではないのですが、10個や20個となってくると1ユーザに複数のテストを割り当てることになり、各テストが均等に割り当てられているかの担保が難しく、時間の掛かる作業となってしまいます。
  • 以前のテストが影響してしまう
    • 前に実施していたABテストがバイアスとなって効いてくる場合があります。
    • あるテストを実施していた100で割った余りのユーザだけ数が少なくなっていたり、アクティビティの高いユーザ群になってしまい、誤った意思決定につながる可能性があります。

新しく導入した手法

現在使用している手法も余りを用いていますが、ハッシュ関数を利用したものになっています。 言葉で説明しても理解しづらいので、まずPythonのコードを載せます(Pythonのバージョンは3系を使用しています)。

import hashlib
ratio = 5
key = 'ab_test_key'
user_id = 4192481
mod = int(hashlib.sha1((key + str(user_id)).encode('utf-8')).hexdigest()[:15], 16) % 100
if mod < ratio:
    print('テスト対象')
elif mod >= 100 - ratio:
    print('比較対象')
else:
    print('どちらでもない')
  • 割り当てアルゴリズムは以下のようになっています。
    • ABテストごとに固有な文字列とユーザIDを連結
    • 連結した文字列をハッシュ化して16進数の文字列に変換
    • ハッシュ化した16進数の文字列の先頭15文字を取得する
      • 先頭15文字を取得しているのはRedshiftで16進数を10進数に変換するときに、文字列が長すぎて桁あふれが起きてしまったため
    • 取得した15文字を10進数に直す
    • 最後に取得した数値の100で割った余りがテスト割合より小さければテスト対象として割り当てる
    • 100からテスト割合を引いた値以上ならば比較対象として割り当てる

メリット

  • テスト対象を選ぶのがラク
    • 対象は自動で選ばれるのでかなりラクになりました。
    • ただ1つ注意するポイントとしては、他のテストがテスト対象と比較対象に同じ割合で含まれている点です。同じ割合で含まれていれば比較する際には問題ないという仮定を置いているのですが、どうしてもこの2つのテストは被らせたくないという条件があれば、ABテストごとに設定する固有な文字列を同じにすると被らなくなります。
  • 以前のテストが影響がなくなる
    • この手法を用いると対象となるユーザはテストごとに変わるので、ずっと同じ対象でテストをし続けるといったことはなくなります。

デメリット

  • 管理が複雑になる
    • 以前のようにスプレッドシートだけで管理することはできなくなりましたが、テスト対象と比較対象をデータベースで管理しています。
    • また、ABテストの拡大・停止のためのツールを作成しているので、管理コストはかなり抑えられていると思います。
  • 集計が複雑になる
    • 以前のようにSQLだけで集計を完結させることは難しくなりましたが、割り当たっているABテストの情報をログに付与したり、テーブルで管理することによって集計も手軽に実施することができます。
    • ヒューマンエラーや手間を考えると集計は自動化することをおすすめします。

まとめ

今回新しく導入した手法で、テスト対象を選ぶ手間だったり、以前のテストが影響してしまうといった従来の手法による課題は解決することができました。 ただし、割り当てがブラックボックス化してしまい、本当に上手く行っているのかがわかりづらくなっているという欠点もあります。 そういったところをきちんと検証したいという場合には、A/AテストやGroup Validationといった仕組みを用いると誤った意思決定を減らすことができるので、興味があれば下記のリンクを参照してみてください。

さくっとトレンド抽出: Pythonのstatsmodelsで時系列分析入門

Python 分析ノウハウ 時系列分析

久しぶりの投稿になってしまいましたが、ニュースパス(現在CM放映中!!)開発部の大曽根です。 作業中はGrover Washington Jr のWinelightを聴くと元気が出ます。参加ミュージシャンが素晴らしいですね。

なぜ時系列分析をするのか

数値を非常に重視している弊社では、数値を知るためのツールとしてRedashやChartioおよびSlackへの通知を活用しています。現在の数値を理解する上では、長期のトレンド(指標が下がっているのか、上がっているのか)を知ることが重要です。しかし、日々変化するデータ(特に売上やKPIと言われる指標)は、ばらつきも大きく、変化を適切に捉えることが難しいこともあります。

特にSlackなどへの通知を行っていると、日々の変化に囚われがちです。例えば、弊社ではニュースを扱っていますが、サッカーや競馬のニュースであれば土日に、技術系のニュースであれば平日に偏る傾向があります。日々の数値の浮き沈みに注目していると全体的なトレンドの変化を見逃す可能性もあります。そこで、時系列分析を行うことで変化を適切に捉えたり、予測をある程度正確にすることができます。 今回は周期的なデータとトレンド成分を分割して表現できる季節調整モデルについて紹介します (1週間であれば移動平均などでもある程度のノイズ除去は可能)。

季節調整

季節調整データをざっくり説明すると、時系列のデータを 観測値 = トレンド成分 + 季節成分 + ノイズ成分 で説明するモデルです。

f:id:dr_paradi:20170201104931p:plain

※詳しい説明は北川源四郎 『時系列解析入門』の12章を参考にいただけると幸いです。

弊社の提供するアプリもですが、人間の生活に密着している以上は、人間の生活リズムに影響を受けます。大きく分けて「月要因」、「曜日要因」、「時間要因」がありますが、今回は曜日に注目してサンプルを実装しようかと思います。

実演

本来ならサービスのデータを使いたいのですが、当然使えないので、東京電力さんのデータを使って実装したいと思います。 まずいつものpandas+jupyter+matplotlibで生データを可視化してみます。

data.gunosy.io

import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline

# csvファイルの読み込み
row_data = pd.read_csv('tokyo2016.csv', index_col=0)

df = row_data.reset_index()
df['date'] = pd.to_datetime(df['date'], format='%Y/%m/%d')
df['power'] = df['power'].values.astype(int)
daily_data = df.groupby('date')['power'].sum().reset_index() #元データが1時間ごとなので日毎のデータにする
daily_data = daily_data.set_index(['date'])
daily_data.plot()

すると下記のようなデータが得られます。1年分のデータを可視化する人間の目でもなんとなく傾向は掴めますが、細かい範囲ではデータが波打っており正確なトレンドを把握することは難しいかと思います。

f:id:dr_paradi:20170201013929p:plain

そこで季節調整モデルの出番です。今回は pythonのstatsmodelsというライブラリを用いてデータをトレンド成分 + 季節成分に分けていきます。 StatsModels: Statistics in Python — statsmodels documentation

追記: Pythonのversionは3.5.1、statsmodelsのバージョンは0.8.0rc1を使用しています。

ソースは下記の通りです。

import statsmodels.api as sm
res = sm.tsa.seasonal_decompose(daily_data.values, freq=7)
res.plot()

freq には周期を入力します。今回はdailyのデータということで1週間を周期として設定します。 出力は以下のようになります。グラフの上から 観測データ(生データ)トレンド成分季節成分残差になります。 それぞれの値の関係として

観測データ(生データ) = トレンド成分 + 季節成分 + 残差 が成り立ちます。

f:id:dr_paradi:20170201014033p:plain

季節成分を除いたトレンドを見てみます。

trend = res.trend
trend = pd.DataFrame({'trend': trend, 'date':daily_data.index})
trend['date'] = pd.to_datetime(trend['date'], format='%Y-%m-%d')
trend = trend.set_index(['date'])
trend = trend.plot()

f:id:dr_paradi:20170201010642p:plain

生データと比較すると滑らかになっているのがわかるかと思います。 5月前半のゴールデンウィークなどの大型連休で凹んでいるので、事前に大きな休日などの変化要因を除くことでよりトレンドを正確に把握できるかもしれません。

また、消費電力外気温と適温との差 と考えることもできるので消費電力から気温のトレンドがより分かるようになります。今回は単年でのトレンドを追っていますが、年度ごとに比較することで年ごとの傾向がよりわかりやすくなるでしょう。

次に曜日成分を見てみます。

seasonal_data = res.seasonal
seasonal = pd.DataFrame({'seasonal': seasonal_data, 'date':daily_data.index})
seasonal['date'] = pd.to_datetime(seasonal['date'], format='%Y-%m-%d')
seasonal['weekday'] = seasonal['date'].dt.dayofweek
seasonal[0:7]

f:id:dr_paradi:20170201011320p:plain

weekdayの0は月曜日なので、土日に消費電力は大きく下がることがわかります。おそらく家庭での消費電力よりも工場やオフィスでの消費電力が多いことが影響していそうです。また、月曜日の消費電力が少なく出ていますが、こちらは月曜日に祝日が多いことも影響してそうです。事前にデータを除去したり、平滑化することでトレンド抽出がより正確になる可能性もあります。

以上のように、ただデータを見るだけではなく、季節成分とトレンド成分を分離することで、トレンドの把握が容易になるだけではなく、いろいろな知見を得ることができます。

おまけ: 時間別に見てみる

上記では、1日単位での変動を分析しましたが、時間帯によっても季節変動の差がでるので時間帯別でも同様にトレンド成分を分解してみてみます。

df = row_data.reset_index()
morning_power = df[df['time'] == '8:00'].reset_index()
morning_power['date'] = pd.to_datetime(morning_power['date'], format='%Y/%m/%d')
morning_power['power'] = morning_power['power'].values.astype(int)
morning_data = morning_power.groupby('date')['power'].sum().reset_index()
morning_data = morning_data.set_index(['date'])
res = sm.tsa.seasonal_decompose(morning_data.values, freq=7)
res.plot()

朝8時台のみ抽出したデータ f:id:dr_paradi:20170131013747p:plain

昼の12時台のみ抽出したデータ f:id:dr_paradi:20170131013738p:plain

非常に当たり前の結論になりますが、 トレンド成分を見ると、朝は冬のほうが消費電力が高く、昼は夏の方が消費電力が高いことがわかります。(生データで見てもわかるといえばわかりますが…) 時間帯別に別にトレンドを見ることで当たり前の結論も可視化で明確にすることができます。

まとめ

  • サービスを開発、改善するにあたってのKPIは多くの要因によって変化することが多いので正確に評価するのはなかなか大変です
  • 今回紹介したPythonのstatsmodelを用いて季節変動などの周期的なノイズを除去することでトレンドの把握がしやすくなります
  • 短期的な変動に左右されずにサービスの状態を把握してサービス改善につなげましょう

今後

ノイズ除去後のデータで変化検知とか異常値検知の自動化もしたいですね。