はじめに

グラフ信号処理に関する日本語の書籍が昨年発売された。

グラフ信号処理の基礎と応用: ネットワーク上データのフーリエ変換,フィルタリング,学習 (次世代信号情報処理シリーズ 5)

作者:田中雄一
コロナ社

Amazon

本記事ではその中で解説されているグラフ信号のサンプリングと部分空間情報を利用した復元について簡単にまとめた上で、実際に試てみた際のコードと結果を紹介する。

グラフ信号処理の諸概念

グラフ信号

グラフ信号は下図のようにグラフの各頂点上に値を持つ信号である。

このような頂点上に値を持つグラフの例としては、空間上に配置された複数のセンサーが挙げられる。これは、近くにあるセンサー同士が辺でつなげば、その計測値はグラフ信号とみなせる。それ以外にも、路線図と各駅の人口、SNSのつながりと各ユーザの特性（年齢などの何らかの数値）等々、グラフ信号としてモデル化できる現実の事象は様々存在する。また、（離散）時系列信号や画像も時刻、画素を頂点とし近傍を辺でつなげばある種のグラフとみなせるので、従来の信号処理で扱っていた対象もグラフ信号処理で扱うことができる。

このようなグラフ信号に対して、従来の信号処理のように周波数解析を行ったりフィルタリングやサンプリングを行ったりするのがグラフ信号処理である。

隣接行列とグラフラプラシアン

グラフ信号処理において重要なのは、まずグラフとしてどの頂点同士がつながっているか、またつながっている場合はどれくらい強くつながっているのかという情報である。今、頂点 $i$ と $j$ が接続されておりその接続の強さ（頂点間の近さ）を $w_{ij}$ とする。そして、この $w_{ij}$ を $i$ 行 $j$ 列に持つ行列 $\mathbf{A}$ を隣接行列という（ $i, j$ 間が接続されていなければその要素は0）。

さて、今 $N$ 個の頂点のグラフ信号 $\mathbf{x} \in \mathbb{R}^N$ が与えられたとき、その信号が滑らかかどうかを測りたい。グラフ信号が滑らかということは、頂点間がつながっていてその接続が強いほど、その信号同士は近い値を持つと考えられる。そのため、

$\displaystyle \Delta_{\mathbf{x}} = \sum_{(i, j) \in E} w_{ij} (x_{i} - x_{j})^2 \tag{1}$

が小さいほど、グラフは滑らかだといえる。ここで、 $d_i = \sum_{j=1}^N w_{ij}$ とし、 $\mathbf{D}$ を各 $d_i$ を対角成分に持つ対角行列とする。すると式（1）は

$\begin{aligned} \sum_{(i, j) \in E} w_{ij} (x_{i} - x_{j})^2 &= \frac{1}{2} \sum_{i = 1}^N \sum_{j = 1}^N w_{ij} (x_{i} - x_{j})^2 \\ &= \frac{1}{2} \left(\sum_{i = 1}^N x_{i}^2 \sum_{j = 1}^N w_{ij} - 2\sum_{i = 1}^N \sum_{j = 1}^N w_{ij} x_i x_j + \sum_{j = 1}^N x_{j}^2 \sum_{i = 1}^N w_{ij} \right) \\ &=\frac{1}{2} \left(2 \sum_{i = 1}^N x_{i}^2 d_i - 2\sum_{i = 1}^N \sum_{j = 1}^N w_{ij} x_i x_j \right) \\ &= \mathbf{x}^T \mathbf{D} \mathbf{x} - \mathbf{x}^T \mathbf{A} \mathbf{x} \\ &= \mathbf{x}^T \mathbf{L} \mathbf{x} \end{aligned} \tag{2}$

という2次形式で表すことができる。ここで $\mathbf{L} = \mathbf{D} - \mathbf{A}$ と置いており、この $\mathbf{L}$ は（組み合わせ）グラフラプラシアンと呼ばれ、グラフ信号処理において重要な役割を果たす。

グラフフーリエ変換とグラフ周波数

グラフに信号においても、通常の時系列信号と同様に各周波数の成分がどれくらいの強度であるかという情報をグラフフーリエ変換によって得ることができる。

