雑記: モデルをアンサンブルしたい話(その2―カステラ本7.11節、8.2節、8.4節、8.8節、10.1~10.4節)

私の誤りは私に帰属します。お気付きの点がありましたらお手数ですがご指摘いただけますと幸いです。
テキスト(カステラ本)関連記事
  1. 雑記: モデルを選択したい話(カステラ本7.4節~7.7節) - クッキーの日記
  2. 雑記: モデルをアンサンブルしたい話(その1―カステラ本7.3節、8.7節) - クッキーの日記(前回)
まとめ(解釈を含む)
  • 並行世界の自分が学習したモデルを取り寄せればバリアンスを小さくできる。が、そんなことはできないので自分で並行世界をつくり出す方法がブートストラップ法である(そのモデルたちを平均するのがバギングである)。
  • カテゴリ値を取るデータの標本 Z からのノンパラメトリックブートストラップ標本は、「標本 Z に対して、Z 内の各カテゴリの割合を生成したディリクレ分布をベイズ推定して、各カテゴリの割合がその事後分布に(ほぼ)したがうように標本をサンプリングする」ことをしたものであるといえる。ことから、バギングしたモデルは近似的にモデルのベイズ事後分布の平均であるとみなせる(ので普通に訓練したモデル=MAP推定より2乗誤差が小さくなることが見込める)。
  • もっともベイズ事後平均を取るのであれば、異なるモデルたちを候補にしてもよい。そのときモデルを足し合わせる重みはBIC(「どのモデルがこの訓練データに対して最も尤もらしいか」という指標だった)を用いてもよいし、直接最適な重みを推定してもよい(スタッキング)。
  • そもそも足し合わせる各モデルが同じ標本を学ばなくてもよい。1つ前のモデルが誤分類したデータを重点的に学んでいくのがブースティング法の最も基本的なアルゴリズムであるアダブーストM1である。アダブーストM1は前向き段階的加法的モデリングにおいて指数損失を採用しても導出される。というより、前向き段階的加法的モデリングを解く手法全般をブースティングという。
キャラクターの原作とは無関係です。
f:id:cookie-box:20180305232608p:plain:w60

前回は「並行世界の自分が学習したモデルを取り寄せればバリアンスを小さくできる。が、そんなことはできないので自分で並行世界をつくり出す方法がブートストラップ法である」というところまででしたね。

f:id:cookie-box:20180305231302p:plain:w60

だいぶブートストラップ法に対する語弊がないかそれ…。そもそもブートストラップ法って何? 自分で並行世界をつくるなんてことしていいの?

f:id:cookie-box:20180305232608p:plain:w60

テキストでブートストラップ法が最初に出てくるのは286ページの7.11節ですね。前回もいった通りここでのブートストラップ標本とは「訓練データ Z = \bigl\{ (x_1, y_1), \cdots, (x_N, y_N)\bigr\} から N 点を復元抽出する」ことを B 回繰り返して B セットの新しい標本 Z^{\ast 1}, Z^{\ast 2}, \cdots, Z^{\ast B} を得たものです。復元抽出なので元の訓練データの点はある Z^{\ast b} に複数回抽出されていることもありますし、1回も抽出されていないこともあります。カギカッコ内は「訓練データの経験分布から N 点生成する」といっても同じですね。なお、これはあくまでノンパラメトリックブートストラップといわれる方法であって、ブートストラップ標本を得る方法は他にも色々あります(テキスト302ページのノイズを添加する方法など)。より一般に、訓練データ自体を用いて並行世界たちの訓練データを複製する方法であればブートストラップ法のようです。そもそもブートストラップとは靴の後部の小さなつまみのことで、"pull oneself up by one's bootstraps" という慣用句で「他人の助けを借りずに自力で達成する」という意味であるようです。ので、この訓練データだけから新たな標本を生成する方法がブートストラップ法と名付けられたのでしょう。初出は1979年とか。

そういう手法なので、期待値や分散や信頼区間モンテカルロ推定に利用されるのがメインなのではないでしょうか。7.11節に出てくるブートストラップ法の用法は以下ですね。汎化誤差の期待値の推定は込み入っています。ブートストラップ標本は元の標本の点を全ては含んでいませんから、そのディスアドバンテージ分をどう補ってやるかということになってくるようですね。
  • 訓練データから計算される何らかの量 S(Z) の分散を推定するのに、B 個のブートストラップ標本上の不偏分散でもって推定値とする。
  • 汎化誤差を推定するのに、B 個のブートストラップ標本で学習したモデル \hat{f}^{\ast 1}, \cdots, \hat{f}^{\ast B} による平均予測誤差の、全訓練データの平均を推定値 \widehat{\rm Err}^{(1)} とする。なお、各 x_i の平均予測誤差を出す際、点 x_i を訓練データとして参照したモデルはとばす。また、この推定値には亜種がある(以下)。
    • 上の推定値では、「個々の \hat{f}^{\ast 1}, \cdots, \hat{f}^{\ast B} は本来の約63.2%のデータしか学べていない(ので汎化誤差を過大評価しているかもしれない)」という問題点がある。これを緩和するために、訓練誤差を36.8%混ぜる(このようになる導出はあるらしい)。これを \widehat{\rm Err}^{(0.632)} とする。
    • 上の推定値でも、「過学習のために訓練誤差が小さくなりすぎているかもしれない(ので訓練誤差を混ぜると逆に過小評価になるかもしれない)」という問題点がある。そこで、「相対過学習率」を、「 \widehat{\rm Err}^{(1)} から訓練誤差を引いたもの」を「でたらめな正解ラベルを付けたときのそのモデルの予測誤差から訓練誤差を引いたもの」で割ったものとして定義する。「相対過学習率」が大きいほど \widehat{\rm Err}^{(1)} を信用する。これを \widehat{\rm Err}^{(0.632+)} とする。
ブートストラップ標本は元の標本の点を全ては含んでいないというのに関連して、元の標本が N 点であったときにあるブートストラップ標本にある点が選ばれる確率は以下ですね。

def arutenga_erabareru_kakuritsu(N):
    p = 1.0 - 1.0 / float(N)  # 1回の抽出で自分以外の点が選ばれる
    p_total = 1.0  # N回の抽出すべてで自分以外の点が選ばれる
    for i in range(N):
        p_total *= p
    return 1.0 - p_total

ある点が選ばれない確率は N = 10N = 20.25 で、N \to \infty の極限でネイピア数の逆数になります。N' = -N とおけば  (1-1/N)^N = [ (1+1/{N'})^{-N'} ]^{-1} \xrightarrow[-N' \to - \infty ]{} e^{-1} ですから。
ただ7.11節にブートストラップ標本によって推定することがどのように正当かという記述はありませんね。そのような話題は8章にあるようです。以下にまとめます。文字の定義はテキストを参照してください。

  1. 真の値のノイズが加法的でガウシアンなときに、パラメトリックブートストラップ標本で学習したモデルの期待値が最尤推定の結果と一致する(8.2節)。
  2. ガウス分布 N(\theta, 1) の期待値 \thetaベイズ推定する状況( \theta の事前分布も事後分布もガウス分布)では、事前分布の分散を無限大にすると(無情報事前分布)、ある点 z を観測したときの事後分布が N(z,1) となるが、パラメトリックブートストラップ標本はこの事後分布から1点 \theta' をサンプリングして N(\theta', 1) を得てそれの期待値を採用したものとみなせる(8.4節)。
    • つまり、パラメトリックブートストラップ標本は、「各点 z に対して、z を生み出したのは N(z,1) なんだろうなとベイズ推測して、その分布から点をサンプリングする」ことを各点に対してしたものであるといえる。
  3. カテゴリ数が L のカテゴリカルなデータの標本を得たとする。各カテゴリが生成される割合 w の事前分布をすべてのカテゴリの集中度母数が0の極限のディリクレ分布とすると、事後分布はすべてのカテゴリの集中度母数が標本内での観測数であるディリクレ分布となるが、ノンパラメトリックブートストラップ標本内の各カテゴリの割合がしたがう分布(標本内の観測割合を母数にもつ多項分布)はこれに近い(平均が等しく分散共分散がほぼ等しい)(8.4節)。
    • つまり、ノンパラメトリックブートストラップ標本は、「標本 Z に対して、Z 内の各カテゴリの割合を生み出したのは {\rm Dir} (N \hat{w}) なんだろうなとベイズ推測して、各カテゴリの割合がその分布にほぼ近くなる方法で標本をサンプリングする」ことをしたものであるといえる。
文字に色を付けた部分の気持ちになれば、天下りでなく演繹的にブートストラップ法にたどり着けるような気がするかもしれません。僕はしませんが。

f:id:cookie-box:20180305231302p:plain:w60

しないのかよ!? あと、「ほぼ等しい」ってなんか中途半端じゃない? どれくらい等しいの??

f:id:cookie-box:20180305232608p:plain:w60

テキストにそうかいてあるものはかいてありますから…まあ確かめますか。{\rm Dir} (N \hat{w}) にしたがう \tilde{w} の分散は \hat{w}_i (1 - \hat{w}_i) / (N + 1) で、N 倍が  {\rm Multi}(N, \hat{w}) にしたがう \tilde{w} の分散は \hat{w}_i (1 - \hat{w}_i) / N ですから近いですね。共分散は -\hat{w}_i \hat{w}_j / (N + 1)-\hat{w}_i \hat{w}_j / N なので近いです。

f:id:cookie-box:20180305231302p:plain:w60

確かに近いな。じゃあこの状況でのブートストラップ法はほぼベイズ推定なんだ。

f:id:cookie-box:20180305232608p:plain:w60

はい、テキスト8.8節「モデルの平均と統合」の冒頭にも「バギングしたモデルは近似的なベイズ事後平均になっていることがわかる」のような文があります。対して、普通に訓練したモデルはベイズ事後分布の最頻値をとっていることに対応すると。そして、2乗誤差に対して最適なのはベイズ事後平均です。ということからも、バギングすることで2乗誤差を小さくできることが見込まれます。これが一応の「自分で並行世界をつくるなんてことをしていいのか」に対する答えになるのではないかと。

f:id:cookie-box:20180305231302p:plain:w60

そういわれるとそうなのか。

f:id:cookie-box:20180305232608p:plain:w60

さらに8.8節では「平均する個々のモデルを異なるモデルにしたらどうだろうか」という議論を巡らせています。どう思いますか?

f:id:cookie-box:20180305231302p:plain:w60

えっ、そんな、別々のモデルにしちゃったら、モデルたちはベイズ事後分布をなさないよな…「並行世界の自分にモデルをもらったけど並行世界の自分が学習したのは違うモデルだった」ってなったら、自分のモデルとどっちが尤もらしいかってわからないわけだし…それじゃ平均したら駄目なんじゃないの?

f:id:cookie-box:20180305232608p:plain:w60

ええ、どちらが尤もらしいかわかりませんね。でも、前々回に、「どのモデルがこの訓練データに対して最も尤もらしい?」という指標を取り扱いませんでしたか?

f:id:cookie-box:20180305231302p:plain:w60

あ、BIC…。

f:id:cookie-box:20180305232608p:plain:w60

はい、そのような場合は BIC が利用できるとあります。また、より直接的に、「どのモデルにどれだけ重みを割り振れば誤差が最も小さくなる?」と重みを最適化するのが331ページ下部で紹介されているスタッキングです。ここで、モデル mi 番目のデータに対する誤差の計算時には i 番目のデータのみを除いた訓練データで学習したモデル m を用意し、それで予測します。明確に記述されていませんが、データ数×モデル数だけ学習することになるので、結構面倒そうな気はします。

