TL;DR

Xを多変量化したガウス過程と、Yを多変量化したガウス過程はYの成分同士の相関性を考えなければ簡単。
多変量ガウス過程において、Yの成分ごとの相関性を考慮することは可能だが、演算の重さが次元の3乗で重くなる上に事前分布の設計次第なようで、使い勝手は悪そうだ(計算ミスの可能性もあるので、似たことが考察されている参考文献が欲しいところ)

多変量の場合のガウス過程

前記事で導出したガウス過程回帰

$$y(X)=\mathcal N(y|m(X), K(X,X))$$

について、xとyがそれぞれ多変数だった場合を考慮するのを忘れていたので、この記事で触れよう。

Xが多変数、yが単変数の回帰

まず、$X$が単変数ではなく多変数で、$y$が単変数の場合を考える。これは意外と単純で、事前分布としてのカーネル関数と平均値関数を若干の修正するだけで済む。

たとえば、

$$m(\mathbf X)=0$$
$$k(\mathbf X_i,\mathbf X_j)=\exp(-(\mathbf X_i-\mathbf X_j)(\mathbf X_i-\mathbf X_j)^T)$$

という関数を事前分布にすれば考慮したことになる。いままで便宜上$\mathbf X$の行方向はサンプル次元を指していたので、列ベクトルに多変数を据えることにした。
平均値関数は定数関数なので引数がベクトルになったところで変化しないし、カーネル関数は一次元上のユークリッド距離$(X_i-X_j)^2$を多次元空間上のユークリッド距離$(\mathbf X_i-\mathbf X_j)(\mathbf X_i-\mathbf X_j)^T$に変更するだけである。距離は内積なのでスカラー量であり、このカーネル関数を使ったところで更新式は変わらず、

