2020-12-31

雑記：モデルを選択したい話（カステラ本7.4節～7.7節）

雑記

私の誤りは私に帰属します。お気付きの点がありましたらお手数ですがご指摘いただけますと幸いです。

テキスト（カステラ本）

統計的学習の基礎 ―データマイニング・推論・予測― | Trevor Hastie, Robert Tibshirani, Jerome Friedman, 杉山将, 井手剛, 神嶌敏弘, 栗田多喜夫, 前田英作, 井尻善久, 岩田具治, 金森敬文, 兼村厚範, 烏山昌幸, 河原吉伸, 木村昭悟, 小西嘉典, 酒井智弥, 鈴木大慈, 竹内一郎, 玉木徹, 出口大輔, 冨岡亮太, 波部斉, 前田新一, 持橋大地, 山田誠 |本 | 通販 | Amazon の261～269ページ。

その他の参考文献

(PDF) Model selection for ecologists: The worldviews of AIC and BIC の633ページ。

Akaike information criterion - Wikipedia から参照されており、AICとBICは目的が異なることを強調している。635ページの表に両者の対比がまとまっている。

まとめ（解釈を含む）

訓練データの被説明変数には独立に同一の分布にしたがう平均 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）とよばれる。
BICとAICはペナルティ項が $\log N$ 倍異なるが、手元のモデル候補のどれが真かを見出したいときはBICが、専ら予測誤差の最小化に興味があるときはAICが適していると考えられる。

キャラクターの原作とは無関係です。

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

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

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

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

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

訓練誤差 $\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乗誤差の場合でこれを示せというのが演習問題になっていますね。ハヤト、試しに示してみますか？

えっ、まずその $\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)$ は訓練データの説明変数は固定した下での被説明変数に関する共分散（ノイズの出方が色々変わったときの共分散）だな。どう、これでいい？

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

ええ解答ないの！？

何にせよ、モデルの ${\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)$ はゼロであってほしいですよね。

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

おそらく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$ より減らないなら増やさない方がいい」などといった判断ができそうですよね。

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

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

はい、実際次の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}$

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

まず $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}$ での尤度）になっていますが、「手元のデータで交差エントロピーを推定したもの」といった方が意味合い的にしっくりくる気がします。同じなんですが。

ふーん…あれでも、最小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$ だけずれていない？

どちらかというと、最小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$ のずれが発生しているわけです。ちなみに（※）はこの記事の最下部でプロデューサーさんが適当に計算しただけなので適当です。

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

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

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

かいてある通りなんですが、例えば…ハヤトはいま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の調整項が面倒だな…そうだ、 $x$ の $m$ 次の項までのうち一番予測に有用な1個だけを使うことにしよう。これならパラメータの次元数は常に3になるから調整項が要らない。10次くらいまで調べて単に尤度が最大になる次数を選べばいい。このアイデアはモテるだろ！」と。しかし、それでは全然駄目なんです。

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

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

へ、行列のトレース？

最小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$ となって成り立っていますね。ここで $d$ は $x$ の次元数です。先ほどの多項式の例だと $S$ がどうなるかというのは示せないんですが、 $S$ が「 $y$ には依存しない」という条件がある以上、適応的に次数を選んだなら調べた次数を全て取り込まないといけない気がします。まずは全ての次数でパラメータを最適化してみなければならないでしょう。その上で、ある特定の次数の項だけを選ぶような強烈な正則化をするのではないかと思うんですが…。

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

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

https://projecteuclid.org/euclid.aos/1176344136

イントロダクション中に、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式とは左辺が既に違います。あちらは「尤度最大のパラメータでの、訓練標本内交差エントロピーは？」と問うているのに対して、こちらは「どのモデルがこの訓練データに対して最も尤もらしい？」なのですから。

ところで、テキストのAICとBICを比べやすくするために、右辺第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は単純なモデルを選択する傾向があるということですね。そして、「モデル選択が目的のとき、AICとBICのどちらを使うべきかがはっきりしているわけではない（269ページ）」と。

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

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

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

(PDF) Model selection for ecologists: The worldviews of AIC and BIC

そしてこの文献は、AICとBICは何を最適としているかが違うと指摘しています。

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 があるようですね。

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