今、 $\mathbf{L}$ の固有値を $\lambda_1, \cdots, \lambda_N \ \ (\lambda_i \leq \lambda_{i + 1})$ とし、それに対応する固有ベクトルを並べた行列を $\mathbf{U} = \left[ \mathbf{u}_1, \cdots, \mathbf{u}_N \right]$ とする。そうすると式（2）は、

$\begin{aligned} \mathbf{x}^T \mathbf{L} \mathbf{x} &= \mathbf{x}^T \mathbf{U} {\rm diag}(\lambda_1, \cdots, \lambda_N )\mathbf{U}^T \mathbf{x} \\ &= \mathbf{\hat{x}}^T {\rm diag}(\lambda_1, \cdots, \lambda_N )\mathbf{\hat{x}} \\ &= \sum_{i = 1}^N \lambda_i \hat{x}^2 \left[ i \right] \end{aligned} \tag{3}$

となり、グラフラプラシアン $\mathbf{L}$ の固有値の重み付き線形和でグラフ信号の滑らかさが表現される。そのため、信号の変動が大きい（高周波）であれば、大きい $\lambda_i$ に対する係数 $\hat{x}\left[i \right]$ が大きくなる。従って、この固有値 $\lambda_1, \cdots, \lambda_N$ はそれぞれグラフにおける周波数（グラフ周波数）とみなすことができ、その係数 $\mathbf{\hat{x}}$ は対応する周波数の強度、つまりはグラフフーリエ変換で得られるグラフフーリエ係数である。改めてグラフフーリエ変換および逆グラフフーリエ変換の操作を以下に記す。

$\begin{align} \mathbf{\hat{x}} &= \mathbf{U}^T \mathbf{x} \\ \mathbf{x} &= \mathbf{U} \mathbf{\hat{x}} \end{align} \tag{4}$

この操作はグラフラプラシアンの固有ベクトル $\mathbf{u}_1, \cdots, \mathbf{u}_N$ を正規直交基底*1とする空間（グラフ周波数領域）への基底変換となっている。通常の時系列信号におけるフーリエ変換も周波数領域の正規直交基底に変換する操作だったことを考えると、式（4）のグラフフーリエ変換の定義は自然だといえる。

サンプリングと復元

信号の部分空間

グラフ信号が頂点数 $N$ よりも低い次元の部分空間に存在するとする。例えば帯域制限されたグラフ信号は、グラフフーリエ変換すると一部の周波数に対応するフーリエ係数が0になる。これはつまり $N$ よりも少ない基底で表現できる=低次元の部分空間にグラフ信号が存在すると言える*2。式で表すと、この信号が存在する部分空間（ $K$ 次元とする）の正規直交基底を列ベクトルとして並べた行列を $\mathcal{A} \in \mathbb{R}^{N \times K}$ とすると、グラフ信号 $\mathbf{x}$ は係数 $\mathbf{d} \in \mathbb{R}^k$ を使って

$\mathbf{x} = \mathcal{A} \mathbf{d} \tag{5}$

となる。

このように低次元部分空間にグラフ信号が存在する場合、その部分空間が既知で、かつ、後述の条件を満たすようにサンプリングすると元のグラフ信号を完全に復元可能になる。

サンプリング

グラフ信号 $\mathbf{x}$ から $M$ 個の頂点をサンプリングする操作は行列 $\mathbf{S} \in \mathbb{R}^{N \times M}$ を用いて、

$\mathbf{c} = \mathbf{S}^T \mathbf{x} \tag{6}$

と表せられる。これは、例えば各列をサンプリングする頂点の要素を1、他を0とした行列 $\mathbf{I}_{\mathcal{TV}}$ を $\mathbf{S}$ とすればよい。こうすると式（6）は単に部分頂点集合を抜き出す操作になる*3。

復元

$\mathbf{c}$ から元の信号 $\mathbf{x}$ を復元したい。この際、部分空間の基底 $\mathbf{\mathcal{A}}$ は既知だとする。 $\mathbf{\mathcal{A}}$ が既知なので $\mathbf{c}$ を適当な変換を施すことにより、式（5）の $\mathbf{d}$ にできれば $\mathbf{x}$ を完全に復元できる。つまり、この変換を行う行列を $\mathbf{H}$ とすると $\mathbf{\tilde{x}} = \mathbf{\mathcal{A}} \mathbf{H} \mathbf{c}$ として復元する。この式を変形すると、

