甲斐性なしのブログ

うどんくらいしか食べる気がしない

グラフ信号の部分空間を推定する

近接勾配法 連続最適化 スパース学習 信号処理

はじめに

本記事は数理最適化 Advent Calendar 2024 23日目の記事です。

以前の記事にてグラフ信号処理におけるサンプリングと部分空間情報を利用した原信号の復元を紹介した。

yamagensakam.hatenablog.com

しかし、この記事の最後でも述べたが、グラフ信号の部分空間が既知であることは稀であり、しかも観測される信号にはノイズが重畳することも多い。 ということで今回は、ノイズが重畳したグラフ信号からの部分空間推定についていくつか思いついた方法を実験してみたので、それらについて紹介する。

なお、実験に使用したコードは適宜記事内でも記載するが全体については以下のgithubリポジトリに公開している。

github.com

また上に貼った前回の記事の実験コードを流用している部分があるので、見慣れない関数などがあればそちらも参照。

問題設定

まずN頂点のグラフにて観測されたグラフ信号\mathbf{x} \in \mathbb{R}^Nは、K次元部分空間（K \lt N）の基底を並べた行列\mathcal{A} \in \mathbb{R}^{N \times K}と係数\mathbf{d} \in \mathbb{R}^Kを用いて以下の式に従うとする。

\mathbf{x} = \mathcal{A} \mathbf{d} + \epsilon \tag{1}

ここで\epsilonは正規分布 N( 0, \sigma^2) に従うノイズ成分である。

このグラフ信号をn回観測を行って得たグラフ信号集合\left\{\mathbf{x}_1, \cdots, \mathbf{x}_n \right\}とグラフの構造情報（隣接行列など）から低次元部分空間の基底\mathcal{A}を推定するというのが今回解きたい問題。なおこの際、\mathcal{A}は固定だが\mathbf{d}は観測ごとに異なるとする。

実験するにあたって真の基底\mathcal{A}を定める必要があるが、いくつか典型的なモデルがあるので、今回はそのうち以下に紹介する帯域制限グラフ信号、区分的一定グラフ信号、周期的グラフスペクトルグラフ信号の3種類で生成し実験を行う。

帯域制限グラフ信号

これは前回の記事でも使った信号ではあるが改めて紹介する。 グラフラプラシアン\mathbf{L}の固有ベクトルを並べた行列を\mathbf{U}=\left[\mathbf{u}_1 \cdots, \mathbf{u}_N \right] とすると、帯域制限グラフ信号は\mathbf{u}_1 \cdots, \mathbf{u}_NからK個だけ選んだ\mathcal{A}=\left[\mathbf{u}_{s_1} \cdots, \mathbf{u}_{s_K} \right] が基底となる。そのため、選択しなかった固有ベクトルに対応する固有値の周波数成分は含まない信号となる。一般的にはどの周波数成分を選んでも帯域制限グラフ信号と言えるが、ここでは低周波成分からK個だけ選んだグラフ信号を扱う。このグラフ信号生成の実装は前回の記事の実験コードを参照。

K=20としたときのグラフ信号とスペクトルおよびそれに対し標準偏差\sigma=1のノイズを加えた信号の例は以下のようになる。

信号（ノイズなし） スペクトル
信号（ノイズあり \sigma=1） スペクトル

区分的一定グラフ信号

頂点のクラスタを考え、そのクラスタ内は同じ信号値を持つとしたグラフ信号。ちゃんとした定義は参考文献[1]を見ていただくとして、ここでは生成の実装を書くにとどめる（generate_random_graphなど前回の実験コードを流用しているので、先にそちらを見ていただいた方が理解しやすいと思う）。

def generate_piecewise_constant_signal(V, n, K=10):
    N = V.shape[0]
    kmeans = KMeans(n_clusters=K).fit(V)
    A = np.zeros((N, K))
    for k in range(K):
        A[kmeans.labels_ == k, k] = 1.0
    D = np.random.randn(n, K)
    X = D @ A.T
    return X, A

# グラフの作成
V, adj = utils.generate_random_graph(N, graph_adj_num)

# 区分的一定グラフ信号の作成
# nは観測数、Kはクラスタ（部分空間次元）数
X, A = generate_piecewise_constant_signal(V, n, K)

このグラフ信号に標準偏差\sigma=2のノイズを加えた例が以下の通り。

信号（ノイズなし） スペクトル
信号（ノイズあり \sigma=2） スペクトル

周期的グラフスペクトルグラフ信号（PGSグラフ信号）

グラフ周波数のスペクトルが周期的になっているグラフ。こちらも早速生成の実装と実際のグラフ信号、スペクトルの例を示す。

