こんにちは。スマートニュースの高橋力矢です。前回のブログでデータ分析+ゲーム理論を題材として、帰納と演繹をまとめる利点をお伝えしました。なんらかの入力 (e.g., ゲーム理論における利得表) があり、特定のアルゴリズム (e.g., 各プレイヤーの戦略的意思決定) を記述することで出力 (e.g., ナッシュ均衡) を得るアプローチは、ほとんどのソフトウェア・エンジニアが慣れ親しんでいるプログラミングそのものです。つまり多くのエンジニアが手がけるプログラミングの実態は演繹的プログラミングです。ではこの対極に位置する帰納プログラミング (Inductive Programming) はどの程度進歩しているでしょうか。

帰納プログラミングの一分野である確率プログラミング (Probabilistic Programming) は統計学や機械学習との関係が密接で、日本でも利用者の多いStanをはじめとして種々の確率プログラミング言語が盛んに開発されています。機械学習分野のトップ国際会議であるNeural Information Processing Systems (NIPS)でも確率プログラミングに関するワークショップがすでに3回開催されていますし、関連技術であるProgram Inductionにも大きな関心が集まっています。

  • Inductive ProgrammingとProgram Inductionは名前が単なる倒置関係でかつ関連がありますが、別の技術エリアですので違いは後述します

帰納プログラミングの発展が産業界・学術界にもたらす恩恵は予想以上に大きいかもしれません。そこで帰納プログラミングが普及し演繹的プログラミングと組み合わさると何が起き得るか考えてみました。将来予測というのはそもそも殆どあたりませんので、確実な内容を求める読者のために背後にある数理についても説明します。キーワードは確率的変分近似法 (Stochastic Variational Inference; SVI) と自動微分変分法 (Automatic Differential Variational Inference; ADVI) です。ADVIを可能にするSVIの変数変換トリックはとても面白く、モンテカルロ法が必要な他の問題でも応用が利きそうです。

要旨

  • 演繹的プログラミングでは、入力(引数)に対して出力を計算する関数をプログラマーが書き下す。そしてプログラムに入力を与えて出力を得る
  • 帰納プログラミングでも同様に、入力を出力に変換する関数を記述する。しかしプログラムの実行時には出力を与えて入力を復元する。コンパイラーがある種の逆問題を解いていることになる
  • 確率プログラミングが発展することで、2つの社会的変化・イノベーションが起きるかもしれない
    • あらゆる種の数理的研究で本質への集中が容易になる。サイエンティストは数学的利便性を気にすることなく好きなモデルを設計し、そのパラメータの推定はコンパイラーに丸投げできる
    • 演繹+帰納が統合されたソフトウェアを、統計学を知らない多くのプログラマーが生み出していく
  • 出力から入力を復元するという自明でない推論のために、多くの確率プログラミング言語はマルコフ連鎖モンテカルロ法 (MCMC)か SVIを用いる。SVIの高速性は実務上の価値を高めている
  • SVIおよびADVIを可能とする変数変換トリックについて、その着眼点の良さを味わってほしい

確率プログラミングによる「推論機」の台頭

計算機 (computer) という名称は、入力から出力を得る機械一般に与えられます。「計算」という単語は、順方向の演繹を意味しているわけです。計算機を使うとき、入力から出力に変換する手順を人が明示的に指定します。たとえば電卓も計算機の一種で、数値(入力)と計算手順(e.g., 足し算, 掛け算)を与えることで計算結果(出力)を得ます。順方向な計算機の性質ゆえに多くのプログラマーは演繹的な思考になれています。彼らは研究者で言えばゲーム理論家やオペレーションズ・リサーチャーに似ています。

一方で確率プログラミングと呼ばれるプログラミング・パラダイムは帰納を実現し、出力から入力が復元されます。出力から入力への復元手順をプログラマーが書くことはありません。従来の演繹的プログラミングと同様に、入力から出力への写像だけを記述します。確率プログラミング言語の代表例であるStanはベイズ推論を通した帰納プログラミングを可能にしています。プログラマーが入力から出力を生成する確率過程を記述するだけで、出力から入力を復元するアルゴリズムをStanが自動生成します。

