雑記

私の誤りは私に帰属します。お気付きの点がありましたらご指摘いただけますと幸いです。
  1. Haoyi Zhou, Shanghang Zhang, Jieqi Peng, Shuai Zhang, Jianxin Li, Hui Xiong, Wancai Zhang. Informer: Beyond Efficient Transformer for Long Sequence Time-Series Forecasting. arXiv preprint arXiv:2012.07436, 2020.

    [2012.07436] Informer: Beyond Efficient Transformer for Long Sequence Time-Series Forecasting
    GitHub - zhouhaoyi/Informer2020: The GitHub repository for the paper "Informer" accepted by AAAI 2021.

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

Transformer は行列 {\rm Softmax}(QK^\top /\sqrt{d}) によって入力の特徴系列を「前後の文脈を考慮した」特徴系列に変えますが、系列長 L が長くなるほどこの行列を用意して適用するのに \mathcal{O}(L^2) の計算量がかかってしまうのですよね? 文脈を考慮する先が増えるわけですから当然ですが。しかし、それでも何とか計算を間引いて Transformer に長い系列を流したいこともあるわけです。そのような手法としてこれまでどのようなものが提案されてきたのでしょう? それぞれの手法にどのような得手不得手があるのでしょうか?

Transformer の1ヘッドがやること
  •  x \in \mathbb{R}^{L \times d_{\rm in}} を受け取る.
  •  Q = x W_Q \in \mathbb{R}^{L \times d} K = x W_K \in \mathbb{R}^{L \times d} V = x W_V \in \mathbb{R}^{L \times d_{\rm out}} を計算する.
  •  y = {\rm Softmax}(QK^\top /\sqrt{d}) V \in \mathbb{R}^{L \times d_{\rm out}} を計算する(※ Softmax は各行を正規化する).

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

まず Informer [1] の間引き方を復習しようか。考え方としては、{\rm Softmax}(QK^\top /\sqrt{d}) の各行の中で一様分布から逸脱している u 行を残したいんだね。逸脱度の指標は一様分布とその行との交差エントロピーだね(論文では KL 情報量といっているけど定数項を省いているから結局交差エントロピーなんだよね)。

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

なるほど…一様分布と i 行目の交差エントロピーは、i 行目の確率分布に最適化した情報量の尺度で一様分布を受け取ったときの平均的な情報量…Qi 行目 q_iKj 行目 k_j を縦ベクトルと考えると、

\displaystyle {\rm CE}_i = \frac{1}{L} \sum_{j=1}^L \left[ \log \frac{\sum_{j'=1}^L \exp(q_i \cdot k_{j'} / \sqrt{d})}{\exp(q_i \cdot k_j / \sqrt{d})} \right] = \log \left[ \sum_{j=1}^L \exp \left( \frac{q_i \cdot k_j}{\sqrt{d}} \right) \right] - \frac{1}{L} \sum_{j=1}^L \frac{q_i \cdot k_j}{\sqrt{d}}
ですね。すべての i についてこの  {\rm CE}_i を計算してこれが大きい u 行を残せば計算量を節約でき…って、これを計算するのにすべての q_i \cdot k_j が必要ですよね? これでは計算量 \mathcal{O}(L^2) じゃないですか…。

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

うんだから厳密に  {\rm CE}_i を計算するんじゃないんだよね。まず、その式から以下が成り立つことはわかるよね。下限の方は交差エントロピーが最小になるのは i 行目の確率分布が一様分布に一致するときということからわかるし、上限はその行の最大値をつかって抑えただけだし。それでこの最右辺第2項と第3項を  \overline{\rm CE}_i とおくみたい。

\displaystyle \log L \leqq {\rm CE}_i \leqq \log L + \max_j \left[ \frac{q_i \cdot k_j}{\sqrt{d}} \right] - \frac{1}{L} \sum_{j=1}^L \frac{q_i \cdot k_j}{\sqrt{d}}
k_j が多変量正規分布にしたがうと仮定すると、 \overline{\rm CE}_i > \overline{\rm CE}_{i'} かつ i 行目の分散の方が i' 行目の分散より大きいなら高確率で  {\rm CE}_i > {\rm CE}_{i'} になるんだって(正確には「となるような q_i, q_{i'} の領域が存在する」といっているけど)。

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