先の文献には「真のモデルが非常に複雑で（自由度が非常に大きく）真のモデルを当てられる見込みがなく、モデル候補の自由度をそれに合わせなくてもよく、モデル候補のパラメータが真のモデルのパラメータを包含していなくてもよいし、余分なパラメータを含んでいてもよい」というときは前者を、「真のモデルがシンプルで、モデル候補のうち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}$

2020-11-23

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 に学習しても性能が向上した。

所感

どのようなデータを出力するべきかの指針を与えてくれる敵対的学習は、損失を上手く設計できない場面で（上手く設計できるとは）広く有効そうにみえる。

関連記事

雑記： NeurIPS 2020 pre-proceedings の「時系列」を含むタイトル - クッキーの日記

GAN
Transformer
問題設定
Sparse
Adversarial
検証結果

GAN

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

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

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

[1406.2661] Generative Adversarial Networks

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$ が訓練データ画像と識別できないようにすることです。
との目標をまとめましょう。
- $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]$

$D$ は $D(G(z))$ を $0$ にしようとし、 $G$ は $D(G(z))$ を $1$ にしようとする敵対が起きていることがわかりますね。ちなみに上の目的関数は、 $D$ が $G$ に対して最適に学習されている場合は訓練データの分布と $G$ が生成する画像の分布のJSダイバージェンスになります（原論文5ページ）。

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

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

Transformer

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

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

[1706.03762] Attention Is All You Need

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 して最終出力をします。これは例えば「すべてのドイツ語の単語上の確率分布」のようなものにすべきでしょう。こうして最初の単語を得て、次は「文頭トークン、最初の単語、＊、＊、＊」をデコーダに入力してデコードを繰り返します。

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

問題設定

時系列予測も「現時点までの数ステップを参照してこれから未来の数ステップを予測する」と考えれば機械翻訳のようなものですよ。問題設定を確認しましょう。第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$	$\cdots$	$t_0$	$t_0 + 1$	$\cdots$	$t_0 + \tau$
説明変数	$X_1$	$\cdots$	$X_{t_0}$	$X_{t_0+1}$	$\cdots$	$X_{t_0+\tau}$
ターゲット時系列 $1$	$y_{1,1}$	$\cdots$	$y_{1,t_0}$	$y_{1,t_0+1}^{(50)}$	$\cdots$	$y_{1,t_0+\tau}^{(50)}$
$\vdots$	$\vdots$	$\ddots$	$\vdots$	$\vdots$	$\ddots$	$\vdots$
ターゲット時系列 $S$	$y_{S,1}$	$\cdots$	$y_{S,t_0}$	$y_{S,t_0+1}^{(50)}$	$\cdots$	$y_{S,t_0+\tau}^{(50)}$

時刻	$1$	$\cdots$	$t_0$	$t_0 + 1$	$\cdots$	$t_0 + \tau$
説明変数	$X_1$	$\cdots$	$X_{t_0}$	$X_{t_0+1}$	$\cdots$	$X_{t_0+\tau}$
ターゲット時系列 $1$	$y_{1,1}$	$\cdots$	$y_{1,t_0}$	$y_{1,t_0+1}^{(90)}$	$\cdots$	$y_{1,t_0+\tau}^{(90)}$
$\vdots$	$\vdots$	$\ddots$	$\vdots$	$\vdots$	$\ddots$	$\vdots$
ターゲット時系列 $S$	$y_{S,1}$	$\cdots$	$y_{S,t_0}$	$y_{S,t_0+1}^{(90)}$	$\cdots$	$y_{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})$ への機械翻訳のようなものでしょう？

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

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

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

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 なんです。

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

Sparse

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

Sparse

時系列を予測するのに、過去の数ステップにしか注目したくなくて、関係ないステップには全く注目したくないという気持ちがあります（過去のどのステップが未来予測に関わってくるかは時系列によると思うんですが）。必要なステップへの注意に集中したいんです。しかし、 $\alpha_m$ はソフトマックスの結果ですから、小さい要素でも完全にゼロにはなりません。そこで α-entmax を採用します。以下で提案されたものですね。

[1905.05702] Sparse Sequence-to-Sequence Models