$\begin{align} \mathbf{\tilde{x}} &= \mathbf{\mathcal{A}} \mathbf{H} \mathbf{c} \\ &= \mathbf{\mathcal{A}} \mathbf{H} \mathbf{S}^T \mathbf{x} \\ &= \mathbf{\mathcal{A}} \mathbf{H} (\mathbf{S}^T \mathbf{\mathcal{A}}) \mathbf{d} \end{align} \tag{7}$

なので、 $\mathbf{H} = (\mathbf{S}^T \mathbf{\mathcal{A}})^{-1}$ とすれば $\mathbf{\tilde{x}} = \mathcal{A}\mathbf{d} = \mathbf{x}$ と完全復元できる。従って、式（6）のサンプリング操作で得た信号から元の信号を復元できる条件は $\mathbf{S}^T \mathbf{\mathcal{A}}$ が逆行列を持つことである。

一方で $\mathbf{S}^T \mathbf{\mathcal{A}}$ の逆行列が存在しない場合は、 $\mathbf{S}^T \mathbf{\tilde{x}} = \mathbf{c}$ かつ $\mathbf{\tilde{x}} \in \mathcal{A}$ を満たす復元信号 $\mathbf{\tilde{x}}$ が一意に決まらないため、完全な復元はできない。この時、 $\mathbf{H}$ を $\mathbf{S}^T \mathbf{\mathcal{A}}$ の疑似逆行列 $(\mathbf{S}^T \mathbf{\mathcal{A}})^{+}$ として復元すれば、条件を満たす中で $\mathbf{\tilde{x}}^T \mathbf{\tilde{x}}$ が最小のものになる。

疑似逆行列で条件を満たす最小ノルム解が得られる理由

一応このことを示しておく。本筋から逸れるのでこの項は読み飛ばしてもらっても支障はない。

$\begin{align} &\mathbf{\tilde{x}} = \arg \min_{\mathbf{z}} \mathbf{z}^T \mathbf{z} \\ &{\rm s.t.} \ \ \mathbf{S}^T \mathbf{z} = \mathbf{c}, \ \ \mathbf{z} \in \mathcal{A} \end{align} \tag{8}$

の解が $\mathbf{\tilde{x}} = \mathbf{\mathcal{A}} (\mathbf{S}^T \mathbf{\mathcal{A}})^{+} \mathbf{S}^T \mathbf{x}$ となることを示す。まず $\mathbf{z} \in \mathcal{A}$ なので適当な係数 $\mathbf{d}$ を使って、 $\mathbf{z} = \mathcal{A} \mathbf{d}$ とできる。これを目的関数に代入すると $\mathbf{z}^T \mathbf{z} = \mathbf{d}^T \mathcal{A}^T \mathcal{A} \mathbf{d}$ となるが、 $\mathcal{A}$ が列ベクトルに正規直交基底を並べた行列であることを思い出すと $\mathcal{A}^T \mathcal{A} = \mathbf{I}$ であるため $\mathbf{z}^T \mathbf{z} = \mathbf{d}^T \mathbf{d}$ となる。さらに、 $\mathbf{c} = \mathbf{S}^T \mathbf{x}$ なので制約条件は $\mathbf{S}^T \mathcal{A} \mathbf{d} = \mathbf{S}^T \mathbf{x}$ と書くことができる。まとめると式（8）の最適化問題は、

$\begin{align} &\mathbf{\tilde{d}} = \arg \min_{\mathbf{d}} \mathbf{d}^T \mathbf{d} \\ &{\rm s.t.} \ \ \mathbf{S}^T \mathcal{A} \mathbf{d} = \mathbf{S}^T \mathbf{x}\\ &\mathbf{\tilde{x}} = \mathcal{A} \mathbf{\tilde{d}} \end{align} \tag{9}$

と変形できる。これはラグランジュの未定乗数法より、