…うーんでも、その主張だと結局  {\rm CE}_i{\rm CE}_{i'} を比較するのに QK^\topi 行目と i' 行目が完全に必要になりませんか? それでは意味が―

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

以下の実装をみると QK^\top を完全に計算しているわけではないね。

https://github.com/zhouhaoyi/Informer2020/blob/main/models/attn.py#L47-L68

QK^\top を計算するんじゃなくて、K のうち sample_k 行だけ残した K' で代用している。いや、正確にいうと、Q の各行に対して sample_k 行の選び方を変えているから、L 通りの {K'}_{1}, {K'}_{2}, \cdots, {K'}_{L} を用意して、以下の QK^\top もどきを計算して、これに基づいて  \overline{\rm CE}_i が大きい n_top 行を選んでいるね(※ この記事では Qi 行目を縦ベクトル q_i と考えているから横ベクトルに戻すために転置しているよ)。
 \left( \begin{array}{c} q_1^\top {K'}_{1}^\top \\ q_2^\top {K'}_{2}^\top \\ \vdots \\ q_L^\top {K'}_{L}^\top \end{array} \right)
それで、一度 n_top 行を選んだら改めてそれらの行に対しては QK^\top をきちんと計算しているよ。実装を抜き出すと以下だね。

Informer のセルフアテンションの間引き方

f:id:cookie-box:20210419234518p:plain:w680
系列長 8 で間引きを発生させるには factor をデフォルト値の 5 より下げないといけない。factor を 2 にすると n_top = 2 * ceil( log(8) ) = 6 になるから、セルフアテンションの計算対象が一様分布から逸脱した 6 行に絞られるよ。上図では 2, 3 行目の計算が間引かれたけど、この 2 行が本当に一様分布に最も近い 2 行だったかは確認していない。

つづいたらつづく

雑記

クロネッカー積と行列積の混合積の公式  (A \otimes B)(C \otimes D) = AC \otimes BD を成分でかいただけです.

まず,クロネッカー積は以下のように定義されると思います.

定義1(クロネッカー積)

A \in \mathbb{R}^{n_a \times m_a}, \; B \in \mathbb{R}^{n_b \times m_b} に対して  A \otimes B \in \mathbb{R}^{n_a n_b \times m_a m_b} を以下のように定義する.

A \otimes B = \left( \begin{array}{cccc} a_{1,1}B & a_{1,2}B & \cdots & a_{1,m_a}B \\ a_{2,1}B & a_{2,2}B & \cdots & a_{2,m_a}B \\ \vdots & \vdots & \ddots & \vdots \\ a_{n_a,1}B & a_{n_a,2}B & \cdots & a_{n_a,m_a}B \end{array} \right)

それで以下の公式が成り立ちます. AC \otimes BD から A, B を追放したいときなどにこれをつかうと思います.

命題1(クロネッカー積と行列積の混合積の公式)

A \in \mathbb{R}^{n_a \times m_a}, \; B \in \mathbb{R}^{n_b \times m_b}, \; C \in \mathbb{R}^{m_a \times m_c}, \; D \in \mathbb{R}^{m_b \times m_d} に対して以下が成り立つ.

 (A \otimes B)(C \otimes D) = AC \otimes BD

この公式は以下のように定義にあてはめれば示せます.

命題1の証明0(ブロックごとに行列積をとる)

\begin{split} (A \otimes B)(C \otimes D) &= \left( \begin{array}{cccc} a_{1,1}B & a_{1,2}B & \cdots & a_{1,m_a}B \\ a_{2,1}B & a_{2,2}B & \cdots & a_{2,m_a}B \\ \vdots & \vdots & \ddots & \vdots \\ a_{n_a,1}B & a_{n_a,2}B & \cdots & a_{n_a,m_a}B \end{array} \right)  \left( \begin{array}{cccc} c_{1,1}D & c_{1,2}D & \cdots & c_{1,m_c}D \\ c_{2,1}D & c_{2,2}D & \cdots & c_{2,m_c}D \\ \vdots & \vdots & \ddots & \vdots \\ c_{m_a,1}D & c_{m_a,2}D & \cdots & c_{m_a,m_c}D \end{array} \right) \\ &= \left( \begin{array}{cccc} \sum_{k=1}^{m_a} a_{1,k} c_{k,1} BD & \sum_{k=1}^{m_a} a_{1,k} c_{k,2} BD & \cdots & \sum_{k=1}^{m_a} a_{1,k} c_{k,m_c} BD \\ \sum_{k=1}^{m_a} a_{2,k} c_{k,1} BD & \sum_{k=1}^{m_a} a_{2,k} c_{k,2} BD & \cdots & \sum_{k=1}^{m_a} a_{2,k} c_{k,m_c} BD \\ \vdots & \vdots & \ddots & \vdots \\ \sum_{k=1}^{m_a} a_{n_a,k} c_{k,1} BD & \sum_{k=1}^{m_a} a_{n_a,k} c_{k,2} BD & \cdots & \sum_{k=1}^{m_a} a_{n_a,k} c_{k,m_c} BD \end{array} \right) \\ &= AC \otimes BD \end{split}


完全に上の証明でいいんですが,より明示的に成分が等しいことに基づく証明もしておきたい気がします.なので,クロネッカー積の定義を成分にかき直します.ちなみに後で気付いたんですが,ここで「これ以降行列のインデックスをゼロ始まりでとることに決める」と宣言しておけばよかったです.以降では右肩に (0) を付けたときにゼロ始まりのインデックスでとるとしていますが証明では全部右肩に (0) が付いています.無駄です.

定義1'(クロネッカー積の成分;通常のインデックスでアクセス)

x を超えない最大の整数を {\rm floor}(x)xy で割った余りを {\rm mod}(x,y) とかくと,A \otimes Bi, \, j 成分は以下のようにかける(  1 \leqq i \leqq n_a n_b 1 \leqq j \leqq m_a m_b ).

(A \otimes B)_{i,j} = a_{\, {\rm floor} ( (i-1) /n_b)+1, \, {\rm floor} ( (j-1) /m_b )+1 \, } b_{\, {\rm mod}( (i-1), n_b) + 1, \, {\rm mod}((j-1), m_b) + 1\, }

以降,表記の便利のため,右肩に (0) を付けたときはゼロ始まりのインデックスで成分をとることにする.すると以下のようになる(  0 \leqq i' \leqq n_a n_b - 1 0 \leqq j' \leqq m_a m_b- 1 ).

(A \otimes B)^{(0)}_{i',j'} = {a}^{(0)}_{\, {\rm floor} (i'/n_b), \, {\rm floor} ( j'/m_b ) \, } {b}^{(0)}_{\, {\rm mod}(i', n_b), \, {\rm mod}(j', m_b) \, }

上の定義でクロネッカー積の成分にアクセスできると思うんですが,{\rm floor}{\rm mod} が多用されていてみづらいという声もあると思うので最初から {\rm floor}{\rm mod} の2重インデックス(行と列があるので4重インデックス)でアクセスする方法も用意します.

定義1''(クロネッカー積の成分;大インデックスと小インデックスでアクセス)

A \otimes Bn_b I' +i, \, m_b J' + j 成分は以下のようにかける(  0 \leqq I' \leqq n_a - 1 0 \leqq J' \leqq m_a - 1 1 \leqq i \leqq n_b 1 \leqq j \leqq m_b ).

(A \otimes B)_{\, n_b I' + i, \, m_b J' + j \,} = a_{\, I' + 1, \, J' + 1 \, } b_{\, i, \, j \, }

ゼロ始まりのインデックスだと以下のようになる(  0 \leqq i' \leqq n_b - 1 0 \leqq j' \leqq m_b- 1 ).

(A \otimes B)^{(0)}_{\, n_b I' + i', \, m_b J' + j' \,} = a^{(0)}_{\, I', \, J' \, } b^{(0)}_{\, i', \, j' \, }

それで定義1' のアクセス方法で命題1を証明すると以下になります.{\rm floor}{\rm mod} のせいで横幅が長くなってしまうことがわかります.

命題1の証明1(定義1' を使用)

 0 \leqq i' \leqq n_a n_b - 1 0 \leqq j' \leqq m_c m_d - 1 について,以下が成り立つ.

 \begin{split} \displaystyle \bigl( (A \otimes B)(C \otimes D)\bigr)^{(0)}_{i',j'} &= \sum_{k'=0}^{m_a m_b - 1}  {a}^{(0)}_{\, {\rm floor} (i'/n_b), \, {\rm floor} ( k'/m_b ) \, } {b}^{(0)}_{\, {\rm mod}(i', n_b), \, {\rm mod}(k', m_b) \, } {c}^{(0)}_{\, {\rm floor} (k'/m_b), \, {\rm floor} ( j'/m_d ) \, } {d}^{(0)}_{\, {\rm mod}(k', m_b), \, {\rm mod}(j', m_d) \, } \\ &= \sum_{k'=0}^{m_a m_b - 1} \Bigl( {a}^{(0)}_{\, {\rm floor} (i'/n_b), \, {\rm floor} ( k'/m_b ) \, } {c}^{(0)}_{\, {\rm floor} (k'/m_b), \, {\rm floor} ( j'/m_d ) \, }  \Bigr) \Bigl( {b}^{(0)}_{\, {\rm mod}(i', n_b), \, {\rm mod}(k', m_b) \, } {d}^{(0)}_{\, {\rm mod}(k', m_b), \, {\rm mod}(j', m_d) \, } \Bigr) \\ &= \sum_{k'=0}^{m_b - 1} \Bigl( {a}^{(0)}_{\, {\rm floor} (i'/n_b), \, 0 \, } {c}^{(0)}_{\, 0, \, {\rm floor} ( j'/m_d ) \, }  \Bigr) \Bigl( {b}^{(0)}_{\, {\rm mod}(i', n_b), \, k' \, } {d}^{(0)}_{\, k', \, {\rm mod}(j', m_d) \, } \Bigr) \\ & \quad \, + \sum_{k'=0}^{m_b - 1} \Bigl( {a}^{(0)}_{\, {\rm floor} (i'/n_b), \, 1 \, } {c}^{(0)}_{\, 1, \, {\rm floor} ( j'/m_d ) \, }  \Bigr) \Bigl( {b}^{(0)}_{\, {\rm mod}(i', n_b), \, k' \, } {d}^{(0)}_{\, k', \, {\rm mod}(j', m_d) \, } \Bigr) \\ & \quad \, + \cdots + \sum_{k'=0}^{m_b - 1} \Bigl( {a}^{(0)}_{\, {\rm floor} (i'/n_b), \, m_a-1 \, } {c}^{(0)}_{\, m_a-1, \, {\rm floor} ( j'/m_d ) \, }  \Bigr) \Bigl( {b}^{(0)}_{\, {\rm mod}(i', n_b), \, k' \, } {d}^{(0)}_{\, k', \, {\rm mod}(j', m_d) \, } \Bigr) \\ &= \left( \sum_{k''=0}^{m_a - 1} {a}^{(0)}_{\, {\rm floor} (i'/n_b) , \, k'' \,} {c}^{(0)}_{\, k'', \, {\rm floor} ( j'/m_d ) \,} \right) \left( \sum_{k'=0}^{m_b - 1} {b}^{(0)}_{\, {\rm mod}(i', n_b), \, k' \,} {d}^{(0)}_{\, k', \, {\rm mod}(j', m_d) \,} \right) \\ &= (AC)^{(0)}_{\, {\rm floor} (i'/n_b), \, {\rm floor} ( j'/m_d ) \, } (BD)^{(0)}_{\, {\rm mod}(i', n_b), \, {\rm mod}(j', m_d) \, } \\ &= (AC \otimes BD)^{(0)}_{i',j'} \end{split}

以上より, (A \otimes B)(C \otimes D) = AC \otimes BD である.□

他方,定義1'' のアクセス方法で命題1を証明すると以下になります.横幅が短いです.また,後から気付いたんですが以下のように右辺から出発した方が(2重ループを1重ループにまとめる方向になるので)すっきりする気がします.

命題1の証明2(定義1'' を使用)

 0 \leqq I' \leqq n_a - 1 0 \leqq J' \leqq m_c - 1 0 \leqq i' \leqq n_b-1 0 \leqq j' \leqq m_d-1 について,以下が成り立つ.

 \begin{split} \displaystyle (AC \otimes BD)^{(0)}_{\, n_b I' +i', \, m_d J' + j' \,} &= \left( \sum_{k=0}^{m_a - 1} a^{(0)}_{\, I', \, k \,} c^{(0)}_{\, k, \, J' \,} \right) \left( \sum_{k'=0}^{m_b - 1} b^{(0)}_{\, i', \, k' \,} d^{(0)}_{\, k', \, j' \,} \right) \\ &= \sum_{k=0}^{m_a - 1} \sum_{k'=0}^{m_b - 1} \Bigl( a^{(0)}_{\, I', \, k \,}  b^{(0)}_{\, i', \, k' \,} \Bigr) \Bigl( c^{(0)}_{\, k, \, J'\,} d^{(0)}_{\, k', \, j' \,} \Bigr) \\ &= \sum_{k=0}^{m_a - 1} \sum_{k'=0}^{m_b - 1}  (A \otimes B)^{(0)}_{\, n_b I' +i', \, m_b k + k' \,} (C \otimes D)^{(0)}_{\, m_b k +k', \, m_d J' + j' \,}  \\ &= \sum_{K=0}^{m_a m_b - 1}   (A \otimes B)^{(0)}_{\, n_b I' +i', \, K \,} (C \otimes D)^{(0)}_{\, K, \, m_d J' + j' \,} \\ &= \bigl( (A \otimes B) (C \otimes D) \bigr)^{(0)}_{\, n_b I' +i', \, m_d J' + j' \,} \end{split}

以上より, (A \otimes B)(C \otimes D) = AC \otimes BD である.□

おまけ

クロネッカー積と行列積の混合積のイメージ図を描きました.クロネッカー積をとることは洗剤を指定の濃さに塗りながらタイルを敷き詰めることに似ていると思ったのですが,いうほど洗剤を塗りながらタイルを敷き詰めることをするかというとしないので完全に駄目です.
f:id:cookie-box:20210329190053p:plain

Time Series Analysis: ノート10章 弱定常ベクトル過程(その2)

以下の本の10章を読みます。私の誤りは私に帰属します。お気付きの点がありましたらご指摘いただけますと幸いです。

Time Series Analysis

Time Series Analysis

その他の参考文献
  1. 経済・ファイナンスデータの計量時系列分析 (統計ライブラリー) | 竜義, 沖本 |本 | 通販 | Amazon
関連記事
f:id:cookie-box:20190101155733p:plain:w60

結局前回は以下のことしかやっていませんね…。

  • p 次ベクトル自己回帰過程 VAR(p) を導入してその定常条件を確認した。
  • 定常 VAR(p) は Vector MA(∞) にかき直せる。
  • 一般に定常ベクトル過程の自己共分散行列 \Gamma_{\nu} = E \bigl[ ({y}_t - \mu)({y}_{t - \nu} - \mu)^\top \bigr] \Gamma_\nu^\top = \Gamma_{-\nu} となる。
では、次は VAR(p) の自己共分散行列  \Gamma_\nu が具体的にどうなるかみていくのでしょうか?

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

いいや、ここで改めて q 次ベクトル移動平均過程 Vector MA(q)  y_t = \mu + \varepsilon_{t} + \Psi_1 \varepsilon_{t-1} + \cdots + \Psi_q \varepsilon_{t-q} を導入しているね。その自己共分散行列をみている。

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

え、ここで Vector MA(q) を導入するんですか? 既に前回 Vector MA(∞) がフライングで出てきていましたが…まあいいです、試しに Vector MA(2) の自己共分散行列を計算してみましょう。

\begin{split} \Gamma_{\nu} &= E \bigl[ ({y}_t - \mu)({y}_{t - \nu} - \mu)^\top \bigr] \\&= E \bigl[ ( \varepsilon_{t} + \Psi_1 \varepsilon_{t-1} + \Psi_2 \varepsilon_{t-2} )(\varepsilon_{t- \nu} + \Psi_1 \varepsilon_{t- \nu-1} + \Psi_2 \varepsilon_{t- \nu-2} )^\top \bigr] \\&= E \bigl[ ( \varepsilon_{t} + \Psi_1 \varepsilon_{t-1} + \Psi_2 \varepsilon_{t-2} )(\varepsilon_{t- \nu}^\top + \varepsilon_{t- \nu-1}^\top \Psi_1^\top + \varepsilon_{t- \nu-2}^\top \Psi_2^\top ) \bigr] \end{split}
これは…ベクトルホワイトノイズ \varepsilon_t は時刻の添字が一致するときに限り共分散が零行列でない \Omega となりますから…結局以下になりますね。|\nu| > 2 では \Gamma_\nu = 0 です。
\begin{split} \Gamma_{-2} &= \Omega \Psi_2^\top \\ \Gamma_{-1} &= \Omega \Psi_1^\top + \Psi_1 \Omega\Psi_2^\top \\ \Gamma_0 \, \, &= \Omega \quad  \, \,+ \Psi_1 \Omega\Psi_1^\top + \Psi_2 \Omega\Psi_2^\top  \\ \Gamma_1 \, \, &= \quad  \, \, \qquad \Psi_1 \Omega \quad  \, \, + \Psi_2 \Omega\Psi_1^\top \\ \Gamma_2 \, \, &= \quad  \, \, \qquad \qquad \qquad \; \; \Psi_2 \Omega \end{split}
\Gamma_\nu\nu ステップ昔の自分との共分散でした。Vector MA(2) は現在のノイズ、1ステップ昔のノイズ、2ステップ昔のノイズからなります。\Gamma_2 = \Psi_2 \Omega は、「現在の自分の2ステップ昔のノイズが、2ステップ昔の自分の現在のノイズと相関する項」ですね。\Gamma_{-2} = \Omega \Psi_2^\top は、「現在の自分の現在のノイズが、2ステップ未来の自分の2ステップ昔のノイズと相関する項」です。そして3ステップ以上ずれた自分とは相関がなくなります。

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

まあそうなるよね。じゃあ Vector MA(∞) の自己共分散行列はどうなる? ただし係数行列 \Psi_{\nu} は絶対総和可能とするよ(前回明示的に「成分ごとに絶対総和可能」とかいたけど、以降テキストの表記にならって「成分ごとに」を省くよ)。ちなみに単変量の MA(∞) は係数が絶対総和可能であるという条件の下で定常過程だったね(52ページ)。

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

Vector MA(∞) とはいきなり次数が増えましたね…そうなると何ステップずれても相関はありますよね…。まあ順に考えると、

  • 0次の自己共分散行列は  \Gamma_0 = \Omega + \sum_{\nu' = 1}^{\infty} \Psi_{\nu'} \Omega\Psi_{\nu'}^\top になりますね。
  • 1次の自己共分散行列は  \Gamma_1 = \Psi_1 \Omega + \sum_{\nu' = 1}^{\infty} \Psi_{\nu' + 1} \Omega\Psi_{\nu'}^\top になります。
  • 2次の自己共分散行列も  \Gamma_2 = \Psi_2 \Omega + \sum_{\nu' = 1}^{\infty} \Psi_{\nu' + 2} \Omega\Psi_{\nu'}^\top になります。
  • では、-1次であったら  \Gamma_{-1} = \Omega \Psi_1 + \sum_{\nu' = 1}^{\infty} \Psi_{\nu'} \Omega\Psi_{\nu' + 1}^\top になります。
  • -2次であったら  \Gamma_{-2} = \Omega \Psi_2 + \sum_{\nu' = 1}^{\infty} \Psi_{\nu'} \Omega\Psi_{\nu' + 2}^\top になります。
…このようになるはずですが、これらの自己共分散行列は値をもつのでしょうか? 係数行列 \Psi_{\nu} が絶対総和可能であるといっても、ただちにこれらが有限の値をもつことにはなりませんよね? もし発散してしまったら Vector MA(∞) は定常過程にならないのでは?

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

\Gamma_s は有限の値をもつよ。もっというと絶対総和可能になる。証明は Appendix にあるけど、いきなり \Gamma_si, j 成分 \gamma_{i,j}^{(s)} が有限であることを考えるのはややこしいから、まず以下の2つの積の期待値を考えよう。あ、以下 \Psi_0 \equiv I_n とするよ(部長がかき下した自己共分散行列もこう定義すればもっとすっきりする)。

  • z_t(i,l) = \sum_{\nu= 0}^{\infty} \psi_{i, l}^{(\nu)} \varepsilon_{t-\nu, l}y_t の第 i 成分の、ノイズの第 l 成分に依存する項。
  • z_{t-s}(j,m) = \sum_{\nu= 0}^{\infty} \psi_{j, m}^{(\nu)} \varepsilon_{t-s-\nu, m}y_{t-s} の第 j 成分の、ノイズの第 m 成分に依存する項。
これらの積の期待値をつかって自己共分散行列の成分が \gamma_{i,j}^{(s)} = \sum_{l=1}^n \sum_{m=1}^n E\Bigl[ z_t(i,l)z_{t-s}(j,m) \Bigr] とかけるからこれが絶対総和可能であることを示そう。

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

それらの積の期待値ですか…やはりベクトルホワイトノイズの定義から異なるステップのノイズの成分どうしの積の期待値は0ですから、\sum_{\nu= 0}^{\infty} \psi_{i, l}^{(\nu + s)} \psi_{j, m}^{(\nu)} \omega_{l, m} でしょうか。ベクトルホワイトノイズの自己共分散行列 \Omega の成分を \omega_{l, m} としました。これが絶対総和可能かといわれると、\sum_{s= 0}^{\infty} \sum_{\nu= 0}^{\infty} | \psi_{i, l}^{(\nu + s)} \psi_{j, m}^{(\nu)} \omega_{l, m} | \leqq |\omega_{l, m}| \sum_{\nu= 0}^{\infty} \sum_{s= 0}^{\infty} |\psi_{i, l}^{(\nu + s)}| |\psi_{j, m}^{(\nu)} | < \infty より絶対総和可能ですね…であれば各項は有限です…では、Vector MA(∞) は係数行列が絶対総和可能であるという条件の下で定常(t に依存せず s のみに依存する有限の自己共分散行列をもつ)というわけですか。

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

ちなみに同じ要領で、もしベクトルホワイトノイズが有限な4次モーメントをもつなら、係数行列が絶対総和可能な Vector MA(∞) の4次モーメントも有限であるといえる(命題 10.2 の (c))。これは y_{i,t} y_{j,t-s} - E(y_{i,t} y_{j,t-s}) が一様可積分といっていることに他ならない(命題 7.7 の (a))。というわけで、193ページの議論と同様に、2次モーメントがエルゴード的であることがいえる(命題 10.2 の (d))。

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

な、なるほど…? いやに Vector MA(∞) の性質を掘り下げますね…今後 Vector MA(∞) をそんなに使っていくんですか? 項数が無限にあったら使いづらいと思うんですが…。

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

思い出してよ、定常 VAR(p) は Vector MA(∞) にかき直せたでしょ? だからこれらは定常 VAR(p) の性質に他ならないよ。そして、ここで2次モーメントのエルゴード性が確認できたことは定常 VAR(p) の標本から色々な推定や検定をすることの土台になる。

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

あっ…いやでも、それをいうには定常 VAR(p) が「係数行列が絶対総和可能な」Vector MA(∞) にかき直せることを示さなければ。

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

示せるよ。前回 VAR(p) の定常条件を示すのに、1次の差分方程式 \xi_t = F \xi_{t-1} + v_t に持ち込んで F を対角化 F=T \Lambda T^{-1} したけど、VAR(p) の Vector MA(∞) 表現における係数行列はこの F\nu 乗の左上 n \times n ブロックだったでしょ? F^\nu =T \Lambda^\nu T^{-1} の任意の成分は263ページの一番下の式の形にかけるから(テキストでは \nu じゃなく s になっているけど)、これの s=1, 2, \cdots についての和は収束するよね。

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

確かに。VAR(p) の定常条件は「どこかの時刻で発生したぶれが爆発しないでゼロに収束してほしい」といったものでしたが、それは「Vector MA(∞) 表現における係数行列が絶対総和可能であってほしい」と同じであったというわけですか…過程の爆発を防ぐために、VAR(p) はある時刻のノイズを再帰してしまうのでその拡大率を1未満にする必要があり、Vector MA(∞) はある時刻のノイズを再帰はしないが過去すべての時刻のノイズを取り込むので係数行列を絶対総和可能にする必要があると…。

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

10.2 節の続きには、n 変量 Vector MA(∞) を変換して r 変量両側(?)Vector MA(∞) をつくれるというような話があるね。後でどう使うのかわからないけど。それでやっと VAR(p) の自己共分散の話になる。さて、どうやって求めればいいだろう?

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

どうやってって…そもそも単変量の AR(p) のときはどうしましたっけ。あ、AR(1) の場合は以下の記事で求めていました。

定義式  Y_t = c + \phi Y_{t-1} + \varepsilon_t の両辺の分散を取るだけですね。過程が定常である以上 V(Y_t) = V(Y_{t-1}) ですから、これについて解けばいいわけです。であれば、VAR(p) の自己共分散も同様に出せるのではないでしょうか?

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

まあその方向でいいんだけどね。まず 0 次の自己共分散行列がどうなるか具体的にやってみよう。項が p 個あると面倒だから式 (10.1.11) の1次の差分方程式形 \xi_t = F \xi_{t-1} + v_t でやってみよう。これの共分散 E(\xi_t \xi_t^\top) を取るとどうなるだろう?

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

簡単です。E(\xi_t \xi_t^\top) = E \bigl( (F \xi_{t-1} + v_t) (F \xi_{t-1} + v_t)^\top \bigr) なので、E(\xi_t \xi_t^\top) = F E(\xi_t \xi_t^\top) F^\top + E(v_t v_t^\top) と整理して、E(\xi_t \xi_t^\top) について解けば…あれ、陽に解けない? 右から F^{-1} をかけても左から (F^\top)^{-1} をかけても E(\xi_t \xi_t^\top) という行列をかきあらわせませんね…。

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

うん、そこで vec 演算子を導入しよう。行列を列ごとにばらして1本の縦ベクトルにつなぐ演算子だ。以下の記事(下から3つ目のセリフ)でも行列の微分の定義に利用したよね(以下の記事でもかいているように沖本本の vech 作用素とも似ているけどあちらは対称行列専用なので少し違うよ)。

この vec 演算子クロネッカー積(定義はテキスト 732~733 ページをみよう)について成り立つ公式 {\rm vec}(ABC)= (C^\top \otimes A){\rm vec}(B) を用いれば―。

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

 {\rm vec} \bigl(E(\xi_t \xi_t^\top)\bigr) = (I_{(np)^2} - F \otimes F)^{-1} {\rm vec} \bigl( E(v_t v_t^\top) \bigr) となるわけですか。I_{(np)^2} - F \otimes F が正則でないと困りますが、定常 VAR(p) であれば必ず正則になるんですね…証明は面倒なのでしませんが。

では 1 次の自己共分散行列はというと、E(\xi_t \xi_{t-1}^\top) = E \bigl( (F \xi_{t-1} + v_t) \xi_{t-1}^\top \bigr) = F E(\xi_{t} \xi_{t}^\top) ですから、0 次の自己共分散行列に左から F をかけるだけですか。

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

ここまで扱ったことをまとめておきましょう。

  • p 次ベクトル自己回帰過程 VAR(p) y_t = c + \Phi_1 y_{t-1} + \cdots + \Phi_p y_{t-p} +\varepsilon_{t} と定義される。
    • 係数行列の成分 \phi_{i,j}^{(\nu)} の意味は「y_t の第 i 成分が \nu ステップ前の自身の第 j 成分にどれだけ依存するか」といえる。
    • 係数行列が定めるある方程式の解がすべて単位円の内側にあるとき VAR(p) は定常であり、定常 VAR(p) は Vector MA(∞) にかき直すことができ、その係数行列は絶対総和可能になる。
  • q 次ベクトル移動平均過程 Vector MA(q) y_t = \mu + \varepsilon_{t} + \Psi_1 \varepsilon_{t-1} + \cdots + \Psi_q \varepsilon_{t-q} と定義される。
    • 係数行列の成分 \psi_{i,j}^{(\nu)} の意味は「y_t の第 i 成分が \nu ステップ前のノイズの第 j 成分にどれだけ依存するか」といえる。
    • Vector MA(∞) は係数行列が絶対総和可能であるという条件の下で定常であり、その自己共分散行列も絶対総和可能になる。
      • さらに、ベクトルホワイトノイズが有限な4次モーメントをもつならば Vector MA(∞) の2次モーメントはエルゴード的である。

つづいたらつづく

Time Series Analysis: ノート10章 弱定常ベクトル過程(その1)

以下の本の10章を読みます。私の誤りは私に帰属します。お気付きの点がありましたらご指摘いただけますと幸いです。

Time Series Analysis

Time Series Analysis

その他の参考文献
  1. 経済・ファイナンスデータの計量時系列分析 (統計ライブラリー) | 竜義, 沖本 |本 | 通販 | Amazon
関連記事
f:id:cookie-box:20190101160814p:plain:w60

10章と11章は多変量時系列への拡張みたいだね。例えば「毎年のGNP」と「毎年の国債金利」を束ねたら多変量時系列になるよね。目次をみるに、10章では単変量時系列のときと同様に改めて多変量時系列のモデルを導入して収束性などの性質を議論して、11章で実際にモデルを使って推定していくのかな。10章冒頭に章の構成について簡単な説明があるね。以下の内容をやっていくみたい。

  • 10.1節 ベクトル自己回帰過程の導入
  • 10.2節 ベクトル過程の自己共分散と収束
  • 10.3節 ベクトル過程の自己共分散生成関数(autocovariance-generating function)
    • この節は 10.4 節の準備になるみたいだね。
  • 10.4節 ベクトル過程のスペクトル
  • 10.5節 ベクトル過程の標本平均
    • この節では命題 7.5 を多変量時系列に拡張するみたい。この結果は後の14章で、自己相関がある場合や不均一分散の場合のOLS推定量を導出するときや、一般化モーメント推定量を考えるときにつかうみたい。あと17章の単位根検定でもつかうって。
ちなみにこの記事のタイトルでは10章のタイトルである Covariance-Stationary Vector Processes を「弱定常ベクトル過程」と訳したけど、これは沖本本の 4.1 節「弱定常ベクトル過程」に倣ったよ。直訳するなら「共分散定常ベクトル過程」だけど、テキスト45ページ(沖本本だと75ページ)にあるように、「弱定常」も「共分散定常」も意味は同じだからね。

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

はあ…ともかく10章は多変量時系列の導入であって、この章で学んだことを後の章でも利用するのですね? それでその命題 7.5 って何でしたっけ…188ページの以下ですか。

命題 7.5(弱定常過程版の大数の弱法則
Y_tE(Y_t)=\mu, \; E(Y_{t} - \mu)(Y_{t-j} - \mu) = \gamma_j であるような弱定常過程とする。但し \gamma_j は絶対総和可能 \sum_{j=0}^\infty |\gamma_j| < \infty とする。このとき標本平均 \overline{Y}_T =(1/T) \sum_{t=1}^T Y_t について以下が成り立つ。
  • \overline{Y}_T\mu に2次平均収束する。\Rightarrow \overline{Y}_T\mu に確率収束する。
  • \lim_{T \to \infty} \bigl[ T E(\overline{Y}_T - \mu)^2 \bigr] = \sum_{j = -\infty}^\infty \gamma_j
これは以下の通常版の大数の弱法則(183ページ)を弱定常過程版にしたものでしたね。
命題(通常版の大数の弱法則
Y_tE(Y_t)=\mu, \; E(Y_{t} - \mu)^2 = \sigma^2 < \infty であるような i.i.d.過程とする。このとき標本平均 \overline{Y}_T =(1/T) \sum_{t=1}^T Y_t について以下が成り立つ。
  • \overline{Y}_T\mu に2次平均収束する。\Rightarrow \overline{Y}_T\mu に確率収束する。
標本サイズを大きくしたときに真のパラメータに確率収束するような推定量を一致推定量とよびますから、これらは「標本平均は母平均の一致推定量である」といっているに他ならないのですね。…それで、多変量時系列版ではどうなるのかを先に覗いてみましょう。10.5節の279ページの最下部ですね。以下の y_t\mu は太字にしませんが n 次元のベクトルとしていますので注意してください。\Gamma_\nun \times n 次元の自己共分散行列です。
命題 10.5(弱定常ベクトル過程版の大数の弱法則
y_tE(y_t)=\mu, \; E(y_{t} - \mu)(y_{t-\nu} - \mu)^\top = \Gamma_\nu であるような弱定常ベクトル過程とする。但し \Gamma_\nu は成分ごとに絶対総和可能 \sum_{\nu=0}^\infty |\gamma_{i,j}^{(\nu)}| < \infty とする。このとき標本平均 \overline{y}_T =(1/T) \sum_{t=1}^T y_t について以下が成り立つ。
  • \overline{y}_T\mu に成分ごとに2次平均収束する。\Rightarrow \overline{y}_T\mu に成分ごとに確率収束する。
  • \lim_{T \to \infty} \bigl[ T E(\overline{y}_T - \mu)(\overline{y}_T - \mu)^\top \bigr] = \sum_{\nu = -\infty}^\infty \Gamma_\nu
こうして見比べると命題 7.5 と似ていますね。多変量になったといっても様変わりしてしまったという感じにはみえないです。しかし、この命題 10.5 が何になるのでしょうか?

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

何って、標本平均がどんな性質をもっているのか押さえておかないと、それをつかった推定や検定がどのような性質をもつのかわからないからね? 後の章で扱うのかと思ったけど、10.5節の続きにも色々な推定量が紹介されているからそれをみればわかるんじゃないかな。…ともかく10.1節の内容に入っていこう。早速257ページで p 次のベクトル自己回帰過程 VAR(p) を以下のように定義しているね。

 y_t = c + \Phi_1 y_{t-1} + \cdots + \Phi_p y_{t-p} +\varepsilon_{t} \tag{10.1.4}
y_t, c,  \varepsilon_{t} はもはや n 次元のベクトルだね。\Phi_\nun \times n 次元の係数行列だね。

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

えっと、自己回帰過程 AR(p) では係数 \phi_jスカラーでしたよね。それがいまや行列になったのですね。…しかし、その行列の各成分は何を意味するのでしょうか。その VAR(p) の第1成分をかき下してみましょう。ここで行列 \Phi_\nuij 列の成分を \phi_{i,j}^{(\nu)} とかき、ベクトル y_ti 番目の成分を y_{t,i} とかきます(ベクトルの添え字をテキストと逆にしたので注意してください)。

\begin{split} y_{t,1} = c_1 &+ \phi_{1,1}^{(1)} y_{t-1,1} + \phi_{1,2}^{(1)} y_{t-1,2} + \cdots + \phi_{1,n}^{(1)} y_{t-1,n} \\ &+ \phi_{1,1}^{(2)} y_{t-2,1} + \phi_{1,2}^{(2)} y_{t-2,2} + \cdots + \phi_{1,n}^{(2)} y_{t-2,n} + \cdots \\ &+ \phi_{1,1}^{(p)} y_{t-p,1} + \phi_{1,2}^{(p)} y_{t-p,2} + \cdots + \phi_{1,n}^{(p)} y_{t-p,n} +\varepsilon_{t,1} \end{split} \tag{10.1.7}
こうですね。となると、\phi_{i,j}^{(\nu)} の意味はy_t の第 i 成分が \nu ステップ前の第 j 成分にどれだけ依存するか」といえますね。

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

そうだね。あと \varepsilon_{t} もいまやホワイトノイズからベクトルホワイトノイズになっているからね。つまり、\varepsilon_{t} は確率ベクトルであって、任意の t で平均ベクトルが E(\varepsilon_{t}) = 0 であって、自己共分散行列 E(\varepsilon_{t} \varepsilon_{t'}^\top) t = t' のときのみ n \times n 正定値行列 \Omega であって、 t \neq t' のときは零行列だ。ちなみに沖本本76ページの記述だけど「\Omega は対角行列である必要はない」よ。例えば、日々の気温とアイスクリームの売り上げを2変量時系列とみなして VAR(p) でモデリングするとしたら、\Omega はきっと正の非対角成分をもつんじゃないかな。だって、ある日に偶発的に気温が高くなったとしたら、その日はアイスクリームも想定より売れそうだもんね。そこは相関があっていい。でも、ある日に偶発的に気温が高くなったのがその次の日以降の気温/アイスクリームの売り上げに影響するのは駄目。 t \neq t' のときは零行列だから。日付をまたいだ影響は係数行列 \Phi_\nu で表現しないといけない。

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

なるほど。異なる日のノイズどうしは無相関だが同じ日のノイズの成分どうしは相関があってもいいんですね。VAR(p) がどのようなモデルなのかは理解した気がします。…しかし、AR(p) がそうであったように VAR(p) も \Phi_\nu に何か条件を課さないと定常にはならないのですよね? 復習すると、過程が(弱)定常であることの定義は以下でした。平均が時刻によらず、自己共分散も時刻によらず時間差だけに依存するということですね。

定義(単変量時系列の弱定常性)
任意の tj に対して E(Y_t) = \mu, \; E(Y_t - \mu)(Y_{t-j} - \mu) = \gamma_j が成立するとき、過程 Y_t は弱定常であるという。
以下、単に定常といった場合弱定常を指すことにします。AR(p) の定常条件は以下でした。
命題(AR(p) の定常条件)
AR(p)  Y_t = c + \phi_1 Y_{t-1} + \cdots + \phi_p Y_{t-p} +\varepsilon_{t} の定常条件は以下である(どちらも同じ)。
  • 方程式 \lambda^p - \phi_1 \lambda^{p-1} - \cdots - \phi_p = 0 の解がすべて単位円より内側にある。
  • 方程式 1 - \phi_1 z - \cdots - \phi_p z^p = 0 の解がすべて単位円より外側にある。
なぜこのような条件になるのかはテキスト7〜13ページですね。AR(p) が定常であると仮定すると、その平均 \mu からの差分を取って、ベクトルに束ねて、式 (1.2.5) の \xi_t = F \xi_{t-1} + v_t の形にすることができるはずです(具体的にかき下すと以下)。
 \displaystyle \left(\begin{array}{c} y_t - \mu \\ y_{t-1} - \mu \\ \vdots  \\  y_{t-p+1} - \mu \end{array}\right) = \left(\begin{array}{ccccc} \phi_1 & \phi_2 & \cdots & \phi_{p-1} & \phi_p \\ 1 & 0 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & 1 & 0 \end{array}\right) \left(\begin{array}{c} y_{t-1} - \mu \\ y_{t-2} - \mu \\ \vdots  \\  y_{t-p} - \mu \end{array}\right) + \left(\begin{array}{c} \varepsilon_t \\ 0 \\ \vdots  \\  0 \end{array}\right) \tag{1.2.5}
この F を対角化 F=T \Lambda T^{-1} すると \Lambda は対角成分に F固有値が並んだ対角行列になりますが、この固有値の中に絶対値が1より大きいものがあったら、どこかの時刻で発生した \mu からの差分が爆発していってしまいますよね。それだと自己共分散が時刻によらないことに矛盾します。逆にすべての固有値の絶対値が1未満ならば、どこかの時刻で発生した \mu からの差分はいつしか消えてなくなるので、この過程の自己共分散は時刻に依存しません。F固有値 \lambda は方程式 \lambda^p - \phi_1 \lambda^{p-1} - \cdots - \phi_p = 0 の解になります(命題 1.1)。なので AR(p) の定常条件は上のようになるわけです。

…単変量の AR(p) のときはこうでしたが、多変量の VAR(p) になると定常条件はどうなるのでしょうか?

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

先に多変量時系列の場合の弱定常性の定義を一応かいておくよ。単変量のときと同様だけど。

定義(多変量時系列の弱定常性)
任意の t\nu に対して E(y_t) = \mu, \; E(y_t - \mu)(y_{t-\nu} - \mu)^\top = \Gamma_\nu が成立するとき、過程 y_t は弱定常であるという。
それで、VAR(p) でも AR(p) のときとやることは同じだよ。VAR(p) をベクトルに束ねて \xi_t = F \xi_{t-1} + v_t の形にする。式 (10.1.11) だね。

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

またしてもベクトルに束ねるのですか。…いえでも、今回は元々ベクトルであるのをさらにベクトルに束ねるということになりますね?? \xi_t = F \xi_{t-1} + v_t をきちんとかき下してみましょう。この Fnp \times np 行列なのですね…。

 \displaystyle \left(\begin{array}{c} y_t - \mu \\ y_{t-1} - \mu \\ \vdots  \\  y_{t-p+1} - \mu \end{array}\right) = \left(\begin{array}{ccccc} \Phi_1 & \Phi_2 & \cdots & \Phi_{p-1} & \Phi_p \\ I_n & 0 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & I_n & 0 \end{array}\right) \left(\begin{array}{c} y_{t-1} - \mu \\ y_{t-2} - \mu \\ \vdots  \\  y_{t-p} - \mu \end{array}\right) + \left(\begin{array}{c} \varepsilon_t \\ 0 \\ \vdots  \\  0 \end{array}\right) \tag{10.1.11}
まあでも、今回も F のすべての固有値の絶対値が1未満であるべきという条件は変わらないですね。多変量時系列のどれか1つの成分でも爆発したら定常でなくなってしまいますから。そうすると、1章のときと同様の操作をして、結局以下のようになりますね。
命題 10.1(VAR(p) の定常条件)
VAR(p)  y_t = c + \Phi_1 y_{t-1} + \cdots + \Phi_p y_{t-p} +\varepsilon_{t} の定常条件は以下である。
  • 方程式  \bigl| \lambda^p I_n - \lambda^{p-1} \Phi_1 - \cdots - \Phi_p \bigr| = 0 の解がすべて単位円より内側にある。
    • この縦棒は行列式の意味である。つまり、これは \lambdanp 次方程式である。
10.1節の続きには、「Vector MA(∞) 表現」というのが出てきますね。AR(p) は定常であれば MA(∞) にかき直すことができましたが、VAR(p) も定常であれば Vector MA(∞) にかき直せるということでしょうか。

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

そうだね。\xi_t = F \xi_{t-1} + v_t を繰り返し適用すれば、

\xi_{t} = v_{t} + F v_{t-1} + F^2 v_{t-2} + F^3 v_{t-3} + \cdots
となって、この両辺のベクトルの最初の n 個の成分を取れば、
y_{t} = \mu + \varepsilon_{t} + F_{1,1} \varepsilon_{t-1} + F_{1,1}^{(2)} \varepsilon_{t-2} + F_{1,1}^{(3)} \varepsilon_{t-3} + \cdots
となるからね。これは Vector MA(∞) だ。ところで、 F_{1,1}^{(\nu)}F^{\nu} の左上 n \times n ブロックを指すよ。ところで、正則な n \times n 行列 H があったとき、以下のように H^{-1} H を挟むことは差し支えないよね。
y_{t} = \mu + H^{-1} H \varepsilon_{t} + F_{1,1} H^{-1} H \varepsilon_{t-1} + F_{1,1}^{(2)} H^{-1} H \varepsilon_{t-2} + F_{1,1}^{(3)} H^{-1} H \varepsilon_{t-3} + \cdots

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

差支えはありませんが、それが何ですか?

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

いまベクトルホワイトノイズは \varepsilon_t だけど、u_t = H \varepsilon_t というように変換されたベクトルホワイトノイズで考えてもいいってことみたい。特に、いま \varepsilon_t の分散共分散行列は \Omega だけど、\Omega を対角化 H\Omega H^\top = D するような H を選んでベクトルホワイトノイズを u_t = H \varepsilon_t と変換すれば、u_t の成分どうしを無相関にできる。

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

ああ、\varepsilon_t確率密度関数f(\varepsilon_t) \propto \exp(-0.5 \varepsilon_t^\top \Omega^{-1} \varepsilon_t) ですから、u_t = H \varepsilon_t と変換したときの確率密度関数f(u_t) \propto \exp(-0.5 u_t^\top (H^{-1})^\top \Omega^{-1} H^{-1} u_t) = \exp(-0.5 u_t^\top (H\Omega H^\top)^{-1} u_t) となって、分散共分散行列が対角行列になり、成分ごとに無相関になりますね。

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

でもその場合 Vector MA(∞) は y_{t} = \mu + H^{-1} u_t + \cdots となっちゃって u_t の項に係数行列がかかっちゃうから、それを許容しないなら成分ごとに無相関にはできないんだけどね。というか Vector MA(∞) の定義がすぐ後に出てくる式 10.2.3 だと思うんだけど、\varepsilon_t に係数行列はかかっていないから、ノイズを変換しちゃうとやっぱり厳密には Vector MA(∞) ではないね。

続く 10.2 節は、まずベクトル過程の自己共分散行列の定義かな。つまり、定常ベクトル過程 y_t\nu 次の自己共分散行列は \Gamma_{\nu} = E \bigl[ ({y}_t - \mu)({y}_{t - \nu} - \mu)^\top \bigr] だ。単変量のときの自己共分散と似ているね。ただ注意してほしいのは、\Gamma_{\nu}\Gamma_{-\nu} が同じになるとは限らない。成り立つのは \Gamma_{\nu}^\top = \Gamma_{-\nu} になる。

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

えっ、何故そのような違いが。

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

もはや積が交換しないしね。({y}_t - \mu)({y}_{t - \nu} - \mu)^\top({y}_{t-\nu} - \mu)({y}_{t} - \mu)^\top って一般に等しくないよね。

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

ああ確かに…\Gamma_\nu の1行目は「y_{t,1}y_{t-\nu,1}, \, y_{t-\nu,2}, \cdots, y_{t-\nu,n}」との共分散ですが、\Gamma_{-\nu} の1行目は「y_{t-\nu,1}y_{t,1}, \, y_{t,2}, \cdots, y_{t,n}」との共分散ですから、意味が違いますね…。

その2につづく

雑記: 単語分散表現の話

word2vecとGloVeが、前者が predicting で後者が counting といわれたり、結果として得られる「 v女王 \fallingdotseq v + v女性 - v男性」という関係式が直感的でめでたしめでたしとなったりしていると思うんですが、どちらも predict しているのではないのかとか、何が直感的なのかとかわからなかったので、word2vecとGloVeが何を目指しているかのメモをかきました。

モデルの定式化は全面的に参考文献1. に依っていますが、記事の最初の方と最後の方には私の適当な解釈を多分に含みます。私の誤りは私に帰属します。お気付きの点がありましたらご指摘いただけますと幸いです。

ゴールとして、「もしword2vecが単語類似度タスクに強くGloVeがニュース分類タスクに強い(要出典)のだとしたら、後者は共起行列をつくってから学習し前者はそうでないため、前者は特定の文脈にフィットする解を目指しやすいが後者はあくまで全体の平均を目指し、それゆえに前者は単語類似度タスク(特定の文脈での意味)に強く後者がニュース分類タスク(ニュースを構成する単語全体の意味)に強いのではないか」というような考察を適当ではなくきちんとできるようになりたかったのですが、そこまでいきませんでした。また参考文献1. の 2-7 節に、「単語の分散表現の性能差は各モデルの理論的な側面よりも、実装上の細かいトリックや実験設定に左右されるとの指摘がある.」との記述があります。

参考文献
  1. 岡崎 直観. 言語処理における分散表現学習のフロンティア(<特集>ニューラルネットワーク研究のフロンティア). 人工知能, 31(2), 189-201, 2016. 言語処理における分散表現学習のフロンティア(<特集>ニューラルネットワーク研究のフロンティア)
    • 人工知能学会の学会誌の2016年3月号の記事です。
  2. Mikolov, T., Sutskever, I., Chen, K., Corrado, G. S., Dean, J. Distributed Representations of Words and Phrases and their Compositionality, Proc. of NIPS, pp. 3111-3119, 2013. [1310.4546] Distributed Representations of Words and Phrases and their Compositionality
    • word2vecの原論文です。
  3. https://code.google.com/archive/p/word2vec/
    • word2vecの実装があり、architecture: skip-gram (slower, better for infrequent words) vs CBOW (fast) などのメモがあります。
  4. Pennington, J., Socher, R., Manning, C.. Glove: Global vectors for word representation, Proc. of EMNLP, pp. 1532-1543, 2014. GloVe: Global Vectors for Word Representation - ACL Anthology
    • GloVeの原論文です。
  5. https://nbviewer.jupyter.org/github/DSKSD/DeepNLP-models-Pytorch/blob/master/notebooks/03.GloVe.ipynb
    • GloVeの実装です。

まとめ(文字の定義は後述)
skip-gram も GloVe も、結局、単語 c の周辺(前後 \delta 語など、適当に決める)に出現する単語 w の分布 p(w|c) を、内積のエクスポネンシャルによるモデル q(w|c) = C_c \tilde{C}_w \exp(v_c^\top \tilde{v}_{w}) でフィッティングすることを通して単語 c のベクトル表現 v_c を得ようとするのは共通している。
  • \tilde{C}_w = 1 のときは q(w|c) = {\rm softmax}(v_c^\top \tilde{v}_{w}) に他ならない。C_c は規格化定数である。
  • q(w|c_a)q(w|c_b) の比に興味があるときは \tilde{C}_w > 0 があってもよい(比をとると消える)。
但し、分布間の距離の指標は様々であるし、そもそも語彙数が膨大であるためにこの分布の規格化定数を直接計算するのも現実的ではない。skip-gram と GloVe は何を最小化しようとしてどのように規格化を回避するかが異なっている。
  • skip-gram は q(w|c) = {\rm softmax}(v_c^\top \tilde{v}_{w}) とし、学習コーパス内に出現するすべての「単語 c - その周辺の単語 w」の組の条件付き確率 q(w|c) の対数の和  \sum_{t=1}^T \sum_{c \in C_{w_t}} \log q(w_t|c) を最大化するように(∴ 交差エントロピーの重み付き和  - \sum_{c \in V} \#(c,\cdot) \sum_{w \in V} p(w|c) \log q(w|c) を最小化するように)学習することを目指す。が、規格化の計算が難しいので  {\rm sigmoid}(v_c^\top \tilde{v}_{w}) による2値分類タスクにしてしまう(正例1個につき負例用の周辺単語を K 個サンプリングする)。
    • CBOW(その名の通り周辺の単語を先に「カバン」に入れて=足してしまう)はこの最適化をさらに簡略化したものとも捉えられる。
  • GloVe は p(w|c) = \#(c, w) / \#(c, \cdot ) であることに着目して、「 p(w|c)q(w|c) = C_c \tilde{C}_w \exp(v_c^\top \tilde{v}_{w}) でフィッティングするタスク」を、「 \log \#(c, w)v_c^\top \tilde{v}_w + b_c + \tilde{b}_w でフィッティングするタスク」に読み替える。損失は2乗誤差を採用するが、あまり発生しない「単語 c - その周辺の単語 w」の組に対しては重みを小さくするようにする。
    • そのため、予め \#(c, w) のテーブルを用意しておくことになる。
いずれにせよ、skip-gram や GloVe で得たベクトルは「内積が対数確率をオフセットしたもの」という性質をもっている。内積は分配則が成り立つ(和の内積内積の和)ので、v_a + v_b は「周辺の単語の対数確率が単語 a の周辺の対数確率と単語 b の周辺の対数確率の和になるような仮想的な単語のベクトル」と捉えられる。となると、v + v女性 - v男性 というベクトルをもつ仮想的な単語が考えられるが、 v女王 がそれに近かったという事実は、「女王」の周りの分布がそのように「王」「女性」「男性」の分布で特徴付けられるということに他ならない。

word2vecのモデル
  • 学習コーパスを単語列 w_1, w_2, \cdots, w_T とする。
  • 単語 w_t を中心語とした前後 \delta 語の文脈窓を C_{w_t} = (w_{t-\delta}, \cdots, w_{t-1}, w_{t+1}, \cdots, w_{t+\delta}) とする。
word2vecでは、文脈窓 C を与えたときの中心語の確率分布 p(w|C) を、C 内の単語のベクトルを用いてモデリングすることを通して、単語のベクトル表現を獲得する。つまり、長さ 2 \delta の単語列を与えたときに、それらの単語のベクトル表現を用いてどんな単語がその中央にはまりそうかの分布を出すことを目指すのだが、2 \delta 個の単語の情報をどう集約するかで流派が分かれる(下図)。
f:id:cookie-box:20210217225124p:plain
  • skip-gramでは p(w|C) \prod_{c \in C} {\rm softmax}(v_c^\top \tilde{v}_{w})モデリングする。
    つまり、「文脈窓の中心にどんな単語がきそうか」を、文脈窓内の各単語の「この単語がいる文脈窓の中心にどんな単語がきそうか」の積で表現する(先に各単語からの確率を出して積をとって集約する派)。
  • CBOWでは p(w|C){\rm softmax} \left( \left( \sum_{c \in C}v_c \right)^\top \tilde{v}_{w} \right)モデリングする。
    つまり、「文脈窓の中心にどんな単語がきそうか」を、文脈窓内の全単語のベクトルの和で表現する(先に全単語のベクトルを和に集約してから確率を出す派)。
但し実際にはソフトマックスではなくシグモイドで活性化し、サンプリングした負例との2値分類を学習する。
skip-gram
skip-gramでは p(w|C) = \prod_{c \in C} p(w|c) と考える。p(w|c) はいわば「単語 c を与えたときに単語 w が単語 c のいる文脈窓の中心語である確率」になる。この p(w|c) をモデル q(w|c) = {\rm softmax}(v_c^\top \tilde{v}_{w}) でフィッティングする。ただし、v_c は単語 c の文脈語としてのベクトル、\tilde{v}_{w} は単語 w の中心語としてのベクトルであり、全単語にこれらのベクトルをどう用意すべきかがフィッティング対象となる。損失関数は以下とする(学習コーパス内のすべての中心語-文脈語の組の確率の積を最大化する=対数確率の和を最小化する)。
 \displaystyle \mathcal{L}^{\rm SG} = \sum_{t=1}^T \sum_{c \in C_{w_t}} \log q(w_t | c) = \sum_{t=1}^T \sum_{c \in C_{w_t}} \log {\rm softmax}(v_c^\top \tilde{v}_{w_t}) = \sum_{t=1}^T \sum_{c \in C_{w_t}} \log \frac{\exp (v_c^\top \tilde{v}_{w_t})}{ \sum_{w' \in V} \exp (v_c^\top \tilde{v}_{w'})}
但し実際には全単語 V 上でソフトマックスをとることは難しいので損失関数を少し変更する。「文脈語 c から中心語 w を当てる  \# V 値分類タスク」ではなく「文脈語 c と中心語 w が正しい文脈語と中心語の組かどうか当てる 2 値分類タスク」と考えることにして、活性化関数をソフトマックスからシグモイドに取り換える。このタスク設定にすると負例が必要になるので、中心語-文脈語の組1つにつき、負例のための中心語 w^{(k)}K 回サンプリングする(単語のユニグラム分布からサンプリングする)。
 \displaystyle \mathcal{L}^{\rm SG} = \sum_{t=1}^T \sum_{c \in C_{w_t}} \left( \log \frac{\exp (v_c^\top \tilde{v}_{w_t})}{1 + \exp (v_c^\top \tilde{v}_{w_t})} + \sum_{k=1}^K \log \left( 1 - \frac{\exp (v_c^\top \tilde{v}_{w^{(k)}})}{1 + \exp (v_c^\top \tilde{v}_{w^{(k)}})} \right) \right)
この損失関数を最小化するように学習し、各単語の文脈語としてのベクトル v_c をその単語のベクトル表現とする。
CBOW(continuous bag-of-words)
文脈窓 C に対して、文脈窓内の単語の文脈語としてのベクトルの和を v_C = \sum_{c \in C} v_c とする。CBOWでは p(w|C)={\rm softmax}(v_C^\top \tilde{v}_{w}) と考える。損失関数は以下となる。
 \displaystyle \mathcal{L}^{\rm CBOW} = \sum_{t=1}^T \log q(w_t | C_{w_t}) = \sum_{t=1}^T \log {\rm softmax}(v_{C_{w_t}}^\top \tilde{v}_{w_t}) = \sum_{t=1}^T \log \frac{\exp (v_{C_{w_t}}^\top \tilde{v}_{w_t})}{ \sum_{w' \in V} \exp (v_{C_{w_t}}^\top \tilde{v}_{w'})}
但し実際には全単語 V 上でソフトマックスをとることは難しいので skip-gram 同様に損失関数を変更する。
 \displaystyle \mathcal{L}^{\rm CBOW} = \sum_{t=1}^T \left( \log \frac{\exp (v_{C_{w_t}}^\top \tilde{v}_{w_t})}{1 + \exp (v_{C_{w_t}}^\top \tilde{v}_{w_t})} + \sum_{k=1}^K \log \left( 1 - \frac{\exp (v_{C_{w_t}}^\top \tilde{v}_{w^{(k)}})}{1 + \exp (v_{C_{w_t}}^\top \tilde{v}_{w^{(k)}})} \right) \right)
GloVeのモデル
GloVeでは、2単語のベクトルの差でそれぞれの単語の周辺に出現する単語の確率分布の比をモデリングすることを通して、単語のベクトル表現を獲得する。但し、「周辺に出現する単語のベクトルと内積をとった上でエクスポネンシャルをとる」というモデルにしたため、結局、個々の単語の確率分布をモデリングすることになっている。
f:id:cookie-box:20210218003922p:plain
具体的に、
  • 単語 i の周辺に単語 j が出現した回数を \#(i, j) とする。
  • 単語 i の周辺に任意の単語が出現した回数を \#(i, \cdot ) とする。
  • 単語 i の周辺に単語 j が出現する確率は p(j|i) = \#(i, j) / \#(i, \cdot ) となる。
が、確率分布 p(j|i) そのものは、どんな単語の周辺にもよく出現するような単語 j で大きくなるため、単語 i のよい特徴になっていない。それよりも、他の単語 k の分布との比を取った p(j|i)/p(j|k) の方が単語 i の単語 k との関係のよい特徴になっている。ので、確率分布の比 p(j|i)/p(j|k)F(v_i, v_k, \tilde{v}_j)モデリングすることを目指す。ここで、F の関数形に以下を仮定する。
  • v_i, v_k に依存する項は v_i, v_k の差 v_i - v_k にのみ依存する。
  • v_i - v_kスカラーに変換するにあたり \tilde{v}_j との内積 (v_i - v_k)^\top \tilde{v}_j をとる。
  • (v_i - v_k)^\top \tilde{v}_j を確率分布の比に変換するにあたり、和を変換したものが変換したものの積になるように \exp を採用する(同じ単語どうしの確率分布の比は常に1にならなければならない)
よって、F(v_i, v_k, \tilde{v}_j) = \exp \bigl( (v_i - v_k)^\top \tilde{v}_j \bigr) となるが、そうなると p(j|i)/p(j|k)\exp(v_i^\top \tilde{v}_j)/\exp(v_k^\top \tilde{v}_j)モデリングすることになるので、p(j|i)\tilde{C}_j \exp(v_i^\top \tilde{v}_j)モデリングすることと同じである(\tilde{C}_jj にのみ依存する係数)。さらに、p(j|i) = \#(i, j) / \#(i, \cdot ) を利用すると、\log \#(i, j)v_i^\top \tilde{v}_j + b_i + \tilde{b}_jモデリングすることになる。損失は2乗誤差を採用する。
 \displaystyle \mathcal{L}^{\rm GloVe} = \sum_{i=1}^{|V|} \sum_{j=1}^{|C|} f(\#(i,j)) \left( v_i^\top \tilde{v}_j + b_i + \tilde{b}_j - \log \#(i, j) \right)^2
但し、f(\#(i,j)) は共起頻度が小さいペアの重みを小さくするための因子で、例えば以下(共起頻度が100回未満のペアの重みを徐々に軽くする例)が用いられる。
 \displaystyle f(x) = \begin{cases} (x / 100)^{0.75}  & (x < 100) \\ 1 & ({\rm otherwise}) \end{cases}

得られたベクトル表現の性質
これらの手法で得られた単語のベクトル表現は、モデルより、別の単語のベクトル表現との内積が意味をもつ(対数確率をオフセットしたものになる)。もし C_a \fallingdotseq C_b であれば、
  • v_a^\top \tilde{v}_x が、単語 x が単語 a の周辺にいそうな度合いになる。
  • v_b^\top \tilde{v}_x が、単語 x が単語 b の周辺にいそうな度合いになる。
  • v_b^\top \tilde{v}_x - v_a^\top \tilde{v}_x = (v_b - v_a)^\top \tilde{v}_x は、それらの差になる。
    • (v女性 - v男性)^\top \tilde{v}スカート は、正になりそう(※ 適当な例)。
    • (v女性 - v男性)^\top \tilde{v}ズボン は、負になりそう(※ 適当な例)。
  •  v_{a'}^\top \tilde{v}_x + v_b^\top \tilde{v}_x - v_a^\top \tilde{v}_x = (v_{a'} + v_b - v_a)^\top \tilde{v}_x は、単語 x が単語 a' の周辺にいそうな度合いに先の差を足したものになる。
    • (v + v女性 - v男性)^\top \tilde{v}スカート は、「王」の文脈に「スカート」がいそうな度合いを大きくしたもの。
    • (v + v女性 - v男性)^\top \tilde{v}ズボン は、「王」の文脈に「ズボン」がいそうな度合いを小さくしたもの。
    • 同様に、任意の単語 x について、(v + v女性 - v男性)^\top \tilde{v}_x は、「王」の周辺に x がいそうな度合いに、「女性」の周辺に x がいそうな度合いを足して、「男性」の周辺に x がいそうな度合いを引いたもの。
word2vecで(GloVeでも)  v女王 \fallingdotseq v + v女性 - v男性 となることがよく知られているが、それは、任意の単語 x について、「女王」の周辺に x がいそうな度合いが、「王」の周辺に x がいそうな度合いに「女性」の周辺に x がいそうな度合いを足して「男性」の周辺に x がいそうな度合いを引いたものに近くなっていることを示唆しているのに他ならない。