$$K_\ast =K(\mathbf X,\mathbf X)-K(\mathbf X,\mathbf X')[K(\mathbf X',\mathbf X')]^{-1}K(\mathbf X',\mathbf X)$$
$$m_\ast =m(\mathbf X)+K(\mathbf X,\mathbf X')[K(\mathbf X',\mathbf X')]^{-1}(y(\mathbf X')-m(\mathbf X'))$$

で定義されるのみである。以降記述簡略化のため、

$$K=K(X,X),m=m(X)$$
$$Q=K(X,X')[K(X',X')]^{-1},R=K(X',X),dy=y(X')-m(X')$$
とすると、

$$K_\ast=K-QR$$
$$m_\ast=m-Qdy$$

で表現できる。以降はこの書き方を採用する。

Xもyも多変数の場合

yが多変数Yだった場合は結構複雑な話である。
一旦ガウス過程を忘れてYが単なるガウス分布だった場合を考えよう。このときYはベクトルであるので、右辺は多変量ガウシアンである必要が出てくる。

$$Y=\mathcal(Y|\mathbf m, \mathbf K)$$

ここで、mは多変量方向のベクトル、Kは多変量方向に成分を持つ分散共分散行列である。

さあガウス過程を思い出そう。Yはサンプル方向にも無限次元の長さのベクトルと考えるのだった。
すなわちYは行方向にサンプル、列方向に多変量成分を展開する行列(あるいはベクトルの関数)として見なせる。

mの形状も同様に行列(あるいはベクトルの関数)として表せる。

問題は分散共分散行列である。これはサンプル方向にグラム行列をなしているため、いわゆる4階のテンソル(あるいは、行列を返却する二変数関数とでもいうべきか)となる。

こうして考えるとひどく複雑そうに見える。
多変量方向を添字s,tとして、一応やれるところまで立式してみよう。

$$Y=\mathcal N(Y|\mathbf m, \mathbb K)$$

$$\mathbf m={m_s(X_i)}$$
$$\mathbb K={k_{st}(X_i,X_j)}$$

この事前分布の置き方は、たとえば以下である。

$$m_s(X_i)=0$$
$$k_{st}(X_i,X_j)=\exp(-(X_i-X_j)(X_i-X_j)^T)$$

つまりデータがなければ$X=X'$のとき多変量方向の全y成分が相関する、ということを事前分布においている。

更新式を導出してみよう。前の記事の導出式

$$p(\hat Y|D')=\frac{p(\hat Y|\hat X)}{p(Y'|X')}$$

によっていったんガウシアンを展開する。

$$\ln p(\hat Y|\mathcal D^\prime)=-(\hat Y_{is}-\hat m_{is})[\hat K_{ijst}]^{-1}(\hat Y_{jt}-\hat m_{jt})+(Y^\prime_{is}-m^\prime_{is})[K^\prime_{ijst}]^{-1}(Y^\prime_{jt}-m^\prime_{jt})$$

ここで、右辺の添字はそれぞれ和を取ることを意味する(縮約という)。
いま、適当なテンソル$V,M$について

$$V_{is}=V_{a}$$
$$V_{jt}=V_{b}$$
$$M_{ijst}=M_{ab}$$

となるようにしよう。すると、

$$K_{ijst}=K_{ab}$$
と見なすことができ、この$K$も$a,b$について対称行列になる。ここで$K_{ab}$は任意の$i,s$と$j,t$の間の分散と共分散を意味する行列だと解釈できるので、$K_{ab}$はサンプル次元と多変量次元の両方を一列に並べたときの分散共分散行列だと解釈することができる。
この式は前回の導出と全く同じ変形で、

$$K_{\ast ab}=K_{ab}-Q_{ab'}R_{b'b}$$
$$m_{\ast a}=m_{a}-Q_{ab'}dy_{b'}$$

とかける。

(この変換のどこかに間違いなどなければ)これはつまりyを多変量化したときの処理はyの値同士の成分の相関性を加味して共分散が更新されるわけではなく、あくまで事前分布で仮定した範囲で成分同士の共分散が考慮され回帰される、ということになる。

もし事前分布として$a,b$の間の共分散を
$s\neq t$のときにすべて$0$であるとした場合は、成分ごとの従属性が完全になくなるので、結局

$$K_{\ast ab}=\delta_{st}(K_{ijss}-Q_{ij'ss}R_{j'jss})$$
$$m_{\ast a}=m_a+\delta_{st}Q_{ij'st}dy_{j't}$$

であり、これはつまり
$$K_{\ast ijst}=0 \quad s\neq t$$
$$K_{\ast s}=K_{s}-Q_sR_s$$
$$m_{\ast s}=m_s+Q_s+dy_s$$
とかける。

ただし、任意の4階テンソル$T$について、$T_{\ast s}=V_{\ast ijss}$を意味するものとする。

ようは、yの成分同士の相関が無相関であるように事前分布を仮定する場合は、単成分出力のガウス過程回帰を成分ごとに行えばいいことになる。

まとめ

X成分が多変数でY成分が単回帰の場合はベクトル同士の距離を使ってカーネル関数を定義すれば、普通のガウス過程回帰と同じ更新式になる。
Y成分が多変数だが互いの成分同士で相関しない事前分布を仮定した場合も、Yの成分ごとに独立してガウス過程回帰を行えばよい。
Y成分が多変数で互いの成分同士が相関する仮定の場合は、その相関の程度を見積もって事前分布を定義する。yの値に応じて成分ごとの相関を加味してくれるわけではなさそうだ。

ついで、計算量についての考察。

演算量はグラム行列を一旦途中計算結果として保存し、かつXの次元がサンプル数に比べて十分小さい前提で、Xの次元はそれほど演算量に効かず、グラム行列の逆行列を求める操作の演算負荷が大きい。
Yが単変数のときグラム行列のサイズは回帰に使うサンプル数$N$について$N\times N$であり、計算オーダーは$\mathcal O(N^3)$である。
Yが多変数($D$次元とする)であるが成分同士の相関を考慮しない場合は上記逆行列をYの次元の数だけ求めることになるので、$\mathcal O(DN^3)$である。一方、事前分布を成分ごとに変えなければ逆行列を使い回せるので、その場合は計算量は増加せず$\mathcal O(N^3)$ですむ。
Yが多変数であり、成分同士の相関を考慮した事前分布で回帰する場合は、グラム行列のサイズは$DN\times DN$となり、逆行列の計算量は$\mathcal O(D^3N^3)$であり、飛躍的に大きくなる。Y成分同士の相互相関の度合いも事前分布で設計しなくてはならないようなので、Y成分同士の従属性まで考えた多変量ガウス過程は、あまり使い勝手のいい回帰とはいえなそうだ。

所感

前記事とあわせて、数式展開するだけで疲れる量の記事になってしまった。次は単回帰記事書いて息抜きしたい。