こんにちは、matyakiです。
大学院で数理最適化に関する研究をしています。

本記事では、「グラフ信号処理に基づいたグラフ学習」について扱います。

グラフ信号処理とは？

まず、 グラフ信号処理(Graph Signal Processing) とは、グラフ（ネットワーク）上の信号を扱うための手法であり、従来の音声信号のような１次元信号や画像のような２次元信号の処理手法を一般化したものです。具体的には、グラフ信号処理では、１次元信号をパスグラフ上の信号として、周期的な１次元信号をサイクルグラフ上の信号として、画像信号を２次元格子グラフ上の信号として表現できます。また、グラフ信号処理では、交通網のような規則的でない構造を持つネットワーク上で定義される信号を扱うことができます。さらに、グラフ信号処理ではグラフ上のデータのフーリエ変換（グラフフーリエ変換）や、フィルタリングなどの操作が可能です。

最近グラフ信号処理に関する和書（グラフ信号処理の基礎と応用）¹が出たので、興味のある方は読んでみてください。

グラフ学習の必要性

交通網や電力網などのグラフ構造は実体があるため、見ることができます。一方で、人間や企業間の関係や脳の機能的ネットワークは実体を伴っていません。そこで、グラフ上の各頂点で観測されたデータを用いて、そのデータの背景にあるグラフの構造を推定するための手法であるグラフ学習 が必要とされてきました。

グラフ学習は以下のような応用が考えられます：

脳の機能的ネットワークなどのグラフ構造を有する現象の理解
グラフ上のデータのノイズ除去（グラフ学習で得られたグラフを基にローパスフィルタリング）
頂点のクラスタリング（グラフ学習で得られたグラフを用いて頂点のクラスタリング）

本記事では、グラフ信号処理に基づいたグラフ学習 について扱いますが、同様の問題設定で、統計的にグラフの構造を推定する手法もあります。（両手法を比較・サーベイした論文[Mateos, et al. 2019]）²

準備

ここでは、グラフ学習の説明に必要な用語の定義を行います。

グラフとグラフ信号

本記事では重み付き無向グラフ$\mathcal{G} = (\mathcal{V}, \mathcal{E}, \boldsymbol{W})$について扱います。ここで、$\mathcal{V}, \mathcal{E}, \boldsymbol{W}$はそれぞれ頂点集合、辺集合、重み付き隣接行列です。また、頂点数$N = |\mathcal{V}|$とします。

また、グラフ信号はグラフ$\mathcal{G}=(\mathcal{V}, \mathcal{E}, \boldsymbol{W})$に対し、頂点集合$\mathcal{V}$を定義域とする写像$x\colon\mathcal{V}\to\mathbb{R}$であり、頂点集合が$\mathcal{V} = \lbrace i\mid i=1,\ldots,N\rbrace$であるとき、$x(i)$を$i$次元目の要素とする$N$次元ベクトル$\boldsymbol{x}=\begin{bmatrix}x(1)&\cdots&x(N)\end{bmatrix}\in\mathbb{R}^N$として表現できます。
また、グラフ$\mathcal{G}$上のグラフ信号$\boldsymbol{x}$の滑らかさは次式で定義されます。（全変動と呼ばれることもあります。）

  \frac{1}{2}\sum_{(i,j)\in\mathcal{E}}W_{ij}|x(i) - x(j)|^2

ここで、$W_{ij}$は隣接行列$\boldsymbol{W}$の$i,j$成分であり、この値が小さいほど、$\boldsymbol{x}$は$\mathcal{G}$上で滑らかであることを意味します。（言葉的には気持ち悪いですが、一般的に使われている用語です。）
グラフ信号の滑らかさは式変形を行うことで、

\frac{1}{2}\sum_{(i,j)\in\mathcal{E}}W_{ij}|x(i) - x(j)|^2 = \boldsymbol{x}^\top\boldsymbol{L}\boldsymbol{x}