def generate_pgs_signal(A, n, K=10):
    L = utils.compute_graph_laplacian(A)
    _, U = np.linalg.eigh(L)
    D = np.random.randn(n, K)
    Ik = np.tile(np.eye(K), U.shape[0]//K)
    X = D @ Ik @ U.T
    return X, U @ Ik.T

# グラフの作成
V, adj = utils.generate_random_graph(N, graph_adj_num)

# PGSグラフ信号の作成
# nは観測数、Kはスペクトルの周期（部分空間次元）
X, A = generate_pgs_signal(adj, n, K)
信号（ノイズなし） スペクトル
信号（ノイズあり \sigma=2） スペクトル

信号の部分空間の推定

主成分分析を使う方法

まず部分空間と言って思いつくのはこれ。主成分分析によって得られる基底のうち第K主成分までを採用する。実装を見た方が速いのでコードを示す。観測したグラフ信号集合の行列\mathbf{X}がn行N列と数式とは転置されていることに注意*1。

def subspace_estimate_with_pca(X, k):
    pca = PCA(k, svd_solver='full')
    pca.fit(X)
    return pca.components_.T

しかしこの方法ではグラフ構造の情報を一切使わない方法なので、ノイズに対するロバスト性が低いと考えられる*2。そこで次はグラフ情報を使った方法を考える。

グラフフーリエ変換を使う方法

帯域制限されたグラフ信号であれば、制限された周波数成分のスペクトルの大きさが0に近くなるはずである。そのため、それぞれの周波数成分の強さを以下の式のように各観測グラフ信号のスペクトルの2乗和平方根とし、s\left[i \right] が上位K件に対応するグラフフーリエの基底を採用する。

\displaystyle s[i ] = \sqrt{\sum_{j=1}^n \hat{x}_j [i ]^2} \tag{2}

実装は以下の通り。

class GraphFourierTransformer:
    def __init__(self, A):
        self.A = np.copy(A)
        self.L = compute_graph_laplacian(A)
        self.lam, self.U = np.linalg.eigh(self.L)

    def transform(self, X):
        return X @ self.U
    
    def inv_transform(self, Xf):
        return Xf @ self.U.T

def subspace_estimate_with_gft(X, adj, k):
    gft = utils.GraphFourierTransformer(adj)
    Xf = gft.transform(X)
    s = np.sqrt(np.sum((Xf * Xf), axis=0))
    remain = np.argsort(-s)[:k]
    return gft.U[:, remain]

ただこの方法では帯域制限された信号以外の部分空間を推定するのは妥当ではなく、推定精度が悪くなると予想される。

グラフ信号に関する仮定や事前知識を正則化として取り入れる（頂点領域）

主成分分析による方法は仮定が弱く、グラフフーリエ変換による方法は仮定が強すぎる。そのためグラフ信号に関する情報を取り入れた正則化項を考えることで適切な仮定を取り入れた部分空間推定を行う。

まず、グラフ信号集合\mathbf{X}をよく近似できるk次元の部分空間の基底\mathcal{\hat{A}}とその係数\mathbf{\hat{D}}を見つける問題を素直に記述すると、

\begin{aligned} & \mathcal{\hat{A}}, \mathbf{\hat{D}} = {\rm arg} \min_{\mathcal{A}, \mathbf{D}} \parallel \mathbf{X} - \mathcal{A} \mathbf{D} \parallel_F^2 \\ & \mathcal{A} \in \mathbb{R}^{N \times k}, \mathbf{D} \in \mathbb{R}^{k \times n} \end{aligned} \tag{3}

となる。この問題に対して、正則化項を加えることによりグラフの信号に関する情報を埋め込むことを考える。例えばグラフ信号が滑らかであるという仮定を最適化問題に埋め込む場合は、滑らかさの指標であるグラフラプラシアンの2次形式\mathbf{x}^T \mathbf{L}\mathbf{x}が小さくなるはずなので、全ての観測した信号に対するグラフラプラシアンの2次形式の和{\rm Tr}( \left(\mathcal{A} \mathbf{D} \right)^T \mathbf{L} \left(\mathcal{A} \mathbf{D}\right) )を正則化項として目的関数に付け加えるとよい。もちろんこの正則化項はあくまで例で、他にも行列\mathbf{L}を他の正定値対称行列にかえれば、別のグラフに関する情報も表現できる。

この一般の正則化項の行列を\mathbf{T}として、式（3）を書き換えると

\begin{aligned} & \mathcal{\hat{A}}, \mathbf{\hat{D}} = {\rm arg} \min_{\mathcal{A}, \mathbf{D}} \parallel \mathbf{X} - \mathcal{A} \mathbf{D} \parallel_F^2 + \alpha {\rm Tr}(\left(\mathcal{A}\mathbf{D}\right)^T \mathbf{T} \left(\mathcal{A} \mathbf{D}\right))\\ & \mathcal{A} \in \mathbb{R}^{N \times k}, \mathbf{D} \in \mathbb{R}^{k \times n} \end{aligned} \tag{4}

となる。さらに\mathbf{W} = \mathcal{A} \mathbf{D}として、

\begin{aligned} & \mathbf{\hat{W}} = {\rm arg} \min_{\mathbf{W}} \parallel \mathbf{X} - \mathbf{W} \parallel_F^2 + \alpha {\rm Tr}(\mathbf{W}^T \mathbf{T} \mathbf{W})\\ & {\rm rank}(\mathbf{W}) \leq k \end{aligned} \tag{5}

と書き換えることができる*3。 式（5）の最適化問題を確率的勾配法などで解く方法も考えられるが、凸でないので実験する上で少々厄介。従って今回は以下のようにランク制約付きの最適化問題をトレースノルム正則化つき最適化問題に緩和して解く。

\begin{aligned} \mathbf{\hat{W}} = {\rm arg} \min_{\mathbf{W}} \parallel \mathbf{X} - \mathbf{W} \parallel_F^2 + \alpha {\rm Tr}(\mathbf{W}^T \mathbf{T} \mathbf{W}) + \beta \parallel \mathbf{W} \parallel_{*} \end{aligned} \tag{6}

この問題は凸な最適化問題なので解きやすい。

この最適化問題を近接勾配法で解く。トレースノルムの正則化項付き最適化問題への近接勾配法適用は参考文献[3]が詳しい。また過去に近接勾配法でトレースノルム正則化項つきの最適化問題を解いた記事も書いてあるのでこちらも貼っておく。

yamagensakam.hatenablog.com

近接勾配法を適用するために、式（6）を少し変形する。

\begin{aligned} \mathbf{\hat{W}} &= {\rm arg} \min_{\mathbf{W}} \parallel \mathbf{X} - \mathbf{W} \parallel_F^2 + \alpha {\rm Tr}(\mathbf{W}^T \mathbf{T} \mathbf{W}) + \beta \parallel \mathbf{W} \parallel_{*} \\ &= {\rm arg} \min_{\mathbf{W}} {\rm Tr} \left(\mathbf{W}^T (\mathbf{I} + \alpha \mathbf{T} ) \mathbf{W} \right) + 2 {\rm Tr} (\mathbf{X}^T \mathbf{W}) + \beta \parallel \mathbf{W} \parallel_{*} + {\rm const} \end{aligned} \tag{7}

ここで、f_1(\mathbf{W}) = \left(\mathbf{W}^T (\mathbf{I} + \alpha \mathbf{T} ) \mathbf{W} \right) + 2 {\rm Tr} (\mathbf{X}^T \mathbf{W}) 、f_2(\mathbf{W}) =\beta \parallel \mathbf{W} \parallel_{*} とすると、近接勾配法の更新式は、

\begin{aligned} \mathbf{W}_{i + 1} = {\rm prox}_{\eta f_2}(\mathbf{W}_i - \eta \nabla f_1(\mathbf{W}_i)) \end{aligned} \tag{8}

となる。{\rm prox}_{\eta f_2}(\mathbf{Z})はトレースノルムの近接作用素であり、具体的な計算方法は以下の記事のトレースノルムの項目などを参照。

yamagensakam.hatenablog.com

なお、近接勾配法では関数f_1(\mathbf{W}) がL-平滑であればステップ幅\etaは\frac{1}{L}にすれば以下にすれば収束する（詳しくは参考文献[4]などを参照）。今回の場合、Lはf_1(\mathbf{W})のヘッセ行列\nabla^2 f_1(\mathbf{W})の最大固有値\lambda_{\max}なので、実装では\eta =\frac{1}{\lambda_{\max}}に固定した。

最後に式（6）解くことによって得られた \mathbf{\hat{W}}から、部分空間の基底を取り出す必要がある。これは \mathbf{\hat{W}}を特異値分解をして、上位K個の特異値に対応する左特異ベクトルを基底とすればよい。

実装は以下の通り。

def proximal_gradient_method(grad, prox, eta, x_init, max_iter=10000, tol=1e-5):
    x = np.copy(x_init)
    for _ in range(max_iter):
        z = x - eta * grad(x)
        x_new = prox(z)
        if np.linalg.norm(x - x_new) <= tol:
            break
        x = np.copy(x_new)
    return x_new

def trace_prox(Y, lam):
    U, S, V = np.linalg.svd(Y)
    return U @ np.diag(np.maximum(0.0, S - lam)) @ V[:S.shape[0], :]

def subspace_estimate_with_regularization(Y, T, alpha, beta):
    Hess = 2.0 * (np.eye(T.shape[0]) + alpha * T)
    grad = lambda W: -2.0 * Y + W @ Hess
    eta = 1.0/np.max(np.linalg.eigvalsh(Hess))
    prox = lambda Y: trace_prox(Y, beta * eta)
    W_init = np.copy(Y)
    W_hat = proximal_gradient_method(grad, prox, eta, W_init)
    return W_hat

def subspace_estimate_with_vertex_regularization(X, adj, k, alpha, beta):
    L = utils.compute_graph_laplacian(adj)
    W = subspace_estimate_with_regularization(X, L, alpha, beta)
    Uk, _, _ = np.linalg.svd(W.T)
    return Uk[:, :k], W

以降、この手法を「頂点領域正則化法」と呼ぶ。

グラフ信号に関する仮定や事前知識を正則化として取り入れる（周波数領域）

グラフ周波数領域における事前知識を入れることも可能である。ここでは、ノイズのない低次元信号の集合 \mathcal{A} \mathbf{D}をグラフフーリエ変換した\bar{\mathbf{W}} = \mathbf{U}^T \mathcal{A} \mathbf{D} に関する最適化問題を考える。つまり、観測信号をグラフフーリエ変換した集合を \bar{\mathbf{X}} として\parallel \bar{\mathbf{X}} - \bar{ \mathbf{W}} \parallel_{F}^2という損失関数、グラフ周波数領域に関する知識を埋め込んだ項{\rm Tr}(\bar{\mathbf{W}}^T \bar{\mathbf{T}} \bar{\mathbf{W}}) 、低次元化のためトレースノルム正則化項で構成される最適化問題となる。あくまで低次元にしたいのは\mathcal{A} \mathbf{D} であり\bar{\mathbf{W}}ではないため、トレースノルム正則化項は\parallel \mathbf{U} \bar{\mathbf{W}} \parallel_{*} とするべきだが、\mathbf{U}が直交行列であり特異値は直交行列をかけても変わらないことより\parallel \mathbf{U} \bar{\mathbf{W}} \parallel_{*} = \parallel \bar{\mathbf{W}} \parallel_{*}となる。従って、グラフ周波数領域に関する知識を埋め込んだ部分空間推定の最適化問題は

\begin{aligned} \mathbf{\hat{\bar{W}}} = {\rm arg} \min_{\bar{\mathbf{W}}} \parallel \bar{ \mathbf{X}} - \bar{\mathbf{W}} \parallel_F^2 + \alpha {\rm Tr}(\bar{\mathbf{W}}^T \bar{\mathbf{T}} \bar{\mathbf{W}}) + \beta \parallel \bar{\mathbf{W}} \parallel_{*} \end{aligned} \tag{9}

となる。頂点領域での議論と同様に求めた \mathbf{\hat{\bar{W}}} から部分空間の基底を取り出す必要があるが、これは \mathbf{U}\mathbf{\hat{\bar{W}}} = \mathbf{\hat{W}} と頂点領域に変換してから特異値分解すればよい。

以降、この手法を「周波数領域正則化法」と呼ぶ。

次にグラフ周波数領域における知識埋め込みの例として、グラフ信号のモデルがPGSグラフ信号であると仮定したときの定式化を示す。得られているグラフ信号がK次元PGS部分空間から生成された信号だと仮定すると、グラフスペクトルの周期はKなので以下のようなK要素ずらす行列

\mathbf{P} = \begin{bmatrix} 0 & 0 & \cdots & 0 & 1 & 0 & \cdots & 0 \\ 0 & 0 & \cdots & 0 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 & 0 & 0 & \cdots & 1 \\ 1 & 0 & \cdots & 0 & 0 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 & 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & & \vdots \\ 0 & 0 & \cdots & 1 & 0 & 0 & \cdots & 0 \\ \end{bmatrix} \tag{10}

を使ってグラフフーリエ係数ベクトル\bar{\mathbf{w}}をK要素分ずらした\mathbf{P} \bar{\mathbf{w}}は、元の\bar{\mathbf{w}}と一致するはずである。これによってPGS部分空間と仮定した場合の正則化項は、

\begin{aligned} {\rm Tr} \left( (\bar{\mathbf{W}} - \mathbf{P} \bar{\mathbf{W}} )^T (\bar{\mathbf{W}} - \mathbf{P} \bar{\mathbf{W}} ) \right) &= {\rm Tr} \left(\bar{\mathbf{W}}^T (\mathbf{I} - \mathbf{P} )^T(\mathbf{I} - \mathbf{P}) \bar{\mathbf{W}}\right) \\ &= {\rm Tr} \left(\bar{\mathbf{W}}^T \bar{\mathbf{T}} \bar{\mathbf{W}}\right) \end{aligned}

とすればよい（\bar{\mathbf{T}} = (\mathbf{I} - \mathbf{P})^T(\mathbf{I} - \mathbf{P}) ）。

以下がこの\bar{\mathbf{T}} を使った周波数領域正則化法の実装（頂点領域正則化法で示した関数を利用していることに注意）。

def subspace_estimate_with_frequency_regularization(X, adj, k, alpha, beta):
    gft = utils.GraphFourierTransformer(adj)
    X_ = gft.transform(X)
    P = (np.eye(X.shape[1], k=-k) + np.eye(X.shape[1], k= X.shape[1] - k))
    P_ = np.eye(X.shape[1]) - P
    W = subspace_estimate_with_regularization(X_, P_ @ P_.T, alpha, beta)
    Uk, _, _ = np.linalg.svd(W.T)
    return  gft.U @ Uk[:, :k], W

実験

実験設定

紹介した4つの部分空間推定手法の性能を測るため問題設定の節で紹介した3種類の生成モデルに基づくグラフ信号を発生させ、その信号から部分空間を推定してみる。またノイズの標準偏差も問題設定の節で述べたものを使う。 グラフの頂点数は全グラフ信号共通で200とし、観測信号数は50とする。また、各観測信号の係数\mathbf{d}は要素ごとに独立して標準正規分布で生成する。

また、本実験では頂点領域正則化法にて\mathbf{T}をグラフラプラシアン\mathbf{L}とする。つまりグラフが滑らかであるという仮定を埋め込む。周波数領域正則化法ではPGS信号であると仮定し、式（10）の\mathbf{P}を用いて\bar{\mathbf{T}} = (\mathbf{I} - \mathbf{P})^T(\mathbf{I} - \mathbf{P}) とする（今回周期は既知とし、\mathbf{P}のずらし量はその周期に基づいて決める）。

加えて、各手法推定する部分空間の次元を与える必要があるが、これは真の部分空間の次元を既知として実験した。

詳しくは実験のコードを見た方がよいと思うで、冒頭で述べたリポジトリにあるgraph_subspace_experiments.pyを参照。

正準角による部分空間の類似度比較

真の部分空間とどれだけ近い部分空間を推定できたかを評価するため正準角の\cosの和を指標とする（以降では\sum_i \cos \theta_iと表記）。正準角は2つの部分空間からなす角が最小となるように選んだベクトル同士のなす角を第1正準角、その選んだベクトルと直交するという制約を加えて同様の条件で選んだベクトル同士のなす角を第2正準角、というように正準角が小さいほど（その\cosが大きいほど）2つの部分空間は類似するとみなせる。正準角については、以下の記事でも触れているのでこちらも参照。

yamagensakam.hatenablog.com

実験結果

早速、各生成信号に対する各手法の\sum_i \cos \theta_iを示す。頂点領域正則化法と周波数領域正則化法についてはハイパーパラメータの設定値も併せて記載している。

帯域制限グラフ信号 区分的一定グラフ信号 PGSグラフ信号
主成分分析 8.18 4.02 4.11
グラフフーリエ変換 18.00 3.94 0.79
頂点領域正則化法 18.97
(\alpha=10, \beta=0) 		4.74
(\alpha=10, \beta=50) 		3.37
(\alpha=10, \beta=50)
周波数領域正則化法 6.39
(\alpha=100, \beta=0) 		0.89
(\alpha=100, \beta=0) 		4.99
(\alpha=100, \beta=0)

この結果からわかることとして、まず主成分分析は最下位になることはなくどの信号に対しても大きくは外していない。 次にグラフフーリエ変換は帯域制限グラフ信号と区分的一定グラフ信号については、信号の変化が滑らかで周波数成分も低周波領域に集中しているためそこそこ良い性能を得ているが、PGSグラフ信号は全周波数領域に満遍なく強い成分が存在するため、かなり悪い結果となっている。 頂点領域正則化法は帯域制限グラフ信号と区分的一定グラフ信号において、比較手法の中でトップの性能であった。やはりグラフ信号が滑らかであるという仮定が実際の生成信号がマッチしているのがポイントだろう。 一方で、周波数領域正則化法は帯域制限グラフ信号と区分的一定グラフ信号に対しては仮定とマッチせず推定能力が低いが、PGSグラフ信号に対しては、それをターゲットにした仮定を埋め込んでいるためか非常に高い推定精度となった。

また、実験をして分かったが帯域制限グラフ信号では頂点領域正則化の\beta=0とするのが最も良い結果となった。これは、グラフラプラシアンの二次形式を小さくすることが帯域制限することがマッチしており、低ランク化はむしろ余計な情報だったためと考えている。同様にPGSグラフ信号でも周波数領域正則化で\beta=0とするのが最も良い結果だった。これも、PGS信号と周期性の正則化は十分にマッチしているからだと考えている。

その一方で、区分的一定グラフ信号においては\betaをある程度大きくした方がよくなる。これは仮定が実態と少しズレていているため、トレースノルム正則化による低ランク化が補助として推定精度向上に寄与しているのだと思われる。これはPGSグラフ信号の時がより顕著で、\beta=0だと\sum_i \cos \theta_i=1.21とグラフフーリエ変換による方法よりはマシだがかなり低い推定性能だった。従って、トレースノルム正則化項は仮定と実態にズレがある時、そのズレを吸収する効果があると考えられる。

サンプリングと復元

この方法で推定した部分空間を使って前回の記事で書いたサンプリングと復元もやってみる。ただし比較のため部分空間はノイズが重畳した信号から推定した基底を使うが、サンプリングと復元はノイズなしのグラフ信号で実験する。サンプリングは200頂点のうち30点をランダムに選択している。

帯域制限グラフ信号
原信号 サンプリング後の信号

以下が各手法で推定した部分空間を使って復元した信号。やはり、部分空間を最もよく推定できている頂点領域正則化法が最も復元誤差が小さい（各結果の図の上部に復元誤差を記載）。

主成分分析 グラフフーリエ変換
頂点領域正則化法 周波数領域正則化法
区分的一定グラフ信号
原信号 サンプリング後の信号

以下が各手法で推定した部分空間を使って復元した信号。こちらも頂点領域正則化法の復元誤差が最も小さい。

主成分分析 グラフフーリエ変換
頂点領域正則化法 周波数領域正則化法
PGSグラフ信号
原信号 サンプリング後の信号

以下が各手法で推定した部分空間を使って復元した信号。こちらは周波数領域正則化法でかなり原信号に近い復元ができている。

主成分分析 グラフフーリエ変換
頂点領域正則化法 周波数領域正則化法

終わりに

本記事ではグラフ信号の部分空間推定について、頂点領域、周波数領域それぞれの仮定に基づく最適化問題に帰着して推定する方法を紹介した。

なお、今回はサンプリングと復元の実験においてはノイズなしのグラフ信号からサンプリングし復元したが、実際はノイズが重畳した信号からしかサンプリングできない。ノイズありのグラフ信号からサンプリングした際の復元方法についても試してみたいことがあるので、次の記事ではそれをやってみようと思う。

参考文献

[1] グラフ信号処理の基礎と応用: ネットワーク上データのフーリエ変換,フィルタリング,学習

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

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

Amazon

[2] 線形代数を基礎とする 応用数理入門: 最適化理論・システム制御理論を中心に

線形代数を基礎とする 応用数理入門: 最適化理論・システム制御理論を中心に (SGCライブラリ 187)

線形代数を基礎とする 応用数理入門: 最適化理論・システム制御理論を中心に (SGCライブラリ 187)

Amazon

[3] スパース性に基づく機械学習

スパース性に基づく機械学習 (機械学習プロフェッショナルシリーズ)

スパース性に基づく機械学習 (機械学習プロフェッショナルシリーズ)

Amazon

[4] 連続最適化アルゴリズム

連続最適化アルゴリズム

連続最適化アルゴリズム

Amazon

*1:以降も同様。

*2:仮定通り低次元部分空間からの線形写像でグラフ信号が生成され、かつ、ノイズがなければ主成分分析で正しい部分空間が得られる。

*3:ちなみに、\alpha=0とした低ランク近似問題は特異値分解により解ける。詳細はこの記事や同著者による書籍[2]が参考なる。また、グラフ信号は平均\mathbf{0}に中心化されている場合、特異値分解よって得られる部分空間は主成分分析で求めたものと一致する。

グラフ信号処理におけるサンプリングと復元

信号処理 線形代数

はじめに

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

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

グラフ信号処理の基礎と応用: ネットワーク上データのフーリエ変換,フィルタリング,学習 (次世代信号情報処理シリーズ 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}とした。

二乗誤差を評価値とする焼きなましの差分計算をEigenで行う（Atcoder Heuristic Contest 030）

競技プログラミング

はじめに

しばらくブログを書いてなかったのでAHC参加ネタを投稿。

先日Atcoder Heuristic Contest 030（AHC030）に参加し最終56位という結果だったが、提出した基本方針に二乗誤差を最小化する焼きなましが含まれていた。この方法では二乗誤差の計算を何度も行うことになるためその計算時間が肝になるのだが、線形代数ライブラリのEigenを使うと下手に差分計算するより愚直に全計算した方が速かったりする*1。とはいえ、少し工夫して差分計算にもEigenを活用すれば、状況によっては愚直計算よりも速くなる。

ということで、本記事では実際にAHC030で使ったEigenによる二乗誤差の差分計算の方法を説明する。このテクニックはAHC030に限らず二乗誤差を評価値として焼きなましする場合にも応用できる（かも？）。

AHC030の問題については以下を参照。 atcoder.jp

AHC030解法基本方針概要

解法基本方針についてはtwitterのリンクおよびそのスレッドを参照。

二乗誤差最小化によるポリオミノの配置最適化

まずは評価値の計算を行列演算の形で表す。占いをn回行ったとしてi回目の占いで得られる値をy_i、y_iをi番目の要素に持つn要素のベクトルを\mathbf{y}、i行目にi回目の占いにおいてどのマスを選んだか？という情報を持つn \times N^2の行列を\mathbf{X}（各列がマスに対応し、選ばれたマスなら1, そうでなければ0）とする*2。そして、\mathbf{w}を現状態のポリオミノの配置において、各マス何個のポリオミノが重なっているかを表すN^2要素のベクトルとし、これが最適化対象の変数である。これらを使って最適化問題として定式化すると、

\displaystyle \begin{aligned} &{\rm arg}\min_{\mathbf{p}} \parallel \mathbf{y} - \mathbf{X} \mathbf{w} \parallel_2^2 \\ &{\rm s.t. \ \ } \mathbf{w}はポリオミノの配置の組み合わせであり得る値 \end{aligned} \tag{1}

となる。これを1つのポリオミノを選択して置く場所を変える焼きなましによって最適化を図る。そのため二乗誤差\parallel \mathbf{y} - \mathbf{X} \mathbf{w} \parallel_2^2の計算を高速化できれば短い時間で焼きなましのループ回数が稼げ、より正解を出せる確率が上がる。

二乗誤差計算の高速化

Eigenの活用

C++においては行列計算を高速に実行できるライブラリEigenを活用すれば、愚直計算でもかなり速く計算できる。コードもこんな感じで1行で書ける。

double cost = (y - X * w).squaredNorm();

EigenではSIMD命令の活用、ループ変形等によりメモリアクセスの効率化などが行われているため演算回数だけ見れば遅くなりそうな処理も高速にできる場合がある。逆に差分計算はメモリアクセスが飛び飛びになりがちで、SIMD化もしにくいため処理時間上不利になることがある。

二乗誤差差分計算の行列表現

差分計算を行列演算で行うことを考える。ポリオミノの配置変更による\mathbf{w}の変化を\Delta \mathbf{w}とすると、二乗誤差の変化量Dは以下のようになる。

\displaystyle \begin{align} D &= \parallel \mathbf{y} - \mathbf{X} \mathbf{w} \parallel_2^2 - \parallel \mathbf{y} - \mathbf{X} \left(\mathbf{w} + \Delta \mathbf{w} \right) \parallel_2^2 \\ &=2 \left(\mathbf{y}^T \mathbf{X} \Delta \mathbf{w} - \mathbf{w}^T \mathbf{X}^T \mathbf{X} \Delta \mathbf{w} \right) - \Delta \mathbf{w}^T \mathbf{X}^T \mathbf{X} \Delta \mathbf{w} \end{align} \tag{2}

この式の中で \mathbf{X}^T \mathbf{X} は焼きなましの中で変化しないので、事前に計算しておくことができる。 加えて、焼きなましの1ループでは1つのポリオミノしか動かさないので、\mathbf{w}の中で変化する要素はごく一部、つまり\Delta \mathbf{w}の大半の要素は0になる。このような多くの要素が0のベクトルを扱うためEigenにはSparseVectorが用意されている。\Delta \mathbf{w}をSparseVector型として扱うことで、式（2）の計算がだいぶ速くなる。コードはこんな感じ。

Eigen::SparseVector<double> sdw(N * N);
vector<Eigen::Triplet<double>> trivec;
rep(i, N * N){
    if(abs(dw(i)) > 1e-5){
        sdw.insert(i) = dw(i);
    }
}
VectorXd XXdw = XX * sdw;
VectorXd Xdw = X * sdw;
double D = 2.0 * (y.dot(Xdw) - old_w.dot(XXdw)) - dw.dot(XXdw);

比較

上記2方法で評価値計算をする焼きなましがatcoderのジャッジサーバー上で400[ms]内に何回回るかをnを変えて比較する。その際のテストデータはseed=1とする（N = 13, M = 6, \epsilon = 0.04）。結果は以下の通り。

n 全計算 差分計算
20 814301 342042
40 598892 338786
60 312957 320429
80 311225 319219
100 175905 303552
120 226957 297810
140 143947 275129
160 176061 250388
180 115353 249101
200 146547 229857

このように\mathbf{X}の行数（=占い回数）nが小さい間は愚直に全計算した方が回数を稼げるが、n \geq 100あたりで逆転し差分計算した方が速い*3。従って、実際提出したコードでは占いした回数nに応じて愚直全計算と差分計算を切り替える処理とした。

ちなみに

今回の手法を使ったseed=0のビジュアライザは以下の通り。

*1:焼きなましでは1回の更新の中で状態が一部しか変化しないことがしばしばある。そのため基本的にはその差分のみを使って評価値を計算できれば評価値計算の時間を短縮につながる、ヒューリスティックコンテストにおいてはこの差分計算は典型的なテクニックの1つとなっている。

*2:実際の提出したコードでは、問題文にある占いで得られる平均値\muの式に基づき\mathbf{v}と\mathbf{X}を補正している。

*3:愚直全計算で顕著だがnが8の倍数の時にLoop回数が盛り返す。これはおそらくEigenの内部でAVX512命令を使っているのが関係していると思われる。

IIRフィルタの設計問題を焼きなまし法で解いてみる

連続最適化 信号処理

これは数理最適化 Advent Calendar 2022 の 23 日目の記事です。

はじめに

ディジタルフィルタとは、離散時間信号に対し所望の周波数帯域の成分だけを通過させて他は阻止する機能のことを言う。中でもIIR (Infinite Impulse Response) フィルタはフィードバック機構を用いることにより無限長の畳み込み計算を行うフィルタのことを言う。 IIRフィルタは特定の制約を満たしつつ所望の周波数特性が得られるよう適切に設計する必要がある。フィルタ設計方法として様々な方法が提案されており、その中には最適化手法を利用したものも数多くある。

本記事は主に参考文献[1]の1章の内容を基に、最適化によるIIRフィルタ設計を実装してその感触を掴んでみたという内容になっている。 文献ではPSO（Particle Swarm Optimization）で最適化問題を解いているが、今回は手始めに（個人的に）実装ハードルが低い焼きなまし法による設計を試している。

ディジタルフィルタの基本

まずはディジタルフィルタについての基本事項を述べる。ここでは最低限のみ触れるので詳しくは参考文献[2][3]等を参照のこと。

フィルタの設計問題

IIRフィルタの処理は、離散時間の入力信号をx [ i ] 、出力信号をy [ i ] \ \ (i = 0, \cdots, \infty)とすると、以下のような差分方程式で表現できる。ただし、i \lt 0の場合はx[i] =0, y[i ] = 0とする。

\displaystyle y [i ] = -\sum_{m = 1}^{M} b_m y [i - m ] + \sum_{n = 0}^{N} a_n x [i - n ] \tag{1}

このb_m \in \mathbb{R}, (m = 1, \cdots, M)とa_n \in \mathbb{R}, (n = 0, \cdots, N)がフィルタ係数で、所望の周波数特性が得られるようにこの係数を決めるのがフィルタの設計にあたる。

式（1）の右辺第1項は過去の出力値を利用するフィードバック機構になっており、第2項が過去の入力値を利用するフィードフォワード機構になっている。このフィードバック機構があることによってフィルタの設計が不適切な場合、出力信号が無限に発散してしまう可能性がある。このことを不安定といい、逆に発散する心配のないフィルタは安定という*1。従って安定かつ所望の周波数特性が得られるようにIIRフィルタを設計する必要がある。

Z変換

連続信号を周期Tでサンプリングした離散信号はデルタ関数を用いて以下のように表すことができる。

\displaystyle x_s(t) = \sum_{n = 0}^{\infty} x(nT) \delta(t - nT) \tag{2}

このサンプリングされた信号をラプラス変換すると、

\displaystyle \begin{align} \mathcal{L}[x_s(t)] &= \int_{0}^{\infty} \sum_{n = 0}^{\infty} x(nT) \delta(t - nT) \exp(-st) dt \\ &= \sum_{n = 0}^{\infty} x(nT) \exp(-nsT) \end{align} \tag{3}

となる。ここでz = \exp(sT)としたのがZ変換であり、次のような式になる*2。

\displaystyle X(z) = \sum_{i = 0}^{\infty} x [ i ] z^{-i} \tag{4}

ラプラス変換ではs = j\omegaを代入することでフーリエ変換を得ていたが、Z変換も同様にz = e^{j\omega}を代入すれば離散時間フーリエ変換を得ることができる（\omegaは角周波数で周波数をfとすると\omega=2\pi fである）。そのため、この操作で周波数領域における特性を知ることができる。ただし、離散時間フーリエ変換ができるのは\sum_{i = 0}^{\infty} | x [ i ]| \lt \inftyという条件を満たす場合のみであるため、これを満たさない信号に対してz = e^{j\omega}を代入する操作を行っても意味はない。

伝達関数

フィルタの周波数特性を見るには伝達関数が便利である。式（1）の両辺をZ変換にすると、

\displaystyle Y (z) = -\sum_{m = 1}^{M} b_m Y(z) z^{-m} + \sum_{n = 0}^{N} a_n X (z) z^{-n}\tag{5}

となる。ここでH(z) = Y(z)/X(z)は伝達関数と呼ばれ、式（5）を整理してIIRフィルタの伝達関数を求めると以下のようになる。

\displaystyle H (z) =\frac{Y(z)}{X(z)} = \frac{\sum_{n = 0}^{N} a_n z^{-n}}{1 + \sum_{m = 1}^{M} b_m z^{-m}} \tag{6}

この伝達関数に対して、 z = e^{j \omega}を代入すると伝達関数を離散時間フーリエ変換した結果が得られる。これを以下のようにH(\omega)と定義する。

\displaystyle H (\omega) = \frac{\sum_{n = 0}^{N} a_n e^{-jn \omega}}{1 + \sum_{m = 1}^{M} b_m e^{-jm\omega}} \tag{7}

H(\omega)はフィルタの周波数特性であり、このフィルタによって各周波数成分がどれくらい減衰するかという振幅特性| H(\omega) |と、どれくらい位相が遅れるかという位相特性\angle H(\omega)がわかる*3。振幅特性はどの周波数成分を通過させてどの周波数成分を阻止するかを左右するのでフィルタにおいて当然重要な特性だが、位相特性も重要で通過させたい帯域に関してはなるべく線形であることが望ましい。位相特性が線形でないと、出力信号の波形が歪んでしまうためである。

極と安定性

伝達関数の分母=0を満たす点を極という。フィルタが安定であるためには全ての極の絶対値が1未満である必要がある。そのため、分母のフィルタ係数b_mは1 + \sum_{m = 1}^{M} b_m z^{-m} = 0の解が全て絶対値1未満となるように定めなければならない。このことは最適化問題における制約条件となる。

焼きなまし法の適用

定式化

まず所望の周波数特性をH_D(\omega)とする。例えば、以下のような周波数特性を設定する。

\displaystyle H_D(\omega) = \begin{cases} e^{-j \omega \tau_d} & \omega \in \Omega_p \\ 0 & \omega \in \Omega_c \end{cases} \tag{8}

ここで、\tau_dは所望の群遅延、\Omega_pは通過させたい角周波数の集合、\Omega_cは阻止したい角周波数の集合である。 この理想的な周波数特性に近づくよう式（6）、（7）の係数を決めたいので、係数を並べたベクトル\mathbf{x} = [a_0, \cdots, a_N, b_1, \cdots, b_M ] を変数として、以下のような最適化問題を解くことになる。

\displaystyle \begin{align} \min_{\mathbf{x}}& \ \ J(\mathbf{x})\\ \rm{s.t.}& \ \ | p_i(\mathbf{x}) | < 1 \ \ (i = 1, \cdots, M) \end{align} \tag{9}

ここでJ(\cdot)は理想周波数特性との誤差関数、p_i(\mathbf{x})はフィルタ係数が\mathbf{x}の時の極である。誤差関数は色々考えられるがよく使われるのは以下の2つ*4。

  • 最大絶対誤差
\displaystyle J(\mathbf{x}) =\max_{\omega \in \Omega} |H_D(\omega) - H(\omega; \mathbf{x})| \tag{10}
  • 平均2乗誤差
\displaystyle J(\mathbf{x}) = \frac{1}{|\Omega|}\sum_{\omega \in \Omega} |H_D(\omega) - H(\omega; \mathbf{x})|^2 \tag{11}

\Omegaは誤差計算に使用する周波数点の集合である。

さて、式（9）には極に関する制約条件があるがこのままでは扱いにくいので、制約条件違反をペナルティ項として組み込んだ新たな目的関数 F(\mathbf{x})を作り、こちらを焼きなましで解くことにする。

\displaystyle \begin{align} \min_{\mathbf{x}}& \ \ F(\mathbf{x}) = \min_{\mathbf{x}} \ \ (J(\mathbf{x}) + c \phi(\mathbf{x})) \\ \phi(\mathbf{x}) &= \begin{cases} \max_{i} |p_i(\mathbf{x})|^2 & \max_i | p_i(\mathbf{x}) | \geq 1 \\ 0 & \max_i | p_i(\mathbf{x}) | < 1 \end{cases} \end{align} \tag{12}

cはペナルティの強さを決めるハイパーパラメータになっている。

極の計算方法

式（12）のままだと更新の度に分母=0のM次方程式を解くなどして安定性を判断する必要がある。これでは効率が悪いため少し工夫する。まず簡単のためフィルタ係数の次数N, Mは偶数とする*5。そうすると式（6）は次のように変形できる。

\begin{align} \displaystyle H (z) &= \frac{\sum_{n = 0}^{N} a_n z^{-n}}{1 + \sum_{m = 1}^{M} b_m z^{-m}} \\ &= a0 \frac{\prod_{n = 1}^{N/2} (1 - q_n z^{-1})(1 - \bar{q}_n z^{-1})}{\prod_{m = 1}^{M/2} (1 - p_m z^{-1})(1 - \bar{p}_m z^{-1})} \\ &= a0 \frac{\prod_{n = 1}^{N/2} (1 - 2\rm{Re} (q_n) z^{-1} + |q_n|^{2} z^{-2}) } {\prod_{m = 1}^{M/2} (1 - 2\rm{Re}(p_m) z^{-1} + |p_m|^{2} z^{-2}) } \end{align} \tag{13}

ここで、a_0 \in \mathbb{R}, p_m, \bar{p}_m, q_n, \bar{q}_n \in \mathbb{C}であり、\bar{p}_m, \bar{q}_nはp_m, q_nの複素共役である。こう変形するとシステムの極はp_m, \bar{p}_mなので直接的にわかる。従って、\mathbf{x} = [a_0, p_1, \cdots, p_{M/2}, q_1, \cdots, q_{N/2} ] *6のように直接極を変数として扱う形で式（12）の最適化問題を解けば更新毎の安定性の判別は簡単にできる。

最適化アルゴリズム

ここでは実際に実装したアルゴリズムの流れを述べる。焼きなまし法自体の解説は割愛するので、詳しくは参考文献[5]等を参照。

まず、理想周波数特性の関数H_D(\cdot)、誤差計算に用いる周波数点集合\Omega、フィルタ次数M, N、ペナルティの強さc、繰り返し回数L、更新量の標準偏差\sigma、初期温度T、冷却率\alpha、使用する誤差関数J(\cdot)はアルゴリズムのパラメータとして入力する。そして以下の処理を行う。

  1. \mathbf{x} = [a_0, p_1, \cdots, p_{M/2}, q_1, \cdots, q_{N/2} ] の各要素を正規分布\mathcal{N}(0, 1)に従う乱数で初期化する。
  2. \mathbf{x}' = \mathbf{x} + \Delta \mathbf{x} を求める。ここで \Delta \mathbf{x}は各要素が正規分布 \mathcal{N}(0, \sigma^2)に従う乱数。
  3. F(\mathbf{x}') \lt F(\mathbf{x} )であれば、 \mathbf{x} \leftarrow \mathbf{x}'に更新する（F(\cdot)は式（12）の目的関数）。 F(\mathbf{x}') \geq F(\mathbf{x})の場合は、\exp ( \frac{F(\mathbf{x}) - F(\mathbf{x}')}{T})の確率で \mathbf{x} \leftarrow \mathbf{x}'に更新する。
  4. 温度を T \leftarrow \alpha Tとして更新する。
  5. 繰り返し回数がL回未満なら2.に戻る。L回に達したら、探索した中で最もF(\mathbf{x})が小さかった\mathbf{x}を出力して終了。

