Multioutput Gaussian Processで遊ぶ（Fun Advent Calendar 2020 #23）

0．はじめに

Fun Advent Calender 2020の23日目です*1。僕はこのAdvent Calendarが作成された大学に所属する大学院生で、現在は博士後期課程3年（D3）に在籍しています。市内にある理学療法士・作業療法士養成校の教員（常勤）も務めており、研究と講義に追われる日々を送っています。研究では片麻痺の患者さんを対象として様々なデータをごにょごにょしています。専門分野はありません。
（※だいぶ前にTwitterの名前を2歳に変えてしまいましたがseibisiです。ちなみに2歳ではなく27歳です）

今回はMulti-output Gaussian Process（MOGP）というもので遊んでみました。MOGPはガウス過程の多出力版であり、正直あまり上手い日本語訳が分からない*2ので、英語のままにしておきます。Gaussian Processはそのままガウス過程であり、機械学習に興味がある人であれば1度は聞いたことのあるものだと思います。

この記事は、ガウス過程……聞いたことはあるけどなんだかよく分からない……難しい……という学部生を対象として、ガウス過程にHello Worldしてみよう！という目的で書いてみました。機械学習つよつよの方にとってはつまらない内容かもしれませんが、この記事を書くに当たって参考にした文献も載せておくので、理論に興味がある方はそちらを参照してください。

1．前準備

ガウス過程の話をする前に多変量正規分布について述べたいと思います。多変量正規分布はガウス過程における大きな立役者*3となります。この記事では終始一貫して多変量正規分布の話しか出てきません。

まずは定義から述べていきます。確率変数ベクトル $\displaystyle {\bf{x}}$ が平均ベクトル $\displaystyle {\boldsymbol{\mu}}$ 、分散共分散行列 $\displaystyle {\boldsymbol{\Sigma}}$ の多変量正規分布に従う時、多変量正規分布は次の式で定義されます。

$\displaystyle N( {\bf{x}} | {\boldsymbol{\mu}}, {\boldsymbol{\Sigma}} ) = \frac{1}{ (2\pi)^D \sqrt{ | {\boldsymbol{\Sigma}} | } } \exp ( -\frac{1}{2} ( {\bf{x}} - {\boldsymbol{\mu}} )^T {\boldsymbol{\Sigma}}^{-1} ( {\bf{x}} - {\boldsymbol{\mu}} ) )$

確率変数が複数存在する場合は確率変数間の対応関係を考える必要があり、この情報は分散共分散行列に存在しています。対応関係については、本学でも統計もしくは統計に関わる講義で共分散という言葉は1度は聞いたことがあるのではないでしょうか。聞いたことが無い人は、共分散と同じ意味を持つ言葉で相関という言葉が存在しているので、イメージとしてはそれでも大丈夫です。分散共分散行列は非常に重要なので、この記事を読む間は頭の片隅に置いておくと良いかもしれません。

2．Single-output Gaussian Process (SOGP)

まずは出力が1つのガウス過程（SOGP）から始めていく事にします。

2-1　SOGPの雰囲気

まずは下の図を見てみます。

図は何やら色とりどりのぐにゃぐにゃがいっぱい並んでいる様子です。横軸は $\displaystyle x$ 、縦軸には $\displaystyle f(x)$ と書いてあります。このぐにゃぐにゃたちの正体は、ガウス過程から生成されたサンプルです。え？グラフを生成？と混乱する人もいるかもしれませんが、次の図を見てみましょう。

図は、5、10、50、100変量の正規分布から乱数を発生させて結果を可視化したものです。横軸は確率変数ベクトル内の要素の番号（確率変数の番号）、縦軸はその実現値となるように可視化しました。なんとなくグラフっぽく見えるような気がします。確かに生成されているような…？

上の図は変量数を5、10、50…と進めていきましたが、これをさらに100、500、1000……さらには変量数が無限の正規分布からデータをサンプルして、同じようにグラフを描いたらどうなるでしょうか？これが1番最初の図になります。1つ1つのグラフは多変量正規分布からの「サンプル」と考える事もできると思います。