となることが分かります。
ここで、$\boldsymbol{D}\in\mathbb{R}^{N\times N}$を$\mathcal{G}$の次数行列として、$\boldsymbol{L}= \boldsymbol{D} - \boldsymbol{W}$は$\mathcal{G}$のグラフラプラシアン行列です。

また、$K$個のグラフ信号をまとめたデータ行列を$\boldsymbol{X} = \begin{bmatrix} \boldsymbol{x}^1&\cdots&\boldsymbol{x}^K \end{bmatrix}\in\mathbb{R}^{N\times K}$とすると、

\sum_{k = 1}^K(\boldsymbol{x}^k)^\top\boldsymbol{L}\boldsymbol{x}^k = \text{tr}(\boldsymbol{X}^\top\boldsymbol{L}\boldsymbol{X}) = \frac{1}{2}\|\boldsymbol{W}\circ\boldsymbol{Z}\|_{1,1}

が成り立ちます。
ここで、行列$\boldsymbol{Z}\in\mathbb{R}^{N\times N}$はrow-wise distances matrixと呼ばれ、各要素は次式で与えられます：

  Z_{ij} = \sum_{k=1}^K|X_{ik} - X_{jk}|^2

また、行列の1,1ノルムは次式で与えられます：

\|\boldsymbol{A}\|_{1,1} = \sum_{i,j}|A_{ij}|

主双対近接分離法

次に、本記事で使用する連続最適化アルゴリズムである主双対近接分離法(Primal Dual Splitting Method) について触れます。

以下では、$\Gamma_0(\mathbb{R}^n)$を$\mathbb{R}^n$上のすべての閉真凸関数の集合とします。

主双対近接分離法では次式で表される最適化問題について適用可能です：

\text{minimize}_{\boldsymbol{x}\in\mathbb{R}^{n}}\qquad f(\boldsymbol{x}) + g(\boldsymbol{A}\boldsymbol{x}) + h(\boldsymbol{x})
\tag{1}

ここで、$\boldsymbol{A}\in\mathbb{R}^{k\times n},f\in\Gamma_0(\mathbb{R}^n), g\in\Gamma_0(\mathbb{R}^k), h\in\Gamma_0(\mathbb{R}^n)$であり，$h$は$\xi$-平滑関数です。

主双対近接分離法の1つの手法がForward-Backward-Forward-based(FBF-based)主双対近接分離法で、次のようなアルゴリズムです：

ここで、任意の$\gamma > 0$に対して、関数$f\in\Gamma_0(\mathbb{R}^n)$の近接写像$\boldsymbol{\text{prox}}_{\gamma f}\colon\mathbb{R}^n\to\mathbb{R}^n$は次式で定義されます：

\boldsymbol{\text{prox}}_{\gamma f}(\boldsymbol{x})= \text{argmin}_{\boldsymbol{y}\in\mathbb{R}^n}\left(\gamma f(\boldsymbol{y}) + \frac{1}{2}\|\boldsymbol{x} - \boldsymbol{y}\|_2^2\right)

このアルゴリズムの収束性について、次の定理が知られています([Komodakis, et al. 2015])³。

以下の３つの仮定

最適化問題(1)が解を持つ。
$\mu = \xi + |\boldsymbol{A}|_S$とする。このとき、ある$\epsilon\in(0,1/(1 + \mu))$が存在して、$i = 0,1,\ldots$について、$\gamma_i\in [\epsilon,(1-\epsilon)/\mu]$が成り立つ。ここで、$|\boldsymbol{A}|_S$は行列$\boldsymbol{A}$のスペクトルノルムである。
$\text{int}(\text{dom}(g))\cap\boldsymbol{A}(\text{dom}(f))\neq\emptyset$ または $\text{dom}(g)\cap\text{int}(\boldsymbol{A}(\text{dom}(f)))\neq\emptyset$が成り立つ。ここで、$\text{int}(\mathcal{S})$は集合$\mathcal{S}$の内部であり、$\text{dom}(f)$は関数$f$の実効定義域である。

が全て成り立つとき、Algorithm 1 によって生成される$\boldsymbol{x}^i$の列は、問題(1)の大域的最適解に収束する。

グラフ学習