それで、バギングの重要例であるランダムフォレストは15章で掘り下げられています(ただ、ランダムフォレストでは決定木の分岐点を決めるときに特徴量もランダムに絞り込んでさらにバリアンスを下げる工夫をしているのでバギングというのかバギングの改造版というのかわかりませんが)。

さておき、ここまでで、「モデルを平均するとよい」「異なるモデルを異なる重みで平均してもよい」というところまできたわけです。であれば、その重みは入力空間内の点によっても異なっていいと思いませんか? いい換えると、ある場所ではモデル1が重く、ある場所ではモデル2が重いといったように、あたかも得意不得意で役割分担するイメージです。

f:id:cookie-box:20180305231302p:plain:w60

その点で予測誤差が小さいモデルを重くするってこと? さすがにそれは後出しじゃんけんじゃない? というか、その点のノイズに引っ張られまくる気が。

f:id:cookie-box:20180305232608p:plain:w60

いえ、ブートストラップ標本で学習したモデルでそれをやるということではないです。あるモデルはある場所から重点的にサンプリングした訓練データで学習し、また別のあるモデルは別の場所から重点的にサンプリングした訓練データで学習し…といったイメージです。それも入力空間を適当に区切るわけではなく、まだ上手く学習できていない箇所に重点を置きにいくんです。もはや個々のモデルの訓練データが異なるということです。…もうネタばらししてしまうんですが、10章のブースティング法がその手法です。その最も基本的なアルゴリズムである、アダブーストM1(10.1節)が以下です。ゴールは2クラス分類問題を解くことで、各分類器 G_m(x)-11 を出力するという状況です。

  1. 各データの重み \{w_i\}_{i=1}^N を均等に初期化する。
  2. 重み \{w_i\}_{i=1}^N で分類器 G_1(x) を学習する。
  3. G_1(x) の重み付き誤分類率 {\rm err}_1 を計算し、正解しやすさの対数オッズ \alpha_1 = \log \bigl( (1 - {\rm err}_1)/{\rm err}_1 \bigr) を計算する。
  4. 各データの重みを更新する。 G_1(x) が正しく予測できたデータの重みは更新しない。 G_1(x) が正しく予測できなかったデータの重みには  e^{\alpha_1} をかける。
  5. 2. に戻る(更新した重み \{w_i\}_{i=1}^N で分類器 G_2(x) を学習する)。
  6. 最終的な予測モデルを G(x) = {\rm sign} \bigl[ \sum_{m=1}^M \alpha_m G_m(x) \bigr] とする。
「正しく予測できなかったデータの重みを増幅しておく(ことで次以降のモデルに託す)」ことをやっているのがわかると思います。また、重みの増幅率やモデルを足し合わせる係数に対数オッズが出てくるのにも理由があります。

f:id:cookie-box:20180305231302p:plain:w60

強引に10章にとんだな…。

f:id:cookie-box:20180305232608p:plain:w60

ブースティングの話に触れておきたかったので。モデルのアンサンブルという文脈でバギングとブースティングが対比される場合が多いと思うんです。というよりは、ランダムフォレストと勾配ブースティング木が対比されるんでしょうか…しかし、テキスト10章の冒頭にもあるように両者は「本質的に異なる手法」なんです。どちらも決定木がたくさんあるくらいで。なので話を上手く10章にもっていけませんでした。なんていうか、一度「モデルを組み合わせよう」という気持ちを捨てた方がいい気がしますね。

f:id:cookie-box:20180305231302p:plain:w60

うん、いきなりアダブーストM1っていわれてもブートストラップ標本が出てきたときより突拍子がないんだけど、ブートストラップ標本がベイズ推定から導出されるなら、そのアダブーストM1は何から導出されるの。

f:id:cookie-box:20180305232608p:plain:w60

それが10.3節の「前向き段階的加法的モデリング」ですね。いま学習したいモデルが10.3式のように他のモデルの線形和でかけるとします。便宜上、足し合わせるモデルを弱学習器とよび、足し合わせたモデルを強学習器とよびます。10.3式のようなモデルは全てのパラメータを一気に最適化すべきですが、計算コストなどの問題でそれが上手くできないとき、弱学習器を1つずつ最適化できないかということになります。それが前向き段階的加法的モデリングに他なりません。391ページのアルゴリズム10.2です。

f:id:cookie-box:20180305231302p:plain:w60

391ページのアルゴリズム10.2ね…って、387ページのアダブーストM1とだいぶ違わない? 重みとか出てこないし。それに、アダブーストM1はステップ m で学習した弱学習器の損失をみていたけど、この式だとステップ m までの強学習器の損失って感じだし。

f:id:cookie-box:20180305232608p:plain:w60

いえ、アルゴリズム10.2で L を指数損失としてみましょう。そうする各ステップですべき最適化が10.9式になります。「ステップ m までの強学習器の損失」が「重み」と「ステップ m で学習した弱学習器の損失」の積になっていますよね。むしろ「重み付き誤分類率」とは「ステップ m までの強学習器の損失」だったんです。だから重み付き誤分類率を最小化する弱学習器 G_m(x) を学習すればいいんです。

f:id:cookie-box:20180305231302p:plain:w60

なるほど。あれでも、弱学習器 G_m(x) をどれだけの重みで強学習器に足せばいいんだ…って、10.12式か。え、これどうやって出すの?

f:id:cookie-box:20180305232608p:plain:w60

演習10.1になっていますね…そうですね、以下で出ます。

  •  \sum_{i=1}^N w_i^{(m)} \exp \bigl( - \beta_m y_i G_m(x_i) \bigr) が最も小さくなればいいんですよね。
  •  \sum_{i \in {\rm OK}_m} w_i^{(m)} \exp \bigl( - \beta_m \bigr) + \sum_{i \in {\rm NG}_m} w_i^{(m)} \exp \bigl( \beta_m \bigr) と書き換えます。{\rm OK}_mG_m(x_i) = y_i であるインデックスの集合、{\rm NG}_mG_m(x_i) \neq y_i であるインデックスの集合とします。
  • \beta_m についての微分がゼロなら  - \sum_{i \in {\rm OK}_m} w_i^{(m)} \exp \bigl( - \beta_m \bigr) + \sum_{i \in {\rm NG}_m} w_i^{(m)} \exp \bigl( \beta_m \bigr) = 0 です。
  • よって  \exp \bigl( 2 \beta_m \bigr) \sum_{i \in {\rm NG}_m} w_i^{(m)} = \sum_{i \in {\rm OK}_m} w_i^{(m)} です。
  • よって  \displaystyle \beta_m  =\frac{1}{2} \log \frac{ \sum_{i \in {\rm OK}_m} w_i^{(m)}}{\sum_{i \in {\rm NG}_m} w_i^{(m)}} =\frac{1}{2} \log \frac{ 1 - {\rm err}_m}{{\rm err}_m} です。

(その3があれば)つづく

雑記: モデルをアンサンブルしたい話(その1―カステラ本7.3節、8.7節)

私の誤りは私に帰属します。お気付きの点がありましたらお手数ですがご指摘いただけますと幸いです。
テキスト(カステラ本)関連記事まとめ(解釈を含む;文字の定義は記事内を参照)
  • 手元の訓練データでモデルを学習しある未知のデータ x_0 に対する2乗誤差を小さくしたいとする。手元の訓練データも x_0 上の真の値も確率的なので、誤差 f(x_0) + \varepsilon - \hat{f}(x_0) は確率変数であり、以下の3つの和で表せる。
    1. x_0 上の真の値がぶれる分(真の値のノイズ成分): \varepsilon
    2. 訓練データのぶれによりモデル予測値が期待値に満たない分:  E\bigl[ \hat{f}(x_0)\bigr] - \hat{f}(x_0)
    3. モデル予測値が期待値であったとしても真の値に満たない分:  f(x_0) - E\bigl[ \hat{f}(x_0)\bigr](これは確率的でない)
  • 上の 1., 2., 3. の和の2乗を小さくしたい。が、確率的な成分があると議論しにくいので確率的な成分について期待値をとることにする。と、クロスタームは消え、1., 2., 3. の2乗の期待値の和だけが残る。このうち 1. の2乗はモデル側で小さくできないので、モデル側で小さくしうる項として 2. の2乗の期待値と 3. の2乗が残る。この前者をバリアンス、後者をバイアスの2乗とよぶ。
  • よって、期待2乗誤差を小さくするにはバリアンスとバイアスを小さくすればよいが、典型的なモデルで両者はトレードオフの関係にあり、一方を小さくすれば他方が大きくなるとわかる(具体的にはリッジ線形回帰など)。よって、手元の訓練データでモデルを学習する分にはこれらのバランスを取るように適当に正則化などするしかない。
  • ここで仮にもし、並行世界(この世界とは異なる訓練データが得られた世界)たちの自分から学習済みのモデルを取り寄せることができるなら、集めたモデルたちの出力の平均を取って新たなモデルとすることでバリアンスを抑えることができる(※ 複数のモデルの出力の平均を新たな出力にできる場合)。
  • しかし、並行世界たちの自分から学習済みのモデルを取り寄せることはできない。
  • そこで、手元で仮想的な並行世界たちをつくり、そこで学習したモデルたちを平均することが考えれる。特に、訓練データからブートストラップ標本たちを生成することによりそれを達成する手法はバギング(bagging: bootstrap aggregating)といわれる。
キャラクターの原作とは無関係です。
f:id:cookie-box:20180305231302p:plain:w60

そういえば、カステラ本の7章でとばした箇所にあった、「バイアスと分散のトレードオフ」って何?

f:id:cookie-box:20180305232608p:plain:w60

以下の問題を考えてみてください。どうなりますか?

  • ある世界では、気温 X の日のアイスクリームの売り上げ YY = f(X) + \varepsilon になるとします。
    • f(\cdot) は決定的な関数で、\varepsilon は確率的な誤差で、{\rm E}(\varepsilon)=0, \; {\rm Var}(\varepsilon)=\sigma^2 とします。
  • さて、この世界で、ある日の気温が x_0 でした。この日の売り上げを予測モデル \hat{f}(\cdot) による予測値 \hat{f}(x_0) で予測したときの期待2乗誤差はいくらでしょうか? なお、予測モデル \hat{f}(\cdot) の確率的な成分は \varepsilon と独立とします。

f:id:cookie-box:20180305231302p:plain:w60

何その世界…ともかく、気温が x_0 の日の真の売り上げは f(x_0) + \varepsilon だよな。それで、いま確率的なのは、\varepsilon と…\hat{f}(x_0) もなのか。とりあえずこれらについての期待値を単に  E [ \cdot ] とかくと、期待2乗誤差は、

\begin{split} E \Bigl[ \bigl( f(x_0) + \varepsilon - \hat{f}(x_0) \bigr)^2 \Bigr] &=  E \Bigl[ f(x_0)^2 + \varepsilon^2 + \hat{f}(x_0)^2 + 2 f(x_0) \varepsilon - 2 f(x_0)\hat{f}(x_0) - 2 \hat{f}(x_0) \varepsilon \Bigr] \\ &= f(x_0)^2 + \sigma^2 + E \Bigl[ \hat{f}(x_0)^2 \Bigr] - 2 f(x_0) E \Bigl[ \hat{f}(x_0) \Bigr]\end{split}

こうなるよな。f(x_0) は確率的じゃないし、独立な確率変数の積の期待値は期待値の積だし。

f:id:cookie-box:20180305232608p:plain:w60

はい、間違っていませんが、それだと解釈しづらいですよね。「真の値の2乗とノイズの分散と予測値の2乗の期待値から真の値と予測値の期待値の積を引いたものです」といわれても。なのでこうしましょう。