2-2　SOGPの定義

なんとなく雰囲気を掴めたところで、SOGPの定義を確認しておくことにします。出力をひとまとめにしたベクトルを $\displaystyle {\bf{f}} = ( f(x_1), f(x_2), ... ,f(x_N) )^T$ とした時、ガウス過程は次の式で定義されます。

$\displaystyle {\bf{f}}～GP( {\boldsymbol{\mu}}, {\bf{K}} )$

ここで、 $\displaystyle {\boldsymbol{\mu}}$ は平均ベクトル、 $\displaystyle {\bf{K}}$ は分散共分散行列であり、出力のベクトルは多変量正規分布に従います。話を簡単にするために、平均ベクトルを $\displaystyle {\bf{0}}$ としておきましょう。

$\displaystyle {\bf{f}}～GP( {\bf{0}}, {\bf{K}} )$

2-3　SOGPの分散共分散行列

さて、ガウス過程では分散共分散行列が大事（主観）という話をしたわけですが、SOGPの分散共分散行列についても見ていきましょう。SOGPでは分散共分散行列の $\displaystyle (i, j)$ 成分は $\displaystyle k(x_i, x_j)$ のように表現され、 $\displaystyle k(x_i, x_j)$ はカーネルと呼ばれています*4。カーネルには以下で述べるもの以外にも色々な種類があります。

ガウシアンカーネル

$\displaystyle k(x_i, x_j) = \theta_1 \exp ( -\frac{ |x_i - x_j| ^2 }{ \theta_2 } )$

指数カーネル

$\displaystyle k(x_i, x_j) = \exp ( -\frac{ |x_i - x_j| }{ \theta } )$

カーネルが何であるかの説明は詳しくはしませんが、SOGPでは分散や共分散を計算するための謎の式くらいでひとまず大丈夫です。試しにガウシアンカーネルと指数カーネルから生成したぐにゃぐにゃを見比べてみる事にします。

上が共分散の計算に指数カーネルを用いたもの、下がガウシアンカーネルを用いた物です。カーネルの選び方によって、だいぶ概形が変わるようですね。ヒートマップは分散共分散行列を可視化したものです。ヒートマップからは、隣接するデータの相関が強く、離れるほど相関が弱くなる事が分かります。この他にもカーネルは数多く存在しており、時に複数のカーネルを組み合わせたりすることもできます。カーネルによってグラフの概形や振る舞いが変わるようですね。

2-4　SOGPまとめ

ここまでで、なんとなくSOGPの雰囲気は掴めたかなと思います。重要なポイントは、以下の2つです。

ガウス過程は関数をサンプルするもので、その正体は多変量正規分布（無限次元）
分散共分散行列の各要素はカーネルなるものを使って計算、カーネルによってグラフの概形や振る舞いが変わる。

次からMOGPを述べていきます。

3．Multi-output Gaussian Process (MOGP)

さて、ここからが本題です。以下のような場面を考えてみます。

3つの色のデータが存在していますが、これらは系列単位で互いに相関しているものとしましょう。系列間の相関を捉えながらどのようにしてガウス過程を適用していきましょうか……。結論、これは圧倒的力技です。何が力技かというと、系列を順番に愚直に並べて1つの多変量正規分布にまとめてしまうという事をします。例えば、系列が $\displaystyle P$ 個あり、それぞれの系列でデータの数が $\displaystyle N$ 個あるとします。ここで、それぞれの系列は多変量正規分布に従う（ガウス過程に従う）ものとします。MOGPでは、この出力を縦に全部並べて長さ $\displaystyle NP$ の確率変数ベクトルを作成します。そうすると、長さ $\displaystyle NP$ の期待値ベクトル、 $\displaystyle NP×NP$ 個の要素を持つ分散共分散行列をパラメータとして多変量正規分布を作ることができます*5。気持ち良いくらいに力技です。

多出力のガウス過程を考える時は色々な手法が存在していますが、この記事では主にLinear Model of Coregionalization (LMC)というものを紹介し、実際にHello Worldしていきたいと思います。LMCはIntrinsic Coregionalization Model (ICM)というものとSemiparametric Latent Factor Model (SLFM)というモデルを拡張したものとなっており、この3つは全て同じ考え方の下で構築されます。以下、順を追って述べていきたいと思います。

