グラフ信号処理におけるローパスフィルタ

グラフでは、頂点と辺を用いて複雑な関係性を表現する。 ローパスフィルタは、低周波成分を強調し、高周波成分を抑制するフィルタである。 グラフにおける周波数成分は、固有値分解や固有ベクトルを利用して定義される。 本記事は、グラフ信号処理におけるローパスフィルタの理論と、その実装手法を解説する。

ソースコード

github

  • jupyter notebook形式のファイルはこちら

google colaboratory

  • google colaboratory で実行する場合はこちら

実行環境

OSはmacOSである。LinuxやUnixのコマンドとはオプションが異なりますので注意していただきたい。

!sw_vers
ProductName:		macOS
ProductVersion:		13.5.1
BuildVersion:		22G90
!python -V
Python 3.9.17

基本的なライブラリをインポートし watermark を利用してそのバージョンを確認しておきます。 ついでに乱数のseedの設定をします。

%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import random

import scipy
import numpy as np

import matplotlib
import matplotlib.pyplot as plt

seed = 123
random_state = 123

random.seed(seed)
np.random.seed(seed)


from watermark import watermark

print(watermark(python=True, watermark=True, iversions=True, globals_=globals()))
Python implementation: CPython
Python version       : 3.9.17
IPython version      : 8.17.2

scipy     : 1.11.2
numpy     : 1.25.2
matplotlib: 3.8.1

Watermark: 2.4.3

1. グラフ信号処理(GSP)の基礎

1.1 GSPの概要

Graph Signal Processingは、グラフ上に定義される信号を対象とする。 頂点を集合$V$、辺を集合$E$として、$G=(V,E)$というグラフを考える。 頂点数を$N$とし、各頂点に対して実数信号を与えると、信号ベクトルは$\mathbf{x} \in \mathbb{R}^N$となる。 この信号をどのように変換し、処理するかを考えるのがGSPの枠組みである。 グラフ構造を考慮するために、隣接行列やラプラシアン行列を活用することが多い。 ネットワーク分析やレコメンド、SNS解析など、多様な応用分野を持つ。 例えば、企業では顧客間の関係をグラフ化し、需要予測や販促戦略に利用している。

1.2 隣接行列とラプラシアン

隣接行列$\mathbf{A} \in \mathbb{R}^{N \times N}$は、辺の有無や重みを表す行列である。 無向グラフならば、$\mathbf{A}$は対称行列となる。 ラプラシアン行列$\mathbf{L}$は、$\mathbf{L} = \mathbf{D} - \mathbf{A}$で定義される。 ただし$\mathbf{D}$は次数行列であり、対角成分に頂点の次数を配置する。 ラプラシアンの固有値と固有ベクトルは、グラフの周波数成分を定義する基盤となる。 ラプラシアンは半正定値行列であり、最小固有値は0である。 これらの行列を用いて、グラフ上の信号をフーリエ変換する手法が確立されている。

1.3 グラフ上の信号の概念

頂点が$N$個あれば、信号ベクトル$\mathbf{x}$は$N$次元となる。 たとえば、各頂点にユーザ、辺に友人関係を設定したSNSを想定する。 ユーザごとの属性や行動量が信号となる。 この信号を空間的(グラフ構造)に処理することで、より洗練された結果を導く。 GSPでは、周波数領域でフィルタをかけたり、平滑化処理を行ったりする。 ローパスフィルタは、その中でも特に低周波成分を強調する手法である。


2. ローパスフィルタの理論的背景

2.1 フィルタとは何か

信号処理におけるフィルタは、特定の周波数帯域を通過させる装置や演算の総称である。 ローパスフィルタは低周波帯域を通し、高周波帯域を遮断または弱める。 時間領域の信号処理では、滑らかな変化を維持し、急峻な変化を抑制する役割がある。 グラフにおいても、辺を跨いで類似度が高い頂点間の信号を平滑化する。 つまり、グラフ上での突然の変化(高周波成分)を抑え、なだらかな構造を強調する。

2.2 グラフフーリエ変換とローパスフィルタ

グラフフーリエ変換(GFT)は、ラプラシアン行列$\mathbf{L}$の固有分解を基盤とする。 固有分解を考えると、$\mathbf{L} = \mathbf{U} \boldsymbol{\Lambda} \mathbf{U}^\top$となる。 ここで$\mathbf{U}$は固有ベクトルを列ベクトルに並べた行列、$\boldsymbol{\Lambda}$は固有値の対角行列である。 グラフ信号$\mathbf{x}$のフーリエ変換は、$\hat{\mathbf{x}} = \mathbf{U}^\top \mathbf{x}$で定義される。 この周波数領域でのローパスフィルタは、低い固有値の成分を通し、高い固有値の成分を抑制する。 たとえば、周波数応答関数$h(\lambda)$を用いて、フィルタ処理結果を$\mathbf{y} = \mathbf{U}h(\boldsymbol{\Lambda})\mathbf{U}^\top \mathbf{x}$と書ける。