\begin{split} E \Bigl[ \bigl( f(x_0) + \varepsilon - \hat{f}(x_0) \bigr)^2 \Bigr] &= f(x_0)^2 + \sigma^2 + E \Bigl[ \hat{f}(x_0)^2 \Bigr] - 2 f(x_0) E \Bigl[ \hat{f}(x_0) \Bigr] \! - \! E \Bigl[ \hat{f}(x_0) \Bigr]^2 \! \! \! \! + \!  E \Bigl[ \hat{f}(x_0) \Bigr]^2 \! \! \\ &= f(x_0)^2 + \sigma^2 + V \Bigl[ \hat{f}(x_0)^2 \Bigr] - 2 f(x_0) E \Bigl[ \hat{f}(x_0) \Bigr] \! + \! E \Bigl[ \hat{f}(x_0) \Bigr]^2 \\ &= \biggl\{f(x_0) -  E \Bigl[ \hat{f}(x_0) \Bigr] \biggr\}^2 + V \Bigl[ \hat{f}(x_0)^2 \Bigr] + \sigma^2 \\ &= E \Bigl[ f(x_0) - \hat{f}(x_0) \Bigr]^2 + V \Bigl[ \hat{f}(x_0)^2 \Bigr] + \sigma^2 \end{split}

こうすると「誤差の期待値(バイアス)の2乗と、予測値の分散(バリアンス)と、ノイズの分散の和」となりますよね。最初から予測値を期待値からのずれにしておけば一発ですけどね(以下)。

\begin{split} E \biggl[ \Bigl( f(x_0) + \varepsilon - \hat{f}(x_0) + E \Bigl[ \hat{f}(x_0) \Bigr] - E \Bigl[ \hat{f}(x_0) \Bigr] \Bigr)^2 \biggr] &= \biggl\{f(x_0) -  E \Bigl[ \hat{f}(x_0) \Bigr] \biggr\}^2 + V \Bigl[ \hat{f}(x_0)^2 \Bigr] + \sigma^2 \\ &= E \Bigl[ f(x_0) - \hat{f}(x_0) \Bigr]^2 + V \Bigl[ \hat{f}(x_0)^2 \Bigr] + \sigma^2  \end{split}

何にせよ、2乗誤差を下図のように区分けして各マスの期待値をとっているだけです。

f:id:cookie-box:20201228135133p:plain:w500

f:id:cookie-box:20180305231302p:plain:w60

俺に計算させた意味! あと図がごちゃごちゃしすぎ!! そもそもさ、期待2乗誤差が「誤差の期待値(バイアス)の2乗と、予測値の分散(バリアンス)と、ノイズの分散の和」だったら何なの? 何かうれしいの?

f:id:cookie-box:20180305232608p:plain:w60

「期待2乗誤差を小さくしたいなら、バイアスとバリアンスが小さいモデルにすべきである」といえるでしょう。ノイズはコントロールできませんから。

f:id:cookie-box:20180305231302p:plain:w60

ああそういうことか。「ある点に対する誤差の期待値も小さいし、予測値自体もぶれないモデル」にすべきってことなんだな…って、いまいちどうすればいいかわからないんだけど?

f:id:cookie-box:20180305232608p:plain:w60

実際にバイアスやバリアンスを小さくするにはどうすればいいのだろうと考えてみましょう。一般的な傾向として、トレードオフがありそうなことに気付くでしょう。

  • まず、バイアス=誤差の期待値を小さくするためには、手元の訓練データ上での誤差を徹底的に小さくするべきです。モデルをいくら複雑にしてもです。だって、精度を上げる余地があるのに手加減してしまったら、その分の誤差を詰められませんから。
  • しかし、バリアンス=予測値のぶれを小さくするには、手元の訓練データを徹底的に学ぶのは得策ではありません。なぜなら、訓練データはたまたま出たノイズを含みます。そのノイズまでしっかり学習してしまうようなモデルは、ノイズの出方によって予測値がぶれてしまいます。ぶれを抑えるには、あえて細かく学習させない工夫が要るでしょう。学習対象パラメータを制約してしまうのがその最たる例です。
実際、カステラ本の256~257ページに載っている、具体的なモデルでのトレードオフの例が以下の表の最初の3つです。k 近傍法では k の大小によってトレードオフが発生し、最小2乗線形回帰ではリッジ正則化をかける大きさ \lambda によってトレードオフが発生すると考えられます。k, \lambda を大きくするほどモデルが滑らかになり、小さくするほどモデルがでこぼこになるのが想像できると思います。また、以下の表の最後2つはウィキペディアにあった例です。
モデル \hat{f} バイアス(誤差の期待値) バリアンス(予測値の分散)
手元の訓練データから予測対象点 x_0 の最近傍 k 点を選び、それらの点の実績値の平均値を予測値とするモデル。 k を大きくするほど x_0 から離れた点が混ざってきてバイアスは大きくなる傾向にある。 k を大きくするほど多くの点の実績値の平均を取るのでバリアンスは小さくなる傾向にある。
手元の訓練データで最小2乗線形回帰したモデル。つまり、\hat{f}(x_0) = x_0^\top (X^\top X)^{-1} X^\top Y とする。横ベクトル h(x_0)=x_0^\top (X^\top X)^{-1} X^\top の各成分は、各データの実績値をどれだけの重みで取り込むべきかを意味する。 バイアスは「真のモデルと最良の線形近似との誤差の期待値」と「最良の線形近似といまの線形近似との誤差の期待値」に分解できる。後者は、何度も訓練データを取って学習したモデルの期待値をとれば、ゼロになる。 \|h(x_0)\|^2 \sigma^2 になる。
1つ上のモデルをリッジ正則化する。つまり、\hat{f}(x_0) = x_0^\top (X^\top X + \lambda I)^{-1} X^\top Y とする。 リッジ正則化したことにより「最良の線形近似といまの線形近似との誤差の期待値」はゼロにならなくなる。 h(x_0)逆行列の箇所に \lambda I が加わった分、バリアンスは小さくなる。
ニューラルネットワーク 例えば隠れ層のニューロン数を増やすとバイアスは小さくなる。 例えば隠れ層のニューロン数を増やすとバリアンスは大きくなる。
決定木。 例えば木を深くするとバイアスは小さくなる。 例えば木を深くするとバリアンスは大きくなる。

f:id:cookie-box:20180305231302p:plain:w60

確かにその表の例だと、バイアスを下げるとバリアンスが上がって、バリアンスを下げるとバイアスが上がる傾向がありそうだな。これがトレードオフか…でも、それだと結局「期待2乗誤差が小さくなるように適宜ハイパーパラメータを調節しましょう」ってことにならない? いまいちバイアスとバリアンスに分解した甲斐がないっていうか…。

f:id:cookie-box:20180305232608p:plain:w60

確かにそうですね。でも、例えば…ハヤトはここ1年の毎日の気温とアイスクリームの売り上げのデータをもとに、アイスクリーム売り上げ予測モデルを学習しました。そこにランプの精が現れていいました。「私は並行世界から来ました。並行世界のあなたが学習したモデルを差し上げます。並行世界の気温とアイスクリームの売り上げを生成する分布はこの世界と同じです。しかしこの世界と実現値はきっと異なるでしょう。予測に用いたモデルと学習方法はあなたと全く同じです」と。ちなみに同様のランプの精が計9人現れました。そしてハヤトは手元に自分が学習したモデル1つと、並行世界の自分が学習したモデル9つの、計10個のモデルを手に入れました。どうしますか?

f:id:cookie-box:20180305231302p:plain:w60

ええ…まあ突っ込みはさておき、同じモデルならバイアスもバリアンスも同じだからどれがよいモデルとかないよな。まあモデルが10個あるならなんか平均してみたくなるけど。

f:id:cookie-box:20180305232608p:plain:w60

ではモデルを平均しましょう。そうすると平均2乗誤差が小さくなることこそあれ、大きくなることはないとテキストの326~327ページにあります。訓練データを生成する分布を \mathcal{P} とし、この世界のハヤトが学習したモデルを \hat{f}^\ast(x) とし、それとあらゆる並行世界のハヤトたちが学習したモデルを平均したモデルを f_{\rm ag}(x) = E_{\mathcal{P}} \bigl[\hat{f}^\ast(x)\bigr] とします。E_{\mathcal{P}} \bigl[\cdot\bigr] は色々な訓練データを得ることについての期待値です(と思っています)。すると、ある (x,y) に対する2乗誤差の期待値について以下の不等式が成り立ちます(テキストでは y が確率変数であるようにかかれていますが、この議論はどちらかというとある y についてのものではないかと思っていますが、間違っていたらすみません)。

 \begin{split} E_{\mathcal{P}} \Bigl[ \bigl( y - \hat{f}^\ast(x) \bigr)^2 \Bigr] &= E_{\mathcal{P}} \Bigl[ \bigl( y- f_{\rm ag}(x) + f_{\rm ag}(x) -\hat{f}^\ast(x)  \bigr)^2 \Bigr] \\ &= E_{\mathcal{P}} \Bigl[ \bigl( y - f_{\rm ag}(x)\bigr)^2 + \bigl( f_{\rm ag}(x) -\hat{f}^\ast(x)  \bigr)^2 + 2\bigl( y - f_{\rm ag}(x)\bigr) \bigl( f_{\rm ag}(x) -\hat{f}^\ast(x)  \bigr) \Bigr] \\ &= \bigl( y - f_{\rm ag}(x)\bigr)^2 + E_{\mathcal{P}} \Bigl[ \bigl( f_{\rm ag}(x) -\hat{f}^\ast(x)  \bigr)^2 \Bigr]  \\ &= \bigl( y - f_{\rm ag}(x)\bigr)^2 + V_{\mathcal{P}} \Bigl[ \hat{f}^\ast(x) \Bigr] \\ &\geqq \bigl( y - f_{\rm ag}(x)\bigr)^2 \end{split}

つまり、1つのモデルの2乗誤差の期待値より、平均したモデルの2乗誤差の方が小さいか同じになります――モデルのバリアンスの分だけ。まあいまは \mathcal{P} から生成されるあらゆる訓練データについて学習したモデルが手元にあるわけではなく、ある10セットの訓練データで学習したモデルが手元にあるだけですが、訓練データの出方の経験分布を \mathcal{P} と考えれば話は同じはずです。

f:id:cookie-box:20180305231302p:plain:w60

本当だ…あれでも、その論理だといつでも複数のモデルを平均した方がいいの? 誤差を小さくしたかったら何でもいいから助っ人のモデルたちを用意すればいいってこと??

f:id:cookie-box:20180305232608p:plain:w60

ではありません。バイアスも平均されるわけですから、助っ人のモデルによってバリアンスが抑えられる以上にバイアスが増加したら意味がないです。何でもいいから寄せ集めればいいということにはなりません。

f:id:cookie-box:20180305231302p:plain:w60

ああ確かに…いまは10個のモデルのバイアスの期待値が同じだから、10個のモデルを平均してもバイアスが増える心配はなかったのか。

f:id:cookie-box:20180305232608p:plain:w60

それに、いまは回帰モデルで2乗誤差を損失とする場合でしたが、損失が異なる場合や、分類が目的の場合もありますよね。例えば個々のモデルが 0 か 1 を出力するような2値分類モデルであって、不正解率を損失とする場合、モデルをどうマージすればいいでしょうか。

f:id:cookie-box:20180305231302p:plain:w60

えっ、うーん、個々のモデルが 0 か 1 を出しているなら、10個のモデルの出力を平均したらきっと 0.3 とか 0.6 とかになるよな。でも、正解か不正解をみたいなら最終的には 0 か 1 かを出さないといけないから…閾値を決めて、0.5 未満なら 0、0.5 以上なら 1 にするとか? ぴったり 0.5 はどう扱うべきか迷うけど…。

