畳み込みニューラルネットワーク(CNN)では、層を深くするほど性能が上がるのが常識でした。ResNetが152層を積んでImageNetを制したように、「深さ=表現力」という直感は画像の世界では正しかったのです。ところが、グラフニューラルネットワーク(GNN)で同じことをしようとすると、奇妙なことが起こります。2層や3層では良い精度が出ていたのに、4層、8層、16層と深くしていくと、精度がむしろガクンと落ちるのです。しかも訓練データに対してすら性能が出ない。これは単なる過学習でも勾配消失でもありません。
層を深くしたGNNのノード表現を取り出して眺めてみると、その理由がわかります。最初はノードごとにバラバラだった特徴ベクトルが、層を重ねるごとに互いに似てきて、最終的にはグラフ上のすべてのノードがほとんど同じベクトルになってしまうのです。空手クラブの会員Aさんと会員Bさんが、本来は別のコミュニティに属しているのに、深い層を通すと区別がつかなくなる。これが過平滑化(over-smoothing)と呼ばれる現象です。ノードを分類したいのに、全ノードが同じ表現に潰れてしまっては分類のしようがありません。
なぜこれを学ぶ価値があるのでしょうか。第一に、過平滑化はGNNを「深くできない」根本的な制約であり、これを理解しないと「なぜ自分のGNNは2層で頭打ちなのか」という実務上の謎が解けません。第二に、過平滑化の数理はグラフラプラシアンの固有値と拡散方程式の定常状態という、物理・信号処理にも通じる美しい構造を持っています。応用先も広く、(1)ソーシャルネットワークの遠距離依存を捉えたい推薦システムや不正検知では深いGNNが必要になり過平滑化対策が不可欠ですし、(2)分子グラフで原子間の長距離相互作用を学ぶ材料・創薬の分野でも、層を深くしながら表現を保つ技術が直接効いてきます。
本記事では、過平滑化を「なんとなく似てくる」という曖昧な話で終わらせず、正規化隣接行列の繰り返し積が定常分布へ収束する過程を固有値分解で厳密に追い、ディリクレエネルギーという量が層を経るごとに指数関数的に減衰することを導出します。そのうえで、Karate Club / Cora データセットでGCNを2/4/8/16層に積み、ノード表現の類似度行列とディリクレエネルギーを可視化して過平滑化を「目で見て」確認します。最後に、PairNorm・初期残差(GCNII風)・Jumping Knowledge という3つの代表的な緩和策をPyTorchで実装し、深層化しても精度が維持されることを定量比較します。
本記事の内容
- 過平滑化現象の直感と、なぜ深層GNNで起こるのか
- 正規化隣接行列の固有値分解による定常状態への収束の導出
- ディリクレエネルギー $E(\bm{H}^l) = \mathrm{tr}(\bm{H}^\top \bm{L} \bm{H})$ の指数減衰の証明
- PyTorchでの過平滑化の可視化(類似度行列・エネルギー減衰)
- PairNorm・初期残差・Jumping Knowledge の実装と精度比較
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が格段に深まります。特にGCNの伝播則とグラフラプラシアンの固有値は、本記事の導出の土台になります。
- グラフニューラルネットワーク(GNN)とは — GNNの全体像とメッセージパッシングの基礎
- GCN(グラフ畳み込みネットワーク)を理解して実装する — 本記事が対象とするGCNの伝播則の導出
- グラフラプラシアンとは — ラプラシアンの定義・固有値・スペクトルの性質
- グラフアテンションネットワーク(GAT)を理解する — 注意機構によるメッセージパッシングの一般化
過平滑化とは — 拡散のたとえで掴む
過平滑化を直感的に理解するには、インクの拡散を思い浮かべるのが一番です。水を張った水槽の各地点に、色とりどりのインクをポツポツと垂らしたとします。最初は赤い点、青い点、緑の点がはっきり分かれています。ところが時間が経つにつれてインクは隣の水へ、さらにその隣へとにじんでいき、十分に時間が経つと水槽全体が均一な「混ぜこぜの色」になってしまいます。どこを見ても同じ色、つまり位置による違いが消えるわけです。
GNNの1層は、まさにこの「にじみ」を1回行う操作です。GCNの伝播則を思い出すと、各ノードの新しい特徴は「自分と隣接ノードの特徴の(重み付き)平均」を取って作られます。隣の情報を混ぜ込むこの操作は、グラフという離散的な空間の上での拡散にほかなりません。1層なら隣接ノードの情報が、2層なら2ホップ先の情報が、$l$層なら$l$ホップ先の情報が混ざり込みます。これは「広い範囲の文脈を見る」という意味では望ましいのですが、混ぜすぎるとインクと同じで、すべてのノードが「グラフ全体の平均」に収束してしまうのです。
ここで重要なのは、これはバグでも実装ミスでもなく、平均化を繰り返すことの数学的必然だという点です。区別したい情報(ノードごとの違い)を、平均化という操作が組織的に破壊していく。CNNの畳み込みは局所的なフィルタで特徴を「抽出」しますが、GCNの畳み込みは隣接の「平均」を取るため、深くするほど高周波成分(ノード間の細かい違い)が削られ、低周波成分(全体の傾向)だけが残ります。これは信号処理でいうローパスフィルタを何度もかけるのと同じ構図です。
では、この「にじみ」がどのくらいの速さで進み、最終的にどこへ収束するのか。それを正確に知るには、GCNの平均化操作の正体である正規化隣接行列を固有値分解して、その繰り返し積の振る舞いを調べる必要があります。次のセクションでその数理に踏み込みましょう。
GCNの伝播則と正規化隣接行列
まず、議論の対象を数式で固定します。$N$個のノードを持つグラフを考え、隣接行列を $\bm{A} \in \mathbb{R}^{N\times N}$、次数行列を $\bm{D} = \mathrm{diag}(d_1, \dots, d_N)$($d_i = \sum_j A_{ij}$)とします。GCN(Kipf & Welling, 2017)では、自己ループを加えた $\tilde{\bm{A}} = \bm{A} + \bm{I}$ と、その次数行列 $\tilde{\bm{D}}$ を使って、対称正規化隣接行列を次のように定義します。
$$ \begin{equation} \hat{\bm{A}} = \tilde{\bm{D}}^{-1/2}\tilde{\bm{A}}\tilde{\bm{D}}^{-1/2} \end{equation} $$
この $\hat{\bm{A}}$ こそが「平均化」の正体です。GCNの第$l$層の伝播則は、ノード特徴行列 $\bm{H}^{(l)} \in \mathbb{R}^{N\times F}$、重み $\bm{W}^{(l)}$、活性化関数 $\sigma$ を使って
$$ \begin{equation} \bm{H}^{(l+1)} = \sigma\!\left(\hat{\bm{A}}\,\bm{H}^{(l)}\bm{W}^{(l)}\right) \end{equation} $$
と書けます。ここで $\bm{H}^{(0)} = \bm{X}$ は入力特徴です。$\hat{\bm{A}}\bm{H}^{(l)}$ という積が、各ノードについて「隣接ノードの特徴を次数で重み付けして集約する」操作に対応します。$\bm{W}^{(l)}$ は特徴を線形変換し、$\sigma$ は非線形性を与えます。
過平滑化の本質を見るために、ここで思い切った単純化をします。重み行列を恒等 $\bm{W}^{(l)} = \bm{I}$、活性化を線形 $\sigma(z)=z$ とおくのです。これは「学習を止めて、純粋な伝播だけを見る」ことに相当します。すると伝播則は単純な行列の繰り返し積になります。
$$ \begin{equation} \bm{H}^{(l)} = \hat{\bm{A}}^{l}\,\bm{X} \end{equation} $$
つまり層を$l$回通すことは、正規化隣接行列を$l$乗して入力にかけることと同じです。「重みや非線形を無視するのは乱暴では?」と思うかもしれませんが、後述するように非線形 ReLU や学習可能な重みを入れても、過平滑化の本質的な傾向は変わらないことが知られています(ReLU は非負性を加えるだけで平均化の性質は保たれ、$\bm{W}$ のスペクトルノルムが1以下なら収束はむしろ加速します)。したがって、$\hat{\bm{A}}^{l}$ が $l\to\infty$ でどう振る舞うかを調べれば、過平滑化の核心がつかめます。
この $\hat{\bm{A}}^{l}$ の極限を知るための鍵は、$\hat{\bm{A}}$ を固有値分解することです。次のセクションで $\hat{\bm{A}}$ の固有値がどんな範囲に収まり、その繰り返し積が何に収束するかを導出します。
正規化隣接行列の固有値分解と定常状態
行列の累乗の振る舞いを調べる定石は、対角化です。$\hat{\bm{A}}$ は対称行列なので(定義式の $\tilde{\bm{D}}^{-1/2}\tilde{\bm{A}}\tilde{\bm{D}}^{-1/2}$ は $\tilde{\bm{A}}$ が対称なら対称)、必ず直交行列で対角化できます。すなわち、実固有値 $\lambda_1, \dots, \lambda_N$ と正規直交固有ベクトル $\bm{u}_1, \dots, \bm{u}_N$ を使って
$$ \begin{equation} \hat{\bm{A}} = \bm{U}\bm{\Lambda}\bm{U}^\top = \sum_{i=1}^{N}\lambda_i \bm{u}_i\bm{u}_i^\top \end{equation} $$
と書けます。ここで $\bm{U} = [\bm{u}_1, \dots, \bm{u}_N]$ は直交行列($\bm{U}^\top\bm{U} = \bm{I}$)、$\bm{\Lambda} = \mathrm{diag}(\lambda_1, \dots, \lambda_N)$ です。
固有値の範囲
過平滑化を理解するうえで決定的に重要なのが、$\hat{\bm{A}}$ の固有値がすべて区間 $(-1, 1]$ に収まり、最大固有値はちょうど1であるという事実です。これを確認しましょう。$\hat{\bm{A}}$ と正規化グラフラプラシアン $\hat{\bm{L}} = \bm{I} – \hat{\bm{A}}$ は同じ固有ベクトルを共有し、固有値の間には次の関係があります。
$$ \begin{equation} \hat{\bm{L}}\bm{u}_i = (1-\lambda_i)\bm{u}_i \quad\Longleftrightarrow\quad \mu_i = 1 – \lambda_i \end{equation} $$
ここで $\mu_i$ は $\hat{\bm{L}}$ の固有値です。正規化ラプラシアンの固有値は $0 \le \mu_i \le 2$ の範囲にあることが知られています(これはラプラシアンが半正定値であることと、隣接性から従う標準的な結果です)。$\mu_i = 1-\lambda_i$ を代入すると、$\lambda_i = 1 – \mu_i$ は $-1 \le \lambda_i \le 1$ の範囲に収まることがわかります。
特に最小固有値 $\mu_1 = 0$ に対応して、最大固有値 $\lambda_1 = 1$ が存在します。$\hat{\bm{L}}$ の固有値0に対応する固有ベクトルは、自己ループ付きグラフが連結なら一意で、各成分が $\sqrt{\tilde{d}_i}$ に比例するベクトルです($\tilde{d}_i$ は自己ループ込みの次数)。規格化すると
$$ \begin{equation} \bm{u}_1 = \frac{1}{\sqrt{\sum_j \tilde{d}_j}}\left(\sqrt{\tilde{d}_1}, \sqrt{\tilde{d}_2}, \dots, \sqrt{\tilde{d}_N}\right)^\top \end{equation} $$
となります。これが「拡散の定常状態」を担う固有ベクトルです。直感的には、次数の平方根に比例した「重み付き一様分布」であり、インクが拡散しきった後の均一な色に相当します。
累乗の極限
固有値の範囲がわかれば、累乗の極限は一瞬で出ます。固有値分解 $\hat{\bm{A}} = \bm{U}\bm{\Lambda}\bm{U}^\top$ を$l$乗すると、直交性 $\bm{U}^\top\bm{U}=\bm{I}$ により中間項がすべて打ち消し合い、
$$ \begin{equation} \hat{\bm{A}}^{l} = \bm{U}\bm{\Lambda}^{l}\bm{U}^\top = \sum_{i=1}^{N}\lambda_i^{l}\,\bm{u}_i\bm{u}_i^\top \end{equation} $$
となります。ここで $\bm{\Lambda}^l = \mathrm{diag}(\lambda_1^l, \dots, \lambda_N^l)$ です。固有値を大きい順に並べ替え($1 = \lambda_1 > \lambda_2 \ge \dots \ge \lambda_N > -1$、連結グラフでは $\lambda_1$ は単純で他は厳密に1未満)、$l\to\infty$ の極限を取ると何が起こるか見てみましょう。$|\lambda_i| < 1$ である $i\ge 2$ の項は $\lambda_i^l \to 0$ となって消え、$\lambda_1 = 1$ の項だけが生き残ります。
$$ \begin{equation} \lim_{l\to\infty}\hat{\bm{A}}^{l} = \lambda_1^l\bm{u}_1\bm{u}_1^\top \Big|_{\lambda_1=1} = \bm{u}_1\bm{u}_1^\top \end{equation} $$
したがって、入力 $\bm{X}$ に十分多くの層をかけると
$$ \begin{equation} \lim_{l\to\infty}\bm{H}^{(l)} = \bm{u}_1\bm{u}_1^\top\bm{X} = \bm{u}_1\,(\bm{u}_1^\top\bm{X}) \end{equation} $$
となります。右辺をよく見ると、$\bm{u}_1^\top\bm{X}$ は$1\times F$の行ベクトル(各特徴次元についてのスカラー)で、それに列ベクトル $\bm{u}_1$ を掛けています。つまり最終的なノード表現は、すべて $\bm{u}_1$ のスカラー倍になります。ノード$i$の表現は $\bm{u}_1$ の第$i$成分 $\sqrt{\tilde{d}_i}$ にだけ依存し、もとの特徴 $\bm{X}$ の情報は $\bm{u}_1^\top\bm{X}$ という1次元への射影に圧縮されてしまう。これが過平滑化の数学的な正体です。すべてのノード表現が、次数で決まる1本の方向 $\bm{u}_1$ に潰れるのです。
「次数だけで決まる表現」では、同じ次数を持つノードはもはや区別できません。ノード分類で大事なのは「どのコミュニティに属するか」であって「次数がいくつか」ではないことが多いので、これは致命的です。では、この収束がどのくらい速いのか、層数とともに「区別のなさ」がどう増えていくのかを定量化したい。そのための物差しが、次に導入するディリクレエネルギーです。
ディリクレエネルギーによる過平滑化の定量化
「全ノードが似てきた」を一つの数値で測りたい。そこで使うのがディリクレエネルギーです。直感的には、これは「隣り合うノードの表現がどれだけ違うか」の総和です。隣接ノード同士が似ていればエネルギーは小さく、バラバラなら大きい。過平滑化とは「隣接ノードが似てくる」ことなので、ディリクレエネルギーが小さくなることと過平滑化が進むことは同じ意味になります。
特徴行列 $\bm{H}\in\mathbb{R}^{N\times F}$ に対するディリクレエネルギーを次のように定義します。
$$ \begin{equation} E(\bm{H}) = \mathrm{tr}\!\left(\bm{H}^\top \hat{\bm{L}}\, \bm{H}\right) = \frac{1}{2}\sum_{i,j}\hat{A}_{ij}\left\|\frac{\bm{h}_i}{\sqrt{\tilde{d}_i}} – \frac{\bm{h}_j}{\sqrt{\tilde{d}_j}}\right\|^2 \end{equation} $$
ここで $\hat{\bm{L}} = \bm{I} – \hat{\bm{A}}$ は正規化ラプラシアン、$\bm{h}_i$ はノード$i$の特徴(行ベクトル)です。右辺の和の形は、ラプラシアン二次形式の標準的な展開で、「正規化した隣接ノードの表現の差の二乗を、エッジ重みで重み付けして足し上げたもの」だとわかります。エネルギーが0になるのは、すべての隣接ノードペアで $\bm{h}_i/\sqrt{\tilde{d}_i} = \bm{h}_j/\sqrt{\tilde{d}_j}$ が成り立つとき、つまり表現が定常状態 $\bm{u}_1$ 方向に潰れたときです。
二次形式表現の導出
定義の左辺(トレース表現)と右辺(和の表現)が等しいことを確認しておきましょう。$\hat{\bm{L}} = \bm{I} – \hat{\bm{A}}$ を代入し、トレースの定義 $\mathrm{tr}(\bm{H}^\top\bm{M}\bm{H}) = \sum_{i,j}M_{ij}\,\bm{h}_i^\top\bm{h}_j$ を使います。
$$ \begin{equation} E(\bm{H}) = \sum_{i}\|\bm{h}_i\|^2 – \sum_{i,j}\hat{A}_{ij}\,\bm{h}_i^\top\bm{h}_j \end{equation} $$
ここで $\hat{A}_{ij} = \tilde{A}_{ij}/\sqrt{\tilde{d}_i\tilde{d}_j}$ であり、対角正規化により $\sum_j \hat{A}_{ij}\cdot(\text{適切な重み})$ が扱いやすくなります。標準的な変形(差の二乗を展開して交差項を集める操作)を施すと、上式は右辺の「差の二乗和」と一致します。具体的には、$\frac{1}{2}\sum_{i,j}\hat{A}_{ij}\|\bm{h}_i/\sqrt{\tilde d_i} – \bm{h}_j/\sqrt{\tilde d_j}\|^2$ を展開すると、二乗項から $\sum_i \|\bm{h}_i\|^2$(対角正規化のもとで$\hat A$の行和が1に対応する成分)、交差項から $-\sum_{i,j}\hat{A}_{ij}\bm{h}_i^\top\bm{h}_j/\sqrt{\tilde d_i\tilde d_j}$ が出てきて、トレース表現に帰着します。要点は、ディリクレエネルギーが半正定値($E(\bm{H})\ge 0$)であり、$\hat{\bm{L}}$ の固有値 $\mu_i\ge 0$ がすべて非負であることと整合する、という点です。
固有モード分解
ディリクレエネルギーがなぜ層を経るごとに減衰するのかを見るには、$\bm{H}$ を $\hat{\bm{L}}$ の固有ベクトルで展開するのが最良です。特徴行列の各列(特徴次元)を固有ベクトル基底で展開し、$\bm{H} = \sum_i \bm{u}_i \bm{c}_i^\top$ と書きます($\bm{c}_i^\top = \bm{u}_i^\top\bm{H}$ は$1\times F$の係数)。これをエネルギーの定義に代入すると、直交性 $\bm{u}_i^\top\bm{u}_j = \delta_{ij}$ と $\hat{\bm{L}}\bm{u}_i = \mu_i\bm{u}_i$ により、
$$ \begin{equation} E(\bm{H}) = \mathrm{tr}\!\left(\sum_{i,j}\bm{c}_i\bm{u}_i^\top\hat{\bm{L}}\bm{u}_j\bm{c}_j^\top\right) = \sum_{i=1}^{N}\mu_i\,\|\bm{c}_i\|^2 \end{equation} $$
と、固有モードごとの寄与の和に分解されます。各モードのエネルギーへの寄与は「ラプラシアン固有値 $\mu_i$ × そのモードの強さ $\|\bm{c}_i\|^2$」です。固有値0の定常モード($i=1$)はエネルギーにまったく寄与せず、高周波モード(大きい $\mu_i$)ほど大きく寄与します。これで「エネルギー=高周波成分の量」という直感が数式で裏付けられました。
次は、伝播 $\bm{H}\mapsto\hat{\bm{A}}\bm{H}$ を一回行うとこのエネルギーがどれだけ減るかを計算し、層数に対する指数減衰を導きます。
ディリクレエネルギーの指数減衰の導出
ゴールを明確にします。「GCNの1層を通すたびに、ディリクレエネルギーは少なくとも一定の比率で減少し、層数$l$に対して指数的に0へ向かう」ことを示します。これが過平滑化が避けられない理由の定量的な核心です。
線形・恒等重みの単純化のもとで、1層の伝播は $\bm{H}^{(l+1)} = \hat{\bm{A}}\bm{H}^{(l)}$ です。$\hat{\bm{A}} = \bm{I} – \hat{\bm{L}}$ なので、伝播後の特徴を固有モードで展開すると、固有値の関係 $\hat{\bm{A}}\bm{u}_i = \lambda_i\bm{u}_i = (1-\mu_i)\bm{u}_i$ から、係数が $\bm{c}_i \mapsto (1-\mu_i)\bm{c}_i$ と変換されます。これをエネルギーの固有モード分解に代入します。
$$ \begin{equation} E(\hat{\bm{A}}\bm{H}^{(l)}) = \sum_{i=1}^{N}\mu_i\,(1-\mu_i)^2\,\|\bm{c}_i^{(l)}\|^2 \end{equation} $$
一方、伝播前のエネルギーは $E(\bm{H}^{(l)}) = \sum_i \mu_i\|\bm{c}_i^{(l)}\|^2$ でした。両者を各モードごとに比べると、伝播後はモード$i$の寄与が $(1-\mu_i)^2$ 倍されていることがわかります。ここで、定常モード以外のすべての非ゼロ固有値について最大の減衰係数を考えます。$\mu_1 = 0$ を除く固有値の中で、$\lambda_i^2 = (1-\mu_i)^2$ の最大値を $s = \max_{i\ge 2}(1-\mu_i)^2 < 1$ と定義しましょう。$s$ は連結グラフでは厳密に1未満です(なぜなら $i\ge 2$ では $\mu_i>0$ かつ $\mu_i<2$ で、$|1-\mu_i|<1$ だから)。
ここがポイントです。エネルギーは定常モード($\mu_1=0$、寄与ゼロ)を含まないので、和は実質 $i\ge 2$ のモードだけで取られます。それらのモードの減衰係数はすべて $s$ 以下なので、
$$ \begin{equation} E(\hat{\bm{A}}\bm{H}^{(l)}) = \sum_{i\ge 2}\mu_i(1-\mu_i)^2\|\bm{c}_i^{(l)}\|^2 \le s\sum_{i\ge 2}\mu_i\|\bm{c}_i^{(l)}\|^2 = s\,E(\bm{H}^{(l)}) \end{equation} $$
という不等式が成り立ちます。1層で少なくとも $s$ 倍に縮むのですから、$l$層を通せば
$$ \begin{equation} E(\bm{H}^{(l)}) \le s^{l}\,E(\bm{H}^{(0)}) \end{equation} $$
となり、$s<1$ ゆえ $l\to\infty$ で $E(\bm{H}^{(l)})\to 0$ が指数的な速さで起こることが示されました。減衰率 $s$ は固有値分布で決まり、$s = \max(\lambda_2^2, \lambda_N^2)$ です。特に $\lambda_2$ はスペクトルギャップ $1-\lambda_2$ と関係し、グラフがよくつながっている(スペクトルギャップが大きい)ほど $s$ が小さく、過平滑化が速く進むことを意味します。これは「混ざりやすいグラフほど早く均一になる」という拡散の直感と完全に一致します。
なお、活性化に ReLU を入れ、重み $\bm{W}^{(l)}$ を学習可能にしても、$\bm{W}^{(l)}$ のスペクトルノルム(最大特異値)を $\beta$ とすると、エネルギーの減衰は概ね $E(\bm{H}^{(l)}) \le (s\beta^2)^l E(\bm{H}^{(0)})$ という形に修正されます(ReLU はエネルギーを増やさないことが示せます)。実用上 $\beta$ が大きすぎると勾配が爆発するため適度に抑えられ、結局 $s\beta^2<1$ となって過平滑化が生じやすい、という構図は変わりません。
理論はこれで尽くせました。次は実際に Karate Club / Cora で GCN を深くし、いま導いた「類似度の上昇」と「ディリクレエネルギーの指数減衰」を自分の目で確かめます。
Pythonで過平滑化を可視化する
理論で予測した2つの現象 — (1)層を重ねるとノード表現の類似度が上がる、(2)ディリクレエネルギーが指数的に減衰する — を実際に確認します。まずは依存ライブラリと、グラフ・特徴の準備から始めます。PyTorch Geometric を使わず、NumPy / PyTorch だけで完結させて中身を透明にします。
import numpy as np
import torch
import torch.nn as nn
import torch.nn.functional as F
import matplotlib.pyplot as plt
torch.manual_seed(0)
np.random.seed(0)
# Zachary の空手クラブグラフ(34ノード)の隣接行列を構築
# networkx があれば利用、なければエッジリストを直接定義
try:
import networkx as nx
G = nx.karate_club_graph()
N = G.number_of_nodes()
A = nx.to_numpy_array(G)
# コミュニティ(2分裂)を正解ラベルに
labels = np.array([G.nodes[i]['club'] == 'Officer' for i in G.nodes()], dtype=int)
except Exception:
# フォールバック: 簡易な2コミュニティのSBMグラフ
N = 34
labels = np.array([0]*17 + [1]*17)
A = np.zeros((N, N))
for i in range(N):
for j in range(i+1, N):
p = 0.4 if labels[i] == labels[j] else 0.05
if np.random.rand() < p:
A[i, j] = A[j, i] = 1
print("ノード数:", N, " エッジ数:", int(A.sum() // 2))
ここでは空手クラブグラフ(34ノード、2つのコミュニティに分裂する有名なソーシャルネットワーク)を使います。labels は各ノードがどちらの派閥に属するかを表し、これが本来GNNで区別したい情報です。エッジ数の出力からグラフの密度感がわかり、コミュニティ間より内部のつながりが強い構造になっていることが、後の過平滑化の進み方に効いてきます。
次に、正規化隣接行列 $\hat{\bm{A}}$ と正規化ラプラシアン $\hat{\bm{L}}$ を作り、ディリクレエネルギーを計算する関数を用意します。
def normalize_adj(A):
"""自己ループ付き対称正規化隣接行列 Â = D̃^{-1/2}(A+I)D̃^{-1/2} を返す"""
A_tilde = A + np.eye(A.shape[0])
d = A_tilde.sum(axis=1)
d_inv_sqrt = 1.0 / np.sqrt(d)
D_inv_sqrt = np.diag(d_inv_sqrt)
return D_inv_sqrt @ A_tilde @ D_inv_sqrt
A_hat = normalize_adj(A)
L_hat = np.eye(N) - A_hat # 正規化ラプラシアン
A_hat_t = torch.tensor(A_hat, dtype=torch.float32)
L_hat_t = torch.tensor(L_hat, dtype=torch.float32)
# Â の固有値を確認(理論: すべて (-1, 1]、最大はちょうど1)
eigvals = np.linalg.eigvalsh(A_hat)
print("Â の固有値の最大:", eigvals.max().round(4),
" 最小:", eigvals.min().round(4))
def dirichlet_energy(H, L):
"""ディリクレエネルギー E(H) = tr(Hᵀ L H) を計算"""
return torch.trace(H.t() @ L @ H).item()
固有値の出力を見ると、最大固有値が理論どおり1.0(数値誤差の範囲で)、最小固有値が$-1$より大きい値になっているはずです。これは前のセクションで導いた「$\hat{\bm{A}}$ の固有値は $(-1,1]$ に収まり、最大はちょうど1」という主張の数値的な裏付けです。最大固有値1に対応する固有モードが、層を重ねたときに生き残る定常状態 $\bm{u}_1$ にあたります。
純粋な伝播による過平滑化
学習を一切せず、正規化隣接行列をひたすら掛けるだけで何が起こるかを見ます。これは式 $\bm{H}^{(l)} = \hat{\bm{A}}^l\bm{X}$ をそのまま実装したものです。各層でディリクレエネルギーと、ノード表現どうしの平均コサイン類似度を記録します。
# 入力特徴: ランダムな16次元ベクトル(ノードごとに独立)
X = torch.randn(N, 16)
def mean_cosine_similarity(H):
"""全ノードペアの平均コサイン類似度(1に近いほど過平滑化)"""
Hn = F.normalize(H, p=2, dim=1)
S = Hn @ Hn.t()
# 対角(自己類似度=1)を除いた平均
off = S - torch.diag(torch.diag(S))
return off.sum().item() / (H.shape[0] * (H.shape[0] - 1))
H = X.clone()
energies, sims = [dirichlet_energy(H, L_hat_t)], [mean_cosine_similarity(H)]
for l in range(1, 21):
H = A_hat_t @ H # 1層分の純粋な伝播
energies.append(dirichlet_energy(H, L_hat_t))
sims.append(mean_cosine_similarity(H))
print("層0のエネルギー:", round(energies[0], 4),
" 層20のエネルギー:", round(energies[20], 6))
print("層0の平均類似度:", round(sims[0], 4),
" 層20の平均類似度:", round(sims[20], 4))
出力では、層0で大きかったディリクレエネルギーが層20ではほぼ0に潰れ、平均コサイン類似度が層0の0付近(ランダムなので無相関)から層20では1に近い値へ上昇しているはずです。これはまさに理論の予言どおりで、「エネルギーの指数減衰」と「全ノード表現が同一方向 $\bm{u}_1$ に揃う(類似度→1)」が同時に起きていることを示しています。学習も非線形も無いのに過平滑化が生じる点が重要で、これが平均化操作そのものの性質であることを物語っています。
エネルギーの減衰が本当に「指数的」なのかを、片対数プロットで確かめましょう。
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
# 左: ディリクレエネルギーの片対数プロット(指数減衰なら直線)
axes[0].semilogy(range(21), energies, 'o-', color='crimson')
axes[0].set_xlabel('Layer (propagation step) l')
axes[0].set_ylabel('Dirichlet energy E(H) [log scale]')
axes[0].set_title('Exponential decay of Dirichlet energy')
axes[0].grid(True, alpha=0.3)
# 右: 平均コサイン類似度(1に近づくほど過平滑化)
axes[1].plot(range(21), sims, 's-', color='navy')
axes[1].set_xlabel('Layer (propagation step) l')
axes[1].set_ylabel('Mean pairwise cosine similarity')
axes[1].set_title('Node representations collapse to one direction')
axes[1].set_ylim(-0.1, 1.05)
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('oversmoothing_decay.png', dpi=150, bbox_inches='tight')
plt.show()
左のグラフが片対数軸でほぼ直線になっていれば、$E(\bm{H}^{(l)})\le s^l E(\bm{H}^{(0)})$ という指数減衰の証拠です。直線の傾きが減衰率 $\log s$ に対応し、空手クラブのようによくつながったグラフでは傾きが急(=過平滑化が速い)になります。右のグラフでは類似度が単調に1へ近づき、十数層で「全ノードがほぼ同じ向き」に潰れる様子が確認できます。理論で導いた定常状態への収束が、わずか10〜15層で実用上完了してしまうことがわかります。
類似度行列のヒートマップ
「潰れていく」様子をもっと直接見るために、特定の層でのノード間コサイン類似度行列をヒートマップにします。初期はコミュニティ構造(ブロック対角)が見えていても、深くすると全面が真っ赤(全ペア類似度≈1)になるはずです。
# 並べ替え: 同じコミュニティのノードを隣接させて構造を見やすく
order = np.argsort(labels)
def sim_matrix(H):
Hn = F.normalize(H, p=2, dim=1)
return (Hn @ Hn.t()).detach().numpy()[order][:, order]
layers_to_show = [0, 2, 8]
fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))
H = X.clone()
mats = {0: sim_matrix(H)}
for l in range(1, 9):
H = A_hat_t @ H
if l in layers_to_show:
mats[l] = sim_matrix(H)
for ax, l in zip(axes, layers_to_show):
im = ax.imshow(mats[l], vmin=-1, vmax=1, cmap='copper'.replace('o','o')
if False else 'RdBu_r')
ax.set_title(f'Cosine similarity matrix at layer {l}')
ax.set_xlabel('node (sorted by community)')
fig.colorbar(im, ax=ax, fraction=0.046)
plt.tight_layout()
plt.savefig('oversmoothing_simmat.png', dpi=150, bbox_inches='tight')
plt.show()
層0のヒートマップでは対角ブロック(同一コミュニティ)がやや明るい程度ですが、層2では同じコミュニティのノードが強く似はじめ、層8ではほぼ全セルが濃い赤(類似度≈1)になって、コミュニティの境界が完全に溶けてしまいます。これがノード分類を不可能にする過平滑化の「見た目」です。2層程度では適度に近傍情報を集約して有益ですが、8層では情報が均一化されて識別力を失う、という深さのトレードオフが一目で理解できます。
ここまでで過平滑化が「起こる」ことと「速い」ことを確認しました。では、深いGNNでも表現を保つにはどうすればよいか。次のセクションで、ディリクレエネルギーを保つよう設計された3つの緩和策を実装し、効果を定量的に比較します。
過平滑化への対策 — 理論と実装
過平滑化の原因は「平均化の繰り返しが高周波(ノード間の差)を削る」ことでした。対策の発想は大きく3つに分かれます。(1)各層で表現を再正規化してエネルギーを補充する(PairNorm)、(2)入力の情報を毎層注入して定常状態への収束を妨げる(初期残差・GCNII)、(3)各層の出力を全部保存して最後に統合し、過平滑化した深い層に依存しすぎないようにする(Jumping Knowledge)。順に理論的な意味づけと実装を見ていきます。
PairNorm
PairNorm(Zhao & Akoglu, 2020)のアイデアは単純で強力です。各層の後で、ノード表現の全体の総ペア距離(=ディリクレエネルギーに比例する量)を一定に保つように再スケールします。過平滑化がエネルギーを減らすなら、人為的にエネルギーを元の水準へ戻してやればよい、という発想です。具体的には、まず中心化(平均を引く)してから、全体のスケールを目標値に正規化します。
$$ \begin{equation} \tilde{\bm{h}}_i = \bm{h}_i – \frac{1}{N}\sum_{j}\bm{h}_j, \qquad \bm{h}_i^{\text{out}} = s\cdot\frac{\tilde{\bm{h}}_i}{\sqrt{\frac{1}{N}\sum_{j}\|\tilde{\bm{h}}_j\|^2}} \end{equation} $$
第1式の中心化は定常モード $\bm{u}_1$(おおむね一様成分)を除去する役割を持ち、第2式で全ノードの二乗ノルムの平均が一定値 $s^2$ になるよう正規化します。これにより、平均化で潰れようとする表現が毎層スケールを取り戻し、ノード間の相対的な差が保たれます。
class PairNorm(nn.Module):
"""総ペア距離(=エネルギー)を一定に保つ正規化"""
def __init__(self, scale=1.0):
super().__init__()
self.scale = scale
def forward(self, x):
x = x - x.mean(dim=0, keepdim=True) # 中心化
norm = x.pow(2).sum(dim=1).mean().sqrt() + 1e-6 # 全体スケール
return self.scale * x / norm
このモジュールは各GCN層の活性化の後ろに挟むだけで機能します。中心化が定常成分を削り、スケール正規化がエネルギーの目減りを打ち消すため、深い層でもノード表現が一点に潰れにくくなります。次に紹介する初期残差は、まったく別の角度から同じ問題を攻めます。
初期残差(GCNII風)
GCNII(Chen et al., 2020)の核心は、毎層で必ず入力特徴 $\bm{H}^{(0)}$ を混ぜ込む初期残差接続です。拡散の比喩で言えば、にじんでいくインクに対して「元の色を毎回少しずつ足し続ける」操作で、定常状態への完全な収束を妨げます。第$l+1$層の伝播を次のように修正します。
$$ \begin{equation} \bm{H}^{(l+1)} = \sigma\!\Big(\big((1-\alpha)\hat{\bm{A}}\bm{H}^{(l)} + \alpha\bm{H}^{(0)}\big)\,\big((1-\beta_l)\bm{I} + \beta_l\bm{W}^{(l)}\big)\Big) \end{equation} $$
第1の括弧の $\alpha\bm{H}^{(0)}$ が初期残差で、$\alpha$(例: 0.1〜0.2)の割合で必ず入力情報を保ちます。これにより、純粋伝播なら $\bm{u}_1$ に潰れるはずの表現が、常に $\bm{H}^{(0)}$ の成分を含み続け、過平滑化が抑えられます。第2の括弧 $((1-\beta_l)\bm{I} + \beta_l\bm{W}^{(l)})$ は恒等写像への接続(アイデンティティマッピング)で、$\beta_l$ を層が深いほど小さくすることで、深層でも重みが恒等に近づき過適合と過平滑化の両方を緩和します。本記事では実装を簡潔にするため $\alpha$ の初期残差を主役にします。
class GCNIILayer(nn.Module):
"""初期残差付きGCN層: (1-α)Â H + α H0 を線形変換"""
def __init__(self, dim, alpha=0.1):
super().__init__()
self.alpha = alpha
self.lin = nn.Linear(dim, dim)
def forward(self, H, H0, A_hat):
# 隣接集約と初期残差の凸結合
H_agg = (1 - self.alpha) * (A_hat @ H) + self.alpha * H0
return F.relu(self.lin(H_agg))
各層で H0(最初に入力を線形射影したもの)を渡し続ける点がポイントです。$\alpha\bm{H}^{(0)}$ の項があることで、たとえ100層積んでも表現が入力情報を保持し、定常状態への収束が原理的に防がれます。3つ目の対策は、層ごとの出力を「捨てずに使う」発想です。
Jumping Knowledge(JK-Net)
JK-Net(Xu et al., 2018)は、最終層だけを使うのをやめ、全層の出力を集めて統合します。浅い層は局所的(過平滑化していない)情報を、深い層は広域の情報を持つので、ノードごとに適切な「受容野」を選べるようにする狙いです。$L$層の出力 $\bm{H}^{(1)}, \dots, \bm{H}^{(L)}$ を集約関数 $\mathrm{AGG}$ でまとめます。
$$ \begin{equation} \bm{H}^{\text{JK}} = \mathrm{AGG}\big(\bm{H}^{(1)}, \bm{H}^{(2)}, \dots, \bm{H}^{(L)}\big) \end{equation} $$
集約には連結(concat)、要素ごと最大(max)、LSTM注意などが使われます。concat は全層を残して下流で重みを学習させ、max は各次元で「最も活性化した層」を選ぶことで、過平滑化していない浅い層の鋭い特徴を救い出します。深い層が潰れていても、浅い層の出力が保存されているため最終表現が壊れない、というのが過平滑化への効き目です。
class JumpingKnowledge(nn.Module):
"""全層出力を concat または max で統合"""
def __init__(self, mode='max'):
super().__init__()
assert mode in ('max', 'cat')
self.mode = mode
def forward(self, layer_outputs):
# layer_outputs: 各層の (N, dim) のリスト
stacked = torch.stack(layer_outputs, dim=0) # (L, N, dim)
if self.mode == 'max':
return stacked.max(dim=0).values # (N, dim)
else: # cat
return torch.cat(layer_outputs, dim=1) # (N, L*dim)
max モードは次元を増やさず、各特徴次元について全層中の最大値を取るので、浅い層の強い信号が最終表現に確実に反映されます。cat モードは情報を一切捨てない代わりに次元が層数倍になります。これら3つの対策をそろえたので、いよいよ深層GCNで精度を比較し、どれだけ過平滑化を防げるかを定量化します。
深層GCNの精度比較実験
最後に、Cora 風のノード分類タスクで、(A)素のGCN、(B)PairNorm入り、(C)初期残差(GCNII風)、(D)Jumping Knowledge を、層数を変えながら比較します。PyTorch Geometric が入っていれば本物の Cora を、なければ空手クラブグラフで簡易検証します。まずモデルを定義します。
class DeepGCN(nn.Module):
"""素のGCN / PairNorm / 初期残差 / JK を切り替えられる深層GCN"""
def __init__(self, in_dim, hid, out_dim, n_layers, mode='vanilla'):
super().__init__()
self.mode = mode
self.n_layers = n_layers
self.in_lin = nn.Linear(in_dim, hid)
self.out_lin = nn.Linear(hid, out_dim)
if mode == 'gcnii':
self.layers = nn.ModuleList([GCNIILayer(hid) for _ in range(n_layers)])
else:
self.layers = nn.ModuleList([nn.Linear(hid, hid) for _ in range(n_layers)])
self.pairnorm = PairNorm() if mode == 'pairnorm' else None
self.jk = JumpingKnowledge('max') if mode == 'jk' else None
def forward(self, X, A_hat):
H0 = F.relu(self.in_lin(X))
H = H0
outs = []
for l in range(self.n_layers):
if self.mode == 'gcnii':
H = self.layers[l](H, H0, A_hat)
else:
H = F.relu(self.layers[l](A_hat @ H)) # 素のGCN伝播
if self.pairnorm is not None:
H = self.pairnorm(H)
outs.append(H)
if self.jk is not None:
H = self.jk(outs)
return self.out_lin(H)
このクラスは mode で4手法を切り替えます。vanilla は何も対策しない素のGCN、pairnorm は各層後にPairNorm、gcnii は初期残差層、jk は全層出力をmax統合します。共通の入出力射影を持たせ、層数 n_layers だけを変えて公平に比較できる設計です。次に学習ループを回します。
def make_data():
"""Cora があれば使い、なければ空手クラブ。特徴・ラベル・分割を返す"""
try:
from torch_geometric.datasets import Planetoid
ds = Planetoid(root='/tmp/Cora', name='Cora')[0]
A = np.zeros((ds.num_nodes, ds.num_nodes))
ei = ds.edge_index.numpy()
A[ei[0], ei[1]] = 1.0
Ah = torch.tensor(normalize_adj(A), dtype=torch.float32)
return (ds.x, ds.y, Ah, ds.train_mask, ds.test_mask, int(ds.y.max())+1)
except Exception:
# 空手クラブ: 特徴は単位行列、ラベルはコミュニティ、半数を訓練に
Xk = torch.eye(N)
y = torch.tensor(labels, dtype=torch.long)
mask = torch.zeros(N, dtype=torch.bool); mask[::2] = True
return (Xk, y, A_hat_t, mask, ~mask, 2)
Xd, yd, Ad, tr, te, C = make_data()
def train_eval(mode, n_layers, epochs=200):
torch.manual_seed(0)
model = DeepGCN(Xd.shape[1], 32, C, n_layers, mode)
opt = torch.optim.Adam(model.parameters(), lr=0.01, weight_decay=5e-4)
for _ in range(epochs):
model.train(); opt.zero_grad()
out = model(Xd, Ad)
loss = F.cross_entropy(out[tr], yd[tr])
loss.backward(); opt.step()
model.eval()
with torch.no_grad():
pred = model(Xd, Ad).argmax(dim=1)
acc = (pred[te] == yd[te]).float().mean().item()
return acc
make_data は環境に応じてデータを切り替え、train_eval は指定の手法・層数で学習し、テスト精度を返します。最適化には Adam と weight decay を使い、すべての手法で同条件にしています。これで層数を振った比較表を作ります。
layer_list = [2, 4, 8, 16]
modes = ['vanilla', 'pairnorm', 'gcnii', 'jk']
results = {m: [] for m in modes}
for m in modes:
for L in layer_list:
acc = train_eval(m, L)
results[m].append(acc)
print(f"{m:9s} L={L:2d} test acc = {acc:.3f}")
# 可視化: 層数に対する精度曲線
plt.figure(figsize=(8, 5))
for m in modes:
plt.plot(layer_list, results[m], 'o-', label=m)
plt.xlabel('Number of layers')
plt.ylabel('Test accuracy')
plt.title('Depth vs accuracy: over-smoothing mitigation')
plt.xticks(layer_list)
plt.legend(); plt.grid(True, alpha=0.3)
plt.savefig('oversmoothing_accuracy.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフが本記事の結論を一枚に凝縮しています。素のGCN(vanilla)は2〜4層では高精度ですが、8層・16層では精度が大きく低下するはずです。これは過平滑化で全ノードが似た表現に潰れ、分類できなくなったためで、まさに理論で導いたエネルギーの指数減衰の帰結です。一方、PairNorm・GCNII・JKは層を深くしても精度の落ち込みが小さく、16層でも実用的な精度を保ちます。とりわけ初期残差(GCNII)は入力情報を毎層注入するため深さに対して頑健で、JKは浅い層の鋭い特徴を救い出すことで深層化の悪影響を吸収します。深さと精度のトレードオフが、対策によって明確に緩和されることが定量的に確認できました。
これらの数値はデータセットや乱数で多少変動しますが、「素のGCNは深くすると劣化し、過平滑化対策はそれを防ぐ」という定性的な傾向は頑健に再現されます。理論(エネルギー減衰)と実験(精度低下)、そして対策(エネルギー保持)が一本の線でつながったことになります。
まとめ
本記事では、深層GNNを悩ませる過平滑化(over-smoothing)を、直感・数理・実装の3面から徹底的に解剖しました。
- 直感: GCNの1層は「隣接ノードの平均化」であり、グラフ上の拡散にほかならない。混ぜすぎるとインクのように全ノードが均一な表現へ潰れる。
- 収束の数理: 正規化隣接行列 $\hat{\bm{A}}$ の固有値は $(-1,1]$ にあり最大はちょうど1。固有値分解 $\hat{\bm{A}}^l=\sum_i\lambda_i^l\bm{u}_i\bm{u}_i^\top$ から、$l\to\infty$ で表現は定常モード $\bm{u}_1$(次数で決まる方向)に潰れる。
- エネルギーの指数減衰: ディリクレエネルギー $E(\bm{H})=\mathrm{tr}(\bm{H}^\top\hat{\bm{L}}\bm{H})$ は固有モード分解で $\sum_i\mu_i\|\bm{c}_i\|^2$ となり、1層ごとに $s=\max_{i\ge2}(1-\mu_i)^2<1$ 倍以下に縮む。よって $E(\bm{H}^{(l)})\le s^l E(\bm{H}^{(0)})$ で指数的に0へ。
- 可視化: 空手クラブグラフで、エネルギーの片対数プロットが直線(指数減衰)になり、類似度行列が深い層で全面赤(全ノード同一化)になることを確認した。
- 対策と比較: PairNorm(エネルギー再正規化)、初期残差/GCNII(入力の毎層注入)、Jumping Knowledge(全層統合)を実装し、深層化しても精度が維持されることを定量比較した。
過平滑化を理解すると、「なぜGNNはCNNのように深くできないのか」「どうすれば深くできるのか」という実務上の核心が見通せます。学習のロードマップとしては、ここから (1)注意機構で集約を適応的にする方向、(2)ラプラシアン固有値そのものを操るスペクトル的アプローチ、(3)残差・正規化による深層化、の3つに枝分かれしていきます。
次のステップとして、以下の記事も参考にしてください。
- グラフアテンションネットワーク(GAT)を理解する — 隣接の重みを学習して適応的に集約する手法
- GCN(グラフ畳み込みネットワーク)を理解して実装する — 本記事の出発点となる伝播則の導出
- グラフラプラシアンとは — 固有値・スペクトルの性質を深掘り