内部のコンパイラーの仕組みを気にしなければ、確率プログラミング言語は計算機ではなく推論機を提供しています。計算機とは演繹を行う機械の総称、推論機とは帰納を行う機械の総称です。推論機がWindows, Mac OS, Linuxが走っているごく普通の計算機=演繹機械である計算機上で動く理由は、統計学が「精確な帰納を実現する方法を確率論から演繹して突き詰めた学問」だからです。つまり、マクロな挙動が推論機であっても一番ミクロな部分は演繹機械の組み合わせで成り立っています。しかしミクロな振る舞いを理解する必要があるのはbias-variance解析を行うような理論統計学者だけです。一般のStanユーザー、統計学をツールとして使うサイエンティストはミクロな構造は忘れてよく、マクロとしての推論機が使えると考えれば十分でしょう。

サイエンティストがStanを使う一番のメリットは、モデルパラメータの推定手順を書かなくてよい点です。サイエンティストは現象を生み出したと予想される背後のモデルを考えればよく、モデルパラメータの推定はStanにまかせられます。現段階では確率プログラミング言語を使うのは「データサイエンティスト」が大半ですが、そのうち物理学や経済学、数理生物学といった多くの分野の専門家が盛んに使うようになるでしょう。彼らの仕事はモデルの設計とそこからのインプリケーション獲得であって、推定アルゴリズム自体を一から書く機械学習エンジニアとは職務が違うからです。例えばStanの場合は以下のようなコードを書きます。これは説明変数$X$が1次元スカラーで従属変数$Y$が$\lbrace 0, 1 \rbrace$の二項選択である単純なロジスティック回帰です。

data {
  int<lower=0> N;
  real x[N];
  int<lower=0, upper=1> y[N];
}
parameters {
  real bias;
  real coefficient;
}
model {
  for (i in 1:N)
    y[i] ~ bernoulli(inv_logit(bias + coefficient * x[i]));
}

dataブロックでは説明変数と従属変数の観測値、つまりデータを指定します。parametersブロックではモデルパラメータを設定し、modelブロックでデータが観測される確率過程を記述しています。inv_logitはシグモイド関数の値を返します。modelブロックにはbias及びcoefficientを生成する事前分布も書くべきですが、事前分布の記述を省略した場合は一様分布が自動的に使われます。

parametersブロックのbias及びcoefficient変数が求めたいパラメータです。通常のプログラミング言語では、データからパラメータを復元するアルゴリズムを統計学者か機械学習エンジニアが書きます。しかしStanではパラメータからデータが生成される確率過程(順問題)を書き下すだけで、逆問題である推定アルゴリズムはコンパイラーが自動生成してくれます。

単なるロジスティック回帰であればPythonのscikit-learnパッケージやRのglmパッケージを用いれば推定アルゴリズムを書く必要はありません。しかしStanの場合は、modelブロックがどんなに複雑でも推定アルゴリズムを書かなくて良い点が大きく違います。モデルの構造ごとに、その都度統計学者に依頼して個別の推定アルゴリズム設計に時間をかけていたプロセスを省けるわけです。

帰納プログラミングがもたらすイノベーション

多くの数理科学は演繹的であって、知見は本質的にはモデルから得られます。科学者は、有用な知見を生み出すのに十分な構造を備えたモデルを最初に考案しなければいけません。ところが知見を重視するモデルはしばしば、統計学者や機械学習エンジニアなどパラメータ推定の専門家からは忌避されます。

統計学者たちは推定における数学的利便性を大切にする傾向があります。例えば尤度が解析的であり、closed-form式による効率的な推定アルゴリズムが存在するモデルは重宝されます。一方でパラメータ推定の困難なモデルは、たとえモデル自体が真実に近くても実務では忌避されます。パラメータ推定に膨大な計算時間が必要なこと、もしくは貧弱な局所最適解しか得られないことがその理由です。