f:id:cookie-box:20180305232608p:plain:w60

要は多数決ですね。しかしこの場合、マージすることが誤差を小さくすることにつながると限りません。テキスト327ページにある例ですが、あるデータの正解が 1 とします。そして、個々のモデルがこのデータに正しく 1 と出力できる確率が 0.4 とします。個々のモデルの正解率の期待値が 0.4 ということです。しかし、あらゆるモデルの多数決をとると必ず 0 を出力してしまいますね。4割のモデルが1、6割のモデルが0と判定するので。なので、多数決によって正解率がかえって 0 になってしまうのです。個々のモデルの正解率の期待値が 0.6 なら多数決によって正解率が 1 になるんですが。多数決を取る場合、多数派がよい予測をしていなければならないんです。こうかいてみると当然ですが。

f:id:cookie-box:20180305231302p:plain:w60

なら、バイアスを小さくする方に振ってバリアンスはモデルの平均で抑えにいくのがいいってことなのか…って、いや、現実には並行世界からランプの精こないじゃん!

f:id:cookie-box:20180305232608p:plain:w60

いかにも。なので、並行世界を自分でつくるしかありません。それがブートストラップ標本です。ブートストラップ標本とは、N 個のデータから N 個のデータを復元抽出(抽出したサンプルを元に戻して次のサンプルを抽出していく)することを B 回繰り返して B セットの標本を得たものです(301ページ;正確にはノンパラメトリックブートストラップですが)。

f:id:cookie-box:20180305231302p:plain:w60

え、うーん、そんなんでいいの? 確かに標本を複数用意できるけど、なんか自作自演みたいな…。

(その2があれば)つづく

雑記: モデルを選択したい話(カステラ本7.4節~7.7節)

私の誤りは私に帰属します。お気付きの点がありましたらお手数ですがご指摘いただけますと幸いです。
テキスト(カステラ本)その他の参考文献まとめ(解釈を含む)
  • 訓練データの被説明変数には独立に同一の分布にしたがう平均 0 のノイズがのっていると考える。と、その訓練データに最適化したモデルはたまたま出たノイズに引っ張られている可能性がある。したがって、訓練データ上の誤差(訓練誤差  \overline{{\rm err}})でモデル間を比較するのは適切ではない。手元に訓練データしかないとしても、あらゆるノイズの出方に対する  \overline{{\rm err}} の期待値をとったもの(訓練標本内誤差  {\rm Err}_{\rm in})で比較するべきである。
  • が、訓練標本内誤差  {\rm Err}_{\rm in} は知り得ないので、あらゆるノイズの出方に対して最適化したときの期待値  E_y \bigl[{\rm Err}_{\rm in}\bigr] を推定することになる。例えば線形回帰モデルで2乗誤差の和を損失としたときは、 E_y \bigl[{\rm Err}_{\rm in}\bigr] = E_y \bigl[\overline{{\rm err}}\bigr] + 2d \sigma_\varepsilon^2/N になることがわかる(d は説明変数の次元数、\sigma_\varepsilon^2 はノイズの分散、N は訓練データの個数)。ので、 E_y \bigl[\overline{{\rm err}}\bigr] を計算可能な  \overline{{\rm err}} で代用した  \overline{{\rm err}} + 2d \sigma_\varepsilon^2/N {\rm Err}_{\rm in} の推定値とすればよい。2d \sigma_\varepsilon^2/N は「手元の訓練データにたまたま出たノイズを学んだせいで誤差が過小評価されている分を補完(?)するペナルティ項」といえる。
  • また、損失を交差エントロピーとしたときはより一般のモデルで漸近的に  E_y \bigl[{\rm Err}_{\rm in}\bigr] \approx E_y \bigl[\overline{{\rm err}}\bigr] + d/N が成り立つ(d はモデルの有効パラメータ数)。この右辺を  \overline{{\rm err}} で代用した  \overline{{\rm err}} + d/N赤池情報量規準AIC)とよばれる。
  • 他方、元より「モデル候補 \mathcal{M}_m \, (m=1, \cdots, M) のうち訓練データ Z に対して最も尤もらしいモデルはどれか」を出発点に、ベイズ的にモデル選択しようとすることもできる。このとき、モデル候補 \mathcal{M}_m のパラメータの最尤推定値を \hat{\theta}_m とすると \log {\rm Pr} (Z|\mathcal{M}_m) = \log {\rm Pr} (Z| \hat{\theta}_m, \mathcal{M}_m) - d_m \log N / 2 + O(1) が成り立つ(d_m はモデル候補 \mathcal{M}_m の自由度)。これに基づいた  \overline{{\rm err}} + (d \log N )/Nベイズ情報量規準(BIC)とよばれる。
  • BICAICはペナルティ項が \log N 倍異なるが、手元のモデル候補のどれが真かを見出したいときはBICが、専ら予測誤差の最小化に興味があるときはAICが適していると考えられる。
キャラクターの原作とは無関係です。
f:id:cookie-box:20180305232608p:plain:w60

モデルを選択するとき、何に注意してどう選択するべきなのかって気になりますよね。そこでカステラ本の7章が「モデルの評価と選択」であるようなんです。面倒なので261ページから読みましょう。

f:id:cookie-box:20180305231302p:plain:w60

いや俺は気になっていないんだけど。とにかくテキストの261ページね…なんか誤差の定義がやたら多くない? 訓練誤差とか汎化誤差とか期待誤差とか…。

f:id:cookie-box:20180305232608p:plain:w60

モデルの誤差には「どの訓練データで学習したのか」「その訓練データにはどんなノイズがのっていたのか」「それをどのテストデータに適用したのか」がすべて関係してきますからね。何を固定して何を固定していないかによって呼び名を変えているようです。以下にまとめてみました。

訓練誤差
 \overline{{\rm err}} = \displaystyle \frac{1}{N} \sum_{i=1}^N L\bigl( y_i, \hat{f}(x_i) \bigr)
訓練データ \mathcal{T} = \bigl\{ (x_1, y_1), \cdots, (x_N, y_N) \bigr\} で学習したモデル \hat{f} で訓練データを予測したときの予測誤差の平均値。\hat{f} は訓練データにたまたま出たノイズまで学習しがちなので、誤差の指標としては「あまりに楽観的(261ページ)」である。
訓練標本内誤差
 {\rm Err}_{\rm in} = \displaystyle \frac{1}{N} \sum_{i=1}^N E_{Y^0} \Bigl[ L\bigl( Y_i^0, \hat{f}(x_i) \bigr)  \, \Big| \, \mathcal{T} \Bigr]
訓練データ \mathcal{T} で学習したモデル \hat{f} で訓練データを予測したときの予測誤差の、各データのノイズの出方に関する期待値の、平均値。訓練データ \mathcal{T} を、たまたま出たノイズは学ばずに、きちんと学べたかという指標になる。興味があるのはこの指標ではなく汎化誤差だが、モデル間を相対比較する分にはこの指標が使用しやすい。訓練誤差を訓練標本内誤差にするための調整分としての最善度 {\rm op} \equiv {\rm Err}_{\rm in} - \overline{{\rm err}} は「ノイズを学んでしまった分」となるがこれは直接推定できない。ふつうは訓練データの被説明変数に関する期待値を取った平均最善度  \omega \equiv E_y \bigl[ {\rm op}\bigr] を推定する(「最善度」と名付けられている割に、大きいとよい指標ではなく、ゼロになるべき指標である)。
汎化誤差
 {\rm Err}_{\mathcal{T}} = E_{X^0, Y^0} \Bigl[ L\bigl( Y^0, \hat{f}(X^0) \bigr) \, \Big| \, \mathcal{T} \Bigr]
訓練データ \mathcal{T} で学習したモデル \hat{f} でテストデータ点 (X^0, Y^0) を予測したときの予測誤差の、テストデータ点 (X^0, Y^0) に関する期待値。(おそらく)この指標に最も興味があるが、訓練データ \mathcal{T} しか手元にない以上これを直接推定することはできない。
期待誤差
 {\rm Err} = E_{\mathcal{T}} E_{X^0, Y^0} \Bigl[ L\bigl( Y^0, \hat{f}(X^0) \bigr) \, \Big| \, \mathcal{T} \Bigr]
汎化誤差の、訓練データ \mathcal{T} の取り方に関する期待値。
そして、2乗損失や0/1損失など色々な損失で  \omega = (2/N) \sum_{i=1}^N {\rm Cov}(\hat{y}_i, y_i) が成り立つそうです。2乗誤差の場合でこれを示せというのが演習問題になっていますね。ハヤト、試しに示してみますか?

f:id:cookie-box:20180305231302p:plain:w60

えっ、まずその \omega が何なのかよくわかっていないんだけど…まあどのみち先に {\rm Err}_{\rm in}\overline{{\rm err}} が要るみたいだから、L を2乗誤差にすると、

\begin{split} {\rm Err}_{\rm in} &= \displaystyle \frac{1}{N} \sum_{i=1}^N E_{Y^0} \Bigl[ \bigl( \hat{f}(x_i) - Y_i^0 \bigr)^2  \, \Big| \, \mathcal{T} \Bigr] \\ &= \frac{1}{N} \sum_{i=1}^N \hat{f}(x_i)^2 + \frac{1}{N} \sum_{i=1}^N E_{Y^0} \Bigl[(Y_i^0)^2 \, \Big| \, \mathcal{T} \Bigr] - \frac{2}{N} \sum_{i=1}^N  E_{Y^0} \Bigl[ Y_i^0 \, \Big| \, \mathcal{T} \Bigr] \hat{f}(x_i) \\ &= \frac{1}{N} \sum_{i=1}^N \hat{f}(x_i)^2 + E_{Y^0} \Bigl[ (Y_i^0)^2 \, \Big| \, \mathcal{T} \Bigr] - \frac{2}{N} E_{Y^0} \Bigl[ Y_i^0 \, \Big| \, \mathcal{T} \Bigr] \sum_{i=1}^N \hat{f}(x_i) \\ \overline{{\rm err}} &= \displaystyle \frac{1}{N} \sum_{i=1}^N \bigl(\hat{f}(x_i) - y_i \bigr)^2 \\ &= \frac{1}{N} \sum_{i=1}^N \hat{f}(x_i)^2 + \frac{1}{N} \sum_{i=1}^N y_i^2 - \frac{2}{N} \sum_{i=1}^N  y_i \hat{f}(x_i) \end{split}

こうなるよな。これを {\rm op} \equiv {\rm Err}_{\rm in} - \overline{{\rm err}} に代入して、

\begin{split} {\rm op} &= E_{Y^0} \Bigl[ (Y_i^0)^2 \, \Big| \, \mathcal{T} \Bigr] - \frac{2}{N} E_{Y^0} \Bigl[ Y_i^0 \, \Big| \, \mathcal{T} \Bigr] \sum_{i=1}^N \hat{f}(x_i) - \frac{1}{N} \sum_{i=1}^N y_i^2 + \frac{2}{N} \sum_{i=1}^N  y_i \hat{f}(x_i) \end{split}

さらにこれを訓練データの被説明変数について期待値をとったのが \omega だから、訓練データのノイズが変わったら変わりうるところに E_{Y^0}[ \, \cdot \, | \, \mathcal{T}] をかぶせればいいわけだから、

