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

Gunosyデータ分析ブログ

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

【これからの強化学習】 Gunosy データマイニング研究会 #118 を実施しました

gunosy-dm.connpass.com

こんにちは。グノシー開発部のアルシャマンです。最近は、KID FRESINOのSalve feat. JJJをよく聴いています。

今日は4/12(水)に開催したGunosy DM #118について紹介します。前回に引き続きこれからの強化学習の1.3~1.5節の輪読と、論文紹介を行いました。

Gunosy DMこれからの強化学習については、以下のブログ記事で紹介しています。 data.gunosy.io

書籍輪読(これからの強化学習)

データ分析部の大曽根と吉田からそれぞれ1.3~1.4節と1.5節についての発表がありました。

1.3節では、MDP(マルコフ決定過程)における価値関数の表現と、それを推定するアルゴリズムについて学びました。具体的には、ある方策πのもとでの行動価値関数について成立する再帰式であるベルマン方程式とSarsaという学習法、最適行動価値関数について成立する再帰式であるベルマン最適方程式とQ-learningという学習法についてです。

1.4節では、方策を行動価値関数とは別のパラメータで表現された確率モデルと考え、そのパラメータについて最適化することで強化学習問題を解く方策勾配に基づく方法について学びました。

1.5節ではPOMDP(部分観測マルコフ決定過程)における強化学習について学びました。「エージェントは状態を部分的にしか観測できない」とした場合の強化学習についてです。どの状態にいるのかを表す確率分布である信念状態を導入して、価値関数を表現します。

論文紹介

次にデータ分析部の関から、Optimizing the Recency-Relevancy Trade-off in Online News Recommendationsという論文の紹介がありました。

ニュースサイトのトップページに掲載するニュースを選択する時に考慮される「Recency(新しさ)」と「Relevancy(ニュースの重要さ)」というトレードオフ関係にある2つの指標についての分析と「Highest Future-Impact(ニュース掲載後にどれだけPVを稼ぐか)」に基づいたニュース選択方法の提案がされています。

最後に

今後もGunosy DMは隔週開催予定です。次回の#119は、4/25(火)を予定しています。

次回の内容は、今回論文発表を担当した関からこれからの強化学習の2.1~2.3節の輪読と、同じくデータ分析部の【Edward】MCMCの数学的基礎からStochastic Gradient Langevin Dynamicsの実装までというブログの著者のmathetakeからA Complete Recipe for Stochastic Gradient MCMCという論文の紹介となっています。

興味ある方は、ぜひご参加下さい。

gunosy-dm.connpass.com

【これからの強化学習】 Gunosy データマイニング研究会を実施しました

gunosy-dm.connpass.com

こんにちは、データ分析部の阿部です。 今回は、先日開催したデータマイニング研究会という勉強会についてご紹介します。

データマイニング研究会とは

本勉強会では書籍の輪読と論文紹介を行い、データマイニングに関する基礎知識の向上及び、先端事例の共有・議論を行うことを目的としています。 2週間に1回のペースで開催されており、社外にも公開し広く知見を共有することを目指しています。 Gunosy創業時から取り組んでいるためこの手の勉強会としては歴史は長く(?)、今回で117回目になりました。

f:id:y-abe-hep:20170407102639j:plain

f:id:y-abe-hep:20170407102646j:plain

これからの強化学習

今回からは「これからの強化学習」を進めていて、1.1と1.2を終わらせました。 内容は強化学習の基礎的なところで、強化学習の構成要素が中心となっています。

これからの強化学習

これからの強化学習

  • 作者: 牧野貴樹,澁谷長史,白川真一,浅田稔,麻生英樹,荒井幸代,飯間等,伊藤真,大倉和博,黒江康明,杉本徳和,坪井祐太,銅谷賢治,前田新一,松井藤五郎,南泰浩,宮崎和光,目黒豊美,森村哲郎,森本淳,保田俊行,吉本潤一郎
  • 出版社/メーカー: 森北出版
  • 発売日: 2016/10/27
  • メディア: 単行本(ソフトカバー)
  • この商品を含むブログ (2件) を見る

本書はこれから研究を始める大学院生や研究者に向けて、最新の研究を理解する手がかりになることを目指しているとあります。 実務でも強化学習は応用されており、例えば広告のクリエイティブ選択に強化学習の一種とも言えるバンディットアルゴリズムが使われていたりします。 この他にも応用できそうな箇所もあるので、将来的に本書の内容を活かしてサービスへの応用もできればと考えています。

ちなみに、前回までは「トピックモデルによる統計的潜在意味解析」を進めていました。 こちらも良書なので、紹介しておきます。

トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ)

トピックモデルによる統計的潜在意味解析 (自然言語処理シリーズ)

論文紹介

論文紹介はこちらです。

Focus on the Long-Term: It’s better for Users and Business(KDD2015)

Focus on the Long-Term: It's better for Users and Business

Web系企業ではA/Bテストによってプロダクトの品質を分析するケースが多いですが、どうしても短期的な指標の評価に重点を置いてしまい長期的な価値が評価しきれていないこともあります。 この論文ではどのようにして長期的価値を評価するのかということを論じており、興味深かったです。

最後に

今回はGunosyデータマイニング研究会について紹介しました。 興味があったらぜひご参加ください!

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

こんにちは。初めまして。
データ分析部新入りの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へ移行する話

はじめに

こんにちは、データ分析部の森本です。主な業務は記事配信アルゴリズムの改善とログ基盤の整備です。
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のような自由度はないため要件に応じた利用が求められます