例えば経済学者は、感度解析を目的として特定の構造を持ったパラメトリック効用モデルを好むかもしれません。しかし機械学習の専門家は、全ての効用関数を近似可能な(=万能近似性を持つ)ノンパラメトリック・モデルを統一的に使う方がシンプルで効率的だと考える傾向があり、両者の視点は食い違っています。特定構造のパラメトリックモデルは複雑な推定アルゴリズムの実装を要求しますが、ノンパラメトリックモデルの多くは凸最適化や二分木の適合など簡潔なアルゴリズムで推定可能なためです。最近使われているDeep Neural Network (DNN)による効用モデルも万能近似性をもつためモデル構造を強制しませんし、バックプロパゲーションによる最適化はシンプルです。つまるところ、知見の取り出しやすさを重視する立場と推定のしやすさを重視する立場との間に、深い河が流れているわけです。

しかし本質に立ち返って考えれば、数学的利便性に過度にこだわってイノベーションを妨げるのは愚かなことであって、モデルを自由に設計できるツールがあればそれを使うべきでしょう。推定の便宜のためにモデル構造を強制するのは、イラスト編集ソフトで圧縮の効率性を保つために一部の編集機能が使えなくなるようなもので、本末転倒です。

Stanのような推論機が存在することで、演繹専門家タイプの自然科学者や社会科学者が特に大きな恩恵を受けます。彼らはモデルのデザインと知見の獲得という本質に集中できるようになります。もはや統計学者からモデルの形を強制されることはありません。

加えて、今後のソフトウェア・エンジニアリングでは、計算機と推論機を自在に組み合わせたソースコードが多く現れるかもしれません。推論部分を精度よく行うブラックボックスが現れることで、エンジニア達も同様にアプリケーションの本質に集中できるようになります。もしそのようなパラダイム転換が起きれば人間と機械との協業が飛躍的に増大して多大な産業的成果が期待できます。

どうやって帰納をブラックボックスするのか

帰納プログラミングのユーザーから見たメリットが明らかになったところで、内部で使われている技術を一つ紹介したいと思います。Stanはベイズ推定を基本とし、パラルータの事後分布サンプルを与えてくれます。事後分布サンプルの平均やヒストグラムを分析することで、パラメータの推定値とその信用区間(=不確実性)がわかるようになっています。

ベイズ推定

ベイズ推定の基本式を見ましょう。例えば入力$X$から出力$Y$を予測する回帰モデルのパラメータを推定するケースを考えます。入力サンプルの集合を${\mathcal X}$、出力サンプルの集合を${\mathcal Y}$とし、推定するパラメータ集合を$\Theta$としましょう。入力${\mathcal X}$とパラメータ$\Theta$から出力${\mathcal Y}$が生成される確率密度関数または確率質量関数を$p({\mathcal Y}\vert {\mathcal X}, \Theta)$、事前分布を$p(\Theta)$とすれば事後分布は一般的に次の式(1)として書けます。以降の利便性のため確率分布$p$上の期待値を演算子${\mathbb E}_p$で記述します。${\mathbb E}_{p(\Theta)}\left[f(\Theta)\right]\triangleq \int_\Theta f(\Theta)p(\Theta)d\Theta$です。

$$ p(\Theta\vert {\mathcal X},{\mathcal Y})= \dfrac{p({\mathcal Y}\vert {\mathcal X}, \Theta)p(\Theta)} {\int_{\Theta}p({\mathcal Y}\vert {\mathcal X}, \Theta)p(\Theta)d\Theta} = \dfrac{p({\mathcal Y}\vert {\mathcal X}, \Theta)p(\Theta)} {{\mathbb E}_{p(\Theta)}\left[p({\mathcal Y}\vert {\mathcal X}, \Theta)\right]} ~~~~~~ (1) $$

式(1)において分母の積分が解析的に得られるケースは稀です。そこでモンテカルロ法を用いて事後分布$p(\Theta\vert {\mathcal X},{\mathcal Y})$に従うサンプルを得るのが一般的です。マルコフ連鎖モンテカルロ法 (Markov Chain Monte Carlo; MCMC)を用いる場合は、定常分布が$p(\Theta\vert {\mathcal X},{\mathcal Y})$となるようなマルコフ過程を設計してシミュレートしますが、マルコフ連鎖の収束には時間がかかります。今回は説明を省きますが、収束性の改善や多峰性のある事後分布への効率改善を目的としてMCMCも著しく発展してきています。