つまり、 $\alpha_{\rm entmax}(h) = [ (\alpha - 1) h - \tau {\bf 1}]^{1/(\alpha - 1)}_{+}$ です。ここで、 $[\cdot]_+$ は ReLU、 ${\bf 1}$ は要素がすべて1のベクトルで、 $\tau$ は $h$ に応じてただ1つ定まる閾値です。 $\alpha > 1$ はスパースさの度合いを定めるパラメータで、 $\alpha=1$ のとき $\alpha_{\rm entmax}(h)$ は softmax と一致するそうです。 $\alpha=2$ のときは sparsemax なるものになるそうです。参照している論文の Figure 3 にいくつかの $\alpha$ に対するグラフがありますね。 $\alpha = 1$ のときは $h$ がどれだけ小さくても小さな正の値になりますが、 $\alpha$ が大きくなるほどゼロに切り捨てられる範囲が広がってくることがわかると思います。ステップ関数に近づいてきますね。この論文では $\alpha = 1.5$ を採用するようです。

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

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

Adversarial

最後に 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パーセンタイルの予測だったら上振れ側に厳しいイメージですね。

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

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

検証結果

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

NeurIPS2018読みメモ： Deep State Space Models for Time Series Forecasting（その1） - クッキーの日記

誤差の指標として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 よりは誤差が大きいようですが。

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

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

おわり

2020-11-21

雑記

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

上の記事の「凸関数の最小点の勾配が満たす必要十分条件」のパートで、雑記： 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")

2020-11-15

Time Series Analysis：ノート5章最尤推定（その2）

本読み

以下の本の5章を読みます。キャラクターは架空のものです。私の誤りは私に帰属します。お気付きの点がありましたらコメントでご指摘いただけますと幸いです。

Time Series Analysis

作者:Hamilton, James D.
発売日: 1994/01/11
メディア: ハードカバー

前回： Time Series Analysis：ノート5章（その1） - クッキーの日記

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

前回の復習ですが、ノイズが正規ホワイトノイズであるような定常 AR(1) 過程 $Y_t = c + \phi Y_{t-1} + \varepsilon_t$ の実現値 $(y_1, \cdots, y_T)$ の同時分布は以下になりました。この分布は「最初に出てくる $y_1$ はどんな分布だろう」「 $y_1$ を観測した上で次の $y_2$ はどんな分布だろう」と積み上げていけば導くことができます。

定常ガウシアン AR(1) 過程の実現値の同時分布

$\displaystyle f_{Y_T, \cdots, Y_1}(y_T, \cdots, y_1 \, ; c, \phi, \sigma^2) = \frac{\sqrt{1 - \phi^2}}{\left( \sqrt{2 \pi \sigma^2} \right)^T} \exp \left( - \sum_{t=2}^T \frac{ \bigl( y_t - c - \phi y_{t-1} \bigr)^2}{2\sigma^2} - \frac{ \bigl( y_1 - c/ (1 - \phi) \bigr)^2}{2\sigma^2/ (1 - \phi^2)}\right)$ 実現値 $(y_1, \cdots, y_T)$ からパラメータ $(c, \phi, \sigma^2)$ を推定するには上の同時分布が最大になるような $(c, \phi, \sigma^2)$ を選べばいいわけです（実際にはこれの対数を取って最大化しますが）。これを最尤推定といいます。

119~122ページでは上の同時分布を求める別解として、「 $y \equiv (y_1, \cdots, y_T)^{\top}$ はどんな多変量正規分布にしたがうだろうか？」というアプローチで同じ同時分布を導きました。

122~123ページでは「 $y_1$ が所与」という考えやすい場合で実際に最尤推定値を求めますが、前回は流してしまったので改めて追ってみます。 $y_1$ が所与なら上の同時分布は以下のように化けますね。

定常ガウシアン AR(1) 過程の実現値の同時分布（ $y_1$ が所与の場合）