ここでは、メインのテーマであるグラフ信号処理に基づいたグラフ学習について扱います。
グラフ学習は、大きく2つに分類されます：

静的グラフ学習(Static Graph Learning; SGL): 時間変化しないグラフの構造を推定する
時変グラフ学習(Time-Varying Graph Learning; TVGL): 時間変化するグラフの構造を推定する

本記事では静的グラフ学習について扱います。
（時変グラフ学習については別記事で。。。）

静的グラフ学習

本記事で扱う静的グラフ学習は[Kalofolias, 2016]⁴を参考にしています。

問題設定

静的グラフ学習では、$K$個の観測されたグラフ信号からなるデータ行列$\boldsymbol{X} = \begin{bmatrix} \boldsymbol{x}^1&\cdots&\boldsymbol{x}^K \end{bmatrix}\in\mathbb{R}^{N\times K}$を入力として受け取り、次の特性P1とP2を満たすようなグラフの隣接行列$\boldsymbol{W}\in\mathbb{R}^{N\times N}$を推定する問題について考えます：

P1: 観測されたグラフ信号は推定されるグラフ上で滑らか
P2: 推定されるグラフは疎

（気持ちとしては、類似度が高い頂点対には重みの大きい辺があるが、全体としては疎なグラフが得られます。）

定式化

この問題は、ハイパーパラメータを$\alpha > 0, \beta\geq 0$として、次の最適化問題で定式化されます：

\text{minimize}_{\boldsymbol{W}\in\mathcal{W}}\qquad\frac{1}{2}\|\boldsymbol{W}\circ\boldsymbol{Z}\|_{1,1} - \alpha\boldsymbol{1}_N^\top\log(\boldsymbol{W}\boldsymbol{1}_N) + \frac{\beta}{2}\|\boldsymbol{W}\|_\text{F}^2\tag{2}

ここで、$\mathcal{W}$はありうる重み付き単純無向グラフの隣接行列の集合であり、次式で定義されます：

\mathcal{W}= \left\{\boldsymbol{W}\in\mathbb{R}_{\geq 0}^{N\times N}\;\middle|\;\boldsymbol{W} = \boldsymbol{W}^\top, \text{diag}(\boldsymbol{W}) = \boldsymbol{0}\right\}

最適化問題(2)の目的関数の気持ちは、次の通りです：

第一項：滑らかさ
第二項と第三項：第一項しかないときに得られる自明な解$\boldsymbol{W} = \boldsymbol{O}$を避けつつ、得られるグラフの疎性を制御する（$\boldsymbol{W}\boldsymbol{1}_N$は各頂点の次数を集めたベクトルである）

ハイパーパラメータについて

この問題について次のような面白い命題が成立します[Kalofolias, 2016, A.2]⁴：

$F(\boldsymbol{Z}, \alpha, \beta)$を問題(2)の最適解とする。
このとき、次式が成立する：

F(\boldsymbol{Z}, \alpha, \beta) = \alpha F(\boldsymbol{Z}, 1, \alpha\beta)

この式は、$\alpha$はグラフの辺全体のスケールを変更させるパラメータであることを意味しています。（辺の重みの定数倍を考えない；辺の重みの大小・グラフの形に興味があるときには、$\alpha = 1$と固定してよくて、$\beta$のみをチューニングすれば良い、ということ。すなわち、この問題の本質的なハイパーパラメータの数は1である。）

最適化

最適化問題(2)を主双対近接分離法で解ける形に落とし込みます。

まず、決定変数である隣接行列$\boldsymbol{W}\in\mathcal{W}$について、実質的な決定変数は$\boldsymbol{W}$の対角成分より上（右）の成分のみ（$\because$対角成分は0かつ対称行列）なので、それらの成分を集めたベクトル$\boldsymbol{w}\in\mathbb{R}_{\geq 0}^{N(N-1)/2}$を決定変数とします。このとき、問題は次のように書き換えることができます：