実装と実験

ここまで述べたことを実装して実験してみた。なお、ここに示す実験用コードはgithubにも上げている。

github.com

焼きなまし法の実装と実験

焼きなまし処理のコード

まずは焼きなまし法でフィルタ係数を求める処理の実装を以下に示す。

import pickle
import numpy as np
import matplotlib.pyplot as plt

def compute_mean_squared_error(H_, ref_h, omega_list):
    h = np.zeros(len(omega_list), dtype=np.complex)
    for i, omg in enumerate(omega_list):
        h[i] = H_(omg)

    e = np.real((ref_h - h) * np.conjugate(ref_h - h))
    return np.mean(e)

def compute_abs_max_error(H_, ref_h, omega_list):
    h = np.zeros(len(omega_list), dtype=np.complex)
    for i, omg in enumerate(omega_list):
        h[i] = H_(omg)
    e = np.abs(ref_h - h)
    return np.max(e)

# フィルタの周波数特性を計算する関数
def frequency_characteristic_func(a0, p_array, q_array, omega):
    denomi = 1.0
    for p in  p_array:
        denomi *= (1 - p * np.exp(-1j * omega)) *\
            (1 - np.conjugate(p) * np.exp(-1j * omega))

    numer = 1.0
    for z in q_array:
        numer *= (1 - z * np.exp(-1j * omega)) *\
            (1 - np.conjugate(z) * np.exp(-1j * omega))
    return a0 * (numer/denomi)