$\displaystyle f_{Y_T, \cdots, Y_1}(y_T, \cdots, y_2 \, ; y_1, c, \phi, \sigma^2) = \frac{1}{\left( \sqrt{2 \pi \sigma^2} \right)^{T-1}} \exp \left( - \sum_{t=2}^T \frac{ \bigl( y_t - c - \phi y_{t-1} \bigr)^2}{2\sigma^2} \right)$ 対数を取るとこうですね。 $\displaystyle \mathcal{L}(c, \phi, \sigma^2) \equiv \log f_{Y_T, \cdots, Y_1}(y_T, \cdots, y_2 \, ; y_1, c, \phi, \sigma^2) = - \frac{T-1}{2} \log \left( 2 \pi \sigma^2 \right) - \sum_{t=2}^T \frac{ \bigl( y_t - c - \phi y_{t-1} \bigr)^2}{2\sigma^2}$ これを $(c, \phi, \sigma^2)$ について最大化したいわけですから、これらで偏微分してみましょう。 $\displaystyle \frac{\partial \mathcal{L}(c, \phi, \sigma^2)}{\partial c} = \sum_{t=2}^T \frac{ 2 \bigl( y_t - c - \phi y_{t-1} \bigr)}{2\sigma^2}$ $\displaystyle \frac{\partial \mathcal{L}(c, \phi, \sigma^2)}{\partial \phi} = \sum_{t=2}^T \frac{ 2 y_{t-1} \bigl( y_t - c - \phi y_{t-1} \bigr)}{2\sigma^2}$ $\displaystyle \frac{\partial \mathcal{L}(c, \phi, \sigma^2)}{\partial \sigma^2} = - \frac{T-1}{2\sigma^2} + \sum_{t=2}^T \frac{ \bigl( y_t - c - \phi y_{t-1} \bigr)^2}{2\sigma^4}$ これらの式をイコールゼロとしたものを満たす $(c, \phi, \sigma^2)$ が最尤推定値 $(\hat{c}, \hat{\phi}, \hat{\sigma}^2)$ になりますが、まず上の2式より $(\hat{c}, \hat{\phi})$ が導かれ、次に最後の式より $\hat{\sigma}^2$ が導かれますね。

$y_1$ が所与の定常ガウシアン AR(1) 過程の最尤推定値 $(\hat{c}, \hat{\phi}, \hat{\sigma}^2)$

$\displaystyle \left(\begin{array}{cc} T-1 & \sum_{t=2}^T y_{t-1} \\ \sum_{t=2}^T y_{t-1} & \sum_{t=2}^T y_{t-1}^2 \end{array}\right) \left(\begin{array}{c} \hat{c} \\ \hat{\phi} \end{array}\right) = \left(\begin{array}{c} \sum_{t=2}^T y_{t-1} \\ \sum_{t=2}^T y_{t-1}y_t \end{array}\right)$ $\displaystyle \hat{\sigma}^2 = \frac{1}{T-1} \sum_{t=2}^T \bigl( y_t - \hat{c} - \hat{\phi} y_{t-1} \bigr)^2$ $y_1$ が所与なら最尤推定値を求めることは容易ですね。ただ、「言い換えると、これは最小2乗回帰（OLS）である」とありますが、どういうことなんでしょう？

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

いま、以下の前者の問からその最尤推定値にたどり着いたよね。でもそれって、以下の後者の問から出発したのと同じになってるんだよね。

$y_1$ が所与の定常ガウシアン AR(1) 過程の、実現値 $(y_2, \cdots, y_T)$ に対する最尤推定値は何だろうか？
$y_1$ が所与の AR(1) 過程の、実現値 $(y_2, \cdots, y_T)$ に対する最小2乗推定値は何だろうか？

だって $(\hat{c}, \hat{\phi})$ を求めるのに $\bigl( y_t - c - \phi y_{t-1} \bigr)^2$ を最小化しているけど、これは各 $y_t$ を $c + \phi y_{t-1}$ で予測したときの2乗誤差に他ならないし、 $\hat{\sigma}^2$ は残差の分散になっているしね。ちなみに、後者は定常性とガウス性を仮定していない。「誤差2乗和を最小化し、残差の分散を求める」という作業は定常性もガウス性も要求しないからね。

あっと、最小2乗推定と等価ということはわかりました。しかし、なら最初から最小2乗推定を紹介していただければよかったのではありませんか？定常性とガウス性を仮定して頑張って尤度関数を求めて、「実は単に最小2乗推定するのと同じだったんだけどね」では私が報われなくありませんか？

