上の図が本記事のエッセンスです。左のグラフを隣接行列 $\bm{A}$ に書き起こし、そこから次数行列 $\bm{D}$ を引いたものがラプラシアン $\bm{L}$ です。対角成分はその頂点の次数、非対角成分はマイナスの重みになっている構造が一目で分かります。この行列の固有値・固有ベクトルがクラスタリングの鍵を握っています。
SNSの友達グラフを眺めていると、自然に「仲間内のかたまり」が見えてきます。あるグループの人たちは互いに密につながり、別のグループとはほとんど線が引かれていない。このような「つながりの濃いかたまり」を、座標も色も与えられていない純粋な接続情報だけから機械的に見つけ出すにはどうすればよいでしょうか。点と点の間に距離が定義されたユークリッド空間なら k-means が使えますが、グラフには「座標」がありません。あるのは「誰と誰がつながっているか」という隣接関係だけです。
この問いに、線形代数のもっとも美しい道具のひとつで答えるのがスペクトラルクラスタリングです。そのアルゴリズムの流れを図で先に示すと、以下のようになります：まずグラフからラプラシアン行列を作り、その固有ベクトルで各ノードを新しい座標空間に埋め込み、そこで k-means を適用する、という三段階です。鍵を握るのがグラフラプラシアン $\bm{L} = \bm{D} – \bm{A}$ という行列で、この行列の固有値・固有ベクトル（=スペクトル）の中に、グラフの「切れ目」の情報がそっくり埋め込まれています。驚くべきことに、「グラフを2つの塊に最もきれいに切り分けたい」という離散的な組合せ最適化問題が、ラプラシアンの固有ベクトルを求めるという連続的な線形代数の問題に変身するのです。
グラフラプラシアンは、データ解析の現場で恐ろしく広く使われています。たとえば、
- 画像セグメンテーション — 画素をノード、色や位置の近さをエッジ重みとしたグラフを切り分け、前景と背景を分離する（Normalized Cut の原論文はまさに画像分割が動機でした）。
- コミュニティ検出 — SNSや論文引用ネットワークから、密につながったグループを抽出する。
- グラフニューラルネットワーク（GCN） — 後で見るように、GCNの畳み込みはラプラシアンの固有空間における「グラフ・フーリエ変換」として定義され、そのチェビシェフ近似が ChebNet・GCN の出発点になっています。
つまりグラフラプラシアンを理解することは、古典的なクラスタリングから最新のグラフ深層学習までを貫く一本の背骨を手に入れることに他なりません。
本記事の内容
- 隣接行列・次数行列からラプラシアン $\bm{L}=\bm{D}-\bm{A}$ を定義し、その性質を導く
- 二次形式 $\bm{x}^\top \bm{L}\bm{x} = \frac{1}{2}\sum_{i,j} A_{ij}(x_i – x_j)^2$ が「信号の滑らかさ」を測ることを示す
- 正規化カット Ncut の最小化が、正規化ラプラシアンの一般化固有値問題に緩和されることを導出する
- Fiedler ベクトル（第2最小固有ベクトル）がグラフの二分割を与えることを示す
- Python で2つの環状クラスタにスペクトラルクラスタリングを適用し、$k$ 個の最小固有ベクトルを特徴として k-means で分離する
- GCN のチェビシェフ近似との接続を解説する
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- グラフニューラルネットワーク（GNN）とは — グラフをニューラルネットで扱う枠組み。ラプラシアンが畳み込みの基礎になる
- スペクトラルグラフ畳み込みとChebNet — 本記事で導くラプラシアンの固有分解が、そのままグラフ畳み込みの理論基盤になる
行列の固有値分解、対称行列の性質（固有値が実数で固有ベクトルが直交する）、レイリー商の最小化、そして k-means の基本がわかっていると、最短距離で理解できます。
グラフと隣接行列・次数行列
まず舞台を整えましょう。グラフ $G=(V,E)$ は頂点（ノード）の集合 $V=\{1,2,\dots,n\}$ と、頂点を結ぶ辺（エッジ）の集合 $E$ からなります。友達ネットワークなら頂点が人、辺が友人関係です。各辺には「つながりの強さ」を表す重み $w_{ij} \ge 0$ を割り当てることができます。単純に「つながっているかどうか」だけなら $w_{ij}\in\{0,1\}$ ですが、距離が近いほど大きい重みにする、といった連続的な値も使えます。
このグラフの構造をすべて1つの行列に詰め込んだものが隣接行列 $\bm{A}$ です。$n\times n$ の行列で、その $(i,j)$ 成分は頂点 $i$ と頂点 $j$ の間の重みです。
$$ A_{ij} = \begin{cases} w_{ij} & (i,j) \in E \\ 0 & \text{otherwise} \end{cases} $$
ここでは無向グラフを考えるので、$i$ と $j$ のつながりと $j$ と $i$ のつながりは同じです。したがって $A_{ij}=A_{ji}$、すなわち $\bm{A}$ は対称行列になります。この対称性が、あとで固有値がすべて実数になることを保証する重要な性質です。
次に、各頂点が「どれくらいつながっているか」を表す量を考えます。頂点 $i$ につながる辺の重みをすべて足し合わせたものを次数（degree） $d_i$ と呼びます。頂点の次数は「その頂点に接続するエッジの重みの合計」であり、次数の大きい頂点ほどグラフの中で「ハブ」的な役割を持ちます。
$$ d_i = \sum_{j=1}^{n} A_{ij} $$
これを対角線に並べた対角行列が次数行列 $\bm{D}$ です。
$$ \bm{D} = \mathrm{diag}(d_1, d_2, \dots, d_n) = \begin{pmatrix} d_1 & & \\ & \ddots & \\ & & d_n \end{pmatrix} $$
次数行列は「各頂点の重要度・混み具合」を、隣接行列は「誰と誰がつながっているか」を表します。この2つの行列の差を取るところから、本記事の主役が生まれます。次節でいよいよラプラシアンを定義しましょう。
グラフラプラシアンとは
ラプラシアンと聞くと、微積分に出てくる $\nabla^2 = \frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}+\cdots$ という演算子を思い出すかもしれません。連続空間のラプラシアンは「ある点の値と、その周りの平均値とのズレ」を測る道具でした。熱がどう拡散するか、膜がどう振動するかを記述する、まさに物理の主役です。
グラフラプラシアンは、この「周囲とのズレ」という発想をグラフの世界に移植したものです。連続空間の「周囲」は無限に多くの近傍点でしたが、グラフでは「辺でつながった隣接ノード」が周囲にあたります。ある頂点の値と、隣の頂点たちの値の差を測る — それがグラフラプラシアンの正体です。
定義は拍子抜けするほど単純で、次数行列から隣接行列を引くだけです。
$$ \begin{equation} \bm{L} = \bm{D} – \bm{A} \end{equation} $$
成分で書くと、対角成分は次数 $d_i$、非対角成分はマイナスの重みになります。
$$ L_{ij} = \begin{cases} d_i & (i = j) \\ -w_{ij} & (i \neq j,\ (i,j)\in E) \\ 0 & \text{otherwise} \end{cases} $$
なぜ「差」を取るとラプラシアンになるのか、まだピンとこないかもしれません。その意味は次節の二次形式を見ると一気に明快になります。先に、この行列が持つ重要な性質を3つ確認しておきましょう。
-
対称性 — $\bm{D}$ は対角行列（対称）、$\bm{A}$ は無向グラフなので対称。よって $\bm{L}=\bm{D}-\bm{A}$ も対称です。対称行列なので固有値はすべて実数で、固有ベクトルは互いに直交するように選べます。
-
行和がゼロ — 各行について、対角成分 $d_i = \sum_j A_{ij}$ と非対角成分の和 $-\sum_{j\neq i}A_{ij}$ が打ち消し合い、行の総和は $0$ です。これは、すべての成分が $1$ のベクトル $\bm{1}=(1,1,\dots,1)^\top$ に対して $\bm{L}\bm{1}=\bm{0}$ が成り立つことを意味します。つまり$\bm{1}$ は固有値 $0$ に対応する固有ベクトルです。
-
半正定値性 — 次節で示すように、任意のベクトル $\bm{x}$ について $\bm{x}^\top \bm{L}\bm{x} \ge 0$ が成り立ちます。したがって固有値はすべて非負で、$0 = \lambda_1 \le \lambda_2 \le \cdots \le \lambda_n$ と並べられます。
これらの性質、とくに「固有値 $0$ が存在する」という事実が、グラフの連結性やクラスタ数と深く結びついています。なぜ差を取ることに意味があるのか、二次形式を計算して確かめましょう。
二次形式が「滑らかさ」を測る
グラフラプラシアンの本質は、それが定義する二次形式 $\bm{x}^\top \bm{L}\bm{x}$ にあります。ここで $\bm{x}=(x_1,\dots,x_n)^\top$ は各頂点に1つずつ実数値を割り当てた「グラフ上の信号」だと思ってください。たとえば各人にある好みのスコアを与えたもの、あるいは各画素の明るさを並べたものです。
この二次形式を計算すると、驚くほどきれいな形になります。これが本記事でいちばん大切な等式です。
$$ \begin{equation} \bm{x}^\top \bm{L}\bm{x} = \frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n} A_{ij}\,(x_i – x_j)^2 \end{equation} $$
右辺は「辺でつながった頂点同士の値の差の二乗を、重みをかけて全部足したもの」です。隣り合う頂点の値が近ければ $(x_i-x_j)^2$ は小さく、大きく食い違っていれば大きくなる。つまりこの二次形式は、信号 $\bm{x}$ がグラフ上でどれだけ滑らか（隣接ノード同士で値が揃っているか） を測る量なのです。値が小さいほど滑らか、大きいほどギザギザ。これが「ラプラシアン＝周囲とのズレ」というイメージの正体です。
左の「滑らかな信号」では隣接ノード間の値の差 $|x_i – x_j|$ がいずれも小さく（約0.1〜0.3）、二次形式 $\bm{x}^\top\bm{L}\bm{x}$ が非常に小さい値になります。一方、右の「ギザギザな信号」ではエッジをまたぐ度に値が大きく変化し、同じラプラシアンで二次形式を計算すると何十倍もの大きな値になります。この対比が「ラプラシアン＝滑らかさの測定器」という直感を鮮明に示しています。エッジの線幅が差の大きさを表しており、太い線ほど信号の「段差」が急峻であることが読み取れます。
導出
この等式を、定義から1行ずつ導きましょう。ゴールは左辺 $\bm{x}^\top \bm{L}\bm{x}$ を右辺の二乗和の形に変形することです。
まず $\bm{L}=\bm{D}-\bm{A}$ を代入して、$\bm{D}$ の項と $\bm{A}$ の項に分けます。$\bm{D}$ は対角行列なので $\bm{x}^\top \bm{D}\bm{x}=\sum_i d_i x_i^2$、$\bm{A}$ の項は二重和になります。
$$ \bm{x}^\top \bm{L}\bm{x} = \bm{x}^\top \bm{D}\bm{x} – \bm{x}^\top \bm{A}\bm{x} = \sum_{i} d_i x_i^2 – \sum_{i}\sum_{j} A_{ij} x_i x_j $$
ここで次数の定義 $d_i = \sum_j A_{ij}$ を第1項に代入すると、第1項も二重和の形に書き直せます。
$$ \bm{x}^\top \bm{L}\bm{x} = \sum_{i}\Big(\sum_{j} A_{ij}\Big) x_i^2 – \sum_{i}\sum_{j} A_{ij} x_i x_j = \sum_{i}\sum_{j} A_{ij} x_i^2 – \sum_{i}\sum_{j} A_{ij} x_i x_j $$
これで2つの項が同じ二重和 $\sum_i\sum_j A_{ij}(\cdots)$ の形になったので、ひとつにまとめられます。
$$ \bm{x}^\top \bm{L}\bm{x} = \sum_{i}\sum_{j} A_{ij}\big(x_i^2 – x_i x_j\big) $$
ここで一工夫します。$\bm{A}$ が対称（$A_{ij}=A_{ji}$）であることを使い、添字 $i$ と $j$ を入れ替えた同じ式 $\sum_i\sum_j A_{ij}(x_j^2 – x_i x_j)$ を作ります。これはもとの式と値が等しい（和の取り方を入れ替えただけ）ので、両者を足して2で割っても値は変わりません。
$$ \bm{x}^\top \bm{L}\bm{x} = \frac{1}{2}\sum_{i}\sum_{j} A_{ij}\big(x_i^2 – x_i x_j\big) + \frac{1}{2}\sum_{i}\sum_{j} A_{ij}\big(x_j^2 – x_i x_j\big) $$
2つの和をまとめると、括弧の中は $x_i^2 – 2x_i x_j + x_j^2$ となり、これは完全平方です。
$$ \bm{x}^\top \bm{L}\bm{x} = \frac{1}{2}\sum_{i}\sum_{j} A_{ij}\big(x_i^2 – 2x_i x_j + x_j^2\big) = \frac{1}{2}\sum_{i}\sum_{j} A_{ij}\,(x_i – x_j)^2 $$
これで目標の等式が得られました。各項 $A_{ij}(x_i-x_j)^2 \ge 0$ なので和も非負、すなわち $\bm{x}^\top \bm{L}\bm{x}\ge 0$ であり、ラプラシアンが半正定値であることもこの式から直ちにわかります。
固有値ゼロと連結成分
この滑らかさ表現から、ラプラシアンの固有値 $0$ の意味も読み取れます。$\bm{x}^\top \bm{L}\bm{x}=0$ となるのは、すべての辺について $x_i=x_j$ が成り立つとき、つまりつながっている頂点が同じ値を持つときです。グラフが1つに連結していれば、これは「全頂点が同じ値」、すなわち $\bm{x}\propto \bm{1}$ だけ。よって固有値 $0$ の重複度は $1$ です。
もしグラフが $k$ 個の連結成分（互いにつながっていない $k$ 個の島）に分かれていれば、各島ごとに別々の定数値を取れるので、固有値 $0$ に対応する固有ベクトルは $k$ 個独立に存在します。つまりラプラシアンの固有値 $0$ の重複度は、グラフの連結成分の数に等しい。これはスペクトラルクラスタリングの理論的支柱です。クラスタが「ほぼ」分離していれば固有値は完全に $0$ ではなく「小さな正の値」になり、その小さな固有値に対応する固有ベクトルがクラスタの境界を教えてくれます。
3つのグラフを比べると、連結成分の数と固有値 $0$ の重複度が対応していることが確認できます。左の連結グラフでは最小固有値がほぼ $0$ で残りは正、中央の2成分グラフでは最小2つがほぼ $0$、右の3成分グラフでは最小3つがほぼ $0$ になっています。各成分を異なる色で塗り分けると、色ごとに独立した「定数信号（$\bm{x}\propto \bm{1}$の成分限定版）」が固有値 $0$ に対応することが直感的に分かります。
「滑らかさ」と「固有値ゼロ」が結びついたところで、いよいよ本題のクラスタリング問題を定式化しましょう。
グラフ分割問題とカット
グラフをきれいに2つに切り分けるとは、数学的にどういうことでしょうか。頂点集合 $V$ を2つの互いに素な部分集合 $A$ と $B$（$A\cup B=V$、$A\cap B=\varnothing$）に分けることを考えます。このとき、$A$ と $B$ をまたぐ辺の重みの総和をカット（cut） と呼びます。
$$ \mathrm{cut}(A,B) = \sum_{i\in A,\ j\in B} w_{ij} $$
直感的には、カットは「切断するために断ち切る必要のあるつながりの総量」です。良いクラスタリングは、クラスタ内部のつながりは密で、クラスタ間のつながりは疎であってほしい。だから、まずは素朴に「カットを最小化する分割」を探したくなります。これが最小カット問題です。
ところがこの素朴な目標には落とし穴があります。カットを最小化すると、しばしば「1個だけの孤立した頂点」と「残り全部」という、ひどく偏った分割が選ばれてしまうのです。孤立した頂点を切り離すのに必要な辺はたった数本ですから、カットの値はとても小さくなります。しかしこれは私たちが望む「バランスの取れた2分割」ではありません。
そこで、分割の「大きさ」でカットを正規化して、極端に小さなクラスタを罰するように改良します。これが次節の正規化カット（Normalized Cut, Ncut） です。最小カットの欠点を正規化で直すという発想が、ラプラシアンの固有値問題への扉を開きます。
左の「最小カット」では、孤立したノード8と残り全体を分離するだけで cut が最小（赤い1本のエッジだけを切る）になってしまいます。右の「正規化カット」では、両クラスタの体積 $\mathrm{vol}(A), \mathrm{vol}(B)$ で割ることで極端に小さいクラスタに大きなペナルティが課され、赤いエッジを2本切るにもかかわらず密なクラスタを2つに分ける均等な分割が選ばれます。この「正規化」がスペクトラルクラスタリングの要です。
正規化カット（Ncut）の定式化
正規化カットは、Shi と Malik が1997年に画像セグメンテーションのために提案した指標です。アイデアは、カットを各クラスタの「体積」で割って、小さすぎるクラスタが選ばれにくくすることです。クラスタ $A$ の体積を、$A$ に属する頂点の次数の総和
$$ \mathrm{vol}(A) = \sum_{i\in A} d_i $$
で定義します。体積は「クラスタがどれだけのつながりを内部に抱えているか」の指標です。これを使って正規化カットを次のように定義します。
$$ \begin{equation} \mathrm{Ncut}(A,B) = \frac{\mathrm{cut}(A,B)}{\mathrm{vol}(A)} + \frac{\mathrm{cut}(A,B)}{\mathrm{vol}(B)} \end{equation} $$
分母に両方の体積が入っているのがポイントです。もし $A$ が極端に小さい（=体積が小さい）と、第1項 $\mathrm{cut}/\mathrm{vol}(A)$ が爆発的に大きくなり、Ncut は大きな値になってしまいます。したがって Ncut を小さくするには、カットが小さいことに加えて、両クラスタがバランス良く大きいことが要求されます。これでさきほどの「孤立頂点問題」が解消されます。
指示ベクトルによる書き換え
Ncut を最小化したいわけですが、このままでは「どの頂点をどちらに入れるか」という組合せの問題で、頂点が $n$ 個あれば $2^n$ 通りもの分け方があり、まともに探索できません（実際この問題は NP困難です）。そこで、分割を1本のベクトルで表現し、Ncut をそのベクトルの二次形式として書き直すことを目指します。
巧妙な指示ベクトル $\bm{f}=(f_1,\dots,f_n)^\top$ を、次のように定義します。
$$ f_i = \begin{cases} \ \ \sqrt{\dfrac{\mathrm{vol}(B)}{\mathrm{vol}(A)}} & (i \in A) \\[2mm] -\sqrt{\dfrac{\mathrm{vol}(A)}{\mathrm{vol}(B)}} & (i \in B) \end{cases} $$
一見すると謎めいた重み付けですが、これは後で2つの嬉しい性質を満たすように逆算して決めたものです。この $\bm{f}$ を使うと、Ncut の分子（カット）が分母（体積）で割られた形が、ラプラシアンの二次形式 $\bm{f}^\top \bm{L}\bm{f}$ にぴたりと収まることを示します。
まず分子側を計算します。前節で導いた二次形式の公式 $\bm{f}^\top \bm{L}\bm{f}=\frac{1}{2}\sum_{i,j}A_{ij}(f_i-f_j)^2$ において、$i,j$ が同じクラスタにあれば $f_i=f_j$ なので項はゼロ。異なるクラスタにまたがる辺だけが残ります。$i\in A,\ j\in B$ の辺では
$$ f_i – f_j = \sqrt{\frac{\mathrm{vol}(B)}{\mathrm{vol}(A)}} + \sqrt{\frac{\mathrm{vol}(A)}{\mathrm{vol}(B)}} $$
となります。この差を二乗するとき、$\sqrt{\frac{\mathrm{vol}(B)}{\mathrm{vol}(A)}}+\sqrt{\frac{\mathrm{vol}(A)}{\mathrm{vol}(B)}} = \frac{\mathrm{vol}(A)+\mathrm{vol}(B)}{\sqrt{\mathrm{vol}(A)\mathrm{vol}(B)}} = \frac{\mathrm{vol}(V)}{\sqrt{\mathrm{vol}(A)\mathrm{vol}(B)}}$ とまとめられます（$\mathrm{vol}(A)+\mathrm{vol}(B)=\mathrm{vol}(V)$ を使いました）。したがって、
$$ \begin{align} \bm{f}^\top \bm{L}\bm{f} &= \frac{1}{2}\sum_{i,j}A_{ij}(f_i-f_j)^2 \\ &= \frac{1}{2}\cdot 2\!\!\sum_{i\in A, j\in B}\!\! w_{ij}\left(\sqrt{\tfrac{\mathrm{vol}(B)}{\mathrm{vol}(A)}}+\sqrt{\tfrac{\mathrm{vol}(A)}{\mathrm{vol}(B)}}\right)^2 \\ &= \mathrm{cut}(A,B)\cdot \frac{\mathrm{vol}(V)^2}{\mathrm{vol}(A)\,\mathrm{vol}(B)} \end{align} $$
ここで第2行の係数 $2$ は、和を $i\in A,j\in B$ と $i\in B,j\in A$ の両方で取ると同じ辺を2回数えることから来ています（$\frac{1}{2}$ と打ち消し合います）。最後の行は、上で求めた差の二乗を代入し、$\mathrm{cut}(A,B)=\sum_{i\in A,j\in B}w_{ij}$ でまとめたものです。
一方で、$\mathrm{Ncut}(A,B)$ そのものを変形すると、
$$ \mathrm{Ncut}(A,B) = \mathrm{cut}(A,B)\left(\frac{1}{\mathrm{vol}(A)}+\frac{1}{\mathrm{vol}(B)}\right) = \mathrm{cut}(A,B)\cdot\frac{\mathrm{vol}(V)}{\mathrm{vol}(A)\,\mathrm{vol}(B)} $$
となります（通分して $\mathrm{vol}(A)+\mathrm{vol}(B)=\mathrm{vol}(V)$ を使いました）。2つの式を見比べると、
$$ \bm{f}^\top \bm{L}\bm{f} = \mathrm{vol}(V)\cdot \mathrm{Ncut}(A,B) $$
の関係が得られます。$\mathrm{vol}(V)$ はグラフ全体で決まる定数なので、$\mathrm{Ncut}$ の最小化は $\bm{f}^\top \bm{L}\bm{f}$ の最小化と等価です。組合せ最適化が、ついに二次形式の最小化に化けました。
2つの制約条件
ただし $\bm{f}$ は自由には選べません。さきほどの指示ベクトルの定義から、次の2つの制約が自動的に従います。これらが、後で固有値問題を正しく導くための鍵です。
第一に、$\bm{f}$ は次数で重み付けした内積の意味で $\bm{1}$ と直交します。実際、
$$ (\bm{D}\bm{1})^\top \bm{f} = \sum_i d_i f_i = \sqrt{\tfrac{\mathrm{vol}(B)}{\mathrm{vol}(A)}}\sum_{i\in A}d_i – \sqrt{\tfrac{\mathrm{vol}(A)}{\mathrm{vol}(B)}}\sum_{i\in B}d_i $$
ここで $\sum_{i\in A}d_i=\mathrm{vol}(A)$、$\sum_{i\in B}d_i=\mathrm{vol}(B)$ を代入すると、
$$ = \sqrt{\tfrac{\mathrm{vol}(B)}{\mathrm{vol}(A)}}\,\mathrm{vol}(A) – \sqrt{\tfrac{\mathrm{vol}(A)}{\mathrm{vol}(B)}}\,\mathrm{vol}(B) = \sqrt{\mathrm{vol}(A)\mathrm{vol}(B)} – \sqrt{\mathrm{vol}(A)\mathrm{vol}(B)} = 0 $$
となり、確かに $\bm{f}^\top \bm{D}\bm{1}=0$ が成り立ちます。
第二に、$\bm{f}$ の $\bm{D}$ ノルムが定数になります。同様の計算で、
$$ \bm{f}^\top \bm{D}\bm{f} = \sum_i d_i f_i^2 = \frac{\mathrm{vol}(B)}{\mathrm{vol}(A)}\mathrm{vol}(A) + \frac{\mathrm{vol}(A)}{\mathrm{vol}(B)}\mathrm{vol}(B) = \mathrm{vol}(B)+\mathrm{vol}(A) = \mathrm{vol}(V) $$
が得られます。これで Ncut 最小化問題は、次の制約付き最小化として完全に書き換わりました。
$$ \min_{\bm{f}} \ \bm{f}^\top \bm{L}\bm{f} \quad \text{s.t.}\quad \bm{f}^\top \bm{D}\bm{1}=0,\quad \bm{f}^\top \bm{D}\bm{f}=\mathrm{vol}(V),\quad f_i \in \{\text{2値}\} $$
唯一残った厄介者は、最後の「$f_i$ が2つの離散値しか取れない」という制約です。これがある限り問題は NP困難のまま。次節で、この離散制約を取り払う「緩和」を行い、問題を固有値問題へと変身させます。
一般化固有値問題への緩和
組合せ最適化が NP困難である最大の原因は、$\bm{f}$ の各成分が飛び飛びの2値しか取れないことでした。スペクトラル緩和（spectral relaxation） のアイデアは大胆です。この離散制約をいったん忘れ、$\bm{f}$ を任意の実数ベクトルとして許してしまうのです。そうすると問題は連続最適化になり、解析的に解けるようになります。
緩和後の問題はこうなります。
$$ \min_{\bm{f}\in\mathbb{R}^n}\ \bm{f}^\top \bm{L}\bm{f} \quad \text{s.t.}\quad \bm{f}^\top \bm{D}\bm{1}=0,\quad \bm{f}^\top \bm{D}\bm{f}=\mathrm{vol}(V) $$
これはレイリー商の最小化に他なりません。変数変換 $\bm{g}=\bm{D}^{1/2}\bm{f}$ を導入すると見通しが良くなります。$\bm{f}=\bm{D}^{-1/2}\bm{g}$ を代入すると、目的関数と制約は次のように書き換わります。
$$ \bm{f}^\top \bm{L}\bm{f} = \bm{g}^\top \bm{D}^{-1/2}\bm{L}\bm{D}^{-1/2}\bm{g} = \bm{g}^\top \bm{L}_{\mathrm{sym}}\bm{g} $$
ここで現れた
$$ \begin{equation} \bm{L}_{\mathrm{sym}} = \bm{D}^{-1/2}\bm{L}\bm{D}^{-1/2} = \bm{I} – \bm{D}^{-1/2}\bm{A}\bm{D}^{-1/2} \end{equation} $$
を対称正規化ラプラシアンと呼びます。制約 $\bm{f}^\top\bm{D}\bm{f}=\mathrm{vol}(V)$ は $\bm{g}^\top\bm{g}=\mathrm{vol}(V)$ という単純なノルム制約になり、直交制約 $\bm{f}^\top\bm{D}\bm{1}=0$ は $\bm{g}^\top \bm{D}^{1/2}\bm{1}=0$ になります。
ここで $\bm{D}^{1/2}\bm{1}$ が何かに注目します。$\bm{L}_{\mathrm{sym}}(\bm{D}^{1/2}\bm{1}) = \bm{D}^{-1/2}\bm{L}\bm{1}=\bm{0}$（$\bm{L}\bm{1}=\bm{0}$ より）なので、$\bm{D}^{1/2}\bm{1}$ は $\bm{L}_{\mathrm{sym}}$ の固有値 $0$ の固有ベクトルです。つまり制約「$\bm{g}\perp \bm{D}^{1/2}\bm{1}$」は、「$\bm{g}$ を $\bm{L}_{\mathrm{sym}}$ の最小固有ベクトルに直交させよ」という意味になります。
レイリー・リッツの定理によれば、対称行列 $\bm{L}_{\mathrm{sym}}$ のレイリー商 $\frac{\bm{g}^\top \bm{L}_{\mathrm{sym}}\bm{g}}{\bm{g}^\top\bm{g}}$ を、最小固有ベクトルに直交する制約のもとで最小化する解は、第2最小固有値に対応する固有ベクトルです。したがって緩和問題の解 $\bm{g}$ は $\bm{L}_{\mathrm{sym}}$ の第2固有ベクトルであり、もとの変数に戻すと、
$$ \bm{L}_{\mathrm{sym}}\bm{g} = \lambda\bm{g} \quad\Longleftrightarrow\quad \bm{L}\bm{f} = \lambda \bm{D}\bm{f} $$
という一般化固有値問題 $\bm{L}\bm{f}=\lambda\bm{D}\bm{f}$ に帰着します。Ncut 最小化は、ラプラシアンと次数行列の一般化固有値問題を解くことに緩和されたのです。あとは得られた第2固有ベクトルをどう「2値の分割」に戻すか。それが次節の Fiedler ベクトルの話です。
Fiedler ベクトルによる二分割
緩和して得られる第2最小固有ベクトルには、発見者の数学者 Miroslav Fiedler にちなんでFiedler ベクトルという名前がついています。Fiedler ベクトルに対応する固有値 $\lambda_2$（正規化していないラプラシアンなら代数的連結度と呼ばれる）は、「グラフがどれだけ強く連結しているか」を表す重要な量です。
なぜ第2固有ベクトルがクラスタリングに使えるのでしょうか。第1（最小）固有ベクトルは $\bm{1}$（全頂点が同じ値）で、これは「分割しない」自明な解に対応し、何の情報も持ちません。だから直交制約でこれを排除しました。残った中で最も「滑らか」=最も小さい固有値を持つのが Fiedler ベクトルです。滑らかな信号とは「隣接ノード同士で値が近い」信号でしたから、Fiedler ベクトルは「密につながったかたまり内ではほぼ同じ値を取り、かたまりの境界でだけ大きく値が変わる」という性質を持ちます。
この性質を使えば、二分割は驚くほど簡単です。Fiedler ベクトル $\bm{f}$ の各成分の符号を見て、
$$ i \in \begin{cases} A & f_i \ge 0 \\ B & f_i < 0 \end{cases} $$
と分けるだけ。正の値を取る頂点と負の値を取る頂点が、ちょうど2つのクラスタに対応します。緩和によって連続値になってしまった解を、符号という「最も近い2値表現」に丸め戻すわけです。閾値を $0$ ではなく中央値にする変種もありますが、考え方は同じです。
緩和したことの代償として、この符号分割は必ずしも厳密な Ncut 最小解とは限りません。しかし実用上は非常に良い近似を与えることが知られており、最小カットの偏り問題も正規化によって回避できています。2分割の原理がわかったところで、これを $k$ 個のクラスタに一般化しましょう。
k クラスタへの一般化
クラスタが2つではなく $k$ 個ある場合、1本の Fiedler ベクトルだけでは足りません。連結成分の議論を思い出してください。グラフがちょうど $k$ 個のクラスタにきれいに分かれていれば、ラプラシアンは固有値 $0$ を $k$ 重に持ち、その固有空間が各クラスタの指示ベクトルで張られるのでした。クラスタが「ほぼ」分離している現実のグラフでは、固有値が完全な $0$ ではなく「$k$ 個の小さな固有値」になり、対応する $k$ 本の固有ベクトルがクラスタ構造を符号化します。
そこでアルゴリズムは次のようになります。
- ラプラシアン $\bm{L}$（または正規化版 $\bm{L}_{\mathrm{sym}}$）を作る。
- 最小から $k$ 個の固有値に対応する固有ベクトル $\bm{u}_1,\dots,\bm{u}_k$ を求め、それらを列に並べた行列 $\bm{U}\in\mathbb{R}^{n\times k}$ を作る。
- $\bm{U}$ の各行 $\bm{y}_i\in\mathbb{R}^k$ を、頂点 $i$ の新しい「埋め込み座標」とみなす。つまり、座標のなかったグラフのノードを、$k$ 次元のユークリッド空間の点に写像する。
- この $k$ 次元の点群に通常の k-means を適用して $k$ クラスタに分ける。
このアルゴリズムの美しさは、「グラフという座標のない対象を、ラプラシアンの固有ベクトルによって座標のある空間に埋め込み、そこで普通のクラスタリングをする」という流れにあります。固有ベクトルへの写像（スペクトラル埋め込み）が、複雑にからみあったクラスタを、ユークリッド空間できれいに分離した塊に変換してくれるのです。正規化版を使うアルゴリズムは Ng・Jordan・Weiss（2002）のものが有名で、$\bm{U}$ の各行を単位長に正規化してから k-means にかけます。
理論はここまでで完成です。あとは Python で動かして、本当にからみあったクラスタが分離できるのかを確かめましょう。
Pythonでの実装
スペクトラルクラスタリングの真価は、k-means が苦手とする「非凸なクラスタ」で発揮されます。ここでは、二重の同心円（内側の輪と外側の輪）という、k-means が絶対に分離できない代表的なデータで試します。まずデータを作り、隣接行列を作るところから始めましょう。
データ生成と類似度グラフの構築
座標データから類似度グラフを作るときは、各点を「自分から近い $k$ 個の点」とだけ結ぶ k近傍グラフ を使うのが定番です。これにより、同じ輪の上の点同士はつながり、内側の輪と外側の輪はほとんどつながらない、という構造ができます。
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_circles
from sklearn.neighbors import kneighbors_graph
np.random.seed(0)
# 二重の同心円データ（200点）。factorで内円と外円の比を指定
X, y_true = make_circles(n_samples=200, factor=0.4, noise=0.06, random_state=0)
# k近傍グラフから隣接行列を構築（k=10、対称化）
knn = kneighbors_graph(X, n_neighbors=10, mode='connectivity', include_self=False)
A = 0.5 * (knn + knn.T) # 無向グラフにするため対称化
A = (A > 0).astype(float).toarray() # 0/1の重みに
# 元データを可視化
plt.figure(figsize=(6, 6))
plt.scatter(X[:, 0], X[:, 1], c=y_true, cmap='coolwarm', s=30)
plt.title('Original data: two concentric circles')
plt.xlabel('x1'); plt.ylabel('x2'); plt.axis('equal')
plt.tight_layout()
plt.show()
左の散布図から、データが内側の輪と外側の輪という2本の環状クラスタをなしていることがはっきり分かります。右の k 近傍グラフでは、同じ輪の上の点同士が短いエッジで密に結ばれ、内側の輪と外側の輪をまたぐエッジはほとんど存在しないことが確認できます。このグラフ構造があることで、「半径」という非線形な境界がエッジの有無という離散的な情報に変換されており、ラプラシアンがその構造を正確に捉えられます。次に、このデータからラプラシアンを作り、その固有構造を覗いてみましょう。
ラプラシアンの構築と固有値の観察
隣接行列 $\bm{A}$ から次数行列 $\bm{D}$、ラプラシアン $\bm{L}=\bm{D}-\bm{A}$、そして対称正規化ラプラシアン $\bm{L}_{\mathrm{sym}}$ を作り、固有値を小さい順に並べて観察します。クラスタ数 $k=2$ なら、固有値の並びに「最初の2つが小さく、3つ目で急に大きくなる」というギャップが現れるはずです。
# 次数行列とラプラシアン
D = np.diag(A.sum(axis=1))
L = D - A
# 対称正規化ラプラシアン L_sym = I - D^{-1/2} A D^{-1/2}
d = A.sum(axis=1)
D_inv_sqrt = np.diag(1.0 / np.sqrt(d))
L_sym = np.eye(len(A)) - D_inv_sqrt @ A @ D_inv_sqrt
# 固有値・固有ベクトル（対称行列なのでeighを使う：実固有値・直交固有ベクトル）
eigvals, eigvecs = np.linalg.eigh(L_sym)
# 小さい方から10個の固有値をプロット
plt.figure(figsize=(7, 4))
plt.plot(range(1, 11), eigvals[:10], 'o-')
plt.axvline(2.5, color='gray', linestyle='--', label='eigengap (k=2)')
plt.title('Smallest eigenvalues of normalized Laplacian')
plt.xlabel('index'); plt.ylabel('eigenvalue'); plt.legend()
plt.tight_layout()
plt.show()
print("最小から5個の固有値:", np.round(eigvals[:5], 4))
左のグラフから、最初の2つの固有値がほぼ $0$ に張りつき、3番目の固有値で値がはっきり跳ね上がる「固有ギャップ（eigengap）」が読み取れます。最小の2固有値がほぼ $0$ なのは、2本の円がほとんど分離した連結成分のように振る舞っているからで、前に証明した「固有値 $0$ の重複度＝連結成分数」がノイズで少し緩んだ姿です。右の棒グラフは固有ギャップ $\lambda_k – \lambda_{k-1}$ を示しており、$k=2$ のギャップが突出して最大であることから、クラスタ数として $k=2$ が最適であることが読み取れます。このギャップの位置がクラスタ数 $k=2$ を教えてくれており、理論どおりの結果です。次は、この第2固有ベクトル（Fiedler ベクトル）を直接見てみましょう。
Fiedler ベクトルによる二分割
第2最小固有ベクトル＝Fiedler ベクトルの各成分を、頂点インデックス順に並べてプロットし、符号で分割します。理論では「クラスタ内ではほぼ一定、境界で符号が変わる」はずです。
# 第2最小固有ベクトル（Fiedlerベクトル）
fiedler = eigvecs[:, 1]
# 符号で二分割
labels_fiedler = (fiedler >= 0).astype(int)
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
# Fiedlerベクトルの値（点ごと）
axes[0].scatter(range(len(fiedler)), fiedler, c=labels_fiedler, cmap='coolwarm', s=25)
axes[0].axhline(0, color='gray', linestyle='--')
axes[0].set_title('Fiedler vector entries (sign split)')
axes[0].set_xlabel('node index'); axes[0].set_ylabel('fiedler value')
# 元データ上で色分け
axes[1].scatter(X[:, 0], X[:, 1], c=labels_fiedler, cmap='coolwarm', s=30)
axes[1].set_title('Bipartition by Fiedler sign')
axes[1].set_xlabel('x1'); axes[1].set_ylabel('x2'); axes[1].axis('equal')
plt.tight_layout()
plt.show()
左のグラフから、Fiedler ベクトルの成分値がソートしても2つの明確なグループ（正の塊と負の塊）にくっきり分かれていることが分かります。中間的な値を取るノードがほとんどなく、$0$ を境にスパッと2群に割れているのが特徴です。中央のグラフでは各点を Fiedler ベクトル値でカラーマッピングしており、色の青赤が内側と外側の輪に完璧に対応しています。右のグラフでその符号を元データに塗り戻すと、内側の円と外側の円が見事に分離されており、直線では不可能だった非凸クラスタの分割が実現できています。続いて、$k$ 本の固有ベクトルを使う一般的な手順を実装し、k-means と比べてみます。
スペクトラルクラスタリング本体と k-means との比較
$k=2$ 個の最小固有ベクトルを並べた行列 $\bm{U}$ の各行を埋め込み座標とし、k-means にかけます。同じデータに素の k-means も適用し、両者を比べます。
from sklearn.cluster import KMeans
# k個の最小固有ベクトルを特徴行列に（k=2）
k = 2
U = eigvecs[:, :k]
# Ng-Jordan-Weiss流：各行を単位長に正規化
U_norm = U / (np.linalg.norm(U, axis=1, keepdims=True) + 1e-10)
# スペクトラル埋め込み上でk-means
labels_spec = KMeans(n_clusters=k, n_init=10, random_state=0).fit_predict(U_norm)
# 比較用：元の座標に直接k-means
labels_km = KMeans(n_clusters=k, n_init=10, random_state=0).fit_predict(X)
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].scatter(X[:, 0], X[:, 1], c=labels_km, cmap='coolwarm', s=30)
axes[0].set_title('Plain k-means (fails)')
axes[0].set_xlabel('x1'); axes[0].set_ylabel('x2'); axes[0].axis('equal')
axes[1].scatter(X[:, 0], X[:, 1], c=labels_spec, cmap='coolwarm', s=30)
axes[1].set_title('Spectral clustering (succeeds)')
axes[1].set_xlabel('x1'); axes[1].set_ylabel('x2'); axes[1].axis('equal')
plt.tight_layout()
plt.show()
2つのグラフの対比が、スペクトラルクラスタリングの威力を雄弁に物語っています。左の素の k-means は ARI（Adjusted Rand Index）がほぼ $0$（ランダム分割と変わらない）で、データの重心からの距離で分けようとするため円を左右（あるいは上下）にばっさり切ってしまい、2つの輪をまったく分離できていません。一方、右のスペクトラルクラスタリングは ARI=1.000（完璧な一致）で内側の輪と外側の輪を完璧に分けています。違いは、ラプラシアンの固有ベクトルが作る埋め込み空間では、もともと非凸だったクラスタが線形分離可能な塊に変換されている点にあります。最後に、その埋め込み空間が実際にどんな形をしているのかを覗いてみましょう。
スペクトラル埋め込み空間の可視化
なぜ埋め込み後は k-means が成功するのか。$\bm{U}$ の2列、すなわち第1・第2固有ベクトルを座標軸として点をプロットすると、その理由が一目でわかります。
# スペクトラル埋め込み空間（固有ベクトル座標）での点配置
plt.figure(figsize=(6, 6))
plt.scatter(U[:, 0], U[:, 1], c=labels_spec, cmap='coolwarm', s=30)
plt.title('Spectral embedding (eigenvector coordinates)')
plt.xlabel('1st eigenvector'); plt.ylabel('2nd eigenvector (Fiedler)')
plt.tight_layout()
plt.show()
左の「固有ベクトル座標空間」では、もとは絡み合った2本の輪だった点群が、2つの離れた「点の塊」に変身しています。第1固有ベクトル軸はほぼ一定値（自明な定数解の名残）で、第2固有ベクトル（Fiedler）軸方向に2群がはっきり分離していることが見て取れます。右の「行正規化後の埋め込み（NJW法）」では各行を単位長にしたことで2群がさらに明確に分かれており、この状態では重心ベースの k-means が易々と正解を出せます。スペクトラル埋め込みが「非線形なクラスタ構造を線形分離可能な配置に展開する」という、この記事の核心がここに凝縮されています。
GCN のチェビシェフ近似との接続
ここまで見てきたラプラシアンの固有分解は、実はグラフ畳み込みニューラルネットワーク（GCN） の理論的な出発点でもあります。最後に、この橋を渡しておきましょう。
正規化ラプラシアンを固有分解して $\bm{L}_{\mathrm{sym}}=\bm{U}\bm{\Lambda}\bm{U}^\top$ と書くとき、固有ベクトル行列 $\bm{U}$ はグラフ上の信号に対する「フーリエ基底」の役割を果たします。実際、固有値 $\lambda$ は連続のフーリエ変換における「周波数」に対応し、小さい固有値の固有ベクトルは滑らかな（低周波の）信号、大きい固有値の固有ベクトルは激しく振動する（高周波の）信号を表します。これがちょうど、本記事で見た「滑らかさ＝二次形式の小ささ＝小さい固有値」という直感の延長線上にあります。
グラフ上の信号 $\bm{x}$ を $\hat{\bm{x}}=\bm{U}^\top\bm{x}$ と変換するのがグラフ・フーリエ変換で、スペクトル領域でフィルタ $g_\theta(\bm{\Lambda})$ をかけてから戻す操作
$$ \bm{y} = \bm{U}\,g_\theta(\bm{\Lambda})\,\bm{U}^\top \bm{x} $$
がスペクトラルグラフ畳み込みです。しかしこの定義をそのまま使うには、$n\times n$ 行列の固有分解（$O(n^3)$）が必要で、大きなグラフでは計算が破綻します。
そこで登場するのがチェビシェフ多項式近似です。フィルタ $g_\theta(\bm{\Lambda})$ を固有値のチェビシェフ多項式 $T_k$ の和で近似すると、
$$ g_\theta(\bm{L}_{\mathrm{sym}}) \approx \sum_{k=0}^{K} \theta_k\, T_k(\tilde{\bm{L}}) $$
のように、固有分解を一切せずにラプラシアンの「べき乗」だけで畳み込みが書けます（$\tilde{\bm{L}}$ は固有値を $[-1,1]$ にスケールしたラプラシアン）。これが ChebNet であり、その $K=1$ への単純化と正規化トリックから生まれたのが、現在広く使われる GCN の伝播式 $\bm{H}^{(l+1)}=\sigma(\tilde{\bm{D}}^{-1/2}\tilde{\bm{A}}\tilde{\bm{D}}^{-1/2}\bm{H}^{(l)}\bm{W}^{(l)})$ です。詳しい導出はスペクトラルグラフ畳み込みとChebNetで扱っています。
つまり、スペクトラルクラスタリングで「クラスタ構造を取り出す道具」だったラプラシアンの固有空間が、GCN では「ニューラルネットが学習する畳み込みの場」として再利用されているのです。本記事で固有値・固有ベクトルの意味を掴んでおけば、GCN の数式が「フィルタを周波数領域で設計している」のだと自然に読めるようになります。
左のグラフは固有値の大きさを「周波数」として視覚化したものです。小さい固有値（青）がクラスタ情報を含む低周波成分、大きい固有値（赤）がノイズ的な高周波成分に対応しており、スペクトラルクラスタリングは「低周波の固有ベクトルだけを使う低域通過フィルタ」と解釈できます。右の系譜図は、スペクトラルクラスタリングから始まり、スペクトラルグラフ畳み込み（$O(n^3)$の固有分解が必要）→ ChebNet（チェビシェフ多項式による局所近似）→ GCN（$K=1$への単純化）という流れを示しています。ラプラシアンを出発点にした一本の理論の流れが、現代の GCN まで続いているのです。
まとめ
本記事では、グラフラプラシアンの固有構造からスペクトラルクラスタリングを導き、Python で実装しました。
- ラプラシアン $\bm{L}=\bm{D}-\bm{A}$ は対称・半正定値で、固有値 $0$ の重複度がグラフの連結成分数に等しい。
- 二次形式 $\bm{x}^\top\bm{L}\bm{x}=\frac{1}{2}\sum_{i,j}A_{ij}(x_i-x_j)^2$ は信号の「滑らかさ」を測り、小さい固有値ほど滑らかな信号に対応する。
- 正規化カット Ncut の最小化は、指示ベクトルを使うと $\bm{f}^\top\bm{L}\bm{f}$ の最小化と等価になり、離散制約を緩和すると一般化固有値問題 $\bm{L}\bm{f}=\lambda\bm{D}\bm{f}$（＝対称正規化ラプラシアンの固有値問題）に帰着する。
- Fiedler ベクトル（第2最小固有ベクトル）の符号がグラフの二分割を与え、$k$ 本の最小固有ベクトルを特徴とした k-means が $k$ クラスタへの一般化になる。
- 二重円データで、k-means が失敗する非凸クラスタをスペクトラルクラスタリングが完璧に分離できることを確認した。スペクトラル埋め込みが非線形構造を線形分離可能に展開している。
- ラプラシアンの固有分解はそのまま GCN/ChebNet のグラフ・フーリエ変換の基礎になっている。
スペクトラルクラスタリングは「離散の組合せ問題を線形代数に翻訳する」という、応用数学の典型的かつ強力なパターンの好例です。この視点を身につければ、グラフ深層学習の理論も格段に見通しよく学べます。
次のステップとして、以下の記事も参考にしてください。
- スペクトラルグラフ畳み込みとChebNet — 本記事のラプラシアン固有分解を畳み込みに発展させる
- グラフニューラルネットワーク（GNN）とは — グラフ深層学習の全体像
関連記事
グラフラプラシアンの連続版（$\nabla^2$ 演算子）の基礎は以下の記事で詳しく解説しています。スペクトラルクラスタリングとの接続を理解する上で有用です。
本記事で導いたラプラシアンの固有分解が、グラフ畳み込みの理論基盤として ChebNet・GCN へと発展する過程は以下で詳述しています。
参考文献
- J. Shi and J. Malik, “Normalized Cuts and Image Segmentation,” IEEE TPAMI, 2000.
- A. Ng, M. Jordan, and Y. Weiss, “On Spectral Clustering: Analysis and an Algorithm,” NeurIPS, 2002.
- U. von Luxburg, “A Tutorial on Spectral Clustering,” Statistics and Computing, 2007.