前書き
グラフを対象とした機械学習について学ぼうと思ったのですが, 日本語の記事だけでは理解しきれなかったので色々調べたことをまとめます. 主にこちらのサーベイ論文を参考にしています（該当箇所で改めて引用します）.
グラフ信号処理
グラフ信号処理では文字通りグラフの頂点の信号（特徴量）に対して, 一般的な信号処理と同様にフーリエ変換やフィルタなどの処理を
行います. 初めはグラフラプラシアンの定義や畳み込みの定義がよくわかりませんでしたが, 色々調べてよくよく考えてみると
通常の信号処理のアナロジーによって自然に導入できることが分かりました. 以下, 自分が戸惑った点や不明だった点についてまとめていきます. また, 本題に入る前に, 記事内で一般的に用いる記号を次の表で定義しておきます.
| 記号 |
意味 |
| |\cdot| |
集合の濃度です. 有限の頂点を持つグラフしか扱わないので要素数と捉えて頂いて構いません. |
| \mathcal{G} |
グラフです. 頂点集合, 辺集合, 隣接行列の3つ組で表します. \mathcal{G}=(\mathcal{V},\mathcal{E},W)
|
| \mathcal{V} |
グラフ\mathcal{G}の頂点集合を表します. i番目の要素はv_iで表します. また, |\mathcal{V}|=Nとします. |
| \mathcal{E} |
グラフ\mathcal{G}の辺集合を表します. k番目の要素e_kが頂点v_iと頂点v_jを繋ぐ場合, e_k=(i,j)のように順序対で表します. |
| W |
グラフ\mathcal{G}の重み付き隣接行列です. W_{ij}=\begin{cases}w_{ij} & \left((i,j)\in\mathcal{E}\right)\\ 0 & \left(\text{otherwise}\right)\end{cases} w_{ij}は辺(i,j)の重みです. 重み付きでないグラフの場合は1とします. N行N列の行列であり, 平易に言えば, 頂点iと頂点jの間に辺が存在すれば, i行j列の要素はその重みに, そうでなければ0であるような行列です. |
| \text{deg}(\cdot) |
頂点の次数です. |
| D |
次数行列です. D_{ij}=\begin{cases}\text{deg}(v_i) & \left(i=j\right)\\ 0 & \left(\text{otherwise}\right)\end{cases} で定義されます. 対角要素に頂点v_iの次数を入れた, N行N列の行列です. |
| \cdot\odot\cdot |
アダマール積（要素積）行列の要素同士の積を計算する二項演算子です. |
| \text{diag}(\cdot) |
\cdotの要素を対角要素に持つ対角行列です. |
グラフラプラシアンとは？
グラフ上での畳み込み演算では基本的にグラフラプラシアンを用います. 以下にその定義を示します.
正規化グラフラプラシアンなどもありますが, ここではこのŁをグラフラプラシアンとします. 例を挙げて確認してみましょう.
上に示したグラフの次数行列は
D=\begin{pmatrix}
4&0&0&0&0\\
0&1&0&0&0\\
0&0&1&0&0\\
0&0&0&1&0\\
0&0&0&0&1
\end{pmatrix}
隣接行列は
W=\begin{pmatrix}
0&1&1&1&1\\
1&0&0&0&0\\
1&0&0&0&0\\
1&0&0&0&0\\
1&0&0&0&0
\end{pmatrix}
従って, グラフラプラシアンは
\textit{Ł}=\begin{pmatrix}
4&-1&-1&-1&-1\\
-1&1&0&0&0\\
-1&0&1&0&0\\
-1&0&0&1&0\\
-1&0&0&0&1
\end{pmatrix}
となります. 計算はとても簡単ですね.
私が初めに参考にした本ではこれを起点に後述する畳み込み演算を展開していったのですが, そもそもこの定義の意味がよくわからず（なんでラプラシアン？）戸惑ってしまいました. 色々と調べる内に前述の英語版Wilipediaの記事に行き着きました. 結論から申し上げると, グラフラプラシアンは我々が普段お世話になっているラプラシアン（のマイナス）を離散化・近似したものです.
一般的な2次元のデカルト座標系で定義された関数f(x,y)について, 点(x,y)における\displaystyle\nabla^2f=\frac{\partial^2 f}{\partial x^2}+\frac{\partial^2 f}{\partial y^2}を中心差分で近似することを考えます. まずはx方向のみについて考えると
\begin{array}{rl}
\displaystyle\frac{\partial^2 f}{\partial x^2}&\displaystyle=\frac{\partial}{\partial x}\left\{\frac{\partial f}{\partial x}\right\}\\
&\displaystyle\simeq\frac{\partial}{\partial x}\left\{\frac{f(x+\frac{h}{2},y)-f(x+\frac{h}{2},y)}{h}\right\}\\
&\displaystyle\simeq\frac{f(x+h,y)+f(x-h,y)-2f(x,y)}{h^2}
\end{array}
となります. 従って
\begin{array}{rl}
\displaystyle\nabla^2 f(x,y) &=\displaystyle\frac{f(x+h,y)+f(x-h,y)+f(x,y+h)+f(x,y-h)-4f(x,y)}{h^2}
\end{array}
を得ます. これをさらに少し変形してやると
\displaystyle\nabla^2 f(x,y) =\displaystyle-\frac{1}{h^2}\begin{pmatrix}4&-1&-1&-1&-1\\\end{pmatrix}\begin{pmatrix}f(x,y)\\ f(x+h,y)\\ f(x-h,y)\\ f(x,y+h)\\ f(x,y-h)\end{pmatrix}
と書けます. 係数の行ベクトルがどこかで見た形ですね. 先ほどの例で挙げたグラフラプラシアンの1行目そのものです. 他の点については, 関数f(x,y)の値を指定しない（即ちNeumann境界条件となる）条件の下でラプラシアンを作用させると
\nabla^2\begin{pmatrix}f(x,y)\\ f(x+h,y)\\ f(x-h,y)\\ f(x,y+h)\\ f(x,y-h)\end{pmatrix}=-\frac{1}{h^2}\begin{pmatrix}
4&-1&-1&-1&-1\\
-1&1&0&0&0\\
-1&0&1&0&0\\
-1&0&0&1&0\\
-1&0&0&0&1
\end{pmatrix}\begin{pmatrix}f(x,y)\\ f(x+h,y)\\ f(x-h,y)\\ f(x,y+h)\\ f(x,y-h)\end{pmatrix}
となり, 先ほどの例のグラフラプラシアンを得ます. このように, グラフラプラシアンは通常のラプラシアンを離散グラフに適用したものだと解釈できます（h=1とし, f(x,y)を先ほどのグラフの頂点v_0, f(x\pm h, y\pm h)をその他の頂点とみなせばそのものです）.
今の例を元に, より一般のグラフについても考えてみましょう. 先ほどの計算途中で現れた
\begin{array}{rl}
\displaystyle\nabla^2 f(x,y) &=\displaystyle\frac{f(x+h,y)+f(x-h,y)+f(x,y+h)+f(x,y-h)-4f(x,y)}{h^2}
\end{array}
より, \nabla^2 f(u)は, 頂点uに隣接したすべての頂点v\in\mathcal{V}\;s.t.\;(u,v)\in\mathcal{E}について, f(v)-f(u)を足し合わせたものといえます. ゆえに
\begin{array}{rl}
\displaystyle\nabla^2 f(v_i)&=\displaystyle\sum_{v_j\in\mathcal{V}\;s.t.\;(v_i,v_j)\in\mathcal{E}}\{f(v_j)-f(v_i)\}\\
&=l\begin{pmatrix}f(v_0)\\ \vdots\\ f(v_N)\end{pmatrix}
\end{array}
と書けば
l_j=\begin{cases}
\text{deg}(v_i)&(i==j)\\
-W_{ij}&(otherwise)\\
\end{cases}
となり, これを他の頂点についても計算して行方向に並べることでグラフラプラシアンを得ます.
グラフにおける微分と勾配
参考文献ではグラフ信号のsmoothnessを議論するために, 有向グラフ上での微分, および勾配を導入している例も紹介しています.
頂点v_iにおける信号f(i)の辺e=(i,j)による微分を
\left.\frac{\partial f}{\partial e}\right|_i=\sqrt{W_{ij}}\{f(j)-f(i)\}
頂点v_iにおける勾配ベクトルを
\nabla f(i) = \left\{\left(\left.\frac{\partial f}{\partial e}\right|_i\right)_{e\in\mathcal{E}\;\text{s.t.}\;e=(i,j)\;\text{for some}\;j\in\mathcal{V}}\right\}
と定義しています（では微分の定義の符号が逆になっていますが, 本記事ではの定義を採用します）. これらはグラフラプラシアンとも関係があります. 次のグラフについて確認してみましょう.
まず, グラフラプラシアンですが
\begin{array}{rl}
\textit{Ł}&=\begin{pmatrix}
2&-1&-1\\
-1&2&-1\\
-1&-1&2\\
\end{pmatrix}
\end{array}
と計算できます. 次に勾配をいくつかの頂点で計算してみましょう. 頂点v_0での勾配は
\begin{array}{rl}
\nabla f(0)&=\begin{pmatrix}
-1&1&0\\
-1&0&1\\
\end{pmatrix}
\begin{pmatrix}
f(0)\\ f(1)\\ f(2)\\
\end{pmatrix}\\
\end{array}
というように記述でき, また頂点v_1における勾配は
\begin{array}{rl}
\nabla f(1)&=\begin{pmatrix}
0&-1&1\\
\end{pmatrix}
\begin{pmatrix}
f(0)\\ f(1)\\ f(2)\\
\end{pmatrix}\\
\end{array}
となります（頂点v_2の勾配は空ベクトルです）. これらの係数の行列が勾配演算子になりますが, ここでこのグラフの接続行列, およびその転置を見てみましょう.
\begin{array}{rl}
B&=\begin{pmatrix}
-1&-1&0\\
0&1&-1\\
1&0&1\\
\end{pmatrix}
\end{array}
\begin{array}{rl}
B^\text{T}&=\begin{pmatrix}
-1&0&1\\
-1&1&0\\
0&-1&1\\
\end{pmatrix}
\end{array}
特に, 接続行列の転置に注目すると, 一部の行が先ほど計算した勾配の係数行列になっていることが分かります. このことから, グラフの接続行列（の転置）が勾配演算子に対応していることが確認できます. さらに, 計算すると確認できますが, ラプラシアン\nabla^2=\nabla\cdot\nablaと同様に
が成立しています.
グラフフーリエ変換の導入
ここまででグラフラプラシアンが一般的なラプラシアンに対応するものだと分かりました. そうなるとグラフ上でのフーリエ変換も自然に導入することができます. フーリエ級数では周期関数をe^{iw_nx}の線形和で表していました. そして, これを周期のない（周期が無限大）の関数に適用できるようにしたものがフーリエ変換で次の式で与えられます.
\braket{f,e^{2\pi ix\xi}}=\int_{-\infty}^{\infty}f(x)e^{-2\pi ix\xi}dx
こちらでは和から積分になっていますが, フーリエ変換では関数をe^{2\pi ix\xi}を基底とする空間に射影していること, 基底となっているe^{2\pi ix\xi}がラプラシアンの固有関数であることに注目します. 実際に計算してみると
\nabla^2(e^{2\pi ix\xi})=-4\pi^2\xi^2e^{2\pi ix\xi}
より, 固有関数であることが分かります. グラフにおいても同様の議論ができそうです. グラフにおけるラプラシアンはもちろんグラフラプラシアンで, 固有関数に対応するものはもちろんグラフラプラシアンの固有ベクトルです. 従って, グラフフーリエ変換, および逆変換は次のように定義するのが自然でしょう.
行列Uをグラフラプラシアン\textit{Ł}の固有ベクトルを並べた行列とします（グラフラプラシアンはその定義により対称行列です. 即ちエルミート行列でもあるので, すべての固有ベクトルは一次独立です. ゆえに, 常に行列Uは存在します）. このとき, グラフフーリエ変換を
逆変換を
と定義します（xはグラフ上の信号です）. ただ固有ベクトルが張る空間に射影するだけですね. 補足にはなりますがグラフフーリエ変換で得られたスペクトルについても周波数スペクトルと似たような性質が見られます.
引用元: David I Shuman, Sunil K. Narang, Pascal Frossard, Antonio Ortega and Pierre Vandergheynst. 2013. "The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains", Figure 2.
上に示した図は, あるグラフの3つの固有ベクトル（添え字は固有値の大きさの昇順で与えられています）について各頂点に対応する値を青（正）と黒（負）のバーで示しています. 最も小さい固有値に対応するu_0では, すべての値が正である（周波数が0に対応）のに対し, u_{50}では正の頂点と負の頂点が入り乱れています（高周波に対応）. より具体的にはzero crossing（信号の符号が反転する辺）の数を見るのが良いでしょう. グラフ\mathcal{G}の信号fに関するzero crossing集合は
\mathcal{Z}_{\mathcal{G}}(f)=\{e=(i,j)\in\mathcal{E}|f(i)f(j)\lt 0\}
と定義します. 上図における固有値の大きさとzero crossing集合の濃度の関係は次の図のようになっており, 固有値が大きくなるほどzero crossing集合の濃度も大きくなることが分かります.
引用元: David I Shuman, Sunil K. Narang, Pascal Frossard, Antonio Ortega and Pierre Vandergheynst. 2013. "The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains", Figure 3.
グラフ畳み込み
アプローチの違いにより, グラフ畳み込みにはSpectralな手法とSpatialな手法の2つがあります. 前者は文字通りスペクトル空間上で処理を行う（という考えに基づいた）方法であるのに対し, 後者は注目した頂点の周囲の局所的な情報を集約して処理するような操作を行う方法です. Spectralな方法ではグラフ全体の情報が集約される（逆に言えば, 工夫をしなければ局所性のある処理ができない）他, 固有値分解に伴う計算量の問題があります. Spatialな方法では固有値分解を用いないので計算量が小さく, また処理を直感的に理解しやすいなどのメリットがあります.
Spectralなグラフ畳み込み
先ほどグラフフーリエ変換を導入したのでSpectralな手法から紹介しましょう. 通常の畳み込み演算は
f*g=\int f(\tau)g(t-\tau)d\tau
で定義できます（1次元の場合）. しかし, グラフにおいては頂点間に順序がないため, これに相当する計算を定義することは難しいです. そこで, 畳み込み演算が周波数空間上での積に相当することを思い出し, これを定義に用います.
f*g=U(U^\text{T}f\odot U^\text{T}g)
あるいは, フィルタ\theta_h\in\mathbb{R}^N=U^\text{T}hを直接スペクトル空間上で定義することで, グラフ信号x\in\mathbb{R}^Nとの畳み込みを
\begin{array}{rl}
x*h&=U(U^\text{T}x\odot \theta_h)\\
&=U(\text{diag}(\theta_h)U^\text{T}x)
\end{array}
と定義しても構いません.
Spatialなグラフ畳み込み
本記事の目的はグラフ信号処理とGNNの関係をまとめることですので, グラフ信号処理とあまり関係がない（関係があるものもある）Spatialなグラフ畳み込みは一旦置いておきます. 気が向いたら書くかもしれません. 文献では様々なGNNが日本語でまとめられておりおすすめです.
Discussion