そんなことないよ。確かに最小2乗法はいつでもできるよ。最小2乗推定値自体はいつでも求まると思う。でも、それがまともな推定値になっていることを支持する材料ってないからね。定常性とガウス性をもつ場合には最小2乗推定値は最尤推定値になっていることがいまわかったけど、逆にいうとそうでない場合は変な値になっているかもしれないってことだからね。まあそう言い出すと「なら最尤推定値ならまともなのか」って話になるけど、単なる最小2乗推定値よりはまともだと信じられる場面が多いってことだと思ってる。

な、なるほど…？それで、123ページからは AR(p) 過程の場合ですね。AR(p) 過程の同時分布は…

まず「最初の $p$ 点 $(y_1, \cdots, y_p)^{\top}$ はどんな多変量正規分布にしたがうだろうか？」を考えた上で、
その後「 $(y_1, \cdots, y_p)$ を観測した上で次の $y_{p+1}$ はどんな分布だろう」を積み上げていっています。

AR(1) 過程でみた2つアプローチの合わせ技をやっていますね。素朴なアプローチではできないんでしょうか？

…5.3.8 式に AR(2) の場合が具体的に書き下してあるから、AR(2) で求めて検算してみればいいんじゃないかな…と思ったけど、素朴にはできない気がするかな。まず、何も観測していない下での最初の1点の分布 $f_{Y_1}(y_1)$ は AR(1) のときと一緒だよね。期待値と分散は AR(2) のそれにしなければならないけど。あと、直近の $p$ 点 $(y_1, y_2)$ を観測した下での $y_3$ の分布 $f_{Y_3 | Y_2, Y_1}(y_3 \, | y_2, y_1)$ も AR(1) のときと同じ要領だね。

$\displaystyle f_{Y_1}(y_1 \, ; c, \phi_1, \phi_2, \sigma^2) = \frac{1}{\sqrt{2 \pi \sigma^2/ (1 - \phi_1^2 - \phi_2^2)}} \exp \left( - \frac{ \bigl( y_1 - c/ (1 - \phi_1 - \phi_2) \bigr)^2}{2\sigma^2/ (1 - \phi_1^2 - \phi_2^2)}\right)$ $\displaystyle f_{Y_3 | Y_2, Y_1}(y_3 \, | y_2, y_1 \, ; c, \phi_1, \phi_2, \sigma^2) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( - \frac{ \bigl( y_3 - c - \phi_1 y_2 - \phi_2 y_1 \bigr)^2}{2\sigma^2}\right)$ でもこれらの中間の $f_{Y_2 | Y_1}(y_2 \, | y_1)$ がこんな感じじゃ求められない気がするんだよね。だって、「何も観測していない下での」でもないし、「直近の $p$ 点を観測した下での」でもないから、どっちのアプローチも使えない。結局 AR(2) の共分散構造を満たすように $(y_1, y_2)$ の同時分布を真面目に求めないといけない気がする。その分散共分散行列が125ページの一番上の式だと思う。

素朴にはいかなさそうなんですね。…あれ、126ページには「ガウス性をもたない場合の最尤推定」とありますね。5章ではガウス性をもたない場合は扱わないはずでは？

…ここにかいてあるのは、ガウス性をもたない過程にガウス性をもつと仮定した最尤推定をすることを「準最尤推定」ということ。生データがガウス性をもたない場合でもシンプルな変換でガウス性をもつようにできることがあること。とかかな。その変換の例が 126ページの下の方にある Box-Cox 変換だね。

この変換は $\lambda$ というパラメータを与えて用いるんですね。どのように決定すればいいんでしょうか。

無論ガウス分布に近くなるように決めるんだけど、scipy.stats.boxcox は $\lambda$ を与えなければ勝手に決めてくれるって。

scipy.stats.boxcox — SciPy v1.5.4 Reference Guide

文明の利器…127ページの一番上には非ガウスの場合の ARMA 過程の推定に関する研究が紹介されていますね。この本自体が17年前の本なのでそれを考慮して参考にしなければなりませんが…。5.4~5.6節は MA 過程、ARMA 過程の話ですが一旦とばします。5.6 節の最後には、ガウシアン ARMA 過程を最尤推定するのに最もシンプルなアプローチはカルマンフィルタであるとありますね。また、他の参考文献も紹介されています。