MCMCのうち、Hamiltonian Monte Carlo (HMC)やLangevin Monte Carlo (LMC)は事後分布上の勾配相当情報を利用しているのが特徴的です。MCMCより高速なベイズ推定法としてはラプラス近似と変分ベイズ法が知られていますが、このうち変分ベイズ法にモンテカルロ法による勾配計算を組み込んだのが確率的変分近似法 (Stochastic Variational Inference; SVI)です。

変分ベイズ法の一般形式

SVIを理解するために、まず変分ベイズ法の本質について説明します。真の事後分布$p(\Theta\vert {\mathcal X}, {\mathcal Y})$が解析的に得られることは希ですので、$p$を解析的にもっと扱いやすい別の分布$q(\Theta\vert {\mathcal X},{\mathcal Y})$で近似することを考えます。$q$を変分事後分布と呼びます。$q$のハイパーパラメータを計算する何らかの最適化アルゴリズムが得られれば目標達成です。

データ尤度を事前分布で積分した項を周辺尤度 (marginal likelihood もしくは evidence) と呼びます。そして周辺尤度の下限を最大化することで最適変分事後分布が得られます。Jensenの不等式を使うことで、真の周辺尤度に対してその下限を式(2)のように導出することができます。式(2)の右辺はEvidence Lower-BOund (ELBO)と呼ばれます。 $$ \log {\mathbb E}_{p(\Theta)}\left[p({\mathcal Y}\vert {\mathcal X}, \Theta)\right] \geq {\mathbb E}_{q(\Theta\vert {\mathcal X},{\mathcal Y})} \left[\log \dfrac{p({\mathcal Y}\vert {\mathcal X}, \Theta)p(\Theta)} {q(\Theta\vert {\mathcal X},{\mathcal Y})} \right]~~~~~(2) $$

式(2)の(左辺)-(右辺)、つまり周辺尤度からELBOを引いた統計量を $\Delta_{pq}({\mathcal X},{\mathcal Y}) \triangleq \log {\mathbb E}_{p(\Theta)}\left[p({\mathcal Y}\vert {\mathcal X}, \Theta)\right] - {\mathbb E}_{q(\Theta\vert {\mathcal X},{\mathcal Y})} \left[\log \dfrac{p({\mathcal Y}\vert {\mathcal X}, \Theta)p(\Theta)} {q(\Theta\vert {\mathcal X},{\mathcal Y})} \right]$と置きます。 ${\mathbb E}_{p(\Theta)}\left[p({\mathcal Y}\vert {\mathcal X}, \Theta)\right] \equiv p({\mathcal Y}\vert {\mathcal X})$ 及び $p(\Theta\vert {\mathcal X},{\mathcal Y})\equiv p({\mathcal Y}\vert {\mathcal X}, \Theta)p(\Theta) / p({\mathcal Y}\vert {\mathcal X})$ を利用すると