3-1　前提

互いに相関のある系列ですが、手元のデータは横軸の値が全ての系列で共有されているとは限りません。具体的には、以下の図で示される様子です。

右側のデータはデータがサンプルされる横軸の値がすべて共通となっていますが、左側はそうはなっていません。このようなデータはそれぞれHeterotopic dataとIsotopic dataという名前で区別されています。HeterotopicとIsotopicの区別をされることで何が変わるのかと言えば、数式の記述が若干変わります。本記事ではどっちで考えても大丈夫です。加えて、それぞれの系列で同じ数だけデータが取れているとも限らないので、そこの仮定も頭に入れる必要があります。とは言え、とりあえずは簡単のために、Isotopic dataで全ての系列で同じ数のデータが取れたことを想定しておきます。

3-2　概要

・イメージ

ICM、SLFMと来て最終的にLMCを説明していくわけですが、この3つには共通点があります。それは、どのモデルも背後に潜在的なガウス過程を仮定する事です。図を見たほうが早いですね。

詳しい記号はすぐ後で説明しますので、こういう事をしてるんだなあ程度で大丈夫です。図の $\displaystyle f$ はモデル化したい過程であり、それが $\displaystyle g$ と書かれるものによって表されることを示しています。ひとまず $\displaystyle g$ が $\displaystyle f$ の上の方にいるんだなあ的な感じでOKです。

・式の展開

ICMとSLFMはなんか似ている事が分かりました。ここから少しだけ式を入れていきます。上で示した図は系列が3つで少し多い感じがするので、 $\displaystyle {\bf{f_1, f_2}}$ の2つに絞る事にします。さらに、2つの系列でそれぞれ2点が得られている事とします。つまり、 $\displaystyle {\bf{f_1}}=(f_1(x_1)　f_1(x_2))^T, {\bf{f_2}}=(f_2(x_2)　f_2(x_2))^T$ です。MOGPの冒頭でMOGPは力技であるという事を述べたので、力技をしてみましょう。

1つの多変量正規分布にまとめてみました。平均は今まで通り全て0にしておくことにします。さて、問題は分散共分散行列になる訳ですが、どのように計算されるのかを見ていく事にします。ひとまず分散共分散行列だけを取り除くと、以下のようになります。

対角成分は分散でそれ以外の成分は共分散である事が分かります。これも最初の方に述べた多変量正規分布の時と話は全く同じです。共分散は期待値を用いて、 $\displaystyle Cov(X, Y)=E[XY] - E[X][Y]$ と計算できる訳なので、1つ1つ地道に計算していけばMOGPの分散共分散行列も計算することができます。ここまではICMからLMCまですべて共通する話です。ここから先、 $\displaystyle f$ が潜在的なガウス過程 $\displaystyle g$ を用いてどのように構築されるかによって、ICMからLMCで分散共分散行列が変わってきます。

3-3　Intrinsic Coregionalization Model (ICM)

ICMは潜在的なガウス過程は1つで、そこからのサンプルによって $\displaystyle f$ を構成します。ここでは簡単のために、サンプルを1つにして次のように $\displaystyle f$ を構築します。

$\displaystyle f_1(x_1) = \omega_1^{(1)} g^{(1)}(x_1)$
$\displaystyle ...$
$\displaystyle f_2(x_2) = \omega_2^{(1)} g^{(1)}(x_2)$

ここで、 $\displaystyle g^{(1)}$ は潜在的なガウス過程 $\displaystyle g$ からのサンプルであるとします*6。この時、共分散を計算してみましょう。共分散の公式は上で示したのでそのまま代入します。例として、 $\displaystyle Cov( f_1(x_1), f_2(x_1))$ を計算してみます。