\begin{split} \omega &= E_{Y^0} \Bigl[ (Y_i^0)^2 \, \Big| \, \mathcal{T} \Bigr] - \frac{2}{N} E_{Y^0} \Bigl[ Y_i^0 \, \Big| \, \mathcal{T} \Bigr] \sum_{i=1}^N E_{Y^0} \Bigl[ \hat{f}(x_i)\, \Big| \, \mathcal{T} \Bigr] \\& \quad \;- \frac{1}{N} \sum_{i=1}^N E_{Y^0} \Bigl[ (Y_i^0)^2 \, \Big| \, \mathcal{T} \Bigr] + \frac{2}{N} \sum_{i=1}^N E_{Y^0} \Bigl[ Y_i^0 \hat{f}(x_i) \, \Big| \, \mathcal{T} \Bigr] \\ &= - \frac{2}{N} E_{Y^0} \Bigl[ Y_i^0 \, \Big| \, \mathcal{T} \Bigr] \sum_{i=1}^N E_{Y^0} \Bigl[ \hat{f}(x_i)\, \Big| \, \mathcal{T} \Bigr] + \frac{2}{N} \sum_{i=1}^N E_{Y^0} \Bigl[ Y_i^0 \hat{f}(x_i) \, \Big| \, \mathcal{T} \Bigr] \\ &= \frac{2}{N} \sum_{i=1}^N \biggl\{ E_{Y^0} \Bigl[ Y_i^0 \hat{f}(x_i) \, \Big| \, \mathcal{T} \Bigr]- E_{Y^0} \Bigl[ Y_i^0 \, \Big| \, \mathcal{T} \Bigr] E_{Y^0} \Bigl[ \hat{f}(x_i)\, \Big| \, \mathcal{T} \Bigr] \biggr\} \end{split}

これって共分散の公式の形だ! だったら  \omega = (2/N) \sum_{i=1}^N {\rm Cov}(\hat{y}_i, y_i) だ(\hat{y}_i = \hat{f}(x_i) とする)。この {\rm Cov}(\cdot, \cdot) は訓練データの説明変数は固定した下での被説明変数に関する共分散(ノイズの出方が色々変わったときの共分散)だな。どう、これでいい?

f:id:cookie-box:20180305232608p:plain:w60

間違っていないとは思いますが、カステラ本には演習の解答がないのでなんとも。

f:id:cookie-box:20180305231302p:plain:w60

ええ解答ないの!?

f:id:cookie-box:20180305232608p:plain:w60

何にせよ、モデルの  {\rm Cov}(\hat{y}_i, y_i) が正である状況とは、「y_i が大きくなる方向にノイズが変化したら \hat{y}_i = \hat{f}(x_i) も大きくなる」といった状況であるわけです。そんなノイズに反応するようなモデルは望ましくありません。 {\rm Cov}(\hat{y}_i, y_i) が負であっても困ります。 {\rm Cov}(\hat{y}_i, y_i) はゼロであってほしいですよね。

f:id:cookie-box:20180305231302p:plain:w60

それはそうだな。じゃあモデルが  {\rm Cov}(\hat{y}_i, y_i) がどうなるかを調べれば、そのモデルがどれくらいノイズを学習してしまうかわかるんだな。それでさっそく263ページに加法的誤差モデルの例が出てきて…いや、 \sum_{i=1}^N {\rm Cov}(\hat{y}_i, y_i) = d \sigma_\varepsilon^2 ってなんで? というか加法的誤差モデルって何?

f:id:cookie-box:20180305232608p:plain:w60

おそらく25ページがここでいう加法的誤差モデルの定義なんじゃないかと思うんですが、25ページの記述は「加法的モデル」なんですよね…まあ、 \sum_{i=1}^N {\rm Cov}(\hat{y}_i, y_i) = d \sigma_\varepsilon^2 になること自体は線形回帰モデルの最小2乗法で確かめるといいのではないかと。257ページの上の方でもその事実をつかっていますし。少なくとも7章に導出はなさそうですが。イメージとしては、最小2乗回帰したモデルは未知の点への予測値を訓練データを重ね合わせてつくり出していると解釈できますから、それを介して訓練データのノイズを取り込んでしまうわけです。

ただ、おそらくここの主張は線形回帰より広いクラスでもこうなんだというものと思います(以下)。

  • 説明変数の各次元に適用したモデルを足し合わせるようなモデル(Ex. 線形回帰)で損失を2乗誤差(2乗誤差以外でも結構成り立つと思われる)とするときは、
    • 平均最善度(「ノイズを学んでしまった分」の期待値)は \omega = 2 d \sigma_\varepsilon^2 / N になる。
    • よって、 E_y \bigl[ {\rm Err}_{\rm in} \bigr] = E_y \bigl[ \overline{{\rm err}} \bigr] + 2 d \sigma_\varepsilon^2 / N になる(7.24式)。
    • よって、説明変数の次元 d が多いほどノイズを学んでしまいやすい。
ここまで平均最善度がはっきりしていると、「説明変数を1次元増やして(自由に調整できる要素を1つ増やして)誤差が 2 \sigma_\varepsilon^2 / N より減らないなら増やさない方がいい」などといった判断ができそうですよね。

ときにハヤト、いま「標本データ上での誤差を、自由度で調整して、モデルのよさを測る」というようなことをしたわけです。「標本データ上での誤差を、自由度で調整して、モデルのよさを測る」と聞いて、何か思い出しませんか?

f:id:cookie-box:20180305231302p:plain:w60

「自由度で調整」って、もしかして赤池情報量規準AIC)のこと? 何次多項式で回帰するべきかって文脈でよくみるよな。意味はよくわかってないけど…。

f:id:cookie-box:20180305232608p:plain:w60

はい、実際次の264ページに、より一般的なモデルに適用できる {\rm Err}_{\rm in} の推定方法として、AICが出てくるんです。まず、天下りですが N \to \infty で7.27式(以下)の関係式が成り立つことがわかっています(いつ成り立つのか曖昧ですが、AICの原論文はフリーアクセスではなさそうなので一旦置いておきます)。

 -2E\Bigl[ \log {\rm Pr}_{\hat{\theta}}(Y)\Bigr] \approx \displaystyle -\frac{2}{N} E\Bigl[ \sum_{i=1}^N \log {\rm Pr}_{\hat{\theta}}(y_i)\Bigr] + \frac{2d}{N} \tag{7.27}

f:id:cookie-box:20180305231302p:plain:w60

確かにめっちゃ降って湧いてきたな。というかこの式どういう意味?

f:id:cookie-box:20180305232608p:plain:w60

まず E[\cdot]Y の真の分布の上での期待値ですね。ここで、真の分布は {\rm Pr}_{\theta^\ast}(Y) で表されることを仮定しています。つまり、真のパラメータ  \theta^\ast さえ突き止めれば真の分布を再現できるという状況です。といっても、手元には有限の N 個のサンプルしかないのでこの理想的な状況でも真の分布の推定が必ずできるわけではありません。限られた手元の標本データで尤度が最大になるようにしたパラメータが  \hat{\theta} です。それで、7.27式の左辺は真の分布と、手元の標本データで最善を尽くした分布の交差エントロピーの2倍ですね。これは一番知りたいものです。しかし、真の分布は知りえませんから、手元の標本データの経験分布で代用すると  -(1/N) \sum_{i=1}^N \log {\rm Pr}_{\hat{\theta}}(y_i) となりますね。この値の期待値(色々標本データを取り直したときの期待値)の2倍が右辺第1項です。しかし、この右辺第1項は、 \hat{\theta} に標本データのクセが反映されている以上、左辺の推定値としては「あまりに楽観的」ですよね。その調整分が右辺第2項であるわけです。なので結局、右辺第1項の期待値を手元のデータでの値で代用した {\rm AIC} = -(2/N) \sum_{i=1}^N \log {\rm Pr}_{\hat{\theta}}(y_i) + 2d/N をモデルのよさの推定値とせよというのがここでの主張ですね。右辺第1項に出てくる \sum_{i=1}^N \log {\rm Pr}_{\hat{\theta}}(y_i) は最大対数尤度(=現在の d で尤度を最大にした \hat{\theta} での尤度)になっていますが、「手元のデータで交差エントロピーを推定したもの」といった方が意味合い的にしっくりくる気がします。同じなんですが。

f:id:cookie-box:20180305231302p:plain:w60

ふーん…あれでも、最小2乗線形回帰モデルだと  E_y \bigl[ {\rm Err}_{\rm in} \bigr] = E_y \bigl[ \overline{{\rm err}} \bigr] + 2 d \sigma_\varepsilon^2 / N だった気がするんだけど、調整分が \sigma_\varepsilon^2 だけずれていない?

f:id:cookie-box:20180305232608p:plain:w60

どちらかというと、最小2乗線形回帰モデルの  E_y \bigl[ {\rm Err}_{\rm in} \bigr] = E_y \bigl[ \overline{{\rm err}} \bigr] + 2 d \sigma_\varepsilon^2 / N が7.29式(7.27式)の両辺に \sigma_\varepsilon^2 をかけたものになっていますね。順を追って説明します。

  • まず、一般にAICを利用するときは、「真の分布からの標本データ \{x_1, \cdots, x_N\} が手元にあって、なるべく真の分布 f(x) に近い予測分布 \hat{f}(x) をつくりたい(近さの基準は交差エントロピー)」ことが大前提になるはずです。
  • 対して7.4節の文脈での当面の目的は「訓練標本内誤差 {\rm Err}_{\rm in} を小さくしたい」です。より詳しくいうと、「訓練データ \bigl\{ (x_1, y_1), \cdots, (x_N, y_N) \bigr\} の説明変数は固定した下で、被説明変数は標本なのだと考えて、真の被説明変数の分布上での期待損失が小さいモデルをつくりたい」です。真の分布 f(x,y) の予測分布 \hat{f}(x,y) がつくりたいのだと考えれば、AICを利用したい状況と同じになるでしょう。ただし、「損失が交差エントロピーであれば」です。
  • 最小2乗線形回帰モデルは損失が2乗誤差なので、AICとは目的が少し違っています。が、ノイズの分散 \sigma_\varepsilon^2 が既知であり、最小2乗線形回帰モデルの予測値を中心に分散 \sigma_\varepsilon^2 をもつガウス分布を予測分布と考えるなら、両者の目的は一致するんです。なぜなら、このときの交差エントロピー  -\sum_i \int f(x_i,y) \log \hat{f}(x_i,y) dy は2乗誤差を 2\sigma_\varepsilon^2 で割ったものに定数バイアスを足したものになるからです(※)。交差エントロピー最小化が2乗誤差最小化と一致するんです。
  • ところでAICは「交差エントロピーの2倍」の推定値として構築されていました。これをいまの状況で「2乗誤差」の推定値にかき換えたいなら、\sigma_\varepsilon^2 をかければいいです。「交差エントロピーの2倍」が2乗誤差を \sigma_\varepsilon^2 で割ったものになっているわけですから。定数バイアスは両辺で打ち消し合いますたぶん。
なので調整項に \sigma_\varepsilon^2 のずれが発生しているわけです。ちなみに(※)はこの記事の最下部でプロデューサーさんが適当に計算しただけなので適当です。

f:id:cookie-box:20180305231302p:plain:w60

あー、何を推定しているのかがずれていたのか。一瞬「AIC を使えば \sigma_\varepsilon^2 が要らないの?」って思ったけど、AIC にしたらしたで右辺第1項の計算時に \sigma_\varepsilon^2 が要るな。

f:id:cookie-box:20180305232608p:plain:w60

ノイズ \sigma_\varepsilon^2 が大きいときほどノイズを学習してしまう危険性は増すわけですからね。極端な話、ノイズが全くないデータなら {\rm Err}_{\rm in}\overline{{\rm err}} は一致して最善度は常にゼロですし。まあそれでも訓練データの隙間や外挿部分の正解は知りえませんが。