$$ \Delta_{pq}({\mathcal X}, {\mathcal Y}) \equiv \log\underbrace{p({\mathcal Y}\vert {\mathcal X})}_{evidence} - {\mathbb E}_{q(\Theta\vert {\mathcal X},{\mathcal Y})} \underbrace{\left[\log p(\Theta\vert {\mathcal X},{\mathcal Y})p({\mathcal Y}\vert {\mathcal X}) \right] }_{model\times prior\equiv posterior\times evidence} -{\mathbb E}_{q(\Theta\vert {\mathcal X},{\mathcal Y})}\left[\log q(\Theta\vert {\mathcal X},{\mathcal Y}) \right] $$ であり、整理すると $$ \Delta_{pq}({\mathcal X}, {\mathcal Y}) \equiv \underbrace{\log p({\mathcal Y}\vert {\mathcal X}) - \log p({\mathcal Y}\vert {\mathcal X})}_{cancelled} + {\mathbb E}_{q(\Theta\vert {\mathcal X},{\mathcal Y})} \left[\log \dfrac{q(\Theta\vert {\mathcal X},{\mathcal Y})} {p(\Theta\vert {\mathcal X},{\mathcal Y})} \right]\nonumber
= D_{KL}( q(\Theta\vert {\mathcal X},{\mathcal Y}) \Vert p(\Theta\vert {\mathcal X},{\mathcal Y}) ) $$ となって$\Delta_{pq}({\mathcal X}, {\mathcal Y})$ は変分事後分布と真の事後分布との間のKL-divergenceであることがわかります。つまり、ELBOを最大化する変分事後分布はKL-divergenceの意味で真の事後分布に最も近いことになります。このことはモデルや変分事後分布の関数形に依存せず成り立ちます。

  • 本来のKL-divergenceは$D_{KL}(q \Vert p)$ではなく$D_{KL}(p \Vert q)$ですが、その違いについては議論を省きます

ELBOは、変分事後分布によってデータ対数尤度+事前分布対数尤度を積分したものです。この積分が解析的に導出可能であれば、その変分事後分布のパラメータに関して微分を取り、勾配を登ることで最適変分事後分布に収束させることができます。

変分ベイズ法が機械学習コミュニティで流行り出した2000年代前半は、主にクラスタリング問題(=混合正規分布の適合)において最適変分事後分布を求める研究が流行りました。これは混合正規分布の推定と変分ベイズ法との相性が良かったためで、Dirichlet-正規逆ガンマ分布を事前分布とした場合の変分事後分布がclosed-formのEMアルゴリズムによって推定できました。

モンテカルロ法による、より自由な変分事後分布

今度はサイエンティストがモデルと事前分布を自由に設定する一般の場合を考えます。ELBOがclosed-formになるような特別な関係がモデル・事前分布・変分事後分布の間に存在する場合は、標準的な変分ベイズ法で最適変分事後分布が得られることに言及しました。しかし一般のケースでELBOがclosed-formになることは殆どありません。この状況でどうやって変分近似事後分布を求めたら良いでしょうか?

この場合も方針はELBO、つまり式(2)右辺の微分を求めて勾配を登ることです。$q$には何らかの関数形を仮定はしますが、ELBOに内在する積分の解析解を得るのは諦めます。代わりに積分をモンテカルロ近似します。変分事後分布をハイパーパラメータ$\Psi$でパラメトライズして$q(\Theta; \Psi)$と書きましょう。

変分事後分布から仮に、$k$個のモンテカルロサンプル$\Theta_1,\ldots,\Theta_k$が得られたとします。これを用いて式(3)により近似勾配を評価します。

$$ \frac{\partial{\mathbb E}_{q(\Theta; \Psi)}}{\partial\Psi}\left[\log {p({\mathcal Y}\vert {\mathcal X}, \Theta)p(\Theta)}\right] \simeq \frac{1}{k}\sum_{\ell=1}^k \left[ \frac{\partial \log p({\mathcal Y}\vert {\mathcal X}, \Theta_\ell)}{\partial\Psi} + \frac{\partial \log p(\Theta_\ell)}{\partial\Psi}\right] ~~~~~~~ (3) $$

しかし式(3)を評価するにあたっては二つのことが致命的で、求めたいハイパーパラメータと事後分布サンプルが鳥と卵の関係になっています (*)。

  • 1. データ尤度$p({\mathcal Y}\vert {\mathcal X}, \Theta_\ell)$及び事前分布尤度$p(\Theta)$がハイパーパラメータ$\Psi$の関数として書けないため微分が取れない
  • 2. 事後分布サンプル$\Theta_1,\ldots,\Theta_k$はハイパーパラメータ$\Psi$が分からないとサンプリングできない

アフィン変換トリックによる確率的更新