$\mathbf{d}^T \mathbf{d} - \mathbf{\lambda}^T (\mathbf{S}^T \mathcal{A} \mathbf{d} - \mathbf{S}^T \mathbf{x}) \tag{10}$

の停留点を求めればよい。よって、

$\mathbf{\tilde{d}} = \frac{1}{2} \mathcal{A}^T \mathbf{S} \mathbf{\lambda} \tag{11}$

となり、これを制約条件式に代入して $\mathbf{\lambda}$ について解くと

$\mathbf{\lambda} = 2 (\mathbf{S}^T \mathcal{A} \mathcal{A}^T \mathbf{S})^{-1} \mathbf{S}^T \mathbf{x} \tag{12}$

が得られる。そして式（9）の第3式、式（11）、式（12）から

$\mathbf{\tilde{x}} =\mathcal{A} \mathcal{A}^T \mathbf{S} (\mathbf{S}^T \mathcal{A} \mathcal{A}^T \mathbf{S})^{-1} \mathbf{S}^T \mathbf{x} = \mathcal{A} (\mathbf{S}^T \mathcal{A})^{+} \mathbf{S}^T \mathbf{x} \tag{13}$

となるため、疑似逆行列 $(\mathbf{S}^T \mathcal{A})^{+}$ により条件を満たす最小ノルム解が得られる。

実験

ということで実際に試してみる。具体的には帯域制限されたグラフ信号に対し適当なサンプリングを行い、そこから元の信号を復元できるかをやってみた。実験用コードは以下の通り。

import numpy as np
from scipy.spatial import distance
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
import argparse

# 組み合わせグラフラプラシアンを返す
def compute_graph_laplacian(A):
    D = np.diag(np.sum(A, axis=1))
    return D - A

# 2次元平面上ランダムに頂点を生成し、その近傍をつないだグラフを生成
# 平面上の座標と隣接行列を返す
def generate_random_graph(N, adj_num=4):
    V = np.random.randn(N, 2)
    dist = distance.cdist(V, V, metric='euclidean')
    connect = np.zeros((N, N))
    for i in range(N):
        sorted = np.argsort(dist[i, :])[1:adj_num + 1]
        connect[i, sorted] = 1
        connect[sorted, i] = 1
    dist[connect == 0] = np.inf
    adj = np.exp(-dist) - np.eye(N)
    return V, adj

# グラフ周波数が低い方からK個までに帯域制限されたグラフ信号を生成
# グラフ信号と部分空間の基底を返す
def generate_band_limiation_signal(A, K=10):
    L = compute_graph_laplacian(A)
    _, U = np.linalg.eigh(L)
    d = np.random.randn(K)
    x_ = d @ U[:, :K].T
    return x_, U[:, :K]

# ランダムに部分集合を選択する行列を返す
def compute_random_sampling_matrix(N, M):
    tau = np.random.choice(range(N), size=M, replace=False)
    T = np.zeros((M, N))
    for i, t in enumerate(tau):
        T[i, t] = 1
    return T

# 指定されたサンプリング方法を用いてM個の点をサンプリング
# サンプリングされた信号とサンプリング作用素を返す
def vertex_sampling(x, M, samplinrg_operator_generator):
    St = samplinrg_operator_generator(x.shape[0], M)
    c = St @ x
    return c, St