$\displaystyle Cov( f_1(x_1), f_2(x_1)) = E[f_1(x_1)f_2(x_1)] - E[f_1(x_1)]E[f_2(x_1)]$
$\displaystyle =\omega_1^{(1)}\omega_2^{(1)} ( E[ g^{(1)}(x_1)g^{(1)}(x_2) ] - E[g^{(1)}(x_1)]E[g^{(1)}(x_2)] ) = \omega_1^{(1)}\omega_2^{(1)} k(x_1, x_2)$

$\displaystyle E[f_1(x_1)]E[f_2(x_1)]$ の項は平均0を仮定しているので現れません。最後の式変形はそのまま共分散の形が出てきたので、ここをカーネルで計算させます。これを真面目にすべて計算すると次のような分散共分散行列が出てきます。

ここで、 $\displaystyle B$ は2行2列の係数の行列で、例えば $\displaystyle b_{11} = (\omega_1^{(1)})^2, b_{21} = \omega_1^{(1)}\omega_2^{(1)}, ...$ というようになっています。 $\displaystyle K$ は2行2列のカーネルを要素に持つ行列です。行列の中身がブロックに区分されているように見える訳ですが、このような物はクロネッカー積というものですっきり書く事が出来ます。最終的に導出されるものを先に見てみましょう。

最右辺の豆電球みたいな演算がクロネッカー積です。クロネッカー積の説明はとりあえず省略しておきます。適当にwikipediaとか貼っておきます。
ja.wikipedia.org

というわけで、力技で1つにまとめた分布をすっきりした形で記述することができました。

3-4　Semiparametric Latent Factor Model (SLFM)

SLFMもICMと同様に背後に潜在的なガウス過程を仮定します。しかし、ICMと違うのは潜在的なガウス過程を複数仮定することです。例えば、次のようにして $\displaystyle f$ を構築することとします。

$\displaystyle f_1(x_1) = \omega_{1, 1} g_1(x_1) + \omega_{1, 2} g_2(x_1)$
$\displaystyle ...$
$\displaystyle f_2(x_2) = \omega_{2, 1} g_1(x_2) + \omega_{2, 2} g_2(x_2)$

ここで、 $\displaystyle g_1, g_2$ はそれぞれ異なるガウス過程です。ICMで計算したのと同じように、共分散を計算してみましょう。例えば、 $\displaystyle Cov( f_1(x_1), f_2(x_1))$ を計算してみましょうか。

$\displaystyle Cov( f_1(x_1), f_2(x_1)) = E[f_1(x_1)f_2(x_1)] - E[f_1(x_1)]E[f_2(x_1)]$
$\displaystyle = \omega_{1, 1}\omega_{2, 1} Cov( g_1(x_1), g_1(x_1) ) + \omega_{1, 2}\omega_{2, 2} Cov( g_2(x_1), g_2(x_1) )$
$\displaystyle = \omega_{1, 1}\omega_{2, 1} k_1(x_1, x_1) + \omega_{1, 2}\omega_{2, 2}k_2(x_1, x_1)$

計算の過程で $\displaystyle E[f_1(x_1)]E[f_2(x_1)]$ の項は平均0を仮定しているので現れません。加えて、 $\displaystyle Cov( g_1(x_1), g_1(x_1) ), Cov( g_2(x_1), g_2(x_1) )$ はカーネルそのものです。以後、こんな形で頑張って計算すると、最終的には以下のようになります。

ICMと比べると別の分散共分散行列が出てきましたね。

3-5　Linear Model of Coregionalization (LMC)

最後にLMCです。まずは以下の図を見てみましょう。

LMCはICMとSLFMを合わせたものとなります。SLFMでまずは潜在的なガウス過程の個数を決め、各ガウス過程から複数のサンプルを取得する…といったような感じです（複数のサンプルを取得する記述は参考文献[3]で記載されていますが、それ以外で僕の見た中でLMCを説明している論文のサンプル数は1です。この辺りどうなんですかね？）。まずは、これまでの流れ通りに $\displaystyle g$ を用いて $\displaystyle f$ を構成してみましょう。ここはいきなり一般化された定義から見ていきます。

$\displaystyle f_n(x) = \sum_{q=1}^{Q} \sum_{r=1}^{r_q} \omega_{n, q}^{(r)} g_{q}^{(r)} (x)$