前節で問題となった鳥と卵の問題を解決するのが事後分布の変数変換トリックです。このテクニックは最近色々な所で類似のものを見かけるようになりましたが、私が(Titsias and Lázaro-Gredilla, 2014)で最初に知った時はそのセンスの良さに感銘を受けました。ここでは推定すべきパラメータはホワイトノイズからのアフィン変換であるという仮定をおきます。

  • まずパラメータ$\Theta$の実態は$d$次元のベクトル $\boldsymbol\theta \in {\mathbb R}^d$ とします
  • $r$次元のホワイトノイズベクトル $\boldsymbol\varepsilon \in \mathbb R^r$、射影行列$\boldsymbol A \in {\mathbb R}^{d \times r}$、バイアス項 $\boldsymbol{b} \in \mathbb R^d$を用いて、$\boldsymbol \theta = \boldsymbol A \boldsymbol \varepsilon + \boldsymbol b$ とおきます
    • $\boldsymbol \varepsilon$には無相関の単位多次元正規分布や多次元ラプラス分布、多次元$t$分布などが使えます
  • パラメータ集合$\lbrace \boldsymbol A, \boldsymbol b \rbrace$こそが最適化すべきハイパーパラメータ$\Psi$の実態です

そして、ホワイトノイズ$\boldsymbol \varepsilon$はその定義上、$\lbrace \boldsymbol A, \boldsymbol b \rbrace$に依存しないことに着目します。式(3)のモンテカルロ評価の実態は次の式(4)(5)により与えられます。

$$ \frac{\partial{\mathbb E}_{q(\Theta; \Psi)}}{\partial\boldsymbol A} \left[\log {p({\mathcal Y}\vert {\mathcal X}, \Theta)p(\Theta)} \right] \simeq \frac{1}{k} \sum_{\ell=1}^k\left[ \frac{\partial \log p({\mathcal Y}\vert {\mathcal X}, \boldsymbol{A}\boldsymbol{\varepsilon}_\ell + \boldsymbol b)}{\partial\boldsymbol{A}} + \frac{\partial \log p(\boldsymbol{A}\boldsymbol{\varepsilon}_\ell + \boldsymbol b)}{\partial\boldsymbol{A}} \right] ~~~~~~ (4) $$ $$ \frac{\partial{\mathbb E}_{q(\Theta; \Psi)}}{\partial\boldsymbol b} \left[\log {p({\mathcal Y}\vert {\mathcal X}, \Theta)p(\Theta)} \right] \simeq \frac{1}{k} \sum_{\ell=1}^k\left[ \frac{\partial \log p({\mathcal Y}\vert {\mathcal X}, \boldsymbol{A}\boldsymbol{\varepsilon}_\ell + \boldsymbol b)}{\partial\boldsymbol{b}} + \frac{\partial \log p(\boldsymbol{A}\boldsymbol{\varepsilon}_\ell + \boldsymbol b)}{\partial\boldsymbol{b}} \right] ~~~~~~ (5) $$

無相関ノイズのアフィン変換として事後サンプルを与えるアプローチの優位性は、勾配評価において明らかになります。今回の前提として、データ対数尤度$f(\boldsymbol \theta) \triangleq \log p({\mathcal Y}\vert {\mathcal X}, \boldsymbol\theta)$ と事前分布対数尤度$g(\boldsymbol \theta) \triangleq \log p(\boldsymbol \theta)$が解析的に与えられていることを思い出してください。すると式(4)(5)の右辺は実際には対数尤度を微分した関数$\boldsymbol{\nabla}f(\boldsymbol{\theta})\triangleq \frac{\partial f(\boldsymbol\theta)}{\partial\boldsymbol\theta}, \boldsymbol{\nabla}g(\boldsymbol{\theta}) \triangleq \frac{\partial g(\boldsymbol\theta)}{\partial\boldsymbol\theta}$を使って式(6)(7)のように変形できます。