f:id:cookie-box:20180305231302p:plain:w60

ところでジュン、265ページの「もし基底関数が適応的に選ばれるなら式 (7.23) はもはや成り立たない」ってどういうこと?

f:id:cookie-box:20180305232608p:plain:w60

かいてある通りなんですが、例えば…ハヤトはいま1次元の入力で1次元の出力を回帰するモデルとして \hat{f}(x) = a_0 + a_1 x + a_2 x^2 + \cdots のような多項式を候補に考えていて、何次多項式にするかを決めたいとします。平均 \hat{f}(x)、分散 \sigma^2 であるようなガウス分布を予測分布としましょう。切片と分散も推定対象とすると、m多項式モデルのパラメータの次元数は d =m + 2 になりますね。これでAICを適用すればいいです。しかしハヤトは考えました。「AICの調整項が面倒だな…そうだ、xm 次の項までのうち一番予測に有用な1個だけを使うことにしよう。これならパラメータの次元数は常に3になるから調整項が要らない。10次くらいまで調べて単に尤度が最大になる次数を選べばいい。このアイデアはモテるだろ!」と。しかし、それでは全然駄目なんです。

f:id:cookie-box:20180305231302p:plain:w60

俺そんなこといわない! 俺のキャラクターの把握が雑! そもそもAICの調整項 2d/N の計算に面倒な要素ないだろ!! …まあともかく、なんか駄目なのはわかるよ。なんか後出しじゃんけんっぽいし。見かけのパラメータ数でAICを計算したら駄目ってことだよな。でも、それならどうやって d を決めればいいんだ?

f:id:cookie-box:20180305232608p:plain:w60

次の7.6節に出てきますが、訓練データへの予測値 \hat{y} を訓練データの被説明変数のベクトル y\hat{y} = Sy と表し、S \in \mathbb{R}^{N \times N} を「訓練データの説明変数のベクトル x には依存するが y には依存しない」行列としたとき、{\rm tr}(S) が有効パラメータ数になるということです。

f:id:cookie-box:20180305231302p:plain:w60

へ、行列のトレース?

f:id:cookie-box:20180305232608p:plain:w60

最小2乗線形回帰 \hat{f}(x) = x^\top (X^\top X)^{-1} X^\top Y なら S = X(X^\top X)^{-1} X^\top になると思います。それならトレースの公式から {\rm tr} \bigl( X(X^\top X)^{-1} X^\top \bigr) = {\rm tr} \bigl( (X^\top X)^{-1} X^\top X \bigr) = {\rm tr} \bigl( I_d \bigr) = d となって成り立っていますね。ここで dx の次元数です。先ほどの多項式の例だと S がどうなるかというのは示せないんですが、S が「y には依存しない」という条件がある以上、適応的に次数を選んだなら調べた次数を全て取り込まないといけない気がします。まずは全ての次数でパラメータを最適化してみなければならないでしょう。その上で、ある特定の次数の項だけを選ぶような強烈な正則化をするのではないかと思うんですが…。

f:id:cookie-box:20180305231302p:plain:w60

歯切れ悪いな…まあそれで次の7.7節は、{\rm BIC} = -2 \sum_{i=1}^N \log {\rm Pr}_{\hat{\theta}}(y_i) + d \log N っていうのが出てきた。これ、AICN 倍になっているというわけでもないな。調整項が 2d だったのが d \log N になっているし。でも、AICの調整項は最小2乗線形回帰の例でも裏付けられていたはずで…どうなってんの??

f:id:cookie-box:20180305232608p:plain:w60

AICとは「測っているもの」が異なるということでしょう。Schwarz の論文(BICの原論文)は以下で閲覧できました。

イントロダクション中に、AICはこうだがBICはこうだとありますね。
  • choosing the model for which \log M_j (X_1, \cdots, X_n) - k_j is largest
  • Choose the model for which \log M_j (X_1, \cdots, X_n) - \frac{1}{2} k_j \log n is largest
イントロダクションに「ある自由度のモデルの最尤推定量は、ベイズ定量のサンプル数を大きくした極限として得られる」とあり、だからベイズ的なアプローチでモデル選択できるということだと思います。というかテキスト7.7節内にも結構かいてありますね。BICの出発点は「モデル候補 \mathcal{M}_m \, (m=1, \cdots, M) のうち最もよいモデルはどれか」であるようです。そして、BICは以下が成り立つことに基づいています。

 \displaystyle \log {\rm Pr} (Z|\mathcal{M}_m) = \log {\rm Pr} (Z| \hat{\theta}_m, \mathcal{M}_m) -\frac{d_m}{2} \log N + O(1) \tag{7.40}

この式を日本語訳するならば、「モデル候補 m の下で訓練データ Z が得られる対数尤度は、そのモデルのパラメータを最尤推定 \hat{\theta}_m にした下での対数尤度から d_m \log N /2 を差し引いたものである(d_m はモデル候補 m の自由度)」といったところでしょうか。そしてこれは、AICのベースとなっている7.27式とは左辺が既に違います。あちらは「尤度最大のパラメータでの、訓練標本内交差エントロピーは?」と問うているのに対して、こちらは「どのモデルがこの訓練データに対して最も尤もらしい?」なのですから。

ところで、テキストのAICBICを比べやすくするために、右辺第1項をそろえてみましょう。

\begin{split} \frac{\rm AIC}{2} &= - \frac{1}{N} \sum_{i=1}^N \log {\rm Pr}_{\hat{\theta}}(y_i) + \frac{d}{N}\\ \frac{\rm BIC}{2N} &= - \frac{1}{N} \sum_{i=1}^N \log {\rm Pr}_{\hat{\theta}}(y_i)  + \frac{d \log N}{N} \end{split}

どちらも N \to \infty で調整項がゼロに収束するのは同じです。ただ、収束のスピードが異なります。相対的に、AICの調整項は甘く、BICの調整項は厳しいことになります。つまり、テキストにもありますが、相対的に、AICは複雑な、BICは単純なモデルを選択する傾向があるということですね。そして、「モデル選択が目的のとき、AICBICのどちらを使うべきかがはっきりしているわけではない(269ページ)」と。

f:id:cookie-box:20180305231302p:plain:w60

ええどっちがいいとかないの? そういえば、269ページに「真のモデルを含むモデル集合が与えられたときにBICが正しいモデルを選ぶ確率は、標本数 N \to \infty のとき1に近づく」(漸近一致性をもつ)ってあるけど、逆にAICはこれが成り立たないの? 成り立ってほしいような気がするんだけど…。

f:id:cookie-box:20180305232608p:plain:w60

そのようですね。つまり、真のモデルより複雑なモデルがノイズを学んでしまう分に釣り合ったペナルティを d/N では課し切れないということでしょう。…しかし、AICの目的は「真のモデルはどれだろうか」だったでしょうか? 確かにBICの目的は「真のモデルはどれだろうか」でした。しかし、AICの目的は「そのモデルを尤度最大にチューニングした下で真の分布とどれだけ離れているか」なんです。真のモデルを外していようと、誤差を小さくすれば正義なんです。

例えば以下の文献にはAICは漸近有効性をもつとあります。つまり、予測誤差を最小にするのはAICであるということです。

そしてこの文献は、AICBICは何を最適としているかが違うと指摘しています。
  • In fact the conflict is easily resolved once it is acknowledged that ‘‘asymptotically optimal’’ can have several meanings. Asymptotic efficiency and (asymptotic) consistency are different kinds of optimality.(633ページ)
この文献によると、McQuarrie and Tsai (1998) は、AICを筆頭とする予測誤差を最小にするための指標を「efficient estimators」、BICを筆頭とするモデルを確認するための指標を「consistent estimators」と区別したそうです。AICの眷属にはテキスト7.5節にも出てきた Mallows’ C_p や交差検証(cross validation)が、BICの眷属には Hannan and Quinn や Geweke and Meese があるようですね。

f:id:cookie-box:20180305231302p:plain:w60

えっと、モデル選択の指標が2つの派閥(?)に分かれるのとか、両者はなんか目的が違うっぽいのはわかったけど、そこまでいうならどっちを使うべきってないの?

f:id:cookie-box:20180305232608p:plain:w60

先の文献には「真のモデルが非常に複雑で(自由度が非常に大きく)真のモデルを当てられる見込みがなく、モデル候補の自由度をそれに合わせなくてもよく、モデル候補のパラメータが真のモデルのパラメータを包含していなくてもよいし、余分なパラメータを含んでいてもよい」というときは前者を、「真のモデルがシンプルで、モデル候補のうち1つが真のモデルであると期待されるとき」は後者を利用するべきだとあります。確かに前者の状況でBICを使う意味は薄いと思うんですよね。モデル候補が真のモデルを含んでいる見込みが薄いなら。でも、ではこの状況でAICは上手く機能するのかはきちんと導出を追っていないのでわからないんですよね。いいんだと思うんですが。

おわり

※ 交差エントロピーが2乗誤差を 2\sigma_\varepsilon^2 で割ったものに定数バイアスを足したものになっている。
 \begin{split} \displaystyle -\int f(x_i,y) \log \hat{f}(x_i,y) dy &=-\frac{1}{\sqrt{2 \pi \sigma_\varepsilon^2}} \int  \exp \left( -\frac{(y - y^\ast)^2}{2 \sigma_\varepsilon^2} \right) \log \left[ \frac{1}{\sqrt{2 \pi \sigma_\varepsilon^2}} \exp \left( -\frac{(y - \hat{y})^2}{2 \sigma_\varepsilon^2} \right) \right] dy \\ &=- \frac{1}{\sqrt{2 \pi \sigma_\varepsilon^2}} \int   \left( -\frac{(y - \hat{y})^2}{2 \sigma_\varepsilon^2} \right) \exp \left( -\frac{(y - y^\ast)^2}{2 \sigma_\varepsilon^2} \right)dy + \log \sqrt{2 \pi \sigma_\varepsilon^2} \\ &= -\frac{1}{\sqrt{2 \pi \sigma_\varepsilon^2}} \int   \left( -\frac{(y - y^\ast)^2 +(y^\ast -\hat{y})^2 - 2(y - y^\ast)(y^\ast -\hat{y})}{2 \sigma_\varepsilon^2} \right) \exp \left( -\frac{(y - y^\ast)^2}{2 \sigma_\varepsilon^2} \right)dy  + \log \sqrt{2 \pi \sigma_\varepsilon^2} \\ &= \frac{\sigma_\varepsilon^2}{2 \sigma_\varepsilon^2} + \frac{(y^\ast -\hat{y})^2}{2 \sigma_\varepsilon^2} + \log \sqrt{2 \pi \sigma_\varepsilon^2} \end{split}

NeurIPS2020読みメモ: Adversarial Sparse Transformer for Time Series Forecasting

以下の論文を読みます。キャラクターの原作とは無関係です。私の誤りは私に帰属します。お気付きの点がありましたらご指摘ください。