それで、尤度 $\mathcal{L}(\theta)$ を出す式が特定できた下で、それを最大化する $\hat{\theta}$ を如何に求めるかという話が残っていますが（ $y_1$ が所与の場合の最小2乗法はやりましたが）、5.7 節は数値的に最大化する方法ですね。紹介されているのは以下でしょうか。

グリッドサーチ
最急降下法
Newton-Raphson 法
Davidon-Fletcher-Powell 法（準 Newton 法）

それで、5.8 節が最尤推定を用いた統計的検定で、5.9 節の「不等式制約」というのは、パラメータ探索は「過程が定常である」「分散は正である」のような制約の下で行わなければならないので、そのために使えるテクニックですかね。

つづいたらつづく

2020-11-11

論文読みメモ： Trellis Networks for Sequence Modeling（その1）

論文読み

以下の論文を読みます。

Shaojie Bai, J. Z. Kolter, V. Koltun. Trellis Networks for Sequence Modeling. In 7th International Conference on Learning Representations, ICLR 2019, 2019. [1810.06682] Trellis Networks for Sequence Modeling

前回：雑記 - クッキーの日記

※ キャラクターは架空のものです。私の誤りは私に帰属します。お気付きの点がありましたらコメント等で指摘ください。

これは TCN を提案された方々の論文ですが、この論文では系列データを学習する新しいモデル構造である Trellis Networks（「格子垣」ネットワーク）を提案するようですが、この TrellisNet はアブストラクトによると

ちょっと特殊な Temporal Convolutional Networks (TCN) ともいえるし、
ちょっと特殊な Recurrent Neural Networks (RNN) ともいえる

とのことですね…しかし、RNN と TCN はどちらも系列データを学習するモデルですがその仕組みは全く異なりますよね？ RNN は現時点の入力と前の時点の出力から現時点の出力をつくりますし、TCN は現時点までの入力たちをどんどん畳込んで現時点の出力をつくりますし、そもそも下図をみても、どちらもお世辞にも「格子」にはみえないと思います。

f:id:cookie-box:20201010183627p:plain:w570

じゃあ RNN のようでもあって、TCN のようでもあって、格子にみえるモデルを考えてみたらいいんじゃないかな。そもそもどんなモデルだったら格子にみえるだろう？

どんなって、「格子」というからには…あれ、格子ってなんでしょう…あ、そうです！複素数平面の格子点です！つまり、格子点のように、何か横方向のインデックスと縦方向のインデックスをもつものがあれば、格子にみえるでしょう。むしろそれが格子の定義かもしれませんが…まあいいです。

それで、RNN や TCN には縦方向と横方向のインデックスをもつものがあるかという話です。この場合、横方向のインデックスは時間ステップでしょう。縦方向のインデックスは何層目の特徴かといったようなものが考えられると思います。入力を0層目、出力を1層目とすれば RNN も TCN も格子とはいえますね。でも、それではちょっとお粗末だと思います。それでは格子というより梯子ですから。

…あっ、でも、RNN も TCN も構造上2つ以上積み重ねることができますね。そうすれば、時刻 $t$ の $i$ 層目の出力が格子にみえると思います。

…いや、よく考えると、RNN に TCN を重ねることも、TCN に RNN を重ねることも妨げられるものではないですね。ステップ毎に出力を出すタイプの系列モデルであれば中身がどうあっても積み重ねらることができますから。であれば、TCN → RNN → TCN → RNN のような積み重ねだって可能ですし、かつ格子にみえます。TrellisNet の正体はもうこれでは？

いい線なんじゃないかな。でも TCN → RNN → TCN → RNN というモデルではないみたいだね。それよりは TCN → TCN → TCN → TCN が近いと思う。3ページ目の Figure 1 をみる限りは。たぶん TrellisNet には純粋な再帰路はない。

うおおおい！私には考えさせておいて自分はちゃっかりカンニングしてるじゃないですか！それに、TCN → TCN → TCN → TCN だと RNN 要素なくないですか？ TrellisNet はちょっと特殊な RNN なんでしょう？

どう RNN なのか解説してくれるのが 4.2 節なんだと思う。