$$ \frac{\partial{\mathbb E}_{q(\Theta; \Psi)}}{\partial\boldsymbol A} \left[\log {p({\mathcal Y}\vert {\mathcal X}, \Theta)p(\Theta)} \right] \simeq \frac{1}{k} \sum_{\ell=1}^k\left[ \boldsymbol{\nabla}f(\boldsymbol{A}\boldsymbol{\varepsilon}_\ell+\boldsymbol b) +\boldsymbol{\nabla}g(\boldsymbol{A}\boldsymbol{\varepsilon}_\ell+\boldsymbol b) \right]\boldsymbol{\varepsilon}_\ell^\top ~~~~~~ (6) $$ $$ \frac{\partial{\mathbb E}_{q(\Theta; \Psi)}}{\partial\boldsymbol b} \left[\log {p({\mathcal Y}\vert {\mathcal X}, \Theta)p(\Theta)} \right] \simeq \frac{1}{k} \sum_{\ell=1}^k\left[ \boldsymbol{\nabla}f(\boldsymbol{A}\boldsymbol{\varepsilon}_\ell+\boldsymbol b) +\boldsymbol{\nabla}g(\boldsymbol{A}\boldsymbol{\varepsilon}_\ell+\boldsymbol b) \right] ~~~~~~ (7) $$

式(6)(7)が実装上優れている点は、変分事後分布からのサンプルが勾配評価に必要ないことです。無相関ノイズベクトルのサンプル$(\boldsymbol\varepsilon_\ell)_{\ell=1}^k$さえ与えればよく、(*)で問題とされたハイパーパラメータと事後分布サンプルとの間の鳥と卵の関係は発生しません。乱数を使って生成したサンプルを元に式(6)(7)を用いて変分事後パラメータ$\lbrace \boldsymbol A, \boldsymbol b \rbrace$の更新を繰り返していくことで、最適変分事後分布を得ることができます。

確率的変分近似法 (SVI)

ビッグデータを用いるここ数年の機械学習コミュニティでは、推定の高速化を目的とした確率的勾配降下法 (Stochastic Gradient Descent; SGD)の適用が盛んです。SGDは対数尤度関数の評価に、観測された全サンプルではなく一部のランダムサブセットを用いて高速なパラメータ更新を行います。変分ベイズ法に対してSGDを適用したものを確率的変分近似法 (Stochastic Variational Inference; SVI)と呼びます。

多くの確率的数理モデルはindependent and identically distributed (i.i.d.)な確率過程によって記述されています。I.i.d.データの場合、式(6)(7)右辺のモンテカルロ勾配におけるデータ対数尤度項は、$\boldsymbol{\nabla}f(\boldsymbol{A}\boldsymbol{\varepsilon}_\ell+\boldsymbol b)=\sum_{j=1}^n \boldsymbol{\nabla} f_j(\boldsymbol{A}\boldsymbol{\varepsilon}_\ell+\boldsymbol b)$ ($n$はサンプルサイズ) と独立した$n$項の和に分解できます。

$n$項の和からなるデータ対数尤度勾配を、そのうちの$m~(\ll n)$部分サンプルの和によって確率的に近似することで、式(6)(7)よりも高速なハイパーパラメータ更新をもたらす式(8)(9)が得られます。式(8)(9)においてインデクス集合$\Pi$は$\Pi=\lbrace j; 1 \leq j \leq n \rbrace$, $\vert \Pi \vert \equiv m$を満たすようにランダムに生成されます。

$$ \frac{\partial{\mathbb E}_{q(\Theta; \Psi)}}{\partial\boldsymbol A} \left[\log {p({\mathcal Y}\vert {\mathcal X}, \Theta)p(\Theta)} \right] \simeq \frac{1}{k} \sum_{\ell=1}^k\left[ \frac{n}{m} \sum_{j \in \Pi} \boldsymbol{\nabla}f_j(\boldsymbol{A}\boldsymbol{\varepsilon}_\ell+\boldsymbol b) +\boldsymbol{\nabla}g(\boldsymbol{A}\boldsymbol{\varepsilon}_\ell+\boldsymbol b) \right]\boldsymbol{\varepsilon}_\ell^\top ~~~~~~ (8) $$ $$ \frac{\partial{\mathbb E}_{q(\Theta; \Psi)}}{\partial\boldsymbol b} \left[\log {p({\mathcal Y}\vert {\mathcal X}, \Theta)p(\Theta)} \right] \simeq \frac{1}{k} \sum_{\ell=1}^k\left[ \frac{n}{m} \sum_{j \in \Pi} \boldsymbol{\nabla}f_j(\boldsymbol{A}\boldsymbol{\varepsilon}_\ell+\boldsymbol b) +\boldsymbol{\nabla}g(\boldsymbol{A}\boldsymbol{\varepsilon}_\ell+\boldsymbol b) \right] ~~~~~~ (9) $$