2.3 周波数応答関数の例

理想的なローパスフィルタは、臨界固有値$\lambda_c$以下を1、それ以上を0とする。 しかし、実運用ではギブズ現象などの問題があり、連続的に減衰させる設計が多い。 たとえば、$$h(\lambda) = e^{-\alpha \lambda}$$という指数関数型の設計もある。 あるいは、$$h(\lambda) = \frac{1}{1 + \beta \lambda}$$という形も用いられる。 このように、特定の帯域を強調または抑制する関数を選択する。 適切な選択は、グラフ構造や信号の特性に依存する。


3. グラフフーリエ変換と固有値分解の関係

3.1 ラプラシアンの固有値の意味

ラプラシアンの固有値は、グラフの「滑らかさ」に関わる指標である。 低い固有値に対応する固有ベクトルは、グラフ上で緩やかな変化を示す。 高い固有値に対応する固有ベクトルは、局所的で急激な変化を示す。 そのため、ローパスフィルタは低固有値の成分を優先的に残す。 スムースな部分を重視し、ノイズや突発的な変化(高周波)を抑える役割である。

3.2 固有ベクトルの直交性

無向グラフのラプラシアンは、対称かつ半正定値である。 よって、固有ベクトルたち$\mathbf{u}_1,\mathbf{u}_2,\dots,\mathbf{u}_N$は直交基底を形成する。 これらを用いて信号空間を完全に展開できる。 フーリエ変換とは、この固有ベクトル基底への変換と対応づけられる。 よって、周波数ドメインでフィルタ処理をすると、空間ドメインでの平滑化が実現できる。

3.3 高次のラプラシアンによる近似

ローパスフィルタを実装する際、全固有ベクトルを求めるのは計算コストが大きい。 大規模グラフでは、固有分解自体が困難なことがある。 そこで、ラプラシアン行列を用いた多項式近似や行列指数演算でローパスを実現する手法がある。 たとえば、$h(\mathbf{L})$をテイラー展開やチェビシェフ多項式で近似する。 このように、直接対角化せずとも、ローカルな演算でフィルタ処理を行う方法が研究されている。


4. Pythonを用いた基本的な実装例

ここでは、Pythonによる実装例を提示する。

4.1 グラフ構造の準備

無向グラフを隣接行列として定義する例を示す。

import numpy as np

N = 5

# 隣接行列の定義(小さな例)
A = np.array([[0, 1, 1, 0, 0], [1, 0, 1, 1, 0], [1, 1, 0, 0, 0], [0, 1, 0, 0, 1], [0, 0, 0, 1, 0]], dtype=float)

# 次数行列D
D = np.diag(A.sum(axis=1))

# ラプラシアン
L = D - A
L
array([[ 2., -1., -1.,  0.,  0.],
       [-1.,  3., -1., -1.,  0.],
       [-1., -1.,  2.,  0.,  0.],
       [ 0., -1.,  0.,  2., -1.],
       [ 0.,  0.,  0., -1.,  1.]])

この例では5頂点のグラフを定義している。 辺のあるところを1とし、ないところを0としている。 ラプラシアン行列は、$ \mathbf{D} - \mathbf{A} $で計算できる。

4.2 固有値分解とローパスフィルタ

ラプラシアンを固有分解し、低い固有値の成分を強調する。

# 固有分解
eigvals, eigvecs = np.linalg.eigh(L)

# ローパスの閾値
lambda_c = 1.5


# 周波数応答関数
def lowpass_func(lam):
    return 1.0 if lam < lambda_c else 0.0


# 対角行列でフィルタを定義
H_diagonal = [lowpass_func(val) for val in eigvals]
H = np.diag(H_diagonal)

# 任意の信号x
x = np.array([1, 0, 2, 3, 1], dtype=float)

# フーリエ変換
xhat = eigvecs.T @ x

# フィルタ適用
yhat = H @ xhat

# 逆変換
y = eigvecs @ yhat

y.round(2)
array([1.21, 1.31, 1.21, 1.55, 1.72])

上記のコードにより、理想的なローパスフィルタを実装できる。 低い固有値成分は残し、高い固有値成分はゼロにしている。

4.3 連続的な減衰関数の例

ギブズ現象を回避したい場合、連続的な減衰でフィルタを定義する。

alpha = 0.5


def exp_lowpass(lam):
    return np.exp(-alpha * lam)


H_diagonal = [exp_lowpass(val) for val in eigvals]
H = np.diag(H_diagonal)