def obj_func(p_array, q_array, a0, ref_h, omega_list, c, error_func):
    H = lambda omega: frequency_characteristic_func(a0, p_array, q_array, omega)
    max_p2 = np.max(np.real(p_array * np.conjugate(p_array)))
    e = error_func(H, ref_h, omega_list)

    return e + (c * max_p2 if max_p2 >= 1.0 else 0.0)

# 焼きなまし法のメイン処理
def annealing(HD, omega_list, M, N, c, L, sigma, T, alpha, error_func,
              init_a0=None, init_p_array=None, init_q_array=None):

    # 理想周波数特性の関数から各周波数点に対する周波数特性を計算する
    h = np.zeros(len(omega_list), dtype=np.complex)
    for i, omg in enumerate(omega_list):
        h[i] = HD(omg)

    # 初期化
    if init_a0 is None:
        a0 = 1.0
    else:
        a0 = init_a0
    if init_p_array is None:
        p_array = np.random.randn(M//2) + 1j * np.random.randn(M//2)
    else:
        p_array = np.copy(init_p_array)
    if init_q_array is None:
        q_array = np.random.randn(N//2) + 1j * np.random.randn(N//2)
    else:
        q_array = np.copy(init_q_array)

    # 目的関数F
    F = lambda p_array_, q_array_, a0_: \
        obj_func(p_array_, q_array_, a0_, h, omega_list, c, error_func)

    now_cost = F(p_array, q_array, a0)

    # 最終的に出力するフィルタ係数を初期化
    best_p_array = np.copy(p_array)
    best_q_array = np.copy(q_array)
    best_a0 = a0
    best_cost = now_cost

    for i in range(L):
        # ロールバック用に保持
        save_a0 = a0
        save_p_array = np.copy(p_array)
        save_q_array = np.copy(q_array)

        a0 += np.random.randn() * sigma
        p_array += (np.random.randn(M//2) + 1j * np.random.randn(M//2)) * sigma
        q_array += (np.random.randn(N//2) + 1j * np.random.randn(N//2)) * sigma

        tmp_cost = F(p_array, q_array, a0)
        if np.exp((now_cost - tmp_cost)/T) > np.random.rand():
            now_cost = tmp_cost
            print("now_cost=", now_cost, "T=", T, "i=", i)
            if now_cost < best_cost:
                best_a0 = a0
                best_p_array = np.copy(p_array)
                best_q_array = np.copy(q_array)
                best_cost = now_cost
                print("best_cost = ", best_cost)
        else:
            # ロールバック
            a0 = save_a0
            p_array = np.copy(save_p_array)
            q_array = np.copy(save_q_array)
        T *= alpha

    return best_a0, best_p_array, best_q_array, best_cost

基本的には先述のアルゴリズムに書いた通りの実装になっている。annealingという関数が実際の焼きなましの処理を行う関数である。F(\cdot)の計算などを工夫すればもっと効率化できるかもしれない。

実験

それでは実際に理想的な周波数特性を設定し、この焼きなまし処理を用いてフィルタ係数を求めてみる。 理想的な周波数特性は式（8）の形式において、通過域の集合を\Omega_p = \{\omega |0.4 \pi \leq \omega \leq 0.6 \pi \} 、阻止域の集合を\Omega_c = \{\omega | 0 \leq \omega \leq 0.3 \pi, 0.7 \pi \leq \omega \leq \pi \} 、群遅延を\tau_d = 12と設定して実験した。なお、\Omega_p, \Omega_cを見てわかる通り、急激な変化を避けるため通過域と阻止域の間に特に周波数特性を指定しない遷移域も設けている。この理想周波数特性のグラフは以下のようになる。

理想周波数特性

焼きなましのパラメータについては、

M = 8, N = 8, c = 5, L = 10000, \sigma = 0.01, T=10, \alpha = 0.998

とした。誤差関数については両方の方法を試している。

また、問題の性質上、得られる結果が初期値や乱数に依存して大きく変わるので、この実装では初期値を変えながら100回最適化処理を行い、最もF(\mathbf{x})が小さくなったものを採用している。

このような実験を行うコードは以下のようになっている。

def visualize_characteristic(H, omega_list, save_file_name):
    h = np.zeros(len(omega_list), dtype=np.complex)
    for i, omg in enumerate(omega_list):
        h[i] = H(omg)

    fig = plt.figure()
    amp_h = np.abs(h)
    ax1 = fig.add_subplot(2, 1, 1)
    ax1.plot(omega_list, amp_h, marker=".")
    ax1.set_ylim(-0.1, 1.1)
    ax1.set_ylabel("amplitude")

    angle_h = np.angle(h)
    ax2 = fig.add_subplot(2, 1, 2)
    ax2.plot(omega_list, angle_h, marker=".")
    ax2.set_ylim(-np.pi - 0.2, np.pi + 0.2)
    ax2.set_xlabel("omega")
    ax2.set_ylabel("phase")
    fig.savefig(save_file_name)
    plt.close()

if __name__ == '__main__':
    np.random.seed(0)

    # 理想的な周波数特性の関数
    HD = lambda omg: np.exp(- 1j * 12 * omg) \
        if omg >= 2 * 0.2 * np.pi and omg <= 2 * 0.3 * np.pi  else 0.0

    # 誤差関数に利用する誤差計算に使用する周波数点の集合
    omega_list1 = np.linspace(0, 2 * 0.15 * np.pi, 200)
    omega_list2 = np.linspace(2 * 0.2 * np.pi, 2 * 0.3 * np.pi, 200)
    omega_list3 = np.linspace(2 * 0.35 * np.pi, np.pi, 200)
    omega_list = np.concatenate([omega_list1, omega_list2, omega_list3])

   # 理想的な周波数特性のグラフを作成
    visualize_characteristic(HD, omega_list, "ideal.png")

    # 100回焼きなましをして最もよかったものを採用
    best_cost = 1e9
    for i in range(100):
        a0_, p_array_, q_array_, cost =\
            annealing(HD, omega_list, 8, 8, 5,
                      10000, 0.01, 10, 0.998,
                      compute_mean_squared_error) # 実験する誤差関数を指定
        if cost < best_cost:
            a0 = a0_
            p_array = np.copy(p_array_)
            q_array = np.copy(q_array_)
            best_cost = cost

    # 求めたフィルタの周波数特性の関数
    H = lambda omega: frequency_characteristic_func(a0, p_array, q_array, omega)

    # 求めたフィルタの周波数特性のグラフを作成
    visualize_characteristic(H, np.arange(0, np.pi, 0.01), "solution.png")

    # 制約条件を満たしていることを確認
    for p in p_array:
        print("|p|=", np.abs(p))

    print("final cost=", best_cost)

    # 求めたフィルタ係数を保存
    save_dict = {"a0":a0, "p_array":p_array, "q_array":q_array}
    with open('filter_coef.pickle', 'wb') as f:
        pickle.dump(save_dict, f)

これを実行して得たIIRフィルタの周波数特性は以下のようになった。

  • 誤差関数が最大絶対誤差の場合

求めたIIRフィルタの周波数特性（最大絶対誤差）

  • 誤差関数が平均2乗誤差の場合

求めたIIRフィルタの周波数特性（平均2乗誤差）

目視の印象だが通過域については最大絶対誤差で求めた方が平均2乗誤差に比べて理想の周波数特性に近いように見える。一方で、阻止域については最大絶対誤差だと振幅特性が0付近まで落ち切っておらず平均2乗誤差の方が優位に見える。

フィルタ処理の実装と実験

フィルタ処理のコード

得られたフィルタ係数を用いて実際の離散信号に対してフィルタをかけてみる。式（13）を展開して式（6）の形式にした上で式（1）の計算でフィルタをかけても良いが、ここでは式（13）の形式のまま2次のフィードバック（もしくはフィードフォワード）処理を縦続に接続することでフィルタ処理を実現している。つまり、2次フィルタの出力をそのまま次の2次フィルタに入力するような形にしている。これについては実装を見た方がわかりやすいだろうと思うので早速実装を貼る。

class IIRFilter:
    def __init__(self, a0_, p_array_, q_array_):
        self.a0 = a0_
        self.p_array = np.copy(p_array_)
        self.q_array = np.copy(q_array_)

    def _feedback_filter(self, x, p):
        b1 = -np.real(p) * 2
        b2 = np.real(p * np.conjugate(p))
        y = np.zeros(x.shape[0])
        for i in range(x.shape[0]):
            s = x[i]
            if i - 1 >= 0:
                s += -b1 * y[i - 1]
            if i - 2 >= 0:
                s += -b2 * y[i - 2]
            y[i] = s
        return y

    def _feedforward_filter(self, x, q):
        a1 = -np.real(q) * 2
        a2 = np.real(q * np.conjugate(q))
        y = np.zeros(x.shape[0])
        for i in range(x.shape[0]):
            s = x[i]
            if i - 1 >= 0:
                s += a1 * x[i - 1]
            if i - 2 >= 0:
                s += a2 * x[i - 2]
            y[i] = s
        return y

    def filtering(self, x):
        y = np.copy(x)
        for q in self.q_array:
            y = self._feedforward_filter(y, q)
        for p in self.p_array:
            y = self._feedback_filter(y, p)
        y *= self.a0
        return y
実験

この実装を使っていくつかの信号にフィルタをかけてみる。この実験用コードについては上述のgithubを参照。

  • \cos (0.1 \pi t)

まずは低周波数帯の遮断域にある周波数の信号をフィルタリングしてみる。以下のように減衰された信号が得られるが、誤差関数が最大絶対誤差の場合は平均2乗誤差の場合に比べて減衰が甘い。この結果は上述の振幅特性通りの結果と言える。

最大絶対誤差によるフィルタリング結果

平均2乗誤差によるフィルタリング結果

  • \sin (0.4 \pi t)

次に通過域ギリギリにある周波数の信号を試す。誤差関数が最大絶対誤差の場合はt=0付近以外ほとんど減衰のない出力が得られているが、平均2乗誤差では少し減衰した結果になっている。

最大絶対誤差によるフィルタリング結果

平均2乗誤差によるフィルタリング結果

  • \sin (1.0 \pi t)

高周波数帯の遮断域にある周波数の信号。低周波数帯の時と同様に、誤差関数が最大絶対誤差の場合は減衰しきれていない結果になっている。

最大絶対誤差によるフィルタリング結果

平均2乗誤差によるフィルタリング結果

  • \cos (0.1 \pi t) + \sin (0.4 \pi t) + \sin (1.0 \pi t)

最後に混合した信号で試す。通過域の成分のみ抽出できているように見える。波形としてきれいなのは、遮断域をきちんと遮断できている平均2乗誤差の方か。

最大絶対誤差によるフィルタリング結果

平均2乗誤差によるフィルタリング結果

所感

上でも述べたが良い解が得られるかどうかは初期値と乱数にかなり依存し、あまり好ましくない局所最適解が得られるケースも多々あった。特に誤差関数として最大絶対誤差を選択した場合では、振幅特性がほとんど0になるような解が得られるケースも多かった。そのため今回は初期値を変えながら何度も焼きなましを行う方法で対処したが、そもそも焼きなましで解くには非線形性が強すぎる問題だった可能性がある。 そう考えると、PSOのように複数の個体を有する最適化手法を使うのは妥当な選択なのかもしれない。この辺は次の記事辺りでも比較してみたいところ。

また、フィルタの次数M, Nを増やすほど最適化が難しくなる。変数が増えるので当然と言えば当然だが、今回の試行ではM=8, N=8（9変数）程度がまともな解が得られる限界だったのは少し残念だった。

加えて、通過域と阻止域の境界となる不連続な領域は遷移域として誤差関数の計算に含めないことも適切な最適化を得るのに思いのほか寄与した。不連続に変化する点は線形フィルタの周波数特性では得ることができず、より最適化が難しい問題になるのであろうと予想する。

最後に

焼きなまし法自体は比較的簡単に様々な最適化問題に適用できるため、今回のように最適化問題の感触を掴むのには良いかもしれない。ただし、そこからより良い解をより速く得たい場合は工夫や調整が必要だろう*7。今回ももう少し工夫すればより良いフィルタが得られた可能性もある。

IIRフィルタの設計には既存のアナログフィルタから変換する等最適化を利用しない方法もある。この方法で得られたフィルタ係数を初期値として最適化するのも面白いかもしれない。

参考文献

[1] 今回の記事を書くきっかけになった書籍。今回は1章の内容を参考にしているが、他の章には他にも凸関数に近似して解く方法や、他の信号処理の問題を最適化手法を用いて解く事例が載っていて面白い。

ディジタル信号処理におけるシステム最適化技術: 基礎技術から音声 ・音響信号処理,産業応用と情報システムまで

ディジタル信号処理におけるシステム最適化技術: 基礎技術から音声 ・音響信号処理,産業応用と情報システムまで

Amazon

[2] Z変換など信号処理における基本事項は以下の書籍を参考にしている*8。

ディジタル信号処理―基礎理論と方法論

ディジタル信号処理―基礎理論と方法論

Amazon

[3] 信号処理の基本に関してはこちらも参考にした。ネットで見れる最も丁寧なディジタル信号処理解説の1つだと思う。

www.ic.is.tohoku.ac.jp

[4] 電子情報通信学会が公開してるディジタルフィルタの解説。こちらの2章も今回の記事を書くにあたって参考にした。

www.ieice-hbkb.org

[5] お世話になっている最適化の書籍。この書籍には焼きなまし法についての解説もある。

しっかり学ぶ数理最適化 モデルからアルゴリズムまで (KS情報科学専門書)

しっかり学ぶ数理最適化 モデルからアルゴリズムまで (KS情報科学専門書)

Amazon

*1:第2項のみのフィルタはFIR（Finite Impulse Response）フィルタと呼ばれる。FIRフィルタはフィードバック機構がないため常に安定だが、一般にIIRフィルタの方が少ない次数で良い周波数特性が得られる。

*2:式（2）、（3）では便宜上xを時刻の関数で表していたが、式（4）では添え字の形式で表現している。つまり、t = 0, T, 2T, \cdotsがそれぞれi = 0, 1, 2, \cdotsに対応している。

*3:このような操作で周波数特性が得られるのはフィルタが安定の場合に限る。フィルタが不安定の場合は伝達関数の絶対値総和が発散し、そもそも離散時間フーリエ変換が成り立たたない。

*4: H(\omega; \mathbf{x})は式（7）のH(\omega)と同じものだが、\mathbf{x}の変数であることを明示的に示すためこのように表現している。

*5:奇数の次数にもできるが少し面倒なので、ここでは偶数にしている。

*6: \mathbf{x} の中に実数の要素と複素数の要素が混在していてあまりよろしくないが、面倒なのでこういう表記をしている。

*7:競プロのマラソンマッチ（ヒューリスティックコンテスト）と呼ばれる形式でも焼きなましが有効になる問題が多々ある。しかし、それでもハイスコアを狙うにはただ単に焼きなましを適用するだけでなく、初期値の工夫、定式化の工夫、高速化等様々な工夫が必要になる。

*8:10年以上前に購入して個人的には割とお世話になっている書籍だし内容も分かりやすいと思うが、どうやら今は中古でしか手に入らないらしい。

カーネル関数は何故『半』正定値でよいのか

関数解析 カーネル法 機械学習

はじめに

タイトルの通りカーネル関数の話。内積のことを考えるとカーネル関数から作られるグラム行列は正定値でないとダメな気がしたがそんなことはなくて半正定値であればよい、それは何故か？という話。

いくつかの書籍に書いてある話なのでわざわざ記事にしなくてもよいかもしれないが、自分が抱いた疑問やその解消に至るまでの経緯などを含めて備忘録として残しておく。

なお、カーネル関数については以前書いた以下の記事でも触れているので参考までに。

yamagensakam.hatenablog.com

カーネル関数と再生性

本題に入る前に前提となるカーネル関数の定義や対象とする空間および再生性という性質について簡単に述べておく。

まず集合\mathcal{X}上の2変数関数であって、対称性と半正定値性を有する関数kをカーネル関数と呼ぶ。 これは、任意 x_1 , \cdots, x_n \in \mathcal{X}に対しi行j列目がk(x_i, x_j)のグラム行列\mathbf{K} を作ると、それが半正定値対称行列になることを意味する。

また、カーネル関数kについて、第2引数をx \in \mathcal{X}で固定した1変数関数をk(\cdot, x) = k_x(\cdot) = k_xと表記する。

そして、この1変数関数を用いて内積空間\mathcal{V}を以下のように定義する。

\displaystyle \mathcal{V} = \left \{ \sum_{i=1}^{n} a_i k_{x_i} \ \middle| \ x_i \in \mathcal{X}, \ a_i \in \mathbb{R}, \ n \in \mathbb{N} \right\} \tag{1}

この空間の内積はf = \sum_{i = 1}^{n} a_i k_{x_i}, \ g = \sum_{j = 1}^{m} b_j k_{y_j} \in \mathcal{V}として、

\displaystyle \langle f, g \rangle = \sum_{i=1}^{n} \sum_{j=1}^{m} a_i b_j k(y_j, x_i) \tag{2}

と定義する。後述の零関数の例にもある通りf \in \mathcal{V}の表現方法は複数存在する可能性がある（つまり、a_iやx_iが異なっていても同じfになることがある）。それでも式（2）で定義する内積はwell-definedであり、f, gの表し方によらず1つに定まる。

このように定義した内積空間では、任意のf \in \mathcal{V}, \ x \in \mathcal{X}に対し以下のような再生性という性質を有する。

\displaystyle f(x) = \langle f, k_x \rangle \tag{3}

これは関数fにxを代入する操作が内積の形で表せることを意味し、この性質から様々な定理が導かれる。

ちなみに、このままでは完備性を有しておらずヒルベルト空間とは呼べない（と思われる*1）。この\mathcal{V}に対し完備化という操作を施し完備性を有した空間は再生核ヒルベルト空間になるが、今回の内容には絡まないのでこれ以上は触れない。

疑問点

このように定めた内積空間\mathcal{V}において

\displaystyle \langle f, f \rangle = \sum_{i = 1}^{n} \sum_{j = 1}^{n} a_i a_j k(x_i, x_j) = \mathbf{a}^T \mathbf{K} \mathbf{a} \tag{4}

を考える（\mathbf{a}はa_iをi番目の要素に持つベクトル）。カーネル関数が半正定値ということはグラム行列\mathbf{K}が半正定値で \rm{rank} (\mathbf{K}) \lt nにもなり得る。ということは、 \rm{rank} (\mathbf{K}) \lt nであれば\langle f, f \rangle = 0となる\mathbf{0}以外の\mathbf{a}が存在することになる。

ところが、内積の定義を考えると\langle f, f \rangle = 0 \Leftrightarrow f = 0のはずで、それにもかかわらず\mathbf{a}が\mathbf{0}でなくとも式（4）を0にできるのはなんか奇妙、というのが今回の論点。

疑問の解消

そもそもの話

考えてみればf = 0とはどういうことだろう？f = \sum_{i = 1}^{n} a_i k_{x_i}でa_i = 0\ ( i = 1, \cdots, n) というなら分かるが、それ以外に意味することはあるのだろうか？

これはベクトル空間に立ち返るとfは零ベクトルであることを意味する。零ベクトルなので任意のg \in \mathcal{V}に対して、f + g = g + f = gとなればよい。 関数なので表現方法が異なっていても入出力関係が変わらなければよい。すなわち、任意のx\in \mathcal{X}に対してf(x) + g(x) = g(x) + f(x) = g(x)であればfは零ベクトルと言える*2。これを達成するfは何を入力しても0を出力する関数、つまり零関数である。fが零関数であれば、a_iが全て0である必要もない。

\langle f, f \rangle = 0 \Leftrightarrow f = 0の証明

f=0 \Rightarrow \langle f, f \rangle = 0は自明なので \langle f, f \rangle = 0 \Rightarrow f = 0のみ示す。 まず、f, g \in \mathcal{V}について以下のコーシー・シュワルツの不等式が成り立つことを示す。

\displaystyle | \langle f, g \rangle |^2 \leq \langle f, f \rangle \langle g, g \rangle \tag{5}

これは\langle tf + g, tf + g \rangle \ \ (t \in \mathbb{R}) を考えると、

\displaystyle \begin{align} \langle tf + g, tf + g \rangle = t^2 \langle f, f \rangle + 2t \langle f, g \rangle + \langle g, g \rangle \geq 0 \end{align}

であり、これが任意の実数tに対して成立するためには判別式が0以下、すなわち4\langle f, g \rangle^2 - 4 \langle f, f \rangle \langle g, g \rangle \leq 0となる必要があるため、これを展開して式（5）が導かれる。

また、任意のx \in \mathcal{X}に対し|f(x)|^2 = 0であればf = 0である。|f(x)|^2を式（3）の再生性と式（5）を用いて展開すると、

\displaystyle | f(x) |^2 = | \langle f, k_x \rangle |^2 \leq \langle f, f \rangle \langle k_x , k_x \rangle \tag{6}

となる。そして、\langle f, f\rangle = 0なので 0 \leq |f(x)|^2 \leq 0となり f = 0が導かれる。

具体例で検証

このことと式（4）より\mathbf{a}^T \mathbf{K} \mathbf{a} = 0となるa_i, \ x_i \ \ (i =1, \cdots, n) の組み合わせではf = \sum_{i = 1}^{n} a_i k_{x_i}は零関数であり、任意のy \in \mathcal{X}についてf(y) = \sum_{i = 1}^{n} a_i k(y, x_i) = 0となる。

これをカーネル関数として2次の多項式カーネル k(\mathbf{y}, \mathbf{x}) = (1 + \mathbf{x}^T \mathbf{y})^2 \ \ (\mathbf{x}, \mathbf{y} \in \mathbb{R}^d) を採用した場合を例にとって考えてみると、\mathbf{a}^T \mathbf{K} \mathbf{a} = 0であれば任意の\mathbf{y}に対してf(\mathbf{y}) = \sum_{i = 1}^{n} a_i (1 +\mathbf{x}_i^T \mathbf{y})^2 = 0が成り立つ。

ここまで見てもまだ直感と反する気がするので、一応このことをスクリプトを書いて確かめてみる。多項式カーネルなのでd \ll nなら \rm{rank} (\mathbf{K}) \lt nとなる。そして、 \mathbf{K} の零固有値に対応する固有ベクトルを\mathbf{a}とすれば\mathbf{a}^T \mathbf{K} \mathbf{a} = 0となる。この\mathbf{a}を用いて作ったfは果たして適当な\mathbf{y}に対してf(\mathbf{y}) = 0となるのか、ということを確かめるスクリプトが以下。

import numpy as np

def func(a, X, y, k):
    res = 0
    for i in range(a.shape[0]):
        res += a[i] * k(y, X[i, :])
    return res

if __name__ == '__main__':
    np.random.seed(3407)
    d = 6
    n = 30
    X = np.random.randn(n, d)

    # 2次多項式カーネル
    k = lambda x, y: (1 + x @ y) ** 2

    K = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            K[i, j] = k(X[i, :], X[j, :])

    print ("rank(K) =", np.linalg.matrix_rank(K))

    # Kの零固有値に対応する固有ベクトルを抽出
    # このベクトルaはa @ K @ a = 0を満たす
    l, U = np.linalg.eigh(K)
    a = U[:, np.abs(l) < 1e-12]

    # 入力を適当に生成
    y = -50 + 100 * np.random.rand(d)

    for i in range(a.shape[1]):
        f = lambda y : func(a[:, i], X, y, k)
        print("f" + str(i) + "(x) =", f(y))

結果としては、

rank(K) = 28
f0(x) = -1.4438228390645236e-11
f1(x) = -1.4743761767022079e-12

演算誤差もあるので完全に0とはならないが、ほぼ0になったといえる。また乱数のseedを変えて何度か試しても同じくらいの値になるので、とりあえず納得。

おわりに

最後に参考文献。冒頭に述べた以前書いた記事でも同様の参考文献を挙げたが今回もお世話なったので改めて記載。

機械学習のための関数解析入門: ヒルベルト空間とカーネル法

機械学習のための関数解析入門: ヒルベルト空間とカーネル法

Amazon

*1:これは自分の中でまだ理解が怪しい。式（1）の定義だと有限個の和なので完備になったりしないか？n \rightarrow \inftyとして解析するときに完備性の議論が必要なる？

*2:関数同士の和は(f + g)(x) = f(x) + g(x)としていいのかという話もある。おそらく一般的にも関数の和はこのように定義するのだと思うが、少なくとも今回の場合は、再生性より(f + g)(x) = \langle f + g, k_x \rangle = \langle f, k_x \rangle + \langle g, k_x \rangle = f(x) + g(x)として得られる。

HTTFで超球面上の最適化を試してみた話

競技プログラミング 連続最適化

はじめに

AtCoder上で開催されたフューチャー社主催のヒューリスティックコンテストHACK TO THE FUTURE2022（HTTF）に参加した*1。問題等詳細は以下を参照。

atcoder.jp

この問題は大きく分けて各作業者のスキル推定と作業スケジューリングという2つの要素があるが、その中でもスキル推定の方に以前から気になってたリーマン多様体上の最適化なるものが使えそうだったので試してみたというのが今回の話。実のところ自分も多様体は何ぞやとか数学的なことはよく分かっていないが、今回行った単位超球面上の最急降下法であればそこまで難しくないし、以下の資料を参考に比較的簡単に実装できる。

http://www.orsj.or.jp/archive2/or60-9/or60_9_549.pdf

結果としては割とうまく推定できるが普通の最急降下法でやった場合とスコア的にはさほど差がなく、最終的にはスコアが高く出た普通の最急降下法を採用した（誤差範囲の気もする）。

以下、具体的行ったことを書いていく。

スキル推定問題の定式化

問題全体は上記のリンク先を見ていただくとして、ここではスキル推定問題について定式化をしておく。スキルは作業者間で独立なので、ここでは1人の作業者のスキルを推定する問題を定式化する。 作業者はK種類のスキルを有しており、そのk番目のスキルをs_k \geq 0とする*2。また、タスクの要求スキルも同様にK種類あり、タスクiの要求スキルをd_{ik} \geq 0とする。そして、実際にその作業者がタスクiを実行したときにかかった時間はt_iであり、この時間は \sum_{k = 1}^{K}\max(0, d_{ik} - s_k) = 0ならt_i = 1、 \sum_{k = 1}^{K}\max(0, d_{ik} - s_k) > 0ならこの値に-3 \sim 3の乱数を加えたもので得られる（この場合も0以下になるならt_i = 1とする）。n個タスクが作業者に割り振られ、\left\{(\mathbf{d}_1, t_1), \cdots, (\mathbf{d}_n, t_n) \right\}が得られているので、この集合から作業者のスキル\mathbf{s}を推定したい。

ということで、このスキル推定問題を2乗誤差最小化問題として以下のように定式化する。

\displaystyle {\rm minimize} \ \ E(\mathbf{s}) = {\rm minimize} \ \ \frac{1}{2n} \sum_{i = 1}^n \left\{(t_i - 1) - \sum_{k = 1}^K \max(0, d_{ik} - |s_k|) \right\}^2 \tag{1}

ここで制約付き最適化になるとややこしいので、s_k \geq 0という制約はつけず負値を取ることを許して、式の中で絶対値を取るようにしている。

通常の最急降下法でスキル推定

最終的に採用したのはこっちの方。式（1）の中にmaxやら絶対値やらが入っているので一部の人から怒られそうだが、気にせず\mathbf{s}で微分して勾配を求める*3。

\displaystyle \nabla E(\mathbf{s}) = \frac{1}{n} \mathbf{h} \odot \sum_{i = 1}^n \left\{(t_i - 1) - \sum_{k = 1}^K \max(0, d_{ik} - |s_k|) \right\} \odot \mathbf{g}_i \tag{2}

ただし、\mathbf{h} と\mathbf{g}_i はK次元のベクトルで、そのk番目の要素をそれぞれh_k, g_{ik}とすると

\displaystyle h_k = \left\{ \begin{array}{ll} 1 & (s_k > 0) \\ 0 & (s_k = 0) \\ -1 & (s_k < 0) \end{array} \right., \ \ \ g_{ik} = \left\{ \begin{array}{ll} 1 & (d_{ik} - |s_k| > 0) \\ 0 & (d_{ik} - |s_k| \leq 0) \end{array} \right.

である。また、\odotはアダマール積を表す。スキル\mathbf{s}はこの勾配を用いて、

\displaystyle \mathbf{s}^{(l + 1)} = \mathbf{s}^{(l)} - \alpha \ \nabla E(\mathbf{s}^{(l)}) \tag{3}

という更新式を適当な回数まわして推定した。本来なら収束判定を勾配のノルムとかで行うべきだが、手を抜いてループ回数を固定している。また、ステップ幅\alphaについても本来なら線形探索などで適切な値を求めるべきだが、これまた手を抜いて適当な値に固定している。

超球面上の最急降下法でスキル推定

ここからが本題。すぐ後に示すように式（1）を少し変形するとノルム1制約の最適化問題が出てくるので、そこに超球面上の最急降下法を使った。

スキルの成分分解

問題文をよく見ると真のスキルの生成方法が書かれている。これを見ると、\mathbf{s} = \frac{{\rm randdouble}(20, 60) | \mathbf{s}' |}{\parallel \mathbf{s}' \parallel_2} \ \ (\mathbf{s}' \sim \mathcal{N}(\mathbf{0}, \mathbf{I}))という式で生成されている。従って、ノルムが1のベクトル\mathbf{w}と実数q \ \ (20 \leq q \leq 60)を用いて、\mathbf{s} = q|\mathbf{w}|と書くことができる。 これを使うと式（1）は

\displaystyle \begin{align} {\rm minimize} \ \ E_{d}(\mathbf{s}) = {\rm minimize} &\ \ \frac{1}{2n} \sum_{i = 1}^n \left\{(t_i - 1) - \sum_{k = 1}^K \max(0, d_{ik} - q|w_k|) \right\}^2\\ {\rm s.t.} &\ \ \sum_{k=1}^K w_k^2 = 1 \\ & 20 \leq q \leq 60 \end{align} \tag{4}

というqと\mathbf{w}に関する最適化に変換される。 この式（4）の最適化のため、\mathbf{w}を固定してqを最適化する処理と、qを固定して\mathbf{w}を最適化する処理を交互に繰り返す方法を取った。

qの最適化

\mathbf{w}を固定すると、式（4）はqについて凸な関数になる。また、qは1次元の実数。ということで三分探索で最適解を求めた。ただ、タイムリミットの関係上三分探索のループは5回にしてるので、かなり粗い探索になっている。

\mathbf{w}の最適化

\mathbf{w}はノルムが1という制約条件があり、これは\mathbf{w}が単位超球面上に存在する制約だといえる。 このような制約条件つき最適化問題に対しては様々なアルゴリズムが提案されているが*4、今回は冒頭で挙げた資料の通り超球面上での制約条件なしの最適化問題とみなす手法を実装した。

まずは一応式（4）の\mathbf{w}に関する通常の勾配を書いておくと、

\displaystyle \begin{align} \nabla E_{d}(\mathbf{w}) = \frac{q}{n} \mathbf{h}' \odot \sum_{i = 1}^n \left\{(t_i - 1) - \sum_{k = 1}^K \max(0, d_{ik} - q|w_k|) \right\} \odot \mathbf{g}_i' \end{align} \tag{5}

ただし、

\displaystyle h_k' = \left\{ \begin{array}{ll} 1 & (w_k > 0) \\ 0 & (w_k = 0) \\ -1 & (w_k < 0) \end{array} \right., \ \ \ g_{ik}' = \left\{ \begin{array}{ll} 1 & (d_{ik} - q|w_k| > 0) \\ 0 & (d_{ik} - q|w_k| \leq 0) \end{array} \right.

となる。この勾配方向に式（3）のように動かしても、超球面のことは全く考慮されていないので明後日の方向に行ってしまう可能性がある。ということで、

  • なるべく超球面から離れないよう現在点における接平面に沿った方向に動かす
  • それでも超球面からはみ出るので超球面上に戻す（レトラクション）

ということを繰り返すのが当該のアルゴリズム。以下もう少し具体的に書く。

動かす方向

今超球面上の点\mathbf{w}^{(l)}にあって、そこから動かす方向を決める。\mathbf{w}^{(l)}における接平面の法線ベクトルは\mathbf{w}^{(l)}である。この法線ベクトルが張る空間の直交補空間は接平面に平行な部分空間（接ベクトル空間）であり、\nabla E_{d}(\mathbf{w}) のその部分空間へ直交射影したベクトルを\mathbf{d}とすると、

\displaystyle \mathbf{d} = \nabla E_d(\mathbf{w}^{(l)}) - \langle \mathbf{w}^{(l)}, \nabla E_d(\mathbf{w}^{(l)}) \rangle \mathbf{w}^{(l)} \tag{6}

となる。これを用いて、

\displaystyle \mathbf{w'}^{(l)} = \mathbf{w}^{(l)} - \alpha \mathbf{d} \tag{7}

にまずは動く。こうすれば、とりあえずは接平面に沿った移動になるので、超球面から離れないという観点で言うと少なくとも-\nabla E_d(\mathbf{w}^{(l)})方向に移動するよりはマシになるはず、というのが直感的な理解。なお、\alphaはステップ幅でこちらも本来は（曲線）探索するべきだが実装では手を抜いて固定値にした。

超球面へのレトラクション

式（7）で得た\mathbf{w'}^{(l)} は、それでもやはり超球面からははみ出ている。従ってこれを超球面上に戻す必要がある。この操作をレトラクションと呼ぶ。これは単純に\mathbf{w'}^{(l)} のノルムで割ればよい。よって、

\displaystyle \mathbf{w}^{(l + 1)} = \frac{\mathbf{w'}^{(l)}}{\parallel \mathbf{w'}^{(l)} \parallel_2} \tag{8}

として\mathbf{w}を更新できる。

結果

seed=0の時のビジュアライザ動画を貼りたかったけどはてなブログに動画がアップロードできなかったので、最後の状態を画面キャプチャした画像を貼る。

うーん、どっちがいいとも悪いとも言いにくい。結局のところそこまで変わらないので、処理時間的にも短い通常の最急降下法を最終的に採用したわけだが。気持ちとしては後者の方がテクくて好きだったので何とかして採用したかった。。。

大半の人はスキル推定に焼きなましやMCMCを使っているので、実際のところそちらの方がいいのかもしれない。

あと、これを書いてて思ったが、振られたタスクの要求スキルが小さいと各スキルは過小評価されるので、最初のうちは満遍なくいろいろなタスクを振ってスキル推定の精度を向上させた方がよかったかもしれない。

ちなみに

作業スケジューリングについては、重み付き二分マッチング問題を誰かがタスクを終える度に実行して作業の割り当てを行った。この際、辺の重みは予測される作業時間に加えて、なるべくそのタスクを終えることで次に実行できるタスクが増えるように調整している。 スキル推定はそこまで重要じゃない気もするので、作業スケジューリングの方にもっと時間を割いて考えるべきだったかもしれない。

所謂ヒューリスティックコンテスト形式のコンテストはAtCoder Heuristic Contestが始まってからまともに参加するようになったが、今回のように連続最適化の手法も使えるコンテストがたまにあって中々楽しい。

*1:最終結果はスコア5,664,033で129位。

*2:スキルとして取りうる値は整数だが、ここでは実数値として考え、後から四捨五入する方法を取った。

*3:一応、最初は絶対値ではなく2乗で0以上を保証するなど微分可能な式を立ててやってたが結果的には今の定式化の方がよかった。今や深層学習なんかでReLU等をガンガン微分してるのでまぁ。。。

*4:以前書いた内点法とかyamagensakam.hatenablog.com

超平面と点の距離の求め方を少し抽象的に書いてみる

関数解析 線形代数

はじめに

以下の書籍を読み関数解析のお気持ちを理解しつつあるので、今回は備忘録がてら超平面と点の距離を抽象ヒルベルト空間の観点で記述してみる。

なお話を簡単にするために、超平面は必ず原点を通るものとする。

機械学習のための関数解析入門: ヒルベルト空間とカーネル法

機械学習のための関数解析入門: ヒルベルト空間とカーネル法

Amazon

n次元実ベクトル空間\mathbb{R}^{n}における超平面と点の距離

まずは通常の線形代数の範疇で超平面と点の距離を記述する。超平面上の点を\mathbf{x} \in \mathbb{R}^{n}とすると、超平面は原点を通るとしているのでその方程式は適当な法線ベクトル\mathbf{v} \in \mathbb{R}^{n}を用いて、

\displaystyle \mathbf{v}^T \mathbf{x} = 0 \tag{1}

と表される。この超平面と任意の\mathbf{x}_0 \in \mathbb{R}^nの距離だが、これは\mathbf{x}_0を超平面上に射影した点\bar{\mathbf{x}}_0と\mathbf{x}_0の距離にあたる。超平面上の部分空間と\mathbf{v}が張る空間は直交しているため\bar{\mathbf{x}}_0と\mathbf{x}_0の関係は\alpha \in \mathbb{R}を用いて

\displaystyle \mathbf{x}_0 = \bar{\mathbf{x}}_0 + \alpha \mathbf{v} \tag{2}

が成り立つ。ということで超平面と\mathbf{x}_0の距離dは

\displaystyle d = \parallel \mathbf{x}_0 - \bar{\mathbf{x}}_0 \parallel_2 = \parallel \bar{\mathbf{x}}_0 + \alpha \mathbf{v} - \bar{\mathbf{x}}_0 \parallel_2 = |\alpha| \parallel \mathbf{v} \parallel_2 \tag{3}

となる。\alphaだが、これは式（2）の両辺に\mathbf{v}^Tをかけて、\bar{\mathbf{x}}_0 は超平面上の点なので\mathbf{v}^T\bar{\mathbf{x}}_0 = 0であることに注意すると、

\displaystyle \alpha = \frac{\mathbf{v}^T \mathbf{x}_0}{\mathbf{v}^T\mathbf{v}} = \frac{\mathbf{v}^T \mathbf{x}_0}{\parallel \mathbf{v} \parallel_2^2} \tag{4}

となるので、後はこれを式（3）に代入して

\displaystyle d =\frac{|\mathbf{v}^T \mathbf{x}_0|}{\parallel \mathbf{v} \parallel_2} = \frac{|\mathbf{v}^T \mathbf{x}_0|}{\sqrt{\mathbf{v}^T \mathbf{v}}} \tag{5}

を得る。

ヒルベルト空間の定義と諸定理

次にこれをヒルベルト空間上の議論へ抽象化してみる。が、そもそもヒルベルト空間とは？その超平面とは？距離とは？という話なので、本題に入る前に今回の議論に必要な定義および定理を紹介する。

ヒルベルト空間の定義

まずヒルベルト空間\mathcal{H}の定義を述べる。

  • ベクトル空間であること
  • 任意の元\mathbf{x}, \mathbf{y} \in \mathcal{H}に対し、内積\langle \mathbf{x}, \mathbf{y} \rangle_{\mathcal{H}}が定義されていること
    • この内積は内積の定義さえ満たしていればよい
    • この内積によってノルムが \parallel \mathbf{x} \parallel_{\mathcal{H}} = \sqrt{\langle \mathbf{x}, \mathbf{x} \rangle_{\mathcal{H}}}として定まる
    • さらにこのノルムより2点\mathbf{x}, \mathbf{y}の距離関数がd(\mathbf{x}, \mathbf{y} ) =\parallel\mathbf{x} - \mathbf{y} \parallel_{\mathcal{H}} として定まる
  • 上述の距離関数に関して完備であること
    • すなわち、任意のコーシー列（\parallel \mathbf{x}_m - \mathbf{x}_n \parallel_{\mathcal{H}} \rightarrow 0 \ (m, n \rightarrow \infty)）が\mathcal{H}の元に収束すること

なお、完備でない内積空間も完備化という操作によって元の追加と内積の拡張を行うことでヒルベルト空間に変換することができる。

射影定理

ヒルベルト空間には完備性があるおかげて射影定理が成り立ち、これにより通常の線形代数と同じような議論を行うことができる。射影定理とは以下のような定理である。

\mathcal{M}を\mathcal{H}の閉部分空間（任意の\alpha, \beta \in \mathbb{R}, \ \ \mathbf{x}_1, \mathbf{x}_2 \in \mathcal{M}に対して\alpha \mathbf{x}_1 + \beta \mathbf{x}_2 \in \mathcal{M}が成り立つ閉集合）、\mathcal{M}^{\perp}は\mathcal{M}と直交するベクトル全体とする。このとき、\mathbf{x} \in \mathcal{H}に対し、\mathbf{x} = \mathbf{y} + \mathbf{z} \ \ (\mathbf{y} \in \mathcal{M}, \mathbf{z} \in \mathcal{M}^{\perp}) となる分解が一意に存在する。また、この時\mathbf{y} は\mathcal{M}上で最も\mathbf{x}との距離が最小になる唯一の点である*1。

この\mathbf{y}を\mathbf{x}の\mathcal{M}上への直交射影と呼びその作用素をP_{\mathcal{M}}と表す。同様に \mathcal{M}^{\perp}上への直交射影作用素をP_{\mathcal{M}^{\perp}}と表す。従って、\mathbf{x} = P_{\mathcal{M}}\mathbf{x}+ P_{\mathcal{M}^{\perp}}\mathbf{x} と表すことができる。また、このような分解を直交分解と呼ぶ。

証明は書籍等に譲るが、完備性が射影定理にどう絡んでくるかだけ見ておくと、d = \inf_{\mathbf{y} \in \mathcal{M}} \parallel \mathbf{x} - \mathbf{y} \parallel_{\mathcal{H}}となる \mathcal{H}の元\mathbf{y} が存在し、d = \min_{\mathbf{y} \in \mathcal{M}} \parallel \mathbf{x} - \mathbf{y} \parallel_{\mathcal{H}}となることを保証するために必要になる。具体的には、任意のn \in \mathbb{N}に対し、 d^{2} \leq \parallel \mathbf{x} - \mathbf{y}_n \parallel_{\mathcal{H}}^2 \lt d^{2} + \frac{1}{n^{2}} となる\mathbf{y}_n \in \mathcal{M}が存在し、また点列(\mathbf{y}_n)_{n=1}^\inftyがコーシー列となることから\parallel \mathbf{y}_n - \mathbf{y}_{*} \parallel = 0 \ \ (n \rightarrow \infty)である\mathbf{y}_{*} \in \mathcal{M}の存在を保証するために完備性が必要となる。

線形汎関数とリースの表現定理

ヒルベルト空間\mathcal{H}上に定義された写像\varphi: \mathcal{H} \rightarrow \mathbb{R}が\varphi(\alpha \mathbf{x} + \beta \mathbf{y}) = \alpha \varphi(\mathbf{x}) + \beta \varphi( \mathbf{y}) \ \ (\alpha, \beta \in \mathbb{R}, \ \ \mathbf{x}, \mathbf{y} \in \mathcal{H})を満たすとき\varphiを\mathcal{H}上の線形汎関数と呼ぶ。また、\parallel \varphi \parallel = \sup_{ \parallel \mathbf{x} \parallel_{\mathcal{H}} \leq 1} |\varphi(\mathbf{x})|と定め、この量が有限である場合\varphiは有界であるという。加えて、\varphiは有界であることと\varphiが連続であることは同値である。

リースの表現定理はこの有界な線形汎関数\varphiについて以下のことを主張する定理である。

\displaystyle \varphi(\mathbf{x}) = \langle \mathbf{x}, \mathbf{v} \rangle_{\mathcal{H}} \tag{6}

となる\mathbf{v} \in \mathcal{H}が一意に存在する。

これまた詳細な証明は書籍に譲り、ここでは存在性についてのみ述べる。存在性の証明には以下の補題が必要になる。

  1. 零空間\mathcal{M} = \ker \varphi=\left\{\mathbf{x} \in \mathcal{H} : \varphi(\mathbf{x}) =0 \right\} は閉部分空間
  2. \dim \mathcal{M}^{\perp} = 1

1.に関しては部分空間であることは\varphiの線形性から明らかで、閉集合であることは\varphiの連続性と\left\{0\right\}が\mathbb{R}の閉集合であることから逆像\mathcal{M} = \varphi^{-1}(\left\{0\right\})によって示される。

2.は\dim \mathcal{M}^{\perp} \gt 1である場合は \mathbf{y} \in \mathcal{M}^{\perp} に対して直交な \mathbf{z} \in \mathcal{M}^{\perp} を取れるが、\varphi(\mathbf{y}) \neq 0, \varphi(\mathbf{z}) \neq 0であるため適当な定数を選んで\varphi(\mathbf{y}) = \lambda \varphi(\mathbf{z})が成り立つ。線形性より\varphi(\mathbf{y}- \lambda \mathbf{z}) = 0となるため、\mathbf{y}- \lambda \mathbf{z} \in \mathcal{M}である。従って、\langle \mathbf{y}, \mathbf{y}- \lambda \mathbf{z} \rangle_{\mathcal{H}} = 0であり、 \mathbf{z} は\mathbf{y}と直交になるように選んでいることから\langle \mathbf{y}, \mathbf{y}- \lambda \mathbf{z} \rangle_{\mathcal{H}} = \langle \mathbf{y}, \mathbf{y} \rangle_{\mathcal{H}} - \langle \mathbf{y}, \lambda \mathbf{z} \rangle_{\mathcal{H}} =\parallel \mathbf{y} \parallel_{\mathcal{H}}^2 = 0が得られるため\mathbf{y} = \mathbf{0}である。よって、\mathcal{M}^{\perp}の中で\mathbf{z} \in \mathcal{M}^{\perp}と直交するのは\mathbf{0}のみということは\dim \mathcal{M}^{\perp} = 1だといえる。

さて、この2によって\mathcal{M}^{\perp}上の任意の点は\mathbf{z} \in \mathcal{M}^{\perp} \ \ (\parallel \mathbf{z} \parallel_{\mathcal{H}}=1)を用いて、\alpha \mathbf{z} \ (\alpha \in \mathbb{R})と表現できる。これと1によって任意の\mathbf{x} \in \mathcal{H}は\mathbf{x} = P_{\mathcal{M}} \mathbf{x} + \alpha \mathbf{z}と直交分解でき、加えて、

\displaystyle \begin{align} \varphi(\mathbf{x}) &= \varphi(P_{\mathcal{M}} \mathbf{x} + \alpha \mathbf{z}) \\ &= \varphi(P_{\mathcal{M}} \mathbf{x}) + \alpha \varphi(\mathbf{z})\\ &= \alpha \varphi(\mathbf{z}) \\ &= \alpha \varphi(\mathbf{z}) \left(\langle \frac{1}{\alpha} P_{\mathcal{M}} \mathbf{x}, \mathbf{z} \rangle_{\mathcal{H}} + \langle \mathbf{z}, \mathbf{z} \rangle_{\mathcal{H}} \right) \\ &= \alpha \varphi(\mathbf{z}) \langle \frac{1}{\alpha} P_{\mathcal{M}} \mathbf{x} + \mathbf{z}, \mathbf{z} \rangle_{\mathcal{H}}\\ &= \langle P_{\mathcal{M}} \mathbf{x} + \alpha \mathbf{z}, \varphi(\mathbf{z}) \mathbf{z} \rangle_{\mathcal{H}} \\ &= \langle \mathbf{x} , \varphi(\mathbf{z}) \mathbf{z} \rangle_{\mathcal{H}} \\ \end{align} \tag{7}

と展開できる。そして、有界性より\varphi(\mathbf{z})が有限値になることが保証されているので、式（7）の最後の式において\mathbf{v} =\varphi(\mathbf{z}) \mathbf{z} とすれば式（6）が得られリースの表現定理における存在性が示される。

ヒルベルト空間における超平面と点の距離

ようやく本題だが、ここまでくると最初に述べた実ベクトル空間上の議論とあまり変わらない。以下、前章の結果を用いてヒルベルト空間上の超平面と点の距離の導出を行い、その後に具体的な元と内積が定義されたヒルベルト空間に対し導出した結果を適用した例を示す。

導出

これまでの議論によりヒルベルト空間における超平面とは\mathcal{H}上の有界線形汎関数\varphiにより、

\displaystyle \varphi(\mathbf{x}) = 0 \ \ (\mathbf{x} \in \mathcal{H}) \tag{8}

となる\mathbf{x}の集合のことだといえる。つまりそれは零空間\mathcal{M} = \ker \varphi であり、それと任意の\mathbf{x}_0 \in \mathcal{H}との距離は

\displaystyle d = \parallel \mathbf{x}_0 - P_{\mathcal{M}} \mathbf{x}_0\parallel_\mathcal{H}=\parallel P_{\mathcal{M}} \mathbf{x}_0 + P_{\mathcal{M}^{\perp}} \mathbf{x}_0 - P_{\mathcal{M}} \mathbf{x}_0\parallel_\mathcal{H}= \parallel P_{\mathcal{M}^{\perp}} \mathbf{x}_0\parallel_\mathcal{H} \tag{9}

だといえる。また、前章の議論より\varphi(\mathbf{x}) = \langle \mathbf{x}, \mathbf{v} \rangle_{\mathcal{H}} と書け、P_{\mathcal{M}^{\perp}} \mathbf{x}_0= \alpha \mathbf{v}となる。よって、

\displaystyle d = \parallel P_{\mathcal{M}^{\perp}} \mathbf{x}_0\parallel_{\mathcal{H}} = | \alpha | \parallel \mathbf{v} \parallel_{\mathcal{H}} \tag{10}

と式（3）とノルム以外同じ形になる。さらに、

\displaystyle \varphi(\mathbf{x}_0) = \varphi(P_{\mathcal{M}} \mathbf{x}_0 + \alpha \mathbf{v}) = \varphi(P_{\mathcal{M}} \mathbf{x}_0 )+ \alpha \varphi(\mathbf{v}) \tag{11}

であるが、\mathcal{M}が\varphiの零空間であることを思い出すと \varphi(P_{\mathcal{M}} \mathbf{x}_0 ) = 0なので、

\displaystyle \alpha = \frac{\varphi(\mathbf{x}_0) }{\varphi(\mathbf{v})} = \frac{\varphi(\mathbf{x}_0) }{\parallel \mathbf{v} \parallel_{\mathcal H}^2} \tag{12}

となる。これを式（10）に代入することで、

\displaystyle d = \frac{|\varphi(\mathbf{x}_0) | }{ \parallel \mathbf{v} \parallel_{\mathcal{H}} } = \frac{|\langle \mathbf{x}_0, \mathbf{v}\rangle_{\mathcal{H}} |}{\parallel \mathbf{v} \parallel_{\mathcal{H}} } = \frac{|\langle \mathbf{x}_0, \mathbf{v}\rangle_{\mathcal{H}} |}{\sqrt{\langle \mathbf{v}, \mathbf{v} \rangle_{\mathcal{H}}} } \tag{13}

が得られる。

具体例

式（13）で超平面と点の距離を求める例として、最初に述べたn次元実ベクトル空間の例とカーネル関数と呼ばれる関数から導かれるヒルベルト空間の例を挙げる。

n次元実ベクトル空間再考

式（13）を使って、n次元実ベクトルにおける超平面と点の距離を求める。まず、\mathbf{x}, \mathbf{y} \in \mathbb{R}^nの内積は\langle \mathbf{x}, \mathbf{y} \rangle_{\mathbb{R}^n} = \mathbf{x}^T \mathbf{y}であり、超平面は\varphi(\mathbf{x}) = \langle \mathbf{x}, \mathbf{v} \rangle_{\mathbb{R}^n} = \mathbf{x}^T \mathbf{v} = \mathbf{v}^T \mathbf{x} = 0なので、そのまま式（13）に当てはめると式（5）が得られる。

カーネル関数によって定義されるヒルベルト空間

まず、2変数x_1, x_2 \in \mathcal{X}上の半正定値対称関数k(x_1, x_2)をカーネル関数と呼ぶ。また、カーネル関数の2変数のうち一方を固定した1変数関数をk(\cdot, x) = k_{x}(\cdot) = k_{x}と表記する。今、この1変数関数の有限個の線形和f = \sum_{i} {a_i} k_{x_i} \ (a_i \in \mathbb{R}) が要素である集合\mathcal{V}を考える。k_{x} も\mathcal{V}の要素である。そして、この集合\mathcal{V}に対して内積を次のように定義する。

f = \sum_{i} {a_i} k_{x_i}, \ \ g = \sum_{j} {b_j} k_{y_j} \in \mathcal{V}に対し、

\displaystyle \langle f, g \rangle_{\mathcal{V}} = \sum_{i} \sum_{j} {a_i} {b_j} k(y_j, x_i) \tag{14}

カーネル関数の半正定値対称性から内積の定義を満たす*2。 このようにして定義した内積空間\mathcal{V}を完備化することによりヒルベルト空間\mathcal{H}_kを得る。これでこの節で考える空間を構築できた。なお、\mathcal{V} \subset \mathcal{H}_kであり、f, g \in \mathcal{V}であれば\langle f, g \rangle_{\mathcal{V}} = \langle f, g \rangle_{\mathcal{H}_k}である。

それでは、\varphi(f) = \langle f, \sum_{i} {c_i} k_{z_i} \rangle_{\mathcal{H}_k} = 0\ \ (c_i \in \mathbb{R}, z_i \in \mathcal{X})となる超平面とf = k_{x} の距離を式（13）を使って計算してみる。そのまま式（13）に代入することで、

\displaystyle d = \frac{|\langle k_x, \sum_{i} {c_i} k_{z_i} \rangle_{\mathcal{H}_k} |}{\parallel \sum_{i} {c_i} k_{z_i} \parallel_{\mathcal{H}_k} } = \frac{ |\sum_{i} {c_i} k( z_i, x) |}{\sqrt{\sum_i \sum_j c_i c_j k(z_i, z_j)}} \tag{15}

が得られる。

もう少しこの空間について詳しく述べると、このように定めた\mathcal{H}_kは再生核ヒルベルト空間と呼ばれ、\mathcal{X}上の関数f \in \mathcal{H}_kにx \in \mathcal{X}を代入する操作がf(x) = \langle f, k_{x} \rangle_{\mathcal{H}_k}と内積の形で表される「再生性」という性質がある。この性質によってリプレゼンター定理などカーネル法における重要な定理が導かれる。また、再生性は完備化によって得られる元についても成り立つので任意のf \in \mathcal{H}_kと上述の超平面の距離は

\displaystyle d = \frac{|\langle f, \sum_{i} {c_i} k_{z_i} \rangle_{\mathcal{H}_k} |}{\parallel \sum_{i} {c_i} k_{z_i} \parallel_{\mathcal{H}_k} } = \frac{|\sum_{i} {c_i}\langle f, k_{z_i} \rangle_{\mathcal{H}_k} |}{\parallel \sum_{i} {c_i} k_{z_i} \parallel_{\mathcal{H}_k} } = \frac{ |\sum_{i} {c_i} f(z_i) |}{\sqrt{\sum_i \sum_j c_i c_j k(z_i, z_j)}} \tag{16}

になる。

最後に

今回、抽象的なヒルベルト空間における超平面と点の距離導くため、必要な定義と諸定理について記述した。もちろん今回挙げた例以外にもヒルベルト空間上であれば同じ議論が成り立つ。例えば、閉区間[a, b]上の連続関数の集合C[a, b]についても内積を定義し、それを完備化すればヒルベルト空間になるため今回の議論が成り立つ*3。

本記事では超平面と点の距離をトピックとして取り上げたが、それに限らず一度抽象的な世界で議論してその中で成り立つ法則を見つけると、具象の世界ではそれを適用するだけでいいのでかなり有難い。関数解析はその有難みを特に感じられる分野な気がする。

参考文献

*1:この最良近似性は閉部分空間よりも広い集合である閉凸集合に対して成り立つ。

*2:『半』正定値だと内積の定義のうち\langle f, f \rangle_{\mathcal{V}} = 0 \Leftrightarrow f = 0 が成り立たない気もしたが、これは後述の再生性とコーシー・シュワルツの不等式を使って成り立つことが示せる（詳しくは書籍を参照）。

*3:これを具体例の1つとして挙げようかとも思ったがルベーグ積分の知識が必要になりそうで今回は力が及ばなかった。