自動微分変分法 (ADVI)

式(8)(9)において、データ対数尤度や事前分布対数尤度が比較的簡単な解析式で与えられている場合は、その微分を手計算で導出してハードコーディングしても良いでしょう。しかし、様々なサイエンティストが異なるモデル式を用いる確率的プログラミング言語においては、その微分は自動的に得られるべきです。ニューラルネットワークのような階層的なモデルの場合、データ対数尤度を式で書き下すととても長くなります。式を明示的に与えるよりも、式を結果的に評価する関数を実装した上でその関数ポインターを渡す方が使いやすそうです。

つまりデータが生成される確率過程をstep-by-stepで書き下すインターフェースがサイエンティストによって使いやすいわけです。その一つの実装形態がStanのmodelブロックになります。コンパイラー側では、このstep-by-stepで書き下された対数尤度関数の微分を自動的に取得する必要があります。その際には自動微分技術が使われます。Pythonであればautogradパッケージが便利です。

データ対数尤度の自動微分を利用する変分ベイズ法を自動微分変分法 (ADVI)といいます。SVIによる高速化とADVIによる汎用性はStanを含む確率的プログラミング言語の必須要素であり続けるでしょう。

まとめと更なる展望

今までソフトウェア・エンジニア達が親しんできた演繹的プログラミングに加えて帰納プログラミング言語が実用に耐える水準に上がってきたこと、それに伴う産業・科学上のイノベーションについて論じました。特に帰納プログラミング言語の一形態である確率的プログラミング言語の発展が著しいこと、その背後には確率的変分近似法 (SVI)と自動微分変分法 (ADVI)の成熟があることを紹介しました。SVI及びADVIの背後には自動微分パッケージを利用可能にする変数変換トリックの美しさがあります。

入力から出力へ変換する関数を書いた上で、その順番通り入力を与えて出力を得るのが演繹的プログラミングです。一方で、同じように入力から出力への変換関数を書くにも関わらず、出力から入力を復元するのが帰納プログラミングです。もう一つ面白いプログラミング・パラダイムが考えられることに気づいたでしょうか。それは入力と出力を与えて変換関数をコンパイラーが自動生成するパラダイムで、プログラム帰納 (Program Induction)と呼ばれています。統計的回帰問題において説明変数・従属変数が数値とは限らず構造化データであるような、一般化された回帰問題がプログラム帰納です。

プログラム帰納の実現においてもその背後では確率的な推論テクニック・最適化テクニックが広範に使われています。この技術が成熟してきた暁には、例えば厳密なフォーマットが分かっていないテキストデータから、重要な情報がどこにあるのかを検知して一覧表にする特殊なパーサーを自動生成できます。演繹的プログラミング、帰納プログラミングに加えて、プログラミング帰納まで統合したソースコードが出回る日も遠くないかもしれません。

統合コードを実装するにあたってその質の鍵となるのは、ソフトウェアのどの部分を人間がプログラムし、どの部分を機械に任せるかのメタレベルでの意思決定です。人間と機械との間の分担または協調について、あるべき姿が今日の人工知能の発展により盛んに議論されています。その議論の具体的試金石となる最初のステップは、演繹・帰納が統合化されたソフトウェア開発の現場かもしれません。つまり、ソフトウェア開発に詳しくない人々の仕事が残る・無くなるという議論よりも先に、ソフトウェア・エンジニア達こそが分業と意思決定に関する最初の試練を受ける可能性があります。自動化される仕事とされない仕事との区分について、現在のメディアで見かける予測とは全く違った統廃合が起きる気がしています。その社会変化に肯定的に付き合う自信を持ったエンジニアの方には、下記の求人は面白いかもしれませんよ。

References