yhat = H @ xhat
y = eigvecs @ yhat

yhat.round(2)
array([-3.13, -0.35,  0.27, -0.16, -0.24])

指数関数型のローパスフィルタにより、よりなめらかな出力が得られる。 このように、実装例のコードは意外とシンプルな構造である。


5. 具体的な計算例

5.1 小さなグラフでの固有値の計算

先ほどのコード例で示した5頂点のグラフを改めて見直す。 隣接行列$ \mathbf{A} $は以下のように定義されていた。

$$ \mathbf{A} = \begin{pmatrix} 0 & 1 & 1 & 0 & 0 \\ 1 & 0 & 1 & 1 & 0 \\ 1 & 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 & 0 \end{pmatrix} $$

次数行列$ \mathbf{D} $は頂点ごとの次数を対角に並べた行列である。 各頂点の次数は、$[2, 3, 2, 2, 1]$である。 よって、

$$ \mathbf{D} = \begin{pmatrix} 2 & 0 & 0 & 0 & 0 \\ 0 & 3 & 0 & 0 & 0 \\ 0 & 0 & 2 & 0 & 0 \\ 0 & 0 & 0 & 2 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{pmatrix} $$

ラプラシアン$ \mathbf{L} = \mathbf{D} - \mathbf{A} $は、

$$ \mathbf{L} = \begin{pmatrix} 2 & -1 & -1 & 0 & 0 \\ -1 & 3 & -1 & -1 & 0 \\ -1 & -1 & 2 & 0 & 0 \\ 0 & -1 & 0 & 2 & -1 \\ 0 & 0 & 0 & -1 & 1 \end{pmatrix} $$

この行列を数値計算ソフトで固有分解すると、固有値は約$[0, 1, 2, 3, 4]$などになる(実際は厳密値とは若干の誤差がある)。 このうち、0に対応する固有ベクトルは定数ベクトル(全て1)に近い方向である。 高い固有値ほど、頂点間の信号差が大きい成分を表す。 ローパスフィルタでは、低い固有値(0や1、2)を優先的に通すことになる。

7.2 信号ベクトルへの適用

たとえば、$\mathbf{x}=[1,0,2,3,1]^\top$を対象とする。 固有ベクトルを$\mathbf{u}_1,\dots,\mathbf{u}_5$とすると、 $$ \hat{\mathbf{x}} = \begin{pmatrix} \mathbf{u}_1^\top \mathbf{x} \\ \mathbf{u}_2^\top \mathbf{x} \\ \mathbf{u}_3^\top \mathbf{x} \\ \mathbf{u}_4^\top \mathbf{x} \\ \mathbf{u}_5^\top \mathbf{x} \end{pmatrix} $$ が周波数領域の表現である。 理想的なローパスなら、たとえば固有値が2以下の成分を通過させる。 すると、固有値が3以上の成分をゼロにして逆変換すれば、平滑化された$\mathbf{y}$が得られる。 実際の数値例では、スパースなグラフなので大きくは変化しない可能性がある。 ただし、エッジ数が多くなると変化が顕著になる。


8. よくある疑問点と対策

8.1 カットオフ値の設定

カットオフ値$\lambda_c$の設定は、利用目的と信号特性に左右される。 グラフの規模、エッジ密度、ノイズ特性を考慮する必要がある。 経験的に適当な値を試すか、最適化問題としてカットオフ値を決定する場合もある。

8.2 フィルタ関数の選択

ステップ状の理想フィルタは、しばしばギブズ現象を起こす。 そのため、指数型や逆数型などのスムーズな関数が好まれることが多い。 応用によって最適な関数形を選び、過剰な平滑化を避ける工夫も必要となる。

8.3 大規模グラフへの適用

頂点数が数百万規模に膨れ上がると、固有分解は非現実的となる。 多項式近似(チェビシェフ多項式など)や行列の分割手法が実用的である。 局所演算を繰り返すモンテカルロ法なども提案されている。 企業レベルの巨大ネットワーク解析では、これらの近似手法が必須となる。


結論

この記事では、グラフ理論におけるローパスフィルタの基本と実装例を概説した。 グラフ構造に基づく信号処理は、ノイズ除去やレコメンドの安定化などに有効である。 一方で、固有分解の計算コストやフィルタ設計の問題などのデメリットもある。

実際には、グラフが大規模化した場合の近似手法やパラメータ選択が鍵となる。 Pythonによる実装例を示したが、ライブラリを駆使して最適化する手法も検討した方がよい。

あくまでも個人的な備忘録であるため、詳細な理論や実装については専門書や論文を参照されたい。

参考文献

  • グラフ信号処理の基礎と応用 - ネットワーク上データのフーリエ変換,フィルタリング,学習 -
    • 素晴らしい本です。とても勉強になりました。