Sifan Wu, Xi Xiao, Qianggang Ding, Peilin Zhao, Ying Wei, Junzhou Huang. Adversarial Sparse Transformer for Time Series Forecasting. In Pre-proceedings of the 33rd International Conference on in Neural Information Processing Systems (NeurIPS 2020), 2020. Paper
まとめ
  • 問題設定:
    • 関連するいくつかの時系列(例. 各家庭の15分ごとの電力消費量など)がある。
    • 各時系列について、少し離れた先の特定のパーセンタイルを予測したい。
      • 例えば、 50 パーセンタイルと 90 パーセンタイル、のように。
  • アプローチ:
    • 予測するパーセンタイルごとに Transformer を用意して学習する(曜日や時刻も入力するのが Positional Encoding の役割を担うと思われる)。損失はパーセンタイルに応じたピンボールロスにする。ただし、
      • 予測に無関係なステップを無視し、予測に関係あるステップに注意を集中したいので、アテンションにはソフトマックスではなく α-entmax を用いる(Sparse)
      • Transformer の学習には、ディスクリミネータも利用する(Adversarial)
  • 結果、electricity, traffic の1日後、7日後予測や wind, solar, M4-Hourly の予測のほとんどでも予測性能が既存手法を上回った(誤差の蓄積を回避できた)。
    • 単純な Transformer や Sparse Transformer よりも Adversarial Sparse Transformer がよかった。
    • なお、DeepAR を Adversatial に学習しても性能が向上した。
所感
  • どのようなデータを出力するべきかの指針を与えてくれる敵対的学習は、損失を上手く設計できない場面で(上手く設計できるとは)広く有効そうにみえる。
関連記事目次
GAN
f:id:cookie-box:20180305232608p:plain:w60

時系列予測モデルとして Adversarial Sparse Transformer (AST) なるモデルを提案しているんですが、GAN の要領で Transformer を学習するのが独特です。これによって誤差の蓄積を回避し、長期予測性能を向上させています。

f:id:cookie-box:20180305231302p:plain:w60

GAN って何だったっけ。てか俺誕生日なのになんで論文読まされてるの…。

f:id:cookie-box:20180305232608p:plain:w60

GAN は第3節の背景の一番最後に簡単に説明がありますが、せっかくですし原論文をみてみましょう。

GAN

GAN: Generative Adversarial Networks(敵対的生成ネットワーク)とは生成モデルなんですが、その学習方法が特徴的です。そうですね、いま、手元に人間の顔の画像の訓練データセットがあって、これを利用して人間の顔っぽい画像をランダムに生成したいとします。以下の手順でそんなモデルを得るのが GAN です。
  • 予め適当な空間上の適当な確率分布  p_z(z) を用意します。この分布から生成した z x = G(z; \theta_g) によって画像空間の元に変換します。これでランダムな画像を生成する機構はできました。問題は、この機構によって生成される画像の分布が、画像空間内の「人間の顔っぽいところ」に広がるようにすることです。そうなるように G を学習しなければなりません。
  • となるとそもそも「人間の顔っぽいとは何か」という話になってきますが、訓練データの画像と容易に識別できる画像ならそれは人間の顔っぽくはないだろう、と考えます。ので、訓練データ画像とランダム画像を識別するモデル  D(x; \theta_d) を用意しましょう。 D(x; \theta_d) は訓練データ画像なら 1、訓練データ画像ではないランダム画像なら 0 を取ってほしいです。そうなるように D を学習します。
  • すると G をどう学習するかも定まります。G の目標は、自分が出すランダム画像を、D が訓練データ画像と識別できないようにすることです。
  • DG の目標をまとめましょう。
    •  D の目標は、訓練データ画像には「訓練データ画像だ」と出力し、G から出てきた画像には「ランダム画像であって訓練データ画像ではない」と出力することです。
    •  G の目標は D を欺き、自分が出した画像を「訓練データ画像だ」と出力させることです。
    そしてこの最適化問題を数式でかくとこうです。
     \displaystyle \underset{G}{\rm min} \; \underset{D}{\rm max} \Bigl[ \mathbb{E}_{x \sim p_{\rm data}(x)}\bigl(\log D(x) \bigr) + \mathbb{E}_{x \sim p_z(z)}\bigl(\log (1 - D(G(z)) ) \bigr) \Bigr]
DD(G(z))0 にしようとし、GD(G(z))1 にしようとする敵対が起きていることがわかりますね。ちなみに上の目的関数は、DG に対して最適に学習されている場合は訓練データの分布と G が生成する画像の分布のJSダイバージェンスになります(原論文5ページ)。

f:id:cookie-box:20180305231302p:plain:w60

敵対的って物騒だな…って思ったけど、それを訊くと本格的に敵対してるな。片や欺こうとして、片や見抜こうとしてるんだから。…でもさ、人間の顔っぽい画像の分布がほしいなら、各訓練データに近いほど密度が大きい分布をつくって足し合わせるとかじゃ駄目なの?

f:id:cookie-box:20180305232608p:plain:w60

その近いとは何かという話だと思います。「この方向にこれくらいずらすのは近いのか」とかわからないでしょう。分散を小さめにしたらしたで、元の訓練データ内にあった画像しか出てこないモデルになってしまいますし。元の画像のどれかとかではなく本当にランダムだが見分けにくいようなものがほしくてこのようなことをするのだと思います。見分ける係(D)が人間ではなくてニューラルネットですから人間の感覚との差異はあるでしょうが、それでもかなり上手くいくんでしょう。

Transformer
f:id:cookie-box:20180305231302p:plain:w60

ふーん…それで、Transformer ってのは超ロボット生命体?

f:id:cookie-box:20180305232608p:plain:w60

ではありません。以下の論文で提案されたネットワークで、入力系列を同じ長さの特徴系列にエンコードし、それを利用して適当な長さの出力系列をデコードするものですね。

Transformer

今回の論文でも第3節の背景の The Transformer の箇所に説明がありますね。つまり、Transformer のエンコーダでは「マルチヘッドセルフアテンション」+「全結合」を N 回繰り返します。具体的には以下を N 回繰り返します。
  • d 次元ベクトルが n 個並んだ入力系列 h \in \mathbb{R}^{n \times d} を、
    • W_m^Q \in \mathbb{R}^{d \times d/M}Q_m = h W_m^Q \in \mathbb{R}^{n \times d/M}写像する(d/M 次元ベクトルが n 個並ぶ)。
    • W_m^K \in \mathbb{R}^{d \times d/M}K_m = h W_m^K \in \mathbb{R}^{n \times d/M}写像する(d/M 次元ベクトルが n 個並ぶ)。
    • W_m^V \in \mathbb{R}^{d \times d_v}V_m = h W_m^V \in \mathbb{R}^{n \times d_v} 次元に写像する(d_v 次元ベクトルが n 個並ぶ)。
  • Q_m K_m^\top \in \mathbb{R}^{n \times n} を計算し、各要素を \sqrt{d/M} で割る(n 次元ベクトルが n 個並ぶ)。その上で各 n 次元ベクトルをソフトマックスする。これで得られた \alpha_m を scaled dot-product attention とよぶ。
  • あとは O_m = \alpha_m V_m を計算する(d_v 次元ベクトルが n 個並ぶ)。
  • 以上のことを M ヘッドやった O_1, \cdots, O_M を concat する(M \times d_v 次元ベクトルが n 個並ぶ)。
  • 各ベクトルを全結合し、ReLU で活性化し、さらに全結合する。つまり {\rm max}(0, O W_1 +b_1)W_2 +b_2 とする。これは入力系列と同じ長さの特徴系列である。
余談ですが、BERT では ReLU ではなく GELU(Gaussian Error Linear Units)で活性化していますよね。ともかく、\alpha_m は、例えば長さ5の系列を入力したら、5×5行列になりますが、この1行目が意味するのは、「1番目の位置は、1~5番目の位置にどれだけずつ注目するか」であり、入力系列は自分に注意している(セルフアテンション)ことになります。エンコーダから出力される特徴系列は、前後の文脈を踏まえたその箇所の単語の特徴とでもいうべきものになっているでしょう。

ところで、最終的にほしいのはそんな特徴系列ではありません。例えば機械翻訳であれば、英語の文章をドイツ語の文章に翻訳したものなどがほしいはずです。なのでそのような出力系列を得るデコーダを用意します。デコーダには「文頭トークン、*、*、*、*」(* は未知)という系列を入れて、エンコーダと同様にマルチヘッドセルフアテンションしてまず O を得ます。が、このとき scaled dot-product attention が * に注意しないように、自分より後ろへの注意を0にします。なのでこの層は Masked Multi-Head Attention とよばれていますね。次に、再度マルチヘッドセルフアテンションしますが、この段では V_m = h_{\rm embed} W_m^V にエンコーダからの出力 h_{\rm embed} を取り込みます。自分への注意ではなく、特徴系列に注意するわけです。その後 全結合-ReLU-全結合 します。これをやはり N 回繰り返した後に、全結合-softmax して最終出力をします。これは例えば「すべてのドイツ語の単語上の確率分布」のようなものにすべきでしょう。こうして最初の単語を得て、次は「文頭トークン、最初の単語、*、*、*」をデコーダに入力してデコードを繰り返します。

f:id:cookie-box:20180305231302p:plain:w60

…処理が込み入ってるけど、Transformer は系列から系列を得るのに使えるんだな。入力文章の各位置に前後の文脈を反映させてエンコード文章にして、デコードするときはそこまででデコードできている範囲で前後の文脈を反映させてからエンコード文章を取り込むって感じか…ってあれ? いまは機械翻訳じゃなくて時系列の予測をしたいんだよな?

問題設定
f:id:cookie-box:20180305232608p:plain:w60

時系列予測も「現時点までの数ステップを参照してこれから未来の数ステップを予測する」と考えれば機械翻訳のようなものですよ。問題設定を確認しましょう。第3節の最初ですね。いま、手元に  \{y_{i, 1:t_0}\}_{i=1}^S という S 本の時系列があるとします。各 y_{i,t}\mathbb{R} の元です。これが興味のある時系列ですね(ターゲット時系列とよぶようです)。これに加えて、X_{i,1:t_0} \in \mathbb{R}^{t_0 \times k} という説明変数の時系列もあるとします。これは時間変化するものでも時間変化しないものでも構いませんが、未来の値も予測に利用しているようなので、未来の値が予めわかるものでないといけませんね。そしていま予測したいのは、S 本の各時系列の今後 \tau ステップ間の \rho パーセンタイルです。モデルは以下になります。

 \hat{Y}_{\rho, t_0+1:t_0+\tau} = f_\rho(Y_{1:t_0}, X_{1:t_0+\tau})
\rho ごとにモデルを用意するようですね。以下にイメージをかきましょう。例えば50パーセンタイルと90パーセンタイルを予測するとして、50パーセンタイルモデルは以下の青字の箇所を入力に緑字の箇所を出力します。90パーセンタイルモデルは以下の青字の箇所を入力にオレンジの字の箇所を出力します。まあ入力はどちらも同じなんですが。
時刻1\cdotst_0t_0 + 1\cdotst_0 + \tau
説明変数X_1\cdotsX_{t_0}X_{t_0+1}\cdotsX_{t_0+\tau}
ターゲット時系列 1y_{1,1}\cdotsy_{1,t_0}y_{1,t_0+1}^{(50)}\cdotsy_{1,t_0+\tau}^{(50)}
\vdots\vdots\ddots\vdots\vdots\ddots\vdots
ターゲット時系列 Sy_{S,1}\cdotsy_{S,t_0}y_{S,t_0+1}^{(50)}\cdotsy_{S,t_0+\tau}^{(50)}
時刻1\cdotst_0t_0 + 1\cdotst_0 + \tau
説明変数X_1\cdotsX_{t_0}X_{t_0+1}\cdotsX_{t_0+\tau}
ターゲット時系列 1y_{1,1}\cdotsy_{1,t_0}y_{1,t_0+1}^{(90)}\cdotsy_{1,t_0+\tau}^{(90)}
\vdots\vdots\ddots\vdots\vdots\ddots\vdots
ターゲット時系列 Sy_{S,1}\cdotsy_{S,t_0}y_{S,t_0+1}^{(90)}\cdotsy_{S,t_0+\tau}^{(90)}
ターゲット時系列 1 のみに注目すれば (y_{1,1}, \cdots, y_{1, t_0}) から (y_{1,t_0 + 1}, \cdots, y_{1, t_0 + \tau}) への機械翻訳のようなものでしょう?

