グラフ信号処理に基づいたグラフ学習(2)

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

本記事は、[グラフ信号処理に基づいたグラフ学習(1)](https://qiita.com/matyaki/items/f1df0b8f3456de79e28b)の続きです。

https://qiita.com/matyaki/items/f1df0b8f3456de79e28b

本記事では、**時変グラフ学習(Time-Varying Graph Learning; TVGL)** について扱います。

脳の機能的ネットワークや人間関係のような、時間変化するグラフの構造を推定するために使用可能な手法です。（教師データなしで）

# グラフ学習
## 時変グラフ学習

本記事で扱う時変グラフ学習は[[Kalofolias, et al. 2017]](https://ieeexplore.ieee.org/abstract/document/7952672)[^1]と[[Yamada, et al. 2019]](https://ieeexplore.ieee.org/abstract/document/8682762)[^2]を参考にしています。

### 問題設定

時変グラフ学習では、各時刻$t = 1, \ldots, T$に観測された$K$個のグラフ信号からなるデータ行列$\boldsymbol{X}_t\in\mathbb{R}^{N\times K}$を入力として受け取り、次の特性**P1**-**P3**を満たすようなグラフの隣接行列の列$(\boldsymbol{W}_1, \ldots, \boldsymbol{W}_T)$を推定する問題について考えます：
- **P1**: 各時刻$t$について、観測されたグラフ信号は推定されるグラフ上で滑らか
- **P2**: 各時刻$t$について、推定されるグラフは疎
- **P3**: 短時間ではグラフの構造（辺の有無とその重み）はあまり変化しない

:warning: 本記事では、オフラインでグラフの構造の推定を行います。

### 定式化

この問題は、ハイパーパラメータを$\alpha > 0, \beta\geq 0, \eta\geq 0$として、次の最適化問題に定式化されます（静的グラフ学習のときと同じく、$\alpha$は辺の重みのスケールを表すパラメータであることを示せます。）：
```math
\text{minimize}_{(\boldsymbol{W}_1, \ldots, \boldsymbol{W}_T)\in\mathcal{W}^T}\qquad\sum_{t=1}^T(\frac{1}{2}\|\boldsymbol{W}_t\circ\boldsymbol{Z}_t\|_{1,1} - \alpha\boldsymbol{1}_N^\top\log(\boldsymbol{W}_t\boldsymbol{1}_N) + \frac{\beta}{2}\|\boldsymbol{W}_t\|_\text{F}^2) + \eta\sum_{t=2}^T\|\boldsymbol{W}_t - \boldsymbol{W}_{t-1}\|_{\text{norm}}\tag{3}
```
ここで、$\boldsymbol{Z}_t$は時刻$t$におけるrow-wise distances matrixです。
(3)式中の最後の項について、[[Kalofolias, et al. 2017]](https://ieeexplore.ieee.org/abstract/document/7952672)[^1]では、normをフロベニウスで、[[Yamada, et al. 2019]](https://ieeexplore.ieee.org/abstract/document/8682762)[^2]では、normを1,1ノルムで定式化しています。$\boldsymbol{W}_t$

### 最適化

この問題を式変形することで、主双対近接分離法で解けるようにしていきます。（クソダル計算）

まず、静的グラフ学習のときと同じように、 $\boldsymbol{W}_t$の対角成分より右上の成分をまとめたベクトルを$\boldsymbol{w}_t$とします。また、$\boldsymbol{z}_t$も同様に定義します。
さらに、$\boldsymbol{w}, \hat{\boldsymbol{w}}, \boldsymbol{z}\in\mathbb{R}^{N(N-1)T/2}$を次のように定義します：
```math
\boldsymbol{w}=\begin{bmatrix}\boldsymbol{w}_{1}^\top&\boldsymbol{w}_{2}^\top&\boldsymbol{w}_{3}^\top&\cdots&\boldsymbol{w}_{T}^\top\end{bmatrix}^\top
```
```math
\hat{\boldsymbol{w}}=\begin{bmatrix}\boldsymbol{w}_{1}^\top&\boldsymbol{w}_{1}^\top&\boldsymbol{w}_{1}^\top&\cdots&\boldsymbol{w}_{T-1}^\top\end{bmatrix}^\top
```
```math
\boldsymbol{z}=\begin{bmatrix}\boldsymbol{z}_{1}^\top&\boldsymbol{z}_{2}^\top&\boldsymbol{z}_{3}^\top&\cdots&\boldsymbol{z}_{T}^\top\end{bmatrix}^\top
```
さらに、$\boldsymbol{v}_{1,t} = \boldsymbol{W}_t\boldsymbol{1}_N\in\mathbb{R}^N$として、
```math
\boldsymbol{v}_1=\begin{bmatrix}\boldsymbol{v}_{1,1}^\top&\boldsymbol{w}_{1,2}^\top&\boldsymbol{w}_{1,3}^\top&\cdots&\boldsymbol{w}_{1,T}^\top\end{bmatrix}^\top
```
```math
\boldsymbol{v}_{2} = \boldsymbol{w} - \hat{\boldsymbol{w}}
```
とします。
加えて、行列$\boldsymbol{S},\boldsymbol{P}$をそれぞれ次式を満たす行列とします：
```math
\boldsymbol{S}\boldsymbol{w} = \boldsymbol{v}_1
```
```math
\boldsymbol{P}\boldsymbol{w} = \boldsymbol{v}_2
```
さらにさらに、行列$\boldsymbol{A}$を次式で定義します：
```math
\boldsymbol{A}=\begin{bmatrix}\boldsymbol{S}\\\boldsymbol{P}\end{bmatrix}
```
このとき、最適化問題(3)は次のように変換されます：
```math
\text{minimize}_{\boldsymbol{w}\geq\boldsymbol{0}}\qquad\boldsymbol{w}^\top\boldsymbol{z}-\alpha\boldsymbol{1}_{NT}^\top\boldsymbol{log}(\boldsymbol{S}\boldsymbol{w}) + \beta\|\boldsymbol{w}\|_2^2 + \eta \|\boldsymbol{P}\boldsymbol{w}\|_{\text{norm}}
```

さらに、指示関数を用いて、制約なし最適化問題に変換します：

```math
\text{minimize}_{\boldsymbol{w}\in\mathbb{R}^{N(N-1)T/2}}\qquad\boldsymbol{w}^\top\boldsymbol{z} + \mathbb{I}(\boldsymbol{w}) - \alpha\boldsymbol{1}_{NT}^\top\boldsymbol{log}(\boldsymbol{S}\boldsymbol{w}) + \beta\|\boldsymbol{w}\|_2^2 + \eta \|\boldsymbol{P}\boldsymbol{w}\|_{\text{norm}}
```

これで、問題(1)の形の最適化問題に落とし込むことができました。
- $f(\boldsymbol{w}) = \boldsymbol{w}^\top\boldsymbol{z} + \mathbb{I}(\boldsymbol{w})$
- $g(\boldsymbol{w}) = - \alpha\boldsymbol{1}_{NT}^\top\boldsymbol{log}(\boldsymbol{S}\boldsymbol{w}) + \eta ||\boldsymbol{P}\boldsymbol{w}||$
- $h(\boldsymbol{w}) = \beta ||\boldsymbol{w}|| _2^2\quad(\text{with }\xi = 2\beta)$
とすればよいです。(関数$g$は行列$\boldsymbol{A}$を用いて表現できています。)

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

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

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

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

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1011247/0e082cba-2217-26e0-e8df-c424c083dd47.png)

# 最後に

本記事では、時変グラフ学習について扱いました。

私は、観測時刻が不等間隔であるような場合に使用可能な手法を提案しています。（ハイパーパラメータの数を増やさずに）[[Takahashi, et al. 2023]](https://www.ieice.org/publications/proceedings/summary.php?expandable=14&iconf=NOLTA&session_num=B2L-3&number=B2L-31&year=2023)[^3]

こちらの手法は、ハイパーパラメータ$\eta$を観測時刻の差$\Delta t$によって時間減衰するようにしています。


[^1]: Kalofolias, Vassilis, et al. "Learning time varying graphs." 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). Ieee, 2017.
[^2]: Yamada, Koki, Yuichi Tanaka, and Antonio Ortega. "Time-varying graph learning based on sparseness of temporal variation." ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2019.
[^3]: Masaki, Takahashi, and Iwasaki Satoru. "Dynamic Graph Learning to Nonuniformly Sampled Sequential Data." IEICE Proceedings Series 76.B2L-31 (2023).