ここから共分散を計算していくわけですが、1から計算する必要はありません。すでにICMとSLFMで分散共分散行列の話は済ませてあるので、それを利用して求めてしまいます。LMCでは最終的に以下のような結果となります。

あれれ？ICMでいっぱいサンプルしたとすればそれはどこで反映されるの？という感じに見えますが大丈夫です。1つのガウス過程から2つ以上サンプルしたとしても $\displaystyle B⊗K$ の形で書くことができます。こんな感じでLMCはあっさりと説明を終えることができました。

4．MOGPのHello World

ここまで非常に長くなってしまいましたが、いよいよ実験してみたいと思います。とりあえずはプログラムを動かして結果を図示してみる所で終わりたいと思います。データはMOGPの冒頭に示したグラフ……ですが、ページ戻るのも大変なので、もう一度貼っておきましょうか。

MOGPはGPyを利用して動かしました。GPyを選んだ理由は特にありません。データセットの作成を含めたプログラムの全体は1番最後に貼っておきます。ひとまず核となるプログラムだけ載せる事とします。いよいよLCMにHello Worldです。まずは以下のコードを動かしてみます。

K1 = GPy.kern.Bias(1)
K2 = GPy.kern.Linear(1)
lcm = GPy.util.multioutput.LCM(input_dim=1, num_outputs=3, kernels_list=[K1, K2])
m = GPy.models.GPCoregionalizedRegression([X1, X2, X3], [Y1, Y2, Y3], kernel=lcm)
m.optimize()

初めてBiasという名前のカーネルを聞いたんですけど、こんなのあったんですね。Biasという名前のカーネルだけだと平均辺りで横一直線しか引っ張られなかったので、Linearも加えてみました。結果は以下のようになりました。

右肩上がりや右肩下がりの傾向は掴んでいます。次はシンプルにガウシアンカーネルだけにしてみましょう。

K = GPy.kern.RBF(1)
lcm = GPy.util.multioutput.LCM(input_dim=1, num_outputs=3, kernels_list=[K])
m = GPy.models.GPCoregionalizedRegression([X1, X2, X3], [Y1, Y2, Y3], kernel=lcm)
m.optimize()

良い感じのぐにゃぐにゃ具合ですね。1つ1つの系列でSOGPを掛けるよりも良い気がします。

5．終わりに

この記事ではSOGPとMOGPについて述べ、MOGPの中でもLMCに着目してHello Worldしてみました。系列全体の相関を考慮してモデル化できるのはすごく魅力的でした。しかも、全部1つにまとめる力技で動いているとは、、、多変量正規分布さまさまですね。ちなみに、この記事では以下の内容をごっそり抜かしました。

(前準備) 2つの確率変数ベクトルの同時確率密度関数と条件付確率密度関数
(前準備) ベイズの定理
(SOGP) ノイズを含んだガウス過程
(SOGP) 予測分布の導出
(SOGP) 各種推論（変分ベイズ、期待値伝播法、MCMC…）
(MOGP) MOGPの根本的な理論

などなど……

あくまでガウス過程の雰囲気をちょこっと掴もう！というHello World目的で記事を書いたのでとりあえず良いかな？という感じがします。なんとなく雰囲気が掴めた所で、より具体的な理論を学びたいと思った時には、何かしらの書籍や論文で上の内容を学ぶと良いと思います。僕も現在進行形で勉強中です！以上、Advent Calendarの23日目でした。皆様良いクリスマスをお過ごしください！それでは！

参考文献

[1] M. van der Wilk et al., "A Framework for Interdomain and Multioutput Gaussian Processes", 2020.
[2] H. Liu, J. Cai and Y-S. Ong, "Remarks on Multi-Output Gaussian Process Remarks", 2017.
[3] M. A. A´ lvarez, L. Rosasco and N. D. Lawrence, "Kernels for Vector-Valued Functions: a Review", 2012.
[4] 持橋大地，大羽成征，”ガウス過程と機械学習"．
[5] C. M. Bishop, ”パターン認識と機械学習（上下）”．