f:id:cookie-box:20180305231302p:plain:w60

まあそうかもしれないけど。

f:id:cookie-box:20180305232608p:plain:w60

あとここにきて気付いたんですが、説明変数 X_t は Transformer でいう Positional Encoding とかトークンタイプエンコーディングの位置付けなんですね。例えば月や曜日などとありましたし。Figure 2 に Positional Encoding はありませんし。あ、Positional Encoding は、機械翻訳ではその単語が何単語目かという特徴ベクトルですね。それを予め各単語ベクトルに足し合わせてからマルチヘッドセルフアテンションに流します。Transformer は再帰や畳込みを行わないので、そうでもしないとその単語が文章の1番目にあったのか、2番目にあったのかが本質的に区別されませんからね。

f:id:cookie-box:20180305231302p:plain:w60

えっと、結局どうやって未来の時系列を予測するの?

f:id:cookie-box:20180305232608p:plain:w60

4節をみてまとめましょう。

  • 現時点までのターゲット系列と説明変数系列 (Y_{1:t_0}, X_{1:t_0}) をエンコーダに入力し、特徴系列 (h_1, \cdots, h_{t_0})エンコードする。
  • 未来の説明変数系列 (X_{t_0 + 1:t_0 + \tau}) と特徴系列 (h_1, \cdots, h_{t_0})デコーダに入力し、順次 \hat{y}_{t_0 + 1}, \cdots, \hat{y}_{t_0 + \tau} をデコードする。
    • 明記されていませんが、このとき最初に機械翻訳における文頭トークンに相当するものとして (Y_{t_0}, X_{t_0}) がフィードされていると思うんです。そうでなければ、特徴系列に注目する主体がいませんから。また、Figure 2 に明記されていませんが Masked Multi-Head Attention であるはずです。
とまあ、これだけなら純粋に Transformer を利用して時系列予測したという話です。しかし、本論文で提案しているのは Adversarial Sparse Transformer (AST) です。ただの Transformer と違って Adversarial で Sparse なんです。

f:id:cookie-box:20180305231302p:plain:w60

あ、敵対的なのとか忘れてた…。

Sparse
f:id:cookie-box:20180305232608p:plain:w60

先に Sparse の説明からいきましょう。

Sparse

時系列を予測するのに、過去の数ステップにしか注目したくなくて、関係ないステップには全く注目したくないという気持ちがあります(過去のどのステップが未来予測に関わってくるかは時系列によると思うんですが)。必要なステップへの注意に集中したいんです。しかし、\alpha_m はソフトマックスの結果ですから、小さい要素でも完全にゼロにはなりません。そこで α-entmax を採用します。以下で提案されたものですね。つまり、\alpha_{\rm entmax}(h) = [ (\alpha - 1) h - \tau {\bf 1}]^{1/(\alpha - 1)}_{+} です。ここで、[\cdot]_+ は ReLU、{\bf 1} は要素がすべて1のベクトルで、\tauh に応じてただ1つ定まる閾値です。\alpha > 1 はスパースさの度合いを定めるパラメータで、\alpha=1 のとき \alpha_{\rm entmax}(h) は softmax と一致するそうです。\alpha=2 のときは sparsemax なるものになるそうです。参照している論文の Figure 3 にいくつかの \alpha に対するグラフがありますね。\alpha = 1 のときは h がどれだけ小さくても小さな正の値になりますが、\alpha が大きくなるほどゼロに切り捨てられる範囲が広がってくることがわかると思います。ステップ関数に近づいてきますね。この論文では \alpha = 1.5 を採用するようです。

f:id:cookie-box:20180305231302p:plain:w60

え、 何その式の形…全然ソフトマックスっぽくないんだけど…。

f:id:cookie-box:20180305232608p:plain:w60

そうですね…元々 α-entmax は参照論文の (10) 式として定義されていて、これだとベクトル p の最適化を含んでしまうので、\tau の最適化に落とし込んだのが今回の式のようなんですが、詳細は参照論文をよく読まなければなさそうです。

Adversarial
f:id:cookie-box:20180305232608p:plain:w60

最後に Adversarial の説明です。

Adversarial

ここまでで現時点までの時系列から向こう何ステップかの予測値を生成する Sparse Transformer を用意しましたが、これを普通に学習するのではなく、ディスクリミネータも導入して学習します。ディスクリミネータには LeakReLu で活性化した3層の全結合層を用いたそうです。以降、Sparse Transformer のことをジェネレータとよびかえます。学習の1イテレーションは以下です。
  • データセットからランダムに  [ X_{1:t_0 + \tau}, Y_{1:t_0 + \tau} ] をサンプリングします。Y_{\rm real} \equiv  Y_{1:t_0 + \tau} とします。
  • 現時点のジェネレータで  \hat{Y}_{t_0 + 11:t_0 + \tau} を予測します。Y_{\rm fake} \equiv  Y_{1:t_0}  \circ \hat{Y}_{t_0 + 11:t_0 + \tau} とします。ディスクリミネータに D(Y_{\rm fake} )=1 と判定させるのがジェネレータの目標です。
  • 以下を小さくするようにジェネレータを更新します。
     \mathcal{L}_\rho (Y_{t_0+1:t_0+\tau}, \hat{Y}_{t_0 + 11:t_0 + \tau}) + \lambda \mathbb{E}[ \log( 1 − D(Y_{\rm fake} )) ]
  • 以下を小さくするようにディスクリミネータを更新します。
     \mathbb{E}[ - \log( D(Y_{\rm real} )) - \log( 1 − D(Y_{\rm fake} )) ]
ここで \mathcal{L}_\rho (Y_{t_0+1:t_0+\tau}, \hat{Y}_{t_0 + 11:t_0 + \tau}) は通常の学習の損失ですが、y_{i, t}\hat{y}_{i, t} と予測したときの損失の定義は、例えば 90 パーセンタイルの予測だったら、「上振れしたときは上振れ幅の 0.1 倍、下振れしたときは下振れ幅の 0.9 倍」です(論文は誤植だと思います)。基本的に絶対誤差なんですが、90パーセンタイルの予測だったら下振れ側に厳しく、10パーセンタイルの予測だったら上振れ側に厳しいイメージですね。

f:id:cookie-box:20180305231302p:plain:w60

…それって90パーセンタイルの予測になるの? だって、仮に \hat{y}_{i, t}y_{i, t} を完璧に予測できたらその損失ってゼロになるけど、それは90パーセンタイルの予測じゃなくない?

f:id:cookie-box:20180305232608p:plain:w60

完璧に予測できたらそうなるでしょうが、まず完璧に予測できないとして、「上振れはしてもほとんど下振れしない予測値」にはなっていると思います。それくらいの意味しかないと思いますよ。本当に90パーセンタイルを予測したかったら、90パーセンタイルのアノテーションデータが必要だと思いますし。

検証結果
f:id:cookie-box:20180305232608p:plain:w60

検証は electricity, traffic, wind, solar, M4-Hourly データセットで行ったようです。electricity, traffic については以下の記事で図示しました。

誤差の指標として6ページの (8) 式のρリスクを導入していますが、これは50パーセンタイルの予測であれば単なる絶対パーセント誤差ですね。90パーセンタイルの予測であれば、上振れに甘く、下振れに厳しくなります。

Table 2, Table 3 は electricity, traffic の1日後、1週間後予測ですが(Table 2 は50パーセンタイル、Table 3 は90パーセンタイル)、AST のρリスクが既存手法に比べて最小になっていますね。AST の左隣2列が T, ST となっていますが、これは単なる Transformer と Sparse Transformer でしょう。T より ST の方がよく、ST より AST の方がさらによいことがはっきりわかりますね。Table 6 には wind, solar, M4-Hourly の結果がありますが、50パーセンタイルは一貫して AST が最良ですが、90パーセンタイルについては DeepState や DeepAR がよい場合もありますね。

面白いのは、Table 5 で DeepAR にも敵対的学習を施してみている点です。DeepAR も敵対的学習によって性能が向上することがわかります。向上後も AST よりは誤差が大きいようですが。

f:id:cookie-box:20180305231302p:plain:w60

敵対的学習がそれだけ有効だったってこと? …でもさ、予測モデルの学習って正解データがあるんだよな。正解に合わせるように学習しさえすればいいんじゃないのか? なんでわざわざ敵(ディスクリミネータ)を用意して敵対させなきゃならないんだ…。

f:id:cookie-box:20180305232608p:plain:w60

今回の特徴は再帰的に何ステップも予測を続けていく点だと思います。例えば一期先予測を大きい側に d 謝るのと、小さい側に d 誤るのとでは、絶対誤差の観点ではどちらも同じだけの誤りです。しかし、「現時点までの時系列に続けてその値が観測されることは尤もらしいだろうか」という観点では同じ誤りではない可能性があります。ディスクリミネータ視点、片や「この系列は訓練データらしい(訓練データにあってもおかしくない)」、片や「この系列は訓練データらしくない」となっている可能性があります。そして、後者に陥ってしまったら、そこから先の予測は途端にくるってくる可能性があるかもしれません。だったら後者を優先して修正しなければならない…とか。

おわり

雑記

以下の記事を投稿しました。

上の記事の「凸関数の最小点の勾配が満たす必要十分条件」のパートで、雑記: KKT条件の話で「証明はここでは割愛します。」といっていた箇所の証明をしています(上の記事では凸関数を仮定しているので必要十分条件になっていますが、KKT条件の話では仮定していないので必要条件です)。

あと上の記事のイラストをGIFにしたものです。
f:id:cookie-box:20201121023037g:plain:w420

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from pylab import rcParams
rcParams['figure.figsize'] = 8, 6
rcParams['font.size'] = 16
rcParams['font.family']='Ume Hy Gothic O5'


fig, ax = plt.subplots()

if True:
    x0 = 1.0  # 記事中のイラストでいう x
    x1 = 4.2  # 記事中のイラストでいう y
    y0 = np.exp(x0)
    y1 = np.exp(x1)
    list_x = np.arange(0.2, 4.5, 0.2)
    ax.plot([x0, x0], [-10, y0], color='#cccccc')
    ax.plot([x1, x1], [-10, y1], color='#cccccc')
    ax.plot([0.2, x1], [y0, y0], color='#cccccc')
    ax.plot([0.2, x1], [y1, y1], color='#cccccc')
    ax.plot(list_x, np.exp(list_x), color='black')
    list_x_ = np.array([x0, x1])
    ax.plot(list_x_, np.exp(list_x_), color='green')

ims = []

for theta in np.arange(0.84, 0.0, -0.02):
    x2 = (1.0 - theta) * x0 + theta * x1
    y2 = np.exp(x2)
    y2_ = y0 + (y2 - y0) / theta
    im = plt.plot([x2, x2], [-10, y2], color='#cccccc') \
         + plt.plot([0.2, x1], [y2_, y2_], color='#cccccc', zorder=-10) \
         + plt.plot([x0, x1], [y0, y2_], color='blue')
    ims.append(im)

for i in range(10):  # 最後の状態でしばらく止めるためにコマ数を取っているだけ
    y2_ = y0 + y0 * (x1 - x0)
    im = plt.plot([0.2, x1], [y2_, y2_], color='#cccccc', zorder=-10) \
         + plt.plot([x0, x1], [y0, y2_], color='blue')
    ims.append(im)

ani = animation.ArtistAnimation(fig, ims, interval=150)
ani.save("animation.gif")