\text{minimize}_{\boldsymbol{w}\in\mathbb{R}_{\geq 0}^{N(N-1)/2}}\qquad\boldsymbol{w}^\top\boldsymbol{z} - \alpha\boldsymbol{1}_{N}^\top\boldsymbol{log}(\boldsymbol{A}\boldsymbol{w})+\beta\|\boldsymbol{w}\|_2^2

ここで、$\boldsymbol{z}\in\mathbb{R}_{\geq 0}^{N(N-1)/2}$は$\boldsymbol{w}$と同様に$\boldsymbol{Z}$の対角成分より上の成分を集めたベクトルで、$\boldsymbol{A}\in\mathbb{R}^{N(N-1)/2 \times N}$は$\boldsymbol{A}\boldsymbol{w} = \boldsymbol{W}\boldsymbol{1}_N$を満たす行列です。

さらに、非負制約を目的関数に取り込むために、次の指示関数$\mathbb{I}$を導入します：

  \mathbb{I}(\boldsymbol{w}) =
  \begin{cases}
    0 & (\boldsymbol{w}\geq\boldsymbol{0})\\
    \infty & (\text{otherwise})
  \end{cases}
  \label{indicator}

このとき、問題(2)はさらに次のように制約なし最適化問題に書き換えることができます：

\text{minimize}_{\boldsymbol{w}\in\mathbb{R}^{N(N-1)/2}}\qquad\boldsymbol{w}^\top\boldsymbol{z} + \mathbb{I}(\boldsymbol{w}) - \alpha\boldsymbol{1}_{N}^\top\boldsymbol{log}(\boldsymbol{A}\boldsymbol{w})+\beta\|\boldsymbol{w}\|_2^2

これで晴れて(?)問題(1)の形の最適化問題に落とし込むことができました。

$f(\boldsymbol{w}) = \boldsymbol{w}^\top\boldsymbol{z} + \mathbb{I}(\boldsymbol{w})$
$g(\boldsymbol{w}) = - \alpha\boldsymbol{1}_{N}^\top\boldsymbol{log}(\boldsymbol{A}\boldsymbol{w})$
$h(\boldsymbol{w}) = \beta ||\boldsymbol{w}|| _2^2\quad(\text{with }\xi = 2\beta)$
とすればよいです。

あとは、Algorithm 1 を用いて、最適化問題を解くことができます。
ここで、反復法の1イテレーションあたりの計算量は$O(N^2)$です。

ソースコードと数値実験結果

静的グラフ学習と時変グラフ学習のPythonプログラムをgithub上で公開しています。
scikit-learnのようなAPIで使用することができます。
↓
https://github.com/matyaki-matyaki/PyGraphLearning

↓は上記のgithubのexampleコードの実行結果です。
パラメータのチューニングはしていないですが、ある程度Ground-truthに近い推定結果が得られています。
教師なし学習ですので、ある程度の誤差は許容ということで。。。

最後に

本記事では、グラフ信号処理に基づいたグラフ学習について扱いました。
この手法は近年、脳の分野に適用されており、良い性能を示しているようです。

かなり汎用性の高い手法ですので、他にも様々な応用が考えられると思います。

田中雄一,田中聡久：グラフ信号処理の基礎と応用,2023 ↩
Mateos, Gonzalo, et al. "Connecting the dots: Identifying network structure via graph signal processing." IEEE Signal Processing Magazine 36.3 (2019): 16-43. ↩
Komodakis, Nikos, and Jean-Christophe Pesquet. "Playing with duality: An overview of recent primal? dual approaches for solving large-scale optimization problems." IEEE Signal Processing Magazine 32.6 (2015): 31-54. ↩
Kalofolias, Vassilis. "How to learn a graph from smooth signals." Artificial intelligence and statistics. PMLR, 2016. ↩ ↩²
Ménoret, Mathilde, et al. "Evaluating graph signal processing for neuroimaging through classification and dimensionality reduction." 2017 IEEE Global Conference on Signal and Information Processing (GlobalSIP). IEEE, 2017. ↩
Gao, Siyuan, et al. "Smooth graph learning for functional connectivity estimation." NeuroImage 239 (2021): 118289. ↩