※GPyを用いて遊んでみましたが、以下を参考にプログラムを作成しました。
nbviewer.jupyter.org
github.com

プログラム

# seibisi(2歳)
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import GPy

def makedata(amp, k, phi, N):
    X = np.random.rand(N)[:,None]
    Y = amp * np.sin(X/k+phi) + np.random.rand(X.size)[:,None]
    return X, Y

def plotall(m):
    fig = plt.figure(figsize=(15, 12))
    # output 1
    ax1 = fig.add_subplot(311)
    m.plot(fixed_inputs=[(1,0)], which_data_rows=slice(0, 50), ax=ax1)
    ax1.set_title("f1(x)")
    # output 2
    ax2 = fig.add_subplot(312)
    m.plot(fixed_inputs=[(1,1)], which_data_rows=slice(50, 100), ax=ax2)
    ax2.set_title("f2(x)")
    # output 3
    ax3 = fig.add_subplot(313)
    m.plot(fixed_inputs=[(1,2)], which_data_rows=slice(100, 150), ax=ax3)
    ax3.set_title("f3(x)")

def main():   
    # 1. データセットの作成
    X1, Y1 = makedata(1, 0.2, 0, 50)
    X2, Y2 = makedata(2, 0.3, 2, 50)
    X3, Y3 = makedata(0.8, 0.15, 1.5, 50)
    # 2. 当てはめ
    K = GPy.kern.RBF(1)
    lcm = GPy.util.multioutput.LCM(input_dim=1, num_outputs=3, kernels_list=[K])
    m = GPy.models.GPCoregionalizedRegression([X1, X2, X3], [Y1, Y2, Y3], kernel=lcm)
    m.optimize()
    # 3. プロット
    plotall(m)
    
if __name__ == "__main__":
    main()

*1:去年も僕はAdvent Calendarを書きました（カーネルロジスティック回帰を自前で実装して遊んだ）

*2:日本語の情報も少ない？良い日本語訳をご存知の方は教えて下さい！

*3:というよりそのもの

*4:カーネルが要素の行列はグラム行列とか名前がついていますが、しばらくは分散共分散行列という言葉のまま進めていきたいと思います。

*5:系列の1つ1つが多変量正規分布に従うならば、それを結合したものもまた多変量正規分布に従うだろうというロジックです。

*6:複数考える場合は独立なサンプルとして扱います

遊具

たくさんのおもちゃ

Multioutput Gaussian Processで遊ぶ（Fun Advent Calendar 2020 #23）

0．はじめに

1．前準備

2．Single-output Gaussian Process (SOGP)

2-1　SOGPの雰囲気

2-2　SOGPの定義

2-3　SOGPの分散共分散行列

2-4　SOGPまとめ

3．Multi-output Gaussian Process (MOGP)

3-1　前提

3-2　概要

・イメージ

・式の展開

3-3　Intrinsic Coregionalization Model (ICM)

3-4　Semiparametric Latent Factor Model (SLFM)

3-5　Linear Model of Coregionalization (LMC)

4．MOGPのHello World

5．終わりに

参考文献

プログラム

0．はじめに

1．前準備

2．Single-output Gaussian Process (SOGP)

2-1 SOGPの雰囲気

2-2 SOGPの定義

2-3 SOGPの分散共分散行列

2-4 SOGPまとめ

3．Multi-output Gaussian Process (MOGP)

3-1 前提

3-2 概要

・イメージ

・式の展開

3-3 Intrinsic Coregionalization Model (ICM)

3-4 Semiparametric Latent Factor Model (SLFM)

3-5 Linear Model of Coregionalization (LMC)

4．MOGPのHello World

5．終わりに

参考文献

プログラム

2-1　SOGPの雰囲気

2-2　SOGPの定義

2-3　SOGPの分散共分散行列

2-4　SOGPまとめ

3-1　前提

3-2　概要

3-3　Intrinsic Coregionalization Model (ICM)

3-4　Semiparametric Latent Factor Model (SLFM)

3-5　Linear Model of Coregionalization (LMC)