…おとなしくイントロダクションから読んでいきます。って、最初の段落、どばっと文献の引用がありますが、個々の内容には何も言及されていないですね…。近年の系列データの学習の研究では、RNN も改良が続けられているし、TCN も高いパフォーマンスを発揮しているし、セルフアテンションに基づくネットワークも快進撃をみせているということですが…セルフアテンションって RNN や TCN の並びに入るんですか？ RNN や TCN は将来のデータをカンニングしてはならないという縛りがあるはずですが、文章に対する品詞分類や翻訳などのタスクでは、ふつう将来のデータ（単語）をカンニングしてもいいですよね？むしろカンニングしろって感じですが…。

セルフアテンションが RNN や TCN の並びに入るかということについては、セルフアテンション層は未来の単語も過去の単語も参照するということはしていないと思う。だから分の途中まででも出力は同じだよ。Positional Encoding とかあるから文章中の座標は単語に貼り付けてある必要はあるし、ネットワークを学習するときに次の単語を説明しようとしたり前の単語の説明しようとしたりということをしないということではないよ。

[1706.03762] Attention Is All You Need

なるほど？ …イントロダクションの3段落目には、アブストラクトにあったように、TrellisNet はちょっと特殊な TCN であるということがもう少し詳しくかかれているようですね。

そこを読む限り、TrellisNet は TCN → TCN → TCN → TCN （多層 TCN）なんだけど、全ての TCN は重みを共有しているし、どの層にも入力データが injection されるんだね（右図）。なるほど格子にみえる。 $z_t^{(0)}$ はゼロベクトルね。

でも「これをみせられてもよい性能をもつとは思われないでしょう」と。まあこれだけぱっとみせられても何がいいのか、全ての TCN が重みを共有している意義は何なのかわからないね。でもこれが記憶をトランケートされた RNN の一般化になっていることをこの論文で示すらしい。それが優れている根拠になっているといいたいのかな。

「RNN とみなせるから優れているのだ」ですか。RNN に一定の信頼があるのは確かでしょうが。イントロダクションの5段落目にあるのは TrellisNet の性能を検証したタスクですね。おそらく先の先行研究にならって選定されていると思われます。個々のタスクの出自やデータサイズなんかの詳細は Appendix C にあるようです。

文字レベル Penn Treebank（PTB）
単語レベル Penn Treebank（PTB）
単語レベル WikiText-103（WT103） ― 大規模データとして採用されているようですね。
sequential MNIST
permuted MNIST
sequential CIFAR-10

そしてこのタスクのどれもで SOTA を達成したようです（厳密には WT103 については同時期に出た論文の transformer モデルに劣っているようですが）。これまでそれぞれ異なるモデルが SOTA を達成していた6タスクを TrellisNet が一気に塗り替えたのだというのが強調されていますね。

ところで sequential MNIST 以下はもう「系列モデルに対するスタンダードなストレステスト」とされていますが、手書き数字や写真って本来全てのピクセルを同時に認識でき、縦横のつながりも明らかなものですよね。明らかに CNN さんに学習させるべきですよね。それを出鱈目な順序でピクセルごとに流すような、無駄な曲芸のようなタスクを系列モデルにさせるのってパワハラに当たらないんですか？

TrellisNet
f:id:cookie-box:20201011140713p:plain:w200

TrellisNet（乙）は私たち（甲）との間に「乙は甲に自然に系列としての構造をもつデータを学習する労務を提供すること」などという雇用契約は結んでいないからパワハラには当たらないね。雇用契約にない仕事をさせるのが問題になるからね。

いや法律的な見解ではなく、本来は2次元としての構造をもつデータがベンチマークとなっているのはどうなのかといいたかったんですが…。2節の「背景」は短いですね。かいてある内容としては、

再帰的ニューラルネットって自然言語以外にも色々なドメインに利用されているよね（という感じでまたどばっと引用だけありますね）。
畳込みニューラルネットも長期記憶が必要なタスクで優れているんだよね（ここでこの方たちの TCN を含んでどばっと引用がありますね）。
self-attention もあるよね（どばっと略）。
今回は RNN と TCN の架け橋となるモデルを提案するよ。

といった感じでまあ何にせよ RNN と TCN を組み合わせるといった感じですが、RNN と TCN を組み合わせた先行研究も引用されていますね。ただどのようなものだったのか詳しくわかりませんね。TrellisNet の方が両者をより深く結びつけたのだということではあるんでしょうが。

