はじめに
グラフニューラルネットワーク(Graph Neural Network : GNN) は文字通りグラフ構造データを扱う深層学習モデルで、大きな注目を集めている。GNNの理解にはspatialとspectralという二つのアプローチがある。理論的なバックグラウンドはspectralの方にあるものの、こちらのアプローチはイメージがしづらく、GNNの説明にはspatialの方が採用されることが多い (以下のようなイメージ図でしばしば説明がなされる)。
本記事では、spectralアプローチの中でも個人的に最もイメージがしづらい"グラフをスペクトル領域へと飛ばす"ということについて実質的にどのような演算がなされているかを説明する。一度スペクトル領域に飛ばすイメージさえできてしまえば、あとは
・周波数成分毎にフィルターをかける
・固有値分解のコストが高いから固有値分解＆フィルタ行列のくだりをグラフラプラシアンの多項式（チェビシェフ：Chebyshev）近似で代用
・ChebNetをさらにK=1近似
という３段階でスペクトラルアプローチを通じたGCN（Graph Convolutional Network）の導出が可能になる。個人的にはこの3段階以上に、グラフをスペクトル領域へと飛ばすことのイメージのしづらさがGNNの理解を妨げているのでは？と感じたため、主としてこの点に絞って解説を行う。本記事はなるべく多く図を用いており、少しでもGNN spectral approachの直感的な理解の手助けとなれば幸いである。
ラプラシアンとグラフラプラシアンのアナロジー
フーリエ変換とはある関数$f(t)$（引数$t$はスカラー、関数$f$の返り値もスカラーとする）を以下の式で周波数$\omega$成分に分解する手法である。
F(\omega) = \int_{-\infty}^{+\infty} f(t) e^{-j\omega t} \,dt
ここで登場する複素正弦波$e^{j\omega t}$は以下で示されるラプラシアン $\Delta \ (= \nabla^2 = \nabla \cdot \nabla = {\rm div \ grad \ }$)の固有関数である。
-\frac{\partial^2}{\partial t^2}e^{j\omega t} = \omega^2e^{j\omega t}
従って、フーリエ変換とはラプラシアンの固有関数$e^{j\omega t}$による$f(t)$の展開である。この展開とは、固有関数が張る波の成分に$f(t)$を分解するという意味である。
これを踏まえ、グラフ構造データをフーリエ変換することを考える。グラフ構造に対してラプラシアン$\Delta$を導入したものを文字通りグラフラプラシアン$L$と呼ぶ。グラフラプラシアンとは、グラフ構造データの接続行列$B$を用いて、$L=BB^T$として定義される行列である。ここではラプラシアンとグラフラプラシアンとでアナロジーを取るため、次の二つの例について考えてみる。
ラプラシアンの例
二次元実数空間$\bf R^2$上で関数$f : {\bf R}^2 \rightarrow {\bf R}$を考える。ある座標$(x, y)$でラプラシアン$\Delta$を作用させると、以下のように表される。
\Delta f(x, y) = \frac{\partial^2f}{\partial x^2}(x, y) + \frac{\partial^2f}{\partial y^2}(x, y)
ここで座標$(x, y)$近傍で関数$f$のテイラー展開を行うと、微小な変数$h$を用いて
\begin{align}
f(x+h, y) = f(x, y) + h\frac{\partial f}{\partial x}(x, y) + \frac{h^2}{2!}\frac{\partial^2f}{\partial x^2}(x, y) + ... \\
f(x, y+h) = f(x, y) + h\frac{\partial f}{\partial y}(x, y) + \frac{h^2}{2!}\frac{\partial^2f}{\partial y^2}(x, y) + ...
\end{align}
と表すことができる。これより$\frac{\partial^2f}{\partial y^2}(x, y)$は
\begin{align}
\frac{\partial^2f}{\partial x^2}(x, y) = \frac{f(x+h ,y) - f(x, y) - h\frac{\partial f}{\partial x}(x, y)}{h^2/2} = \frac{f(x+h ,y) - f(x, y) - h\times (\frac{f(x, y) - f(x-h, y)}{h})}{h^2/2}
\end{align}
なる変換から
\frac{\partial^2f}{\partial x^2}(x, y) = \frac{f(x+h,y) + f(x-h,y) - 2f(x,y)}{h^2/2}
となる。同様にy成分についても求めると、ラプラシアンは
\frac{\partial^2f}{\partial x^2}(x, y) + \frac{\partial^2f}{\partial y^2}(x, y) = \frac{f(x+h,y) + f(x-h,y) + f(x,y+h) + f(x,y-h)- 4f(x,y)}{h^2/2}
と表されることが確認できる。この式（の分子部分）を図示化すると以下のようになり、隣接する４点それぞれについて、座標$(x, y)$との勾配(grad)$\left (f(x\pm h, y) - f(x, y)\right ) / h$ or $\left (f(x, y\pm h) - f(x, y)\right ) / h$を求め、４つの勾配を足し合わせて(=発散)$h(/2)$で割っていることが分かる。
このように、ラプラシアンは、以下の図の四角で囲われた超微小2次元空間からの正味の関数$f$の勾配の湧き出しを意味していることが分かる。この"正味の湧き出し"はしばしば発散(div)と呼ばれ、$\Delta \ (= \nabla^2 = \nabla \cdot \nabla = {\rm div \ grad \ })$のようにラプラシアンは勾配の発散として定義される。
グラフラプラシアンの例
以下の図のようなグラフ構造データを扱う。この例では、各ノード(v1 ~ v5)に6要素からなる特徴ベクトルが付与されており、ノード間に有向のエッジ(e1 ~ e4)が張られているグラフを扱う。
ここで、このグラフを表現するために以下のような変数を定義する。$B$は接続行列、$\bf f$は特徴量行列である。グラフラプラシアンとは、この接続行列$B$を用いて$L=BB^T$なる$L$で表される作用素のことである。
特徴量行列$\bf f$にグラフラプラシアン$L=BB^T$を作用させることが何に対応するのかを以降説明していく。まず$B^T\bf f$により以下の図のようにエッジ毎に特徴量の勾配が求まる。
ここで１つのエッジに対して求まる勾配とは、以下の図に示す特徴量の差分である。
このように各エッジに対して勾配が求まった"エッジ数×特徴量の次元"(今回だと4×6)の行列$B^T\bf f$に対して、$B$を作用させて$BB^T{\bf f} = L\bf f$ を求める。
この処理により各ノードで、"自身から出ていくエッジの勾配の合計"から"自身に入ってくるエッジの勾配の合計"を差し引いた値が得られる。これより、グラフラプラシアンを作用させることは、以下の図のように注目ノード（v3）近傍の四角で囲われた領域からのエッジの勾配の正味の湧き出し(=発散)を求めることを意味しており、通常のラプラシアン同様勾配の発散を求めているということがわかる。
グラフフーリエ変換
ここまで説明したグラフラプラシアンを作用させることの意味を踏まえ、グラフフーリエ変換について説明を行う。通常のフーリエ変換では、入力として与えた関数$f(t)$に対して、ラプラシアン$\Delta=\frac{\partial^2 f}{\partial t^2}$の固有関数$e^{j\omega t}$を乗じて変数$t$方向に足し合わせることで、周波数成分への分解を行った。これは、変数$t$方向の関数$f$の値の分布と固有関数との類似度を算出して周波数成分ごとの重みを算出する。
一方、グラフフーリエ変換ではグラフラプラシアンの固有ベクトルを乗じてノード方向に足し合わせることで、周波数成分への分解を行う。これは、ノード方向の特徴量の分布と固有ベクトルとの類似度を算出して周波数成分ごとの重みを算出する。
以上より、通常のラプラシアンの固有関数のイメージは以下の図（左：矢印下側）、グラフラプラシアンの固有ベクトルのイメージは以下の図（右：矢印下側）のようにまとめることができる。
グラフラプラシアンの固有ベクトル$U$は、グラフラプラシアン$L$の固有値分解$L=U\Lambda U^T$を通じて得られる。固有ベクトルは（制約はあるものの基本的には）ノード数×ノード数の行列となる。グラフは先ほどと同じ例を考えるとすると、$U$, $\Lambda$, $\bf f$は以下のようになる。
(※この図の$U$や$\Lambda$の値はテキトーに入れている。実際に計算して求めた値ではない点ご了承頂きたい。行列の形と、$\Lambda$に関して対角成分以外が0という点だけは真実を反映している)
$U$は各列に固有ベクトル$\bf u$が入っている。今回の例だと全部で５つの固有ベクトルが入っているイメージ。これらについて、$\Lambda$の$U$上での対応する列${\bf u}_{\lambda_i}$の対角成分が対応する固有値$\lambda_i$となっている。グラフフーリエ変換は$U^T \bf f$で求まる (以下の図)。
この図に示すように、固有ベクトル(水色枠部分)とノード方向での特徴の分布(緑色枠部分)との類似度を特徴量の各要素で算出(赤色枠部分)することで、周波数成分に分解することが可能になり、待ち望んでいたグラフをスペクトル領域へと飛ばすことに成功した事になる。
まとめ
グラフをスペクトル領域へと飛ばすことに関して、以下のようにまとめられる。
・グラフラプラシアンもラプラシアンも局所的な勾配の発散を求める作用素であり、データ構造は違うものの同じ概念を表す量である。
・通常の信号データはラプラシアンの固有関数を作用させ、（信号データの）時間方向の分布と固有関数の分布との類似度を該当周波数成分の大きさとする変換を通じて、スペクトル領域に飛ばすことができる。
・グラフ構造データはグラフラプラシアンの固有ベクトルを作用させ、（ノード特徴量の）ノード方向の分布と固有ベクトルの分布との類似度を該当周波数成分の大きさとする変換を通じて、スペクトル領域に飛ばすことができる。
残りの部分（おまけ）
本記事はあくまでスペクトル領域へ飛ばすイメージをまとめたものだが、ついでにGCN（Graph Convolutional Network）までの道のりをざっくり示しておく。グラフ畳み込み（Graph Convolution）では、特徴量$\bf f$にグラフラプラシアン$L=U\Lambda U^T$を作用させ、各ノードでの勾配の発散を求め、特徴量を更新する。この時各周波数成分をどのくらいのバランスで混ぜ合わせるかを学習するのが原始的なスペクトルアプローチでのグラフ畳み込みである。このバランスを求めるフィルター$H$は元の固有値分解して得られた対角行列$\Lambda$を引数として受ける関数である。
一回のグラフ畳み込みでのノード特徴量更新は、$UH(\Lambda)U^T\bf f$なる式で表される。$U^T\bf f$は既に確認してある通り、スペクトル領域に飛ばした後の特徴量行列を意味しており、"固有値の数"（$\lambda_1 \sim \lambda_5$） × "ノード特徴量の次元"（6）の行列である。これに対して周波数成分毎にフィルターをかけ、逆フーリエ変換することで特徴量更新が行われる。これがスペクトル領域での原始的なグラフ畳み込みの説明となる。
その後固有値分解の計算コストの高さを解消するために、固有値分解を回避してグラフラプラシアンの多項式（K次の項まで）としてノード特徴量を更新するChebNetが登場し、さらにこの次数の上限Kを1まで下げることで精度を向上させたGraph Convolutional Network (GCN)が誕生する。
参考
以下の記事がいずれも大変わかりやすく参考にさせて頂いた。
・https://www.jstage.jst.go.jp/article/essfr/8/1/8_15/_pdf
・https://hackmd.io/@kkume/rkK3tmpHd
・https://www.slideshare.net/tm1966/graph-neural-networks
・https://qiita.com/shionhonda/items/d27b8f13f7e9232a4ae5
・https://qiita.com/suepiiin/items/4815be134c5424cc9493
・https://buildersbox.corp-sansan.com/entry/2021/05/27/110000
・https://ogyahogya.hatenablog.com/entry/2018/01/26/%E3%82%B0%E3%83%A9%E3%83%95%E3%83%A9%E3%83%97%E3%83%A9%E3%82%B7%E3%82%A2%E3%83%B3
・https://tm23forest.com/contents/graph-laplacian-incidence-matrix
・https://qiita.com/silva0215/items/0d1d25ef51b6865a6e15#%E6%8E%A5%E7%B6%9A%E8%A1%8C%E5%88%97
・https://rikeinotame.com/fourier2/