# グラフの辺を描画
def draw_graph_edge(X, adj, ax):
    lines = []
    for i in range(X.shape[0]):
        for j in range(i + 1, X.shape[0]):
            if(adj[i, j] > 0):
                lines.append([X[i, :], X[j, :]])
    lc = LineCollection(lines, linewidths=0.5)
    ax.add_collection(lc)

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("N", type=int)
    parser.add_argument("K", type=int)
    parser.add_argument("M", type=int)
    parser.add_argument("seed", type=int)
    args = parser.parse_args()

    N = args.N #グラフの頂点数
    K = args.K #部分空間の次元
    M = args.M #サンプリング点数
    np.random.seed(args.seed)

    # グラフの生成
    V, adj = generate_random_graph(N)

    # 帯域制限されたグラフ信号の生成
    x, A = generate_band_limiation_signal(adj, K)

    # グラフと源信号を描画
    _, ax = plt.subplots()
    draw_graph_edge(V, adj, ax)
    plt.scatter(V[:, 0], V[:, 1], c=x, cmap='cool', vmax=1.0, vmin=-1.0)
    plt.colorbar()
    plt.title("original graph N={}".format(N))
    plt.savefig("origin_graph.png")
    plt.close()
    
    # ランダムに頂点をサンプリング
    c, St = vertex_sampling(x, M, compute_random_sampling_matrix)

    # サンプリング後のグラフ信号を描画
    _, ax = plt.subplots()
    draw_graph_edge(V, adj, ax)
    plt.scatter(V[:, 0], V[:, 1], c='white', linewidths=0.3, edgecolors='gray')
    plt.scatter(V[np.where(St==1)[1], 0], V[np.where(St==1)[1], 1],
                c=x[np.where(St==1)[1]], cmap='cool',  vmax=1.0, vmin=-1.0)
    plt.colorbar()
    plt.title("sampled graph N={} M={}".format(N, M))
    plt.savefig("sampled_graph.png")
    plt.close()

    # 補正行列の計算
    H = np.linalg.pinv(St @ A)

    # 復元
    x_hat = A @ H @ c

    # 復元後の信号と源信号の誤差計算
    err = np.linalg.norm(x_hat - x)

    # 復元後のグラフ信号を描画
    _, ax = plt.subplots()
    draw_graph_edge(V, adj, ax)
    plt.scatter(V[:, 0], V[:, 1], c=x_hat, cmap='cool', vmax=1.0, vmin=-1.0)
    plt.colorbar()
    plt.title("recovered graph M={} error={:.3e}".format(M, err))
    plt.savefig("recovered_graph.png")
    plt.close()

if __name__ == '__main__':
    main()

このコードでseed=10, N=200, K=40とし、サンプリング数Mは40, 30, 20と変えて評価した。

原信号

M=40のサンプリング後の信号と復元信号

原信号との2乗和誤差の平方根もグラフのタイトルに記述しているがほとんど0であり、ほぼ完全に復元できている。グラフの右側はほとんどサンプリングしていないのにもかかわらず復元できてるのが面白い。

M=30のサンプリング後の信号と復元信号

M < Kだとやはり誤差は0にならず原信号の復元はできていない。また原信号に比べて信号の変化が滑らかになっている。

M=20のサンプリング後の信号と復元信号

さらにサンプリング点数を減らすと誤差はより大きくなる。とはいえ、パッと見たところの印象としては原信号から大きく乖離しているというわけではなく、原信号の大まかな傾向は保てている。

終わりに

本記事ではグラフ信号からのサンプリングと部分空間情報を利用した復元方法を紹介し、実際に条件さえ整えばうまく復元できること確認した。一方で今回のように信号の部分空間が既知であることは稀であり、さらに観測される信号にはノイズが重畳することも多い。従って、今回紹介した復元を行うためには、ノイズを含む信号から部分空間を推定しなければならない。そのため、次回の記事ではいくつかの部分空間推定方法を試し、それにより元の信号を復元可能かを検討する。

なお、本記事では部分空間に関する情報により信号を復元する方法を紹介したが、冒頭に述べた書籍では信号の滑らかさに関する情報を利用して復元する方法等も紹介されているので、興味ある方は読んでみるとよいだろう。

*1: $\mathbf{L}$ は対称行列なので、異なる固有値に対応する固有ベクトル同士は直交する。

*2:信号が低次元部分空間に存在する例は帯域制限されたグラフ信号に限らない。例えば、グラフスペクトルが周期的であるPGS部分空間のグラフ信号は帯域は制限されていないが低次元の部分空間にグラフ信号が存在する。

*3:より一般には $\mathbf{x}$ にグラフフィルタ $\mathbf{G} \in \mathbb{R}^{N \times N}$ をかけて部分頂点集合を抜き出す行列 $\mathbf{I}_{\mathcal{TV}}\mathbf{G}$ を $\mathbf{S}$ とするが、今回は話を簡単にするために $\mathbf{G} = \mathbf{I}$ とした。