「convolutional LSTMs combine convolutional and recurrent units」というのは RNN ユニットと CNN ユニットを組み合わせたというようにみえるけどね。なら RNN の結果を CNN したり、CNN の結果を RNN したり、それこそ TCN → RNN → TCN → RNN のようなものなのかな。知らないけど。

3節の SEQUENCE MODELING AND TRELLIS NETWORKS に進みます。冒頭では「系列モデルとは長さ T の系列を受け取って長さ T の系列を出すものであって t 番目の出力は未来の入力に依存してはいけないものとする」という TCN の論文同様の定義が示されています。その後にはもう TrellisNet が以下のように定式化されていますね。

系列モデルの入力を $x_{1:T} = (x_1, \cdots, x_T) \in \mathbb{R}^{T \times p}$ とする。
時刻 $t$ の $i$ 層目の出力を $z_{t}^{(i)} \in \mathbb{R}^{q}$ とする（つまり何層目の出力であっても $q$ 次元である）。 $z_{t}^{(0)}$ はゼロベクトルとする。
TrellisNet の $z_{t+1}^{(i + 1)} \in \mathbb{R}^{q}$ のつくり方としては、まず $\tilde{z}_{t+1}^{(i + 1)} \in \mathbb{R}^{r}$ を下図のように線形結合で生成し（Conv1D とも解釈できますね）、 $\tilde{z}_{t+1}^{(i + 1)}$ と $z_{t}^{(i)}$ に適当な非線形関数 $f$ を適用して $z_{t+1}^{(i + 1)}$ とする。

f:id:cookie-box:20201011173741p:plain:w360

$W_1$ と $W_2$ はすべての $i$ 層目で共通ですがら、 $x_t$ と $x_{t+1}$ を変換した $\tilde{x}_{t+1} = W_1^{x} x_t + W_2^{x} x_{t+1}$ は使い回せるとありますね（ $W_j^x \in \mathbb{R}^{r \times p}$ ）。

…いっていることはわかりますが、腑に落ちませんね。特に最終出力 $z_{t+1}^{(i + 1)} = f \bigl( \tilde{z}_{t+1}^{(i + 1)}, z_{t}^{(i)} \bigr)$ について、なぜ $z_{t}^{(i)}$ を使うのか、 $f$ は具体的にどういった関数なのかという点です。 $f$ がわからなければ TrellisNet が何をするモデルなのかイメージすることは難しいと思います。もっといえばすべての層でウェイトを共有する理由だってわかりません。本文中には Vogel & Pock (2017) が画像処理で層方向のウェイトの共有をやっているとありますがだから何やら。

3節の最後に LSTM セルに基づく活性化をするとあるよ。4節にその理論的根拠があるとあるけど、先に活性化の式がある 5.1 節をみてみようか。気になるしね。

LSTM で活性化？ LSTM は活性化関数ではないと思うんですが…5.1 節をみると、ゲート付き活性化というのがそもそもあるんですね？ 5.1 節を読む限り、その LSTM に基づくゲート付き活性化というのは以下になるはずです。

f:id:cookie-box:20201011230518p:plain:w560

えっと、これは何なんでしょう。何やら LSTM めいてはいますが…。

…もし「前回の入力」と「前回の前の層の出力」を使用しないなら、これは時刻 i から時刻 i+1 に進む LSTM だよね。ただし、実際には i とは何層目かであって時刻を表すインデックスではないけど。それに、実際には「前回の入力」と「前回の前の層の出力」が混ぜ込まれるし、いわば「入力とゲートに混ぜ物があり、他人の記憶を植え付けられた LSTM」かな。微塵も LSTM ではないな…。

読み飛ばした 4 節に戻りますね。4.1 節は、TrellisNet と TCN の共通点と相違点について、「各層で未来のデータを参照せず畳込みを行うこと」「多層化するほど参照範囲が広がること」の2点は TCN と同様ですが、「各層の重みが共有されていること」「各層に入力データがインジェクションされること」は TCN と異なるといっていますね。それはまあそうですが。「各層の重みが共有されていること」は学習の安定化とモデルサイズの節約に寄与する正則化になっているとのことです。

その2があればつづく