甲斐性のない男が機械学習とか最適化とかを綴るブログ

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

近接勾配法応用編その3~トレースノルム正則化項付きロジスティック回帰による行列データの分類~

はじめに

前回の記事では、ベクトルデータの分類問題に対するスパースなロジスティック回帰を説明しましたが、今回はその拡張となる行列データに対するロジスティック回帰を見ていきます。この行列データ分類問題においても正則化項は重要な役割を持ちますが、中でもトレースノルムの正則化項を用いることで低ランクな解が得られ、多次元の時系列データの分類などに有効であることを見ていきます。

行列データの分類

まず、学習データ集合\left\{{\bf X}_i, y_i \right\}(i=1⋯N)があるとします。ここで、{\bf X}D \times Tの行列、 yy \in \left\{0,1 \right\}のラベルです。この学習データ集合から係数行列{\bf W}とバイアス項bを学習し、以下に示すような線形モデルで分類することが今回のテーマです。

\displaystyle a({\bf X}) = {\rm Tr}({\bf X}^T {\bf W}) + b \tag{1}

行列データ分類問題の特徴

では、行列データを直接扱うことにどのようなメリットがあるのでしょうか。確かに、行列データ{\bf X}をベクトルに変換して(列ベクトルを縦に連結するなどして)、学習・分類するという方法も考えられます。しかしベクトルにしてしまうと行列構造がなくなり、行と列の意味合いが異なる場合に重要な情報が失われる可能性があります。
例えば、複数のセンサーから取得される多次元の時系列データの場合、行がセンサー、列が時間に対応する行列データとしてとらえることができます。このようなデータを連結してベクトルとして扱ってしまうと、空間方向と時間方向の特徴を一緒くたにしてしまい、重要な情報が失われてしまう可能性があります。特にこのような複数センサーのデータの場合、実際に観測できるデータは{\bf X} = {\bf A}{\bf Z}のように、潜在変数(信号源){\bf Z}の空間的な線形結合で表されることがしばしばあり、クラス分類において有用な情報はこの潜在変数に隠れていることも多いです(下図参照)。


f:id:YamagenSakam:20180613220313p:plain

一方で、行列データとして直接扱うことの意味は、{\bf W} =\sum_{i = 1}^{r} \sigma_i {\bf u}_i {{\bf v}_i}^T 特異値分解して式(1)を以下のように変形すると見えてきます。

\displaystyle a({\bf X}) = \sum_{i=1}^{r}(\sigma_i {{\bf u}_i}^T {\bf X} {\bf v}_i) + b \tag{2}

このようにすると、{\bf u}は空間方向の特徴を抽出する空間フィルタ、{\bf v}は時間方向の特徴を抽出する時間フィルタとして考えることができ、時空間両方の特徴を捉えられることがわかります。上述のように観測されるデータが潜在変数の空間的な線形結合である場合にも、観測データに空間フィルタをかけることでクラス分類に有用な潜在的な情報を抽出した上で時間方向に係数をかけるので、時空間両方の特徴をとらえることができます。

行列データロジスティック回帰の目的関数

行列データの分類問題においてもロジスティック回帰は使えます。具体的には下記の誤差関数を最小化する{\bf W}bを学習データ集合から求めればよいです。
\displaystyle E({\bf W}, b) =-\sum_{i=1}^{N} \left\{y_i \ln \sigma({\rm Tr}({{\bf X}_i}^T {\bf W}) + b ) + (1 - y_i) \ln (1 - \sigma({\rm Tr}({{\bf X}_i}^T {\bf W}) + b )) \right\} \tag{3}

トレースノルム正則化による低ランク化

さて、式(2)にあるように係数行列{\bf W}のランクrが高いほど、かけ合わせる時空間フィルタの個数は増えます。もし、{\bf W}がフルランクならr = \min(D,T)個の時空間フィルタをかけて重み付き和を取ることになります。しかし、上述のような潜在変数の中でもクラス分類に有用な情報は、(空間的に)ごく一部の成分でありそれ以外はノイズというケースも多くあります。例えば、上の図の例だと、3つの信号源のうちクラス分類に有用な成分は1つです。この場合に{\bf W}がフルランクだと必要以上に複雑なモデルとなってしまい過学習となる恐れがあります。
ということで、{\bf W}を低ランクにしたいのですが、そのためにはこちらの記事で少し述べたトレースノルムの正則化 \|{\bf W} \|_{\rm trace}= \sqrt{{\rm Tr}({\bf W}^T {\bf W})}をつけると特異値がスパースになり、結果低ランクな解が得られます。具体的な式としては、
\displaystyle F({\bf W}, b) =E({\bf W}, b)  +  \|{\bf W} \|_{\rm trace} \tag{4}
であり、これを近接勾配法で解いていきます。なお、今回はバイアス項b正則化の中に含めていないことに注意してください*1b正則化項に含まれないので、{\bf W}は近接勾配法の更新式、bは通常の勾配法の更新式で更新していきます。それぞれの更新式で必要な式(3)の{\bf W}による微分および式(4)のbによる微分は、
\displaystyle \nabla E({\bf W}) = \sum_{i=1}^{N} (\sigma({\rm Tr}({{\bf X}_i}^T {\bf W})  + b) - y_i) {\bf X}_i \tag{5}
\displaystyle \frac{ dF}{db} = \sum_{i=1}^{N} (\sigma({\rm Tr}({{\bf X}_i}^T {\bf W})  + b) - y_i)  \tag{6}
となるので、後はこれを元に{\bf W}bの更新を交互に行えばよいです。

実装と実験

アルゴリズムの実装

今回の行列データロジスティック回帰は下記のように実装しました。基本的に今までと同じですが、proximal operatorはトレースノルム正則化のproximal operator計算を行うようにしており、近接勾配法の関数はバイアス項の更新(勾配法計算)も行うように変更しています。

import numpy as np
from numpy.random import *
import matplotlib.pyplot as plt
import sys

#ソフトしきい値作用素の計算
def soft_thresh(b, lam):
    x_hat = np.zeros(b.shape[0])
    x_hat[b >= lam] = b[b >= lam] - lam
    x_hat[b <= -lam] = b[b <= -lam] + lam
    return x_hat

#トレースノルム正則化のproximal operatorの計算
def prox_trace_norm(V, gamma, lam):
    [L, sig, R] = np.linalg.svd(V)
    sig_ = soft_thresh(sig, lam * gamma)
    if L.shape[1] > sig.shape[0]:
        L = L[:,:sig.shape[0]]
    if R.shape[1] > sig.shape[0]:
        R = R[:sig.shape[0], :]
    return np.dot(np.dot(L, np.diag(sig_)), R)

#近接勾配法とバイアス項の勾配法計算
def proximal_gradient_with_bias(grad_f, grad_b, prox, gamma, objective, init_x, init_b, tol = 1e-9):
    x = init_x
    b = init_b
    result = sys.maxint
    while 1:
        x_new = prox(x - gamma * grad_f(x, b), gamma)
        b_new = b - gamma * grad_b(x, b)
        result_new = objective(x_new, b_new)
        if result_new < tol or (np.abs(result - result_new)/np.abs(result) < tol) == True :
            break;
        x = x_new
        b = b_new
        result = result_new
    return x_new, b_new, result_new

#トレースノルム正則化項付きロジスティック回帰
#X:学習用データベクトルを列ベクトルとして並べた行列
#y:ラベルを並べたベクトル
def trace_norm_LogisticRegression(X, y, lam):
    sigma = lambda a : 1.0/(1.0+np.exp(-a))
    p=lambda Z, W, b:sigma(np.trace(np.dot(Z.T, W).T) + b) #メモ化したい・・・
    objective = lambda W, b: - np.sum( y * np.log(p(X, W, b)) +
                        (1 - y) *np.log(( 1 - p(X, W, b)))) + lam * np.linalg.norm(W, 'nuc')
    grad_E = lambda W, b: np.dot(X, p(X, W, b) - y) #勾配
    grad_b = lambda W, b: np.sum(p(X, W, b) - y)

    (u, l, v) = np.linalg.svd(X)
    Gamma =1.0/max(np.average(l.real*l.real,0)) 
    prox = lambda V, gamma:prox_trace_norm(V, gamma, lam) #proximal operator
    W_init = np.zeros((X.shape[0], X.shape[1]))
    b_init = 0
    (W, b, result) = proximal_gradient_with_bias(grad_E, grad_b, prox,
                                                    Gamma, objective, W_init, b_init, 1e-4)
    return W, b, result

なお、停止判定の方法については収束性を考慮して元の方法に戻しています。また、p(Z, W, b)の計算がかなり重いですが1回の更新で4回計算しているので、pをメモ化するなどの工夫を行えばもう少し速くできます(pの引数ZとWはnumpyのarray型なので簡単に辞書のキーにできず断念)。
あと、学習のステップ幅は今回かなりテキトーに決めた実装になっているので、その辺りはもう少し検討が必要です。

実験用疑似データ

実験用疑似のデータは下記のコードで生成しました。

    #学習データ
    #30×100の行列データ負例500個、正例500個
    #信号源のデータZ
    Z_train1=np.random.randn(30,100,500)
    Z_train2=np.random.randn(30,100,500)
    #30個の信号源のうち2個のみクラス分類に有用とする
    Z_train2[0, 20:50, :] += 0.2
    Z_train2[29, 80:90, :] -= 0.5
    Z_train = concatenate((Z_train1,Z_train2),2)

    #信号源データZを観測データに変換する行列A
    A = np.exp(-np.abs(np.tile(np.linspace(0, 1,30), (30,1)) -
            np.tile(np.linspace(0, 1,30), (30,1)).T)/0.1) + np.random.randn(30, 30)

    #実際に観測されるデータX
    X_train = np.dot(Z_train.T, A).T

    y_train1 = ones((500,1))
    y_train2 = zeros((500,1))
    y_train = concatenate((y_train1,y_train2),0).squeeze()

    #テストデータ
    #学習データと同じく30×100の行列データ負例500個、正例500個
    Z_test1=np.random.randn(30,100,500)
    Z_test2=np.random.randn(30,100,500)
    Z_test2[0, 20:50, :] += 0.2
    Z_test2[29, 80:90, :] -= 0.5
    Z_test = concatenate((Z_test1,Z_test2),2)

    X_test = np.dot(Z_test.T, A).T

    y_test1=ones((500,1))
    y_test2=zeros((500,1))
    y_test=concatenate((y_test1,y_test2),0).squeeze()

学習、テストに使用するデータはX_train、X_testですが、上述の信号源からの線形結合を模擬するため、まず潜在変数(信号源)Z_train、Z_testの生成を行いました。この信号源のデータは30成分ありますが、0番目の成分と29番目の成分のみがクラス分類に有用となるようなデータとしています。これを行列Aで線形結合して観測データXを生成しています。なお、行列Aは近くの信号程大きな重み、離れるほど小さな重みで足し合わせるような線形結合としており、更にこのAにもランダム成分を加えています。以下、正例・負例におけるXとZの平均波形を表示します。

  • y=1に対応する観測データ集合X_train1の平均波形

f:id:YamagenSakam:20180613230346p:plain

  • y=0に対応する観測データ集合X_train2の平均波形

f:id:YamagenSakam:20180613230616p:plain

  • y=1に対応する信号源データ集合Z_train1の平均波形

f:id:YamagenSakam:20180613230428p:plain

  • y=0に対応する信号源データ集合Z_train2の平均波形

f:id:YamagenSakam:20180613230713p:plain


やはり、観測データはy=0クラスとy=1クラスの平均データを見比べても、ノイズにまみれて差異があまり分からなくなっていますが、信号源データの場合は、2信号のみy=0クラスとy=1クラスで差異があることがはっきりわかります。

実験結果

このように生成した疑似データに対し、以下のコードのように\lambda=0, 50, 100 \cdots, 2000と変化させて、その時の推定精度とWのランクを求めてみました。

    num = 41
    scale = 50
    rank_of_W  = np.zeros(num)
    acc_arr = np.zeros(num)
    for i in range(0, num):
        (W, b, r) = trace_norm_LogisticRegression(X_train, y_train.squeeze(), i * scale)
        y_result = np.trace(np.dot(X_test.T, W).T) + b
        y_result[y_result>=0] = 1
        y_result[y_result<0] = 0
        rank_of_W[i] = np.linalg.matrix_rank(W)
        acc_arr[i] = np.sum(y_result == y_test)/float(y_test.shape[0])


その結果が以下の図です(赤が推定精度、青がランクです)。

f:id:YamagenSakam:20180613231839p:plain

L1ノルムロジスティック回帰と同様、\lambdaを大きくしすぎても、制約がきつくなりすぎるのか精度が徐々に下がっていきます。実際今回の条件で推定精度が最大になったのは\lambda = 1250のときで、精度が68.0%、{\bf W}のランクが4でした。以下には、この\lambda = 1250での{\bf W}特異値分解して得られる空間フィルタ{\bf U}を、y=0クラスのX_trainとy=1クラスのX_trainの平均行列に施した波形を示します。

  • y=0に対応する観測データ集合X_train1の平均に{\bf U}をかけて得られる波形

f:id:YamagenSakam:20180613232012p:plain

  • y=1に対応する観測データ集合X_train2の平均に{\bf U}をかけて得られる波形

f:id:YamagenSakam:20180613231951p:plain

この図の中で青の波形は最大特異値に対応する空間フィルタで抽出された波形、緑の波形は2番目に大きい特異値に対応する空間フィルタで抽出された波形です。この両者の注目すると、空間フィルタを施すことによって信号源に潜在しているクラス分類に有用な成分を抽出できていることがわかります。

まとめ

今回は行列データの分類アルゴリズムについて、ロジスティック回帰を例にとってトレースノルム正則化項によって低ランクな解が得られること、多次元の時系列データの分類問題などに有効であることを見てきました。トレースノルムの正則化は信号源の中でも一部の成分のみがクラス分類に有用であるとき特に威力を発揮します。
なお、時系列データの分類には今回紹介した方法以外にも隠れマルコフモデルを使った手法やリカレントニューラルネットワークの一種であるLong Short-Term Memory(LSTM)を使った手法もあります。特に最近は深層学習の流れもあってLSTMが主流になってきているとのことなので、そろそろニューラルネット辺りにも手を出していきたいところです。

参考文献

以前も紹介したこちらの本の第8章が非常に参考になります。

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

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

あと、同著者によるこちらのスライドもわかりやすいです。
行列およびテンソルデータに対する機械学習

*1:どうやら、前回のベクトルデータの分類においても、基本的にバイアス項は正則化には含めない方が良いようです

近接勾配法応用編その2 ~L1ノルム正則化項によるスパースなロジスティック回帰~

今回は前々回の記事で書いた近接勾配法の応用例第2段ということで、分類問題などで使われるロジスティック回帰の目的関数にL1ノルムの正則化項をつけて、スパースな解を導いてみます。スパースな解を導出するメリットとしては、クラス分類に寄与しない無駄な属性の係数が自動的に0となり、クラス間で有意差がある属性を見つけ出せ、かつ、分類モデルの複雑度も下げられるので過学習の防止も期待できるというメリットがあります。なお、ロジスティック回帰は多クラス分類にも利用できますが、今回は2クラス分類の例で説明します。

ロジスティック回帰

ロジスティック回帰について詳しいことはたくさん文献が出ているのでそちらを参照いただくとして、ここでは概略と最終的な目的関数を示します。
まず、学習データ集合 \left\{ {\bf x}_i, y_i \right\} (i = 1 \cdots N)があるとします。ここで、{\bf x}d次元のデータベクトル、yy \in \left\{0, 1\right\}のラベルです。
ここで、{\bf x}が与えられてy = 1となる事後確率は
\displaystyle p(y = 1| {\bf x}) = \frac{p(y=1, {\bf x})}{p(y=1, {\bf x}) + p(y = 0, {\bf x})} =  \frac{1}{1 + \exp(-a) }= \sigma(a) \tag{1}
と表せます。ここで、a
\displaystyle a= \ln{\frac{p(y=1, {\bf x})}{p(y=0, {\bf x})}} \tag{2}
という式であり、ロジット関数とも呼ばれます。{\bf x}がクラスy = 1に所属する確率が大きければこのロジット関数は大きな値を取るので事後確率\sigma(a)は1に近づきますし、逆にクラスy=0に所属する確率が大きければロジット関数は小さな値を取るので\sigma(a)は0(つまり、クラスy=0となる事後確率1-\sigma(a)は1)に近づきます。
そしてこのa{\bf x}の線形関数としてモデル化したのがロジスティック回帰です。具体的には、


 \displaystyle a = {\bf w}^T {\bf x} + w_0 = \tilde{{\bf w}}^T \tilde {{\bf x}} \tag{3}
ただし、 \tilde{{\bf w}}= \left[ w_0 , {\bf w} \right],  \tilde{{ \bf x}}= \left[1,  {\bf x}\right] *1
となります。
さて、これまでの内容を踏まえて、学習データ集合に対する尤度関数を考えます。上記の議論より、{\bf x}_i\tilde{{\bf w}}が与えられて、そのクラスがy_i = 1である尤度は\sigma(a) = \sigma(\tilde{{\bf w}}^T \tilde{{\bf x}}_i)です。一方で、y_i = 0である尤度は (1-\sigma(\tilde{{\bf w}}^T \tilde{{\bf x}}_i))です。このように考えると尤度関数は、
 \displaystyle p({\bf y}| \tilde{{\bf w}})= \prod_{i=1}^{N}  \sigma(\tilde{{\bf w}}^T \tilde{{\bf x}}_i)^{y_i}  (1 - \sigma(\tilde{{\bf w}}^T \tilde{{\bf x}}_i))^{(1- y_i)} \tag{4}
となります。後は式(4)を最大化する\tilde{{\bf w}}を求めればよいです。ただ、この式(4)のままだと計算しにくいので、以下のような式(4)の負の対数を取った誤差関数を最小化することを考えます。
\displaystyle E(\tilde{{\bf w}}) = - \ln { p({\bf y}| \tilde{{\bf w}}) } = -\sum_{i=1}^{N} \left\{  {y_i} \ln{\sigma(\tilde{{\bf w}}^T \tilde{{\bf x}}_i)} + (1- y_i)\ln{(1 - \sigma(\tilde{{\bf w}}^T \tilde{{\bf x}}_i))} \right\} \tag{5}
これがロジスティック回帰の目的関数です。

正則化項によるスパース化と近接勾配法の適用

ここまでで普通のロジスティック回帰の目的関数が得られました。このロジスティック回帰をスパースな解にするためには、L1ノルム正則化項をつけた式
 \displaystyle F( \tilde{{\bf w}})= -\sum_{i=1}^{N} \left\{  {y_i} \ln{\sigma(\tilde{{\bf w}}^T \tilde{{\bf x}}_i)} + (1- y_i)\ln{(1 - \sigma(\tilde{{\bf w}}^T \tilde{{\bf x}}_i))} \right\} + \| \tilde{{\bf w}} \|_1 \tag{6}
を近接勾配法により最小化すればよいです。近接勾配法を解くためには、元の目的関数E(\tilde{{\bf w}})の勾配と正則化\| \tilde{{\bf w}} \|_1のproximal operatorが必要でした。
正則化項のproximal operatorは以前の記事で導出したとおりです。また、E(\tilde{{\bf w}})の勾配は、
\displaystyle \nabla E(\tilde{{\bf w}}) = \sum_{i=1}^{N} (\sigma(\tilde{{\bf w}}^T \tilde{{\bf x}}_i) - y_i) \tilde{{\bf x}}_i \tag{7}
となります。後は、これらに基づき近接勾配法を適用するだけです。

実装と実験

これまでの内容pyhonコードで実装すると下記のようになります。なお、近接勾配法のアルゴリズム部分についても、前々回記事で記したコードから停止条件部分を修正したので、改めて近接勾配法の関数も書きます(以前の記事から変更していない関数については下記に記していません)。

import numpy as np
from numpy.random import *
import matplotlib.pyplot as plt

#近接勾配法
def proximal_gradient(grad_f, prox, gamma, objective, init_x, tol = 1e-9):
    x = init_x
    result = objective(x)
    while 1:
        x_new = prox(x - gamma * grad_f(x), gamma)
        result_new = objective(x_new)
        #停止条件を説明変数差分のノルムに変更
        #if result_new < tol or (np.abs(result - result_new)/np.abs(result) < tol) == True :
        if np.linalg.norm(x_new - x, 2) < tol:
            break;
        x = x_new
        result = result_new
    return x_new, result_new

#l1ノルム正則化ロジスティック回帰
#X:学習用データベクトルを列ベクトルとして並べた行列
#y:ラベルを並べたベクトル
def l1norm_LogisticRegression(X, y, lam):
    sigma = lambda a : 1.0/(1.0+np.exp(-a))
    p=lambda Z, w:sigma(np.dot(Z.T, w))  # ※1下記に補足あり
    X=concatenate((ones((1,size(X,1))), X), 0) #バイアス項の追加
    objective = lambda w: - np.sum( y * np.log(p(X, w)) +
                        (1 - y) *np.log(( 1 - p(X, w)))) + lam * np.sum(np.abs(w))
    grad_E = lambda w: np.dot(X, p(X, w) - y) #※2下記補足あり
    (u, l, v) = np.linalg.svd(X)
    gamma = 1.0/max(l.real*l.real)
    prox = lambda v, gamma:prox_norm1(v, gamma, lam) #proximal operator
    w_init = np.zeros(X.shape[0])
    (w_hat, result) = proximal_gradient(grad_E, prox, gamma, objective, w_init, 1e-12)
    return w_hat, result

繰り返し構文を避けるため、行列・ベクトル演算で勾配計算などを行っています。これを少し補足しておくと、まずソースコード※1に当たる部分は、p({\bf Z}, {\bf w}) = \left[\sigma({\bf w}^T {\bf z}_1) , \cdots, \sigma({\bf w}^T {\bf z}_N) \right]^Tというベクトルを返す関数として定義しています。そして、※2の部分では、
\displaystyle{\bf X}(p({\bf X}, {\bf w}) - {\bf y}) = \left[ {\bf x}_1, \cdots, {\bf x}_N \right] \left[\sigma({\bf w}^T {\bf x}_1)  - y_1, \cdots, \sigma({\bf w}^T {\bf x}_N)  - y_N  \right]^T= \sum_{i=1}^{N} {\bf x}_i (\sigma({\bf w}^T {\bf x}_i)  - y_i)
というように、行列・ベクトルの演算で式(7)の勾配計算を行う関数を定義しています。
このL1ノルム正則化つきロジスティック回帰を疑似データで実験してみました。疑似データは以下のようなコードで生成しました。

def main():
    #学習データ
    #最初の2次元だけ分類に寄与する属性
    X_train1=np.random.randn(2,200)
    X_train2=np.random.randn(2,200) + ones((2,200))
    X_train = concatenate((X_train1,X_train2),1)
    #残り123次元は分類に寄与しないランダムな属性
    X_train = concatenate((X_train, np.random.randn(123, 400)), 0)
    y_train1 = ones((200,1))
    y_train2 = zeros((200,1))
    y_train = concatenate((y_train1,y_train2),0).squeeze()

    #テストデータ
    X_test1=np.random.randn(2,200)
    X_test2=np.random.randn(2,200)+ones((2,200))
    X_test = concatenate((X_test1,X_test2),1)
    X_test = concatenate((X_test, np.random.randn(123, 400)), 0)
    y_test1=ones((200,1))
    y_test2=zeros((200,1))
    y_test=concatenate((y_test1,y_test2),0).squeeze()

疑似データの次元は125次元ですが、その内最初の2次元だけクラス分類において意味があり、それ以外の123次元は全くのランダムでノイズとなる属性です。また、y=0のクラスにのみオフセットをはかせている形なので、バイアス成分もクラス分類に寄与します。そのため、今回のアルゴリズムを適用すると、最初の2次元の係数およびバイアス項のみが非0となり、それ以外の係数は0になることが期待されます。
このように生成した疑似データに対し、以下のように\lambda=0, 1, \cdots, 39と変化させて、その時の推定精度とwの非0要素の数を求めてみました。

    non_zero_cnt_arr = np.zeros(40)
    acc_arr = np.zeros(40)
    for lam in range(0, 40):
        (w, r) = l1norm_LogisticRegression(X_train, y_train.squeeze(), lam)
        y_result = np.dot(w[1:126], X_test) + w[0]
        y_result[y_result>=0] = 1
        y_result[y_result<0] = 0
        non_zero_cnt_arr[lam] = np.sum(np.abs(w)>1e-12)
        acc_arr[lam] = np.sum(y_result == y_test)/float(y_test.shape[0])

    fig, ax1 = plt.subplots()
    ax1.plot(acc_arr, "r.-")
    ax1.set_ylabel("accuracy")
    ax2=plt.twinx()
    ax2.plot(non_zero_cnt_arr, "bx-")
    ax2.set_ylabel("non zero element count")
    plt.legend()
    plt.show()

その結果は下記のようになりました(赤線が推定精度、青線が係数の非0要素数です)。
f:id:YamagenSakam:20180603015719p:plain

\lambda = 0のときは、すべての要素が非0になりノイズである属性も判別に用いてしまっているので精度が60.5%とあまり良くありません。そこから\lambdaを大きくしていくと非0要素数は減少していくことがわかります。一方で推定精度については、適切な\lambdaであれば精度が上がっていますが、逆に\lambdaが大きすぎると制約がきつくなりすぎるためか精度は下がっていきます。面白いのは、非0要素数が3のとき(クラス分類に寄与する2次元 + バイアス項のみが非0のとき)、推定精度がピークとなると予想していましたが、予想に反して\lambda=15で非0要素数8のときに推定精度が最大(75.5%)になる結果が得られました*2。これは、適切にクラス分類に寄与する属性を見つけられるまで\lambdaを上げると、ロジスティック回帰の制約条件としては厳しくなりすぎてしまい結果、精度の低下を招くのだと考えられます。そのため、クラス分類問題とどの属性に有意差があるか見つける問題とは一緒くたにして解くのではなく、分けて解くべきなのだと思います。

まとめ

今回はロジスティック回帰にL1ノルム正則化項をつけることによりスパースな解が得られ、分類器としての精度も向上させられることを見てきました。また、分類に寄与する属性もこのL1ノルム正則化を使うことで見つけられますが、だからと言って制約条件が厳しすぎる\lambdaを用いると分類精度低下を招くので注意が必要です。一旦、変数選択の目的でL1ノルム正則化項付きのロジスティック回帰を行って、そこで係数が非0となった要素だけを使って今度は制約なしのロジスティック回帰を再学習させる方法がよいのかも知れません。なお、今回はL1ノルムの正則化を試しましたが、もちろんL1ノルム以外にもグループL1ノルムやL1/L2の混合ノルムも適用できます。
今後の予定としては、もうひとつ近接勾配法の応用として行列データの分類について紹介したり、Alternating Direction Method of Multipliers (ADMM)を調べて書いたり、ICML2018の論文(そろそろ公開?)を読んだりなんかをしたいと思います。

*1:バイアス項w_0{\bf w}に統合し、d + 1次元のベクトルとして扱います。

*2:ランダムな疑似データなので試行毎に結果は変わりますが、傾向としては毎回同じで非0要素数が3よりも多い時点で推定精度のピークを迎えます。

近接勾配法応用編その1 ~スパースコーディング、辞書学習からの超解像~

はじめに

前回の記事で近接勾配法の基本と実装を見てきました。今回はこの近接勾配法を利用して、スパースコーディングの応用の1つである超解像技術を調べPythonで実装してみました*1。画像処理関係の実装はほぼ未経験だったので、その分野の人が見れば前処理等が足りなかったり、アンチパターン的なことをやってたりするかもしれませんが、その点は指摘頂けると幸いです。

スパースコーディングとは

スパースコーディングは観測した信号{\bf y}を、「辞書」と呼ばれるm個の基底ベクトルの集合  {\bf D} = \left[{\bf d}_1, {\bf d}_2, \cdots, {\bf d}_m \right]のスパースな線形結合で表現する符号化手法です。
今、この線形結合の係数を{\bf w}とすると、スパースコーディングは下記のように定式化されます。
\displaystyle {\bf y }\simeq {\bf D} {\bf w} \tag{1} 
ただし、{\bf w}の要素の多くは0(つまり{\bf w}はスパース)です。

このようなスパースコーディングを使うことによりデータの圧縮を行えるだけでなく、画像からのノイズ除去や今回取り上げる超解像など様々な技術へ応用が可能です。
スパースコーディングを成功させる上で重要となるのは辞書{\bf D}です。{\bf D}がテキトーだとうまくいきません*2。ということで、辞書{\bf D}は実データからの学習により構築します。ここで、学習データの個数をn、学習データをまとめた行列{\bf Y} = \left[{\bf y}_1 \cdots, {\bf y}_n\right] 、 各線形結合の係数をまとめた行列{\bf W} = \left[{\bf w}_1, \cdots, {\bf w}_n\right]を定義すると、辞書学習の目的関数は下記のようになります。
\displaystyle {\rm arg} \min_{{\bf D}, {\bf W}} \| {\bf Y} - {\bf D}{\bf W} \|^2_2 \tag{2} 
ただし、{\bf W}はスパース。

この方法で辞書{\bf D}が得られたら、次はその辞書を用いて入力画像の符号化します。そして、その係数と辞書から画像を再構成し超解像画像を得るという流れになります。
ということで、スパースコーディングによる超解像生成は(1)辞書を学習するフェーズと(2)その辞書を元に画像を再構成するフェーズの2フェーズで実現できます。

画像再構成のための辞書学習

実装した手法の説明

では、今回実装した学習用の画像データを元に辞書を学習する方法を説明します。
まず、最終的に得たい解像度の画像データ(高解像画像データ)複数枚を学習データとして用意します。この時、超解像を行いたい画像のカテゴリが決まっているなら、そのカテゴリに合った画像を用意するのがよさそうです。
この学習データを用いて、辞書{\bf D}を学習する流れを下記に示します。

  1. 学習用画像を小領域を取り出した「画像パッチ」をn個用意し、それぞれ縦に並べたベクトル{\bf y}_1, {\bf y}_2, \cdots, {\bf y}_nおよび、これらのベクトルを横に並べた行列を{\bf Y} = \left[{\bf y}_1 \cdots, {\bf y}_n\right] を作る。
  2. {\bf Y}の各行が平均0、分散1になるように標準化する。
  3. {\hat {\bf D}}をランダム行列で初期化する。
  4. {\hat {\bf D}}を固定して、次の最適化問題i = 1, 2, \cdots, nについて解き、行列{\hat {\bf W}} = \left[{\hat {\bf w}_1}, \cdots, {\hat {\bf w}_n} \right]を得る。
  5. \displaystyle {{\hat {\bf w}_i}} = {\rm arg} \min_{{\bf w}_i} \| {\bf y}_i - {\hat {\bf D}}{\bf w}_i \|^2_2  + \lambda \| {\bf w}_i \|_1\tag{3}
  6. \hat{{\bf W}}を固定して、次の最適化問題を解き\hat{{\bf D}}を得る。
  7. \displaystyle \hat{{\bf D}}=  \arg \min_{{\bf D}} \| {\bf Y} - {\bf D} {\hat{ \bf W}} \|^2_2 \tag{4}
  8. 手順4,手順5の操作を収束するまで、もしくは一定回数繰り返す。

ここで、式(3)はlassoの式であるため解がスパースになることが期待され、更に近接勾配法で解くことができます。なお、今回、近接勾配法の応用編ということで手順4の方法でスパースな係数\bf{w}_iを求めていますが、この方法以外にもOrthogonal Matching PursuitやIterative Reweighted Least Squaresなどのアルゴリズムがあります。
また、手順5で式(4)の解を得る方法としてはK-SVDという手法がよく使われますが、今回は簡単のために下記のように解析的に解いた結果を使用します。
\displaystyle {\hat {\bf D}} = {\bf Y} {\hat {\bf W}}^T ( {\hat {\bf W}} {\hat {\bf W}}^T  )^{-1} \tag{5}

この手順4、手順5を繰り返して{\bf W}, {\bf D}の同時最適化を行うわけですが、この両変数に対しては凸関数ではないため(式(3)、式(4)はそれぞれ凸関数ですが)、求まる結果は局所最適解です。
なお、この学習を収束まで待っているとかなり時間がかかるため、今回は適当な回数ループさせたら学習を完了させています。

ちなみに、最初は手順2の標準化を省いた実装をしていましたが、これだとうまく超解像ができませんでした。やはり、学習データの標準化は重要のようです。

実装したコード

ここまでの手順を実現するpythonコードです。

import numpy as np
import scipy as sci
from numpy.random import *
from PIL import Image
import os
import copy
import sys

#dorectory以下のファイルを再帰的に検索する関数
def find_all_files(directory):
    for root, dirs, files in os.walk(directory):
        yield root
        for file in files:
            yield os.path.join(root, file)

#おれおれ画像・ベクトル・行列変換関数群
def color2gray_img(img):
    return img.convert('L')

def im2vec(img):
    img_arr = np.array(img.convert('L'),'float')
    return img_arr.reshape(-1)

def gray_im2float_vec(gray_img):
    return im_scale2float((im2vec(gray_img)))

def color_im2float_vec(color_img):
    return gray_im2float_vec(color2gray_img(color_img))

def vec2image(vec, H, W):
    tmp_vec = copy.deepcopy(vec)
    tmp_vec = tmp_vec.reshape((H, W))
    return Image.fromarray(np.uint8(tmp_vec))

def float_vec2image(vec, H, W):
    return vec2image(float2im_scale(vec), H, W)

def arr2image(arr):
    arr = copy.deepcopy(arr)
    return Image.fromarray(np.uint8(arr))

def float_arr2image(arr):
    return arr2image(float2im_scale(arr))

def im_scale2float(vec):
    return vec.astype(np.float)

def float2im_scale(vec):
    vec[vec > 255.0] = 255.0
    vec[vec < 0.0] = 0.0
    return vec

#画像データからパッチを切り出す関数
#shift_numで切り出す際のずらし量を指定
def gray_im2patch_vec(im, patch_hight, patch_width, shift_num):
    patchs = np.zeros((patch_hight*patch_width,
                        (1 + (im.size[0] - patch_hight)/shift_num + (1 if (im.size[0] - patch_hight) % shift_num  != 0 else 0)) *
                        (1 + (im.size[1] - patch_width)/shift_num + (1 if (im.size[1] - patch_width) % shift_num  != 0 else 0))))
    im_arr = im_scale2float(np.array(im))
    k = 0
    for i in range(1 + (im.size[0] - patch_hight)/shift_num + (1 if (im.size[0] - patch_hight) % shift_num  != 0 else 0)):
        h_over_shift = min(0, im.size[0] - ((i* shift_num) + patch_hight)) #ずらしたことによって画像からはみ出る量(縦方向)
        for j in range(1 + (im.size[1] - patch_width)/shift_num + (1 if (im.size[1] - patch_width) % shift_num  != 0 else 0)):
            w_over_shift = min(0, im.size[1] - ((j* shift_num) + patch_width)) #ずらしたことによって画像からはみ出る量(横方向)
            #パッチの切り出し 画像からはみ出たら、そのはみ出た分を戻して切り出し
            patchs[:, k] = im_arr[(i* shift_num) + h_over_shift:((i* shift_num) + patch_hight) + h_over_shift,
                                    (j* shift_num) + w_over_shift:((j* shift_num) + patch_width) + w_over_shift].reshape(-1)
            k += 1
    return patchs


#式(5)を計算する関数
def dictionary_learning_with_mod(W, Y):
    D = np.dot(np.dot(Y, W.T), np.linalg.inv(np.dot(W, W.T)))
    result = np.linalg.norm(Y - np.dot(D, W)) ** 2
    return D, result

#辞書クラス
class dictionary_with_learning():
    def __init__(self, weight_learner, dictionary_learner, base_num, max_loop_num = 8):
        self.Dic = np.array([])
        self.weight_learner = weight_learner
        self.dictionary_learner = dictionary_learner
        self.base_num = base_num

    #学習データ行列Yから辞書を学習するメソッド
    def learn(self, Y, eps = 1e-8):
        #手順3:ランダム行列で辞書の初期化
        D = np.random.randn(Y.shape[0], self.base_num)

        result = -1
        loop_num = 0
        while 1:
            loop_num += 1
            W = np.zeros((self.base_num, Y.shape[1]))
            #手順4:Dを固定して係数を求める
            for i in range(Y.shape[1]):
                [W[:, i], tmp] = self.weight_learner(D, Y[:, i])

            #手順5:Wを固定してDを求める
            [D, new_result] = self.dictionary_learner(W, Y)

            if result < 0 or np.abs(new_result - result) > eps or loop_num < max_loop_num:
                result = new_result
            else: #手順4、5を収束するまで繰り返す。
                break;
        self.Dic = D
        return D

    def get_dic_array(self):
        return self.Dic



#画像パッチから学習される辞書クラス
class image_patch_dictionary():
    def __init__(self, dictionary):
        self.dictonary = dictionary
        self.patch_hight = 0
        self.patch_width = 0

    #学習画像の集合から辞書を学習するメソッド
    def learn_images(self, image_list, patch_hight, patch_width, eps=2e-2, shift_num=1):
        self.patch_hight = patch_hight
        self.patch_width = patch_width
        Y = np.empty((patch_hight * patch_width, 0))
        #手順1:画像パッチの切り出し
        for img in image_list:
            gray_img = color2gray_img(img)
            Ytmp = gray_im2patch_vec(gray_img, patch_hight, patch_width,shift_num) 
            Y = np.concatenate((Y, Ytmp), axis=1)

        #手順2:学習データ行列の標準化
        Y_ = (Y - np.mean(Y, axis=1, keepdims = True))/np.std(Y, axis = 1, keepdims=True)
        self.dictonary.learn(Y_, eps)

    def get_dic(self):
        return self.dictonary.get_dic_array()

    def get_patch_hight(self):
        return self.patch_hight

    def get_patch_width(self):
        return self.patch_width


def main():
    lam = 5
    base_num = 144
    patch_hight = 12
    patch_width = 12
    train_img_list = []
    for file in find_all_files('C:\\data'):
        if os.path.isfile(file):
            try:
                train_img_list.append(Image.open(file))
                if len(train_img_list) >= 50:
                    break;
            except:
                print file

    weight_learning_method = lambda W, y:ISTA(W, y, lam) #近接勾配法を行う関数を使用
    dic_learning_method = lambda W, Y:dictionary_learning_with_mod(W, Y)#解析的に辞書Dを求める方法を使用
    im_dic = image_patch_dictionary(dictionary_with_learning(weight_learning_method, dic_learning_method, base_num))
    im_dic.learn_images(train_img_list, patch_hight, patch_width, shift_num = patch_hight)

    np.save('dic_array.npy', im_dic.dictonary.Dic)

本実装は、画像⇔行列・ベクトル変換を行う関数群、パッチ切り出しを行う関数、画像パッチ辞書クラスimage_patch_dictionary、辞書クラスdictionary_with_learningから構成されています。
image_patch_dictionaryは実際の辞書学習オブジェクトであるdictionary_with_learningをインスタンス変数として持ち、与えられた学習画像群からパッチの切り出しを行い、辞書を学習させる機能を持ちます。
そのため、辞書学習を行うのはdictionary_with_learningになりますが、係数の学習を行うweight_learnerと辞書の学習を行うdictionary_learnerをインスタンスとして持ち、生成時にどのアルゴリズムで係数推定、辞書推定を行うかを指定します。今回、係数の計算には近接勾配法を使用しますので、前回の記事記載のISTA関数をweight_learnerとして指定し、dictionary_learnerには新たに作成したdictionary_learning_with_mod関数を指定します(main関数参照)。もし係数推定として別のアルゴリズムを使いたければweight_learnerを変えればよいし、辞書学習もK-SVDなどを使いたければdictionary_learnerを変更すればよいです。

画像再構成による超解像

実装した手法の説明

ここからは、スパースコーディングによる超解像の手法の説明です。手法の大まかな考え方は「学習した辞書{\bf D}をダウンサンプリングした辞書{\bf D}_lで入力画像を符号化する係数を求めた後、その係数と元の辞書\bf Dから画像を再構成する」というものです。以下、順を追って説明します。

  1. 画像パッチの切り出し
  2. 今、入力画像(低解像画像)を{\bar{\bf Y}}とし、そこから辞書学習時と同様に画像パッチを切り出します。その切り出した画像パッチを{\bar{\bf y}_i}とします。切り出す画像パッチのサイズは、辞書学習に使用した画像(高解像画像)のパッチサイズをh \times w、低解像画像から高解像画像への倍率をsとすると、\frac{h}{s} \times \frac{w}{s}になります。ということで、{\bar {\bf y}_i}\frac{hw}{s^2}次元のベクトルになります。ここで画像パッチは、もれなくかつ領域を重複させて取得します。領域を重複させないと、画像再構成後パッチ間のつなぎ目が滑らかにならずブロック状のアーティファクト発生します。なるべく重複度が大きい方が精度は良くなりますが、取得する画像パッチも多くなるので、その分計算に時間がかかります。
  3. 辞書のダウンサンプリング
  4. 次に辞書のダウンサンプリングです。ここの方法も色々考えられますが、今回は単純に学習した{\bf D}にダウンサンプリングの作用素をかけて{\bf D}_lを得る方法をとりました。 {\bf D}は学習画像パッチのサイズhw \times基底数mの行列です。各基底{\bf d}h \times wの画像を縦に並べたベクトルと見なせるので、これをh \times wの画像に戻し、s近傍の画素を平均することでダウンサンプリング後の辞書{\bf D}_lを得ます({\bf D}_l\frac{hw}{s^2} \times mの行列になります)。
  5. 符号化係数の推定
  6. 画像パッチが切り出せたら次は、それぞれ符号化する係数を求めます。これは、ダウンサンプリング後の辞書{\bf D}_lを使い、下記のように式(3)と同様の最適化問題を解けばよいです。 \displaystyle {{\hat {\bf w}_i}} = {\rm arg} \min_{{\bf w}_i} \| {\bar {\bf y}_i} - {\bf D}_l {\bf w}_i \|^2_2  + \lambda \| {\bf w}_i \|_1\tag{6}
  7. 高解像画像パッチの生成
  8. ここまでで、各パッチを符号化するスパースな係数が求まりました。次はこの係数を用いて高解像画像を構成する画像パッチ{\bf x}_iを求めます。やることは簡単で、以下のように元々の辞書(高解像画像の辞書){\bf D}と低解像辞書で推定した係数{\hat{\bf w}_i}をかけるだけです。 \displaystyle {\bf x}_i = {\bf D}{\hat {\bf w}_i} \tag{7}
  9. パッチをつなぎ合わせて高解像画像を再構成
  10. 得られた高解像画像パッチから高解像画像{\bf X}_0を再構成します。手順1で画像パッチを領域を重複させて切り出したことに注意が必要です。今、手順1で画像を切り出す際pピクセルずらして切り出したとすると、再構成時はpsピクセルずらしてパッチを貼り合わせます。全てのパッチを貼り合わせた後、各ピクセルを加算回数分で割り平均をとります。この画像パッチのつなぎ合わせについても色々な工夫が提案されているようです。
  11. バックプロジェクションによる画像の補正
  12. 最後に仕上げです。パッチを貼り合わせ高解像画像{\bf X}_0は、パッチごとに最適化されているため、画像全体の仮定が考慮されていません。その仮定というのは、「入力画像は{\bar {\bf Y}}、高解像画像{\bf X}に、劣化やダウンサンプリングを施す作用素{\bf L}をかけることによって得られたものである。」というものです。つまり、{\bar {\bf Y}} = {\bf L} {\bf X}という仮定が前提にあり、それと推定した{\bf X}_0とのバランスをとることで、解の補正を行います。具体的には下記の最適化問題を解きます。 \displaystyle {\bf X}^* = {\rm arg} \min_{{\bf X}} \|{\bar {\bf Y}} - {\bf L} {\bf X}\|^2_2 + \mu \|{\bf X} - {\bf X_0}  \|\tag{8} ここで、\bf{X}_0, {\bar{\bf Y}}, {\bf X}は画像全体の各要素を縦に並べたベクトルです。また、今回の実験では、{\bf L}はダウンサンプリングのみを施す行列として定義しています。 この最適化は解析的に求められそうですが、多くの論文では勾配法を使って求めています。逆行列計算が大変だからでしょうか?とりあえず、今回の実装でも論文に倣って勾配法で解きました。

実装したコード

以上の内容の実装が下記です。モジュールのimportやおれおれ画像変換関数群は辞書学習時のコードと同様のものを使っています。

#バックプロジェクションで使うダウンサンプラー行列を求める関数
#ダウンサンプラーはスパースかつサイズが大きい行列なのでscipy.sparseのスパース行列型として返す
def get_downsampler_arr(from_H, from_W, to_H, to_W):
    S = sci.sparse.lil_matrix((to_H * to_W, from_H * from_W))
    k = 0
    for i in range(0, from_H, from_H/to_H):
        for j in range(0, from_W, from_W/to_W):
            s = sci.sparse.lil_matrix((from_H, from_W))
            s[i : i + from_H/to_H, j:j + from_W/to_W] = (float(to_W)/float(from_W)) * (float(to_H)/float(from_H))
            S[k, :] = s.reshape((1, from_H*from_W))
            k += 1
    return  S.tocsr()

#ベクトルデータのダウンサンプリング関数
def vec_downsampling(vec, from_H, from_W, to_H, to_W):
    arr = vec.reshape((from_H, from_W))
    ret_arr = arr_downsampling(arr, from_H, from_W, to_H, to_W)
    return ret_arr.reshape(-1)

#行列データのダウンサンプリング関数
def arr_downsampling(arr, from_H, from_W, to_H, to_W):
    ret_arr = np.zeros((to_H, to_W))
    k = 0
    for i in range(0, from_H, from_H/to_H):
        for j in range(0, from_W, from_W/to_W):
            ret_arr[i/(from_H/to_H), j/(from_W/to_W)] = (np.average(arr[i : i + from_H/to_H, j:j + from_W/to_W]))
            k += 1
    return ret_arr

#画像パッチをつなぎ合わせて画像を再構成する関数
#shift_numがずらし量、パッチがオーバーラップしている個所は足し合わせて平均化
def patch_vecs2gray_im(patchs, im_hight, im_width, patch_hight, patch_width, shift_num):
    im_arr = np.zeros((im_hight, im_width))
    sum_count_arr = np.zeros((im_hight, im_width)) #各領域何回パッチを足したか表す行列
    i = 0
    j = 0
    w_over_shift = 0
    h_over_shift = 0
    for k in range(patchs.shape[1]):
        P = patchs[:, k].reshape((patch_hight, patch_width))
        im_arr[(i* shift_num) + h_over_shift:((i* shift_num) + patch_hight) + h_over_shift,
            (j* shift_num) + w_over_shift:((j* shift_num) + patch_width) + w_over_shift] += P #指定の領域にパッチの足しこみ
        sum_count_arr[(i* shift_num) + h_over_shift :(i* shift_num) + patch_hight + h_over_shift,
            (j* shift_num) + w_over_shift:(j* shift_num) + patch_width + w_over_shift] += 1 #当該領域のパッチを足しこんだ回数をカウントアップ
        if j < ((im_width - patch_width)/shift_num + (1 if (im_width - patch_width) % shift_num  != 0 else 0)):
            j += 1
            w_over_shift = min(0, im_width - ((j* shift_num) + patch_width)) #パッチを足しこむ領域が画像からはみ出た時、そのはみ出た分を戻すための変数(横方向)
        else:
            j = 0
            i += 1
            h_over_shift = min(0, im_hight - ((i* shift_num) + patch_hight)) #パッチを足しこむ領域が画像からはみ出た時、そのはみ出た分を戻すための変数(縦方向)
            w_over_shift = 0
    im_arr /= sum_count_arr
    return float_arr2image(im_arr), im_arr

#バックプロジェクションを行う関数
def back_projection(x0, y, L, c, mu, x_init = [], eps = 1e-7):
    if x_init == []:
        xt = sci.matrix(np.random.rand(x0.shape[0], x0.shape[1])) *255.0
    else:
        xt = x_init
    result = -1
    while 1:
        x_new = xt - mu * (L.T.tocsr()*(L * xt - y) + c* (xt - x0))
        result_new = sci.linalg.norm(L * x_new - y, 2) + c * sci.linalg.norm(x_new - x0, 2)
        if np.abs(result_new - result) < eps:
            break;
        xt = x_new
        result = result_new
    return x_new, result_new

#超解像計算クラス
class super_resolution:
    def __init__(self, dictionary, weight_learner, back_projection):
        self.Dh = dictionary.get_dic()
        self.patch_hight = dictionary.get_patch_hight()
        self.patch_width = dictionary.get_patch_width()
        self.weight_learner = weight_learner
        self.back_projection_method = back_projection

    def get_super_resolution(self, image, scaling, shift_num = 1):
        #手順1:画像パッチの切り出し
        image = color2gray_img(image)
        Y = gray_im2patch_vec(image, self.patch_hight/scaling, self.patch_width/scaling, shift_num)

        #手順2:辞書のダウンサンプリング
        Dl = np.zeros((self.Dh.shape[0]/(scaling**2), self.Dh.shape[1]))
        for i in range(self.Dh.shape[1]):
            Dl[:, i] = vec_downsampling(self.Dh[:, i], self.patch_hight, self.patch_width,
                                        self.patch_hight/scaling, self.patch_width/scaling)
       
        im_arr = np.zeros(image.size)
        i = 0
        j = 0
        Patchs = np.zeros(((self.patch_hight) * (self.patch_width), Y.shape[1]))
        W = np.zeros((self.Dh.shape[1], Y.shape[1]))
        for k in range(Y.shape[1]):
            #手順3:符号化係数の推定
            [w, result] = self.weight_learner(Dl, Y[:, k])

            #手順4:高解像画像パッチの生成
            Patchs[:, k] = np.dot(self.Dh, w)

        #手順5:パッチをつなぎ合わせて高解像画像を再構成
        [reconstructed_im, X0] = patch_vecs2gray_im(Patchs, image.size[0] * scaling , image.size[1] * scaling , self.patch_hight, self.patch_width, shift_num* scaling)
        L = get_downsampler_arr(image.size[0] * scaling, image.size[1] * scaling, image.size[0], image.size[1])

        #手順6:バックプロジェクションによる画像の補正
        x = np.array(self.back_projection_method(sci.matrix(X0.reshape(-1)).T, sci.matrix(np.array(image).reshape(-1)).T, L)[0])
        return float_arr2image(x.reshape((image.size[0] * scaling, image.size[1] * scaling)))


def main():
    test_img = Image.open('Lenna.png')
    scale = 3
    lam = 50
    base_num = 144
    patch_hight = 12
    patch_width = 12
    shift_num = 1
    c = 0.8
    mu = 0.9
    weight_learning_method = lambda W, y:ISTA(W, y, lam)
    back_projection_method = lambda x0, y, L:back_projection(x0, y, L, c, mu)
    dic_learning_method = lambda W, Y:dictionary_learning_with_mod(W, Y)
    im_dic = image_patch_dictionary(dictionary_with_learning(weight_learning_method, dic_learning_method, base_num))
    im_dic.dictonary.Dic = np.load('dic_array.npy')
    im_dic.patch_hight = 12
    im_dic.patch_width = 12
    SR = super_resolution(im_dic, weight_learning_method, back_projection_method)
    result_img = SR.get_super_resolution(test_img, scale, shift_num)
    result_img.save("Lenna_SR.png")

超解像計算クラスsuper_resolutionは、画像パッチ辞書オブジェクト、係数推定関数オブジェクト、バックプロジェクション関数オブジェクトを生成時に指定し、get_super_resolutionメソッドで指定した辞書、関数に基づいた超解像画像の生成を行います。推定方法は上記で解説したとおりです。
画像パッチをつなぎ合わせて画像の再構成を行う関数patch_vecs2gray_imは、パッチのオーバーラップを許し、オーバーラップした部分は足し合わせた回数分割ることで平均をとっています。また、バックプロジェクションを行う関数back_projectionは、サイズが大きくスパースな行列Lを扱うので、scipy.matrixクラスで各引数を受け取る関数としています。ダウンサンプラー行列はget_downsampler_arr関数で求めています。

実験

まず辞書学習ですが、サイズ250×250の顔画像50枚から12×12のパッチ22000枚程度切り出して学習を行いました。また、辞書の基底数はパッチベクトルの次元と同じ144、学習時の\lambda\lambda = 5にしています。なお、辞書の学習繰り返し回数は8で打ち切っています(それ以上繰り返した学習辞書も作りましたが、生成される超解像画像の精度が逆に悪くなるという結果に・・・)。
超解像画像生成時は、\lambda=50、バックプロジェクション時のパラメータc = 0.8\mu = 0.9、パッチのずらし量1で実験しています。実験に使用する画像データはおなじみのレナさんです。

83×83の画像を249×249の画像へ(拡大率3倍)

  • 入力画像

f:id:YamagenSakam:20180411231300p:plain

f:id:YamagenSakam:20180411232337p:plain f:id:YamagenSakam:20180411231626p:plain

左の通常の拡大はPILのresizeで得られる画像です。少しわかりにくいですが、通常の拡大に比べて右の超解像画像は滑らかになっているのがわかると思います。

62×62の画像を248×248の画像へ(拡大率4倍)

  • 入力画像

f:id:YamagenSakam:20180411231846p:plain

f:id:YamagenSakam:20180411231342p:plain f:id:YamagenSakam:20180411231925p:plain

もう少しわかりやすくするため、更に拡大率を上げた画像を生成してみました。通常の拡大に比べると滑らかになってはいますが、ブロック的なノイズが乗ってしまっています。アルゴリズム的な限界なのか、辞書やパラメータを洗練すればもっとうまくいくのか、その辺りの検証は今後の課題です。

まとめ

今回は近接勾配法の応用例として辞書学習および更にその応用である超解像画像生成をやってみました。
実際やってみると、パッチのつなぎ目が汚くなったりノイズまみれになったりと結構苦労しました。その過程で色々わかったこと・思ったことを下記にまとめます(その分野では当たり前のことばかりかも知れませんが・・・)。また、やり残したこと、今後やりたいことも併せて記します。

  • やってみてわかったこと、思ったこと
    • 辞書学習時はパッチベクトルを平均{\bf 0}、分散1へ標準化した方がよい。
    • 画像再構成時はパッチ間の境界を滑らかにするため、元画像から切り出すパッチは領域を重複させて切り出し、再構成後のパッチを重ね合わせることで画像全体を再構成する方がよい。
    • 学習させすぎると精度が落ちるので、途中で学習を打ち切った方が超解像の結果がよい。
    • \lambdaによって結構結果が変わる。係数推定結果の非零要素数がパッチベクトルの要素数と大体同じになる程度が妥当か。
    • バックプロジェクションはあまり効かない?(ダウンサンプリング以外にボケなどを考慮した{\bf L}すると有効か?)
  • やり残したこと、今後やりたいこと
    • 学習用画像データのバリエーションを増やして実験。
    • 学習回数が増えると精度が落ちるのはなぜか数理的に考察。
    • 推定した画像の一部の要素が0~255の範囲を超えることがある。制約条件をつけられないか検討。
    • K-SVDなどの他の辞書学習手法の適用。
    • 独立成分分析やウェーブレットを用いた辞書行列の適用。
    • 他のダウンサンプリング辞書生成方法の適用。
    • 辞書学習以外の超解像手法との比較。


参考文献

*1:実務等で利用するときはscikit-learn辺りのライブラリを使うべきですが、今回は理解を深めるためnumpy、scipy、PILを使った実装を行っています。

*2:実は{\bf D}正規分布に従ったランダムな行列でも、符号化としてはある程度うまくいくという研究もあります。ただ超解像やその他画像再構成への応用という観点ではどうなんでしょう?

近接勾配法とproximal operator

はじめに

前回前々回とl1ノルム正則化項をつけた離散フーリエ変換により、スパースな解が得られること、そして基底ベクトルを並べた行列{\bf W}がユニタリ行列(直交行列)のため解析的に解けることを見てきました。それでは、{\bf W}が一般の行列の場合どうするか?今度こそ解析的に解けないので、反復法により最適解を求めます。しかし、目的関数にはl1ノルムという微分できない項があるため、最急降下法などのアルゴリズムは使えません。
このような微分できない項を含む最適化問題を効率的に解く手法の1つに近接勾配法(Proximal Gradient Method)があります。今回はこの近接勾配法と近接勾配法において重要な役割を果たすproximal operatorについて書いていきます。また、pythonで近接勾配法の実装を行いましたので、そのソースコードも載せておきます。

近接勾配法、proximal operatorとは

まずは、近接勾配法で解ける最適化問題の定義です。

\displaystyle\min_{{\bf x}} f({\bf x} )+ g({\bf x}) \tag{1}

ただし、f({\bf x})微分可能な凸関数、g({\bf x})微分可能とは限らない凸関数*1です。
l1ノルム正則化項付きの問題を式(1)に当てはめると、g({\bf x}) = \lambda \|{\bf x}\|_1となります。

次にproximal operatorの説明です。とりあえず定義を見てみましょう。関数g({\bf x})のproximal operatorは以下のように定義されます。

\displaystyle {\rm prox}_{\gamma g}({\bf v}) = \arg \min_{{\bf x}} g({\bf x}) + \frac{1}{2 \gamma} \|{\bf x} - {\bf v}  \|^2_2 \tag{2}

ここで\gammaは後述の更新式のステップ幅です。
このproximal operatorは何者か?というのは、g({\bf x})として以下のような指示関数を考えると見えてきます。

\begin{eqnarray}  
g({\bf x}) = I_c(\bf{x}) &=& \left\{  
\begin{array}{ll}
0 & ({\bf x} \in C) \\
\infty & ({\bf x} \notin C) \\ 
\end{array}
 \right.
\end{eqnarray}  \tag{3}

ここでCは閉凸集合です。この関数は{\bf x}が集合Cの要素なら0、Cの要素でないなら\inftyに吹っ飛ばすという関数です。そして、この指示関数のproximal operatorは以下のようになります。

\displaystyle \begin{eqnarray} {\rm prox}_{\gamma g}({\bf v})
&=& \arg \min_{{\bf x}}  I_c(\bf{x})  + \frac{1}{2 \gamma} \|{\bf x} - {\bf v}  \|^2_2 \\
&=&  \arg \min_{{\bf x} \in C} \|{\bf x} - {\bf v}  \|^2_2 \\
&=&   \Pi_C({\bf v})\\
\end{eqnarray}  \tag{4}

{\bf v}が集合Cの要素の場合、そのまま{\bf x} = {\bf v}とすれば\|{\bf x} - {\bf v}  \|^2_2を最小化できます。一方で、{\bf v}が集合Cの要素ではない場合、集合Cの要素で{\bf v}とのユークリッド距離を最小とするものが選ばれます。すなわち、{\bf v}の集合Cへの射影です(このような射影作用素\Pi_C({\bf v})として定義しています)。
よって、proximal operatorは射影作用素の関数gへの一般化と考えることができ、{\bf v}gによって定められる領域からはみ出していたら、その領域に入るように射影する作用素と見なせます。
proximal operatorの具体例を挙げると、g({\bf x})g({\bf x}) = \lambda\| \bf{x}\|_1というl1ノルムであれば、前回述べたソフトしきい値作用素がproximal operatorになります。つまり、

 \begin{eqnarray}  
{\rm prox}_{\gamma g}({\bf v}) = S_{\gamma \lambda} ({\bf v}) &=& \left\{
\begin{array}{ll}
 v_i - \gamma \lambda & (v_i \geq \gamma \lambda) \\
0 & (- \gamma \lambda  < v_i < \gamma \lambda) \\
v_i + \gamma \lambda & (v_i \leq  -\gamma  \lambda) \\
\end{array}
 \right.
\end{eqnarray}  \tag{5}

です。そして、近接勾配法の更新式はこのproximal operatorを使って以下のようになります。

\displaystyle{\bf x}_{t+1} = {\rm prox}_{\gamma g}({\bf x}_t - \gamma \nabla f({\bf x}_t)) \tag{6}

この更新式の収束性や正しさなどの説明は他の文献を参照いただくとして、直観的には次のような解釈ができると思います。
まず、{\bf x}_t - \gamma \nabla f({\bf x}_t)最急降下法の更新式です。また、上述のようにproximal operatorはgによって定められる領域に射影する作用素です。従って、最急降下法の1ステップで得られる値がgの領域からはみ出してたら、その領域へ射影するという操作を繰り返して最適解を探すアルゴリズムと見ることができます。
このように、関数f({\bf x})の勾配と関数g({\bf x})のproximal operatorがわかれば、近接勾配法により最適解が得られます。
なお、\gammaの決め方も収束のためには重要ですが、こちらについても別の文献を参照してください。

Moreau decompositionと双対ノルム

さて、ここまで見たように近接勾配法を使うためにはproximal operatorを求める必要があります。しかし、関数g({\bf x})によっては、式(2)から直接的にproximal operatorを求めるのが難しいことがあります。例えば、g({\bf x})=\lambda \| {\bf x} \|_2とl2ノルムとした場合(2乗が付いていないことに注意)、式(2)を解こうとして{\bf x}微分しても、1/\| {\bf x} \|_2がかかる項が出てきてしまうので厄介です。
そこで便利なのが以下に示すMoreau decompositionが成り立つというproximal operatorの性質です。

\displaystyle {\bf v} = {\rm prox}_{g}({\bf v}) +{\rm prox}_{g^*}({\bf v})  \tag{7}

ここで、{\rm prox}_{g^*}({\bf v})g({\bf x})の共役関数g^*({\bf x})のproximal operatorです。この共役関数とは下記の式で表される関数のことを言います。

\displaystyle g^*({\bf z}) = \sup_{{\bf x}} {\bf z}^T {\bf x} - g({\bf x})  \tag{8}

よって、この共役関数のproximal operatorさえ分かれば、元の関数のproximal operatorも求められます。
では、g({\bf x}) = \lambda\| \bf{x}\|というノルム関数の場合どうなるかを見ていきます。まずノルム関数の共役関数は以下の指示関数となります。

 \begin{eqnarray}  
g^*({\bf z}) = I_{\it B} \left({\bf z} \right) &=& \left\{  
\begin{array}{ll}
0 & (\|z \|_{*} \leq \lambda) \\
\infty & (\|z \|_{*} > \lambda) \\ 
\end{array}
 \right.
\end{eqnarray}\tag{9}

ここで、\| \bf{x}\|_{*}は双対ノルムと呼ばれるもので、以下の式で定義されます。

\displaystyle\| {\bf z} \|_{*} = \sup_{{\bf x}} {\bf z}^T {\bf x} \quad  ({\rm s.t.} \quad \| {\bf x} \| \leq 1)  \tag{10}

ということで、共役関数のproximal operatorは式(4)でみたように下記の射影作用素になります。

 {\rm prox}_{g^*}({\bf v}) = \Pi_{\it B} ({\bf v}) \tag{11}

ただし、{\it B}{\it B}=\left\{ {\bf x} \, | \,   \| \bf{x} \|_* \leq \lambda \right\}で表される集合です。
つまり、{\bf v}\| {\bf v} \|_* \leq \lambdaの範囲に射影する作用素がノルム関数の共役proximal operatorとなります。
おそらく式だけをザッと書いているので何故ノルムの共役関数が指示関数なのか?双対ノルムってそもそも何?どうやって求めるの?などの疑問が湧くかと思います。その辺りを少し補足しておきます。

双対ノルムについての補足

説明が前後しますが、後の説明のためにまずここで双対ノルムについて補足を書きます。そのための準備としてヘルダーの不等式という不等式を紹介します。
 1 \leq p,q \leq \inftyという2つの整数に、1/p + 1/q = 1という関係があれば、

{\bf z}^T{\bf x} \leq \| {\bf x}\|_p  \| {\bf z} \|_q \tag{12}

が成り立つというのがヘルダーの不等式です。実は、 \| \cdot \|_p の双対ノルムは \| \cdot \|_qであり、その逆も成り立ちます。そのため、l2ノルムの双対ノルムはそのままl2ノルムですし、l1ノルムの双対ノルムはl\inftyノルムです。
このことはヘルダーの不等式から示すことができます。まずヘルダーの不等式より、任意の\|{\bf x}\|_p \leq 1に対して\| {\bf z} \|_q \geq {\bf z}^T {\bf x}が成り立ちます。lpノルムの双対ノルム\|{\bf z} \|_*\|{\bf x}\|_p \leq 1という条件下での{\bf z}^T {\bf x}の最大値です。ここで{\bf x}をうまくとることにより、等式を成立させることができます。つまり、

\displaystyle \begin{eqnarray}
\| {\bf z} \|_{*} &=& \sup_{{\bf x}} {\bf z}^T {\bf x} \quad ({\rm s.t.} \quad \| {\bf x} \|_p \leq 1) \\
& =& \| {\bf z} \|_q \\
\end{eqnarray}
\tag{13}

が成立します。

何故、ノルム関数の共役関数が指示関数なのか?

まず、\| {\bf z} \|_* > \lambdaの時を考えます。式(8)の定義よりg^*({\bf z}) = \sup_{{\bf x}} {\bf z}^T {\bf x} - \lambda \| {\bf x} \| なわけですが、双対ノルムの定義式(10)と前述の前提より、 {\bf z}^T {\bf x} > \lambdaかつ\|{\bf x}\| \leq 1となる{\bf x}が存在することになります。これを満たす{\bf x}{\hat{\bf x}}とした場合、その{\hat{\bf x}}に定数tをかけたt {\hat{\bf x}}も式(8)の最適解の候補として入ってきます。これはもちろん、t \rightarrow \inftyとした場合も同様です。これらのことから、

 t{\bf z}^T {\hat {\bf x}} - \lambda \| t{\hat {\bf x}} \| =  t({\bf z}^T {\hat {\bf x}} - \lambda \| {\hat {\bf x}} \| ) \rightarrow \infty \tag{14}

となるため、\| {\bf z} \|_* > \lambdaのときはg^*({\bf z}) = \inftyとなります。
一方で、\| {\bf z} \|_* \leq \lambdaの時は、式(12)のヘルダーの不等式より{\bf z}^T{\bf x} \leq \| {\bf x}\| \| {\bf z} \|_* なので、

{\bf z}^T{\bf x}  - \lambda \| {\bf x}\| \leq \| {\bf x}\| \| {\bf z} \|_* -  \lambda  \| {\bf x} \| = \| {\bf x}\| ( \| {\bf z} \|_* - \lambda) \tag{15}

です。前提より、\| {\bf z} \|_* \leq  \lambdaなので、左辺(つまりはg^*({\bf z}) )が0以上になることはありません。ここで、{\bf x}= {\bf 0}という取り方もできるので、式(15)の左辺を0にすることができます。よって、\| {\bf z} \|_* \leq \lambdaのときはg^*({\bf z}) = 0となります。

結局のところproimal operatorはどうやって求めればよいか?

2つ方法があると思います。1つ目は式(2)から直接求める方法です。例えば、前回記したようの方法でl1ノルム正則化項のproximal operator式(5)を求めるというものです。
もう1つの方法としては、式(7)のMoreau decompositionの性質を使う方法です。特に正則化項がノルムなら、共役関数のproximal operatorは{\it B}=\left\{ {\bf x} \, | \,   \| \bf{x} \|_* \leq \lambda \right\}という超球への射影なので簡単に求められます。例えば、l2ノルムの双対ノルムはl2ノルムなので、

\displaystyle\begin{eqnarray}
{\rm prox}_{g^{*}} ({\bf v})= \Pi_{\it B}({\bf v}) = \left\{
\begin{array}{ll}
\lambda \frac{\bf{v}}{\| \bf {v} \|_2} \quad (\|  {\bf v} \|_2 > \lambda) \\
{\bf v}  \quad (\|  {\bf v} \|_2  \leq \lambda) \\
\end{array}
\right. 
\end{eqnarray}
\tag{15}

であり、求めたいproximal operatorは、

\displaystyle\begin{eqnarray}
{\rm prox}_{g} ({\bf v})= {\bf v} - \rm {prox}_{g^*}  ({\bf v})= \left\{
\begin{array}{ll}
{\bf v} -\lambda  \frac{{\bf v}}{\| {\bf v} \|_2} \quad (\|  {\bf v} \|_2 > \lambda) \\
0 \quad (\|  {\bf v} \|_2  \leq \lambda) \\
\end{array}
\right. 
\end{eqnarray}
\tag{16}

となります。

なお、式(2)にはステップ幅の\gammaが付いています。\gammaが付いていても基本的には上式の\lambda\gamma \lambdaに置き換えればよいだけですが、\gammaを含めたMoreau decompositionの一般形を書いておきましょう。

\displaystyle {\bf v} = {\rm prox}_{\gamma g}({\bf v}) +\gamma{\rm prox}_{g^*/\gamma} \left(\frac{{\bf v}}{\gamma} \right)  \tag{17}

という形式になります。したがって、gがノルム関数の場合は、

\displaystyle {\rm prox}_{\gamma g}( {\bf v}) = {\bf v} - \gamma{\rm prox}_{g^*/\gamma} \left(\frac{{\bf v}}{\gamma} \right)  =  {\bf v} - \gamma \Pi_{\it B} \left(\frac{{\bf v}}{\gamma} \right) \tag{18}

になります。ちなみに、多くの文献ではg({\bf x })=\| {\bf x}\|としているので、{\it B}=\left\{ {\bf x} \, | \,   \| \bf{x} \|_* \leq 1 \right\}として理論を展開している点に注意してください。

様々な正則化項のproximal operator

ここまでノルム関数が正則化項の時のproximal operator導出方法を見てきました。ここで、いくつかの正則化項に対するproximal operatorを記載しておきます。

l2ノルムの2乗正則化g({\bf v}) = \lambda \| {\bf v} \|^2_2

リッジ回帰等でよく使われる正則化項です。普通にproximal operatorの式を微分して0とおけば求められます。2乗しているので、以下のように式(16)とは異なる結果が得られます。

\displaystyle prox_{\gamma g} ({\bf v})= \frac{1}{1+2\gamma\lambda} {\bf v} \tag{19}

Elastic Net:g({\bf v}) = \lambda_1 \| {\bf v} \|_1 + \lambda_2 \| {\bf v} \|^2_2

l1ノルム正則化とl2ノルム2乗正則化の混合です。l1ノルムは2つ変数が同程度の寄与度だったとしてもどちらかしか選ばれない(一方が0になってしまう)ケースが多いですが、Elastic Netではこの点を改善できるようです。このproximal operatorはl1ノルム正則化の時と同様、変数ごとに分離して求めることができます。

\displaystyle prox_{\gamma g} ({\bf v})= \frac{1}{1+2\gamma\lambda_2} S_{\gamma \lambda_1} ({\bf v}) \tag{20}

グループl1ノルム(グループに重複がない場合):\displaystyle g({\bf v}) = \lambda \sum_{g \in G}  \| {\bf v}_g \|_p

変数がいくつかのグループに分けられるとき、それぞれのグループのlpノルムを取って、更にその和を取るノルムです( ベクトル{\bf v}のうち、グループgに所属する変数を並べたのが{\bf v}_g)。これを正則化項に用いるとグループごとにスパースになる、ならないが決まります。このグループl1ノルムを正則化項として使った回帰をgroup lassoと言います。
これはグループに重複がなく完全に分離ができるため、グループごとにlpノルムのproximal operatorを求めればよいです。p=2の時のグループgのproximal operatorは下記の通りです。(式(16)とほぼ同じです。)

\displaystyle\begin{eqnarray}
{\rm prox}_\gamma{g} ({\bf v}_g)= \left\{
\begin{array}{ll}
{\bf v}_g  -\gamma \lambda  \frac{{\bf v}_g}{\| {\bf  v}_g \|_2} \quad (\|  {\bf v}_g \|_2 > \gamma \lambda) \\
0 \quad (\|  {\bf v}_g \|_2  \leq \gamma \lambda) \\
\end{array}
\right. 
\end{eqnarray}
\tag{21}

ちなみに、グループに重複がある場合も提案されていますが、ここでは割愛します(まだ追いきれてない)。

トレースノルム :g({\bf Z}) = \lambda \| {\bf Z} \|_{{\rm trace}} = {\rm Tr}\left(\sqrt{{\bf Z}^T {\bf Z}}\right)

最後は行列関数に関する正則化項です。トレースノルムは\sqrt{{\bf Z}^T {\bf Z}}のトレースを正則化項とするものなので、{\bf Z}の特異値を\sigma({\bf Z})=\left[\sigma_1, \sigma_2, \cdots, \sigma_m \right]とすると、 \| {\bf Z} \|_{{\rm trace}} = \sum_{i=1}^{m} \sigma_iです。そのため、このトレースノルムを正則化項として使うと特異値が0になりやすく、結果低ランクな解が得られるようになります。
それではproximal operatorを見ていきます。{\bf Z}特異値分解{\bf Z} = {\bf U} {\rm diag}(\sigma ({\bf Z})) {\bf V}^Tとすると、proximal operatorは以下のようになります。

\displaystyle {\rm prox}_{\gamma g} ({\bf Z}) =  {\bf U} {\rm diag}(S_{\gamma \lambda}(\sigma ({\bf Z}))) {\bf V}^T \tag{22}

特異値分解して特異値に対してソフトしきい値作要素を施した後に再構築すればよいわけです。

近接勾配法を実装

ということで式(1)を解く近接勾配法を実装しました。

近接勾配法を行う関数

ということで、まずは近接勾配法そのものの実装です。

def proximal_gradient(grad_f, prox, gamma, objective, init_x, tol = 1e-9):
    x = init_x
    result = objective(x)
    while 1:
        x_new = prox(x - gamma * grad_f(x), gamma)
        result_new = objective(x_new)
        if (np.abs(result - result_new)/np.abs(result) < tol) == True :
            break;
        x = x_new
        result = result_new
    return x_new, result_new

引数について説明します。
grad_fは\nabla f({\bf x})にあたります。つまり、grad_fはxを引数として、その勾配を返す関数です。
proxはproximal operatorです。proxも関数で、入力変数vとステップ幅のgammaの2引数をとり、proximal operatorの計算結果を返します。
objectiveは目的関数です。objectiveはxを引数として、目的関数の計算結果を返す関数です。
init_xは初期値、tolは収束の許容誤差です。
やっていることは、式(6)の更新式をそのまま素直に実装しているだけです。

l1ノルム正則化項のproximal operator

次にproximal operatorの実装です。

def prox_norm1(v, gamma, lam):
    return soft_thresh(v, lam * gamma)

def soft_thresh(b, lam):
    x_hat = np.zeros(b.shape[0])
    x_hat[b >= lam] = b[b >= lam] - lam
    x_hat[b <= -lam] = b[b <= -lam] + lam
    return x_hat

prox_norm1がl1ノルムのproximal operatorを計算する関数です。補助としてsoft_threshという関数を作っています。

lasso回帰問題を解く関数

以下は、\min_{{\bf x}}\| {\bf y} - {\bf W}{\bf x} \|_2^2 + \lambda \| {\bf x}\|_1というlasso回帰問題を解く関数です。なお、l1ノルム正則化項に対する近接勾配法はISTA(Iterative Shrinkage Thresholding Algorithm)と呼ばれるため、関数名をISTAとしています。

def ISTA(W, y, lam):
    objective = lambda x:np.sum(np.power(y - np.dot(W, x), 2))/2.0 + lam * np.sum(np.abs(x)) #(1)
    grad_f = lambda x: np.dot(W.T, np.dot(W, x) - y)  #(2)
    (u, l, v) = np.linalg.svd(W)   #(3)
    gamma = 1/max(l.real*l.real)   #(4)
    prox = lambda v, gamma:prox_norm1(v, gamma, lam)  #(5)
    x_init = randn(W.shape[1])  #(6)
    (x_hat, result) = proximal_gradient(grad_f, prox, gamma, objective, x_init, 1e-5)  #(7)
    return x_hat, result

(1)ではラムダ式を使ってobjectiveをxを引数としてとり、上記の目的関数を計算する関数として定義しています。
また、\nabla f({\bf x})={\bf W}^T ({\bf W} {\bf x} - {\bf y})であり、(2)はこの計算を行うgrad_fをxに引数とる関数として定義しています。
(3)、(4)は\gammaを定めています。最大特異値の逆数より小さい値を\gammaとして定めると収束するようですが、その辺りはまだ理解ができておりません。とりあえず、文献の通りに実装します。
(5)がproximal operatorです。prox_norm1関数はv,gamma,lamの3引数を受け取りますが、proximal_gradient関数に渡すproxはvとgammaの2引数のみ取る関数なので、ラムダ式を使って再定義しています。
(6)はランダムに初期値を決めて、(7)でproximal_gradient関数にて問題を解きます。

proximal operatorを別なものにする

上記ISTA関数のproximal operatorとobjectiveの部分だけ変えれば他の正則化にも対応可能です。例としてグループl1ノルム正則化(group lasso)の場合を示します。
まずは、グループl1ノルムのproximal operatorの実装です。

def prox_group_norm(v, groups, gamma, lam):
    u = np.zeros(v.shape[0])
    for i in range(1, max(groups) + 1):
        u[groups == i] = block_soft_thresh(v[groups == i] , gamma * lam)
    return u

def block_soft_thresh(b, lam):
    return max(0, 1 - lam/np.linalg.norm(b)) * b

groupsは1~G(Gはグループ数)を要素に持つリストで、各変数がどのグループに属しているかを表します。後はグループごとにblock_soft_thresh関数にて、l2ノルムのproximal operatorを計算しています。
そして、上記のISTA関数内の(5)でproxとして、prox_group_normを指定すればgroup lassoになります。

    prox = lambda v, gamma:prox_group_norm(v, groups, gamma, lam) #(5)'

もし、Elastic Netを使いたかったらここにElastic Netのproximal operatorを持ってくればよいし、他のノルムの場合も同様です(本当はobjectiveも変更する必要がありますが上記では省略しています)。
また、f({\bf x})を別なものにしたかったら\nabla fを別途定義して、grad_fとすればよいです。
行列関数の場合も同様で、grad_fとproximal operator (あと、objective)さえ定義してあげたら機能します。そのため、proximal operatorとしてトレースノルムを指定してあげると、低ランク解を得ることができます。

まとめ

今回は近接勾配法とproximal operatorについて見てきました。上記の実装については次の機会に応用事例で実験したいと思います(ちなみに、疑似データで実験したところそれっぽく動いていました)。
また、今回は取り上げませんでしたが、この界隈ではISTAの高速版FISTAや拡張ラグランジュ法、その改良であるAlternating Direction Method of Multipliers (ADMM) など様々なアルゴリズムが提案されており、非常に奥が深いです。

参考文献

最後にこの記事を書くにあたって参考にした文献です(基本的に前回とほぼ同じです)。

また日本語の書籍で、proximal operator界隈が詳しく書かれているのは以下2冊かなと思います(他にもあったら教えてください)。

*1:もうひとつ、「単純」という条件も付きます。

l1ノルム正則化フーリエ変換を複素数のまま解く

はじめに

前回の記事では、離散フーリエ変換を係数をLasso回帰として求めスパースな解を得る方法を書きましたが、この際に複素数を取り扱いたくないため実数のみの問題に変換して、ソルバに解かせるということを行いました。今回は、複素数のまま問題を解くことを考え、複素数の最適解を得る方針で考えます。この方針で考察することで、以下のことが見えてきたので記事に残そうと思った次第です。

  • 前回の記事では「この問題は解析的には解けず」などと書いているが、実は{\bf W}が離散フーリエの基底ベクトルを並べた行列であれば解析的に解けること
  • 結局は、離散フーリエ変換してしきい値未満の要素を0にしていることと(ほぼ)等価であること
  • 複素数としてのl1ノルムを考えると、別な解が得られること

ちなみに、自分は「この方法だとなくうまくいきそうだな」と思ったことはまず試してみて、後から理屈を考える(特にうまくいかなかったとき)というスタンスなので、理論的、数学的な厳密性はあまり考慮していません。したがって、本来理論的にはあり得ない式展開、考え方をしている場合もあるかもしれませんがご了承ください。

定式化

前回と同様、正則化項をつけた回帰問題として定式化します。

\displaystyle \min_{\bf x} \| {\bf y} - {\bf W}{\bf x}\|^2_2 + \lambda \Omega({\bf x}) \tag{1}

{\bf y}は元の時間軸の信号で、{\bf y} = \left[ y(0), y(1), \cdots, y(T - 1) \right]^T{\bf x}が求めたい離散フーリエ係数であり、{\bf x} =  \left[ x(0), x(1), \cdots, x(F - 1) \right]^Tです。ここで{\bf x}の各要素は複素数です。{\bf y}は実信号なので実数ベクトルですが、虚部がすべて0の複素数と見なしてもよいでしょう。また、行列{\bf W}の各列ベクトル{\bf w}_k ( k = 0, 1, \cdots, T-1)
{\bf w}_k = \frac{1}{T}\left[\exp(i\frac{2\pi k \times (0)}{T}) ,\exp(i\frac{2\pi k \times (1)}{T}), \cdots, \exp(i\frac{2\pi k \times (T - 1)}{T} ) \right]^T
とします。前回の定義とほぼ一緒ですが、各要素を要素数Tで割っています。こうすると、{\bf W}はユニタリ行列になります。つまり、{\bf W}^* {\bf W} = {\bar {\bf W}}^T {\bf W}  = {\bf I}です({\bf W}^*は随伴行列、 {\bar {\bf W}}複素共役を表します)。
後は正則化\lambda \Omega({\bf x})です。普通に考えると
\Omega({\bf x}) = \| {\bf x} \|_1 = | x(0)| + |x(1)|+ \cdots + |x(T-1)| \tag{2}

(ただし、|x|xの絶対値を表し |x|  = \sqrt{{\bar x} x})としたいところですが、実はこの定義だと前回と同じ結果にはなりません。このことは後ほど説明するとして、今は
\Omega({\bf x}) = | {\rm  Re}( x(0))| +  | {\rm  Im}( x(0))| + \cdots +  | {\rm  Re}( x(T - 1))| +  | {\rm  Im}( x(T - 1))|  \tag{3}

とします({\rm Re(x)}xの実部、{\rm Im(x)}xの虚部を表す)。式(2)と式(3)の大きな違いは、実部と虚部が関連し合っているか、完全に独立なものとして扱うか、です。そのため、式(3)を正則化項として用いることは、実部と虚部を分けて解を求める前回の考え方と一致します。

解析解の導出と実装

当初は「スパースモデリングのための凸最適化」の解説論文を読んで近接勾配法で解こうなんて考えていましたが、解説を読んでいくと実は、{\bf W}がユニタリ行列であるため解析的に解が得られることがわかりました。
ということでその導出です。まず、解説に従うと、式(1)は以下のように式を変形できます。
\begin{eqnarray} 
\displaystyle (1) &=& \min_{\bf x} \| {\bf y} - {\bf W}{\bf x}\|^2_2 + \lambda \Omega({\bf x})  \\
\displaystyle  &=& \min_{\bf x} \| {\bf W}^*{\bf y} -  {\bf W}^*{\bf W}{\bf x}\|^2_2 + \lambda \Omega({\bf x}) \\
\displaystyle  &=&  \min_{\bf x} \| {\bf b} -  {\bf x}\|^2_2 + \lambda \Omega({\bf x}) \ \ \  (\because {\bf b} =  {\bf W}^*{\bf y})\\
\displaystyle  &=&  \min_{\bf x}  \sum_{i = 0}^{T - 1} \left\{ |(b(i) -  x(i))|^2+ \lambda( | {\rm  Re}( x(i))| +  | {\rm  Im}( x(i))| )  \right\} \\
\displaystyle  &=&  \min_{\bf x}  \sum_{i = 0}^{T - 1} \left[  \left\{ {\rm Re} (b(i)) -  {\rm Re} ( x(i) ) \right\} ^2 + \lambda | {\rm  Re}( x(i))|   +  \left\{ {\rm Im} (b(i)) -  {\rm Im} ( x(i) ) \right\} ^2 + \lambda | {\rm  Im}( x(i))|   \right]   \\ 
\displaystyle  &=&\sum_{i = 0}^{T - 1}   \min_{x(i)}   \left[  \left\{ {\rm Re} (b(i)) -  {\rm Re} ( x(i) ) \right\} ^2 + \lambda | {\rm  Re}( x(i))|   +  \left\{ {\rm Im} (b(i)) -  {\rm Im} ( x(i) ) \right\} ^2 + \lambda | {\rm  Im}( x(i))|   \right]    
\end{eqnarray}
このように式(1)の問題は各要素独立に最適解を選ぶ問題に変形されます。しかも、実部と虚部も分離しているので、こちらについても独立して最適解を選ぶことができます。ようするに、以下の最適化問題を各要素の実部、虚部それぞれ解けばよいわけです。
\displaystyle\min_{x}    (b-  x)^2 + \lambda | x|  \tag{4}
この解は次のような式で得られます。
\begin{eqnarray} 
{\hat x} = S_{\lambda} (b) &=& \left\{
\begin{array}{ll}
 b - \frac{\lambda}{2} & (b \geq \frac{\lambda}{2}) \\
 0 & (- \frac{\lambda}{2}  < b < \frac{\lambda}{2}) \\
b + \frac{\lambda}{2} & (b \leq  - \frac{\lambda}{2}) \\
\end{array}
 \right.
\end{eqnarray} \tag{5}
解説論文と異なり式(1)および式(4)の第1項に1/2がかかっていないため、最適解も\lambda1/2かけていることに注意してください。ここで式(5)のS_{\lambda}(b)はソフトしきい値作用素と呼ばれる作用素です。
このソフトしきい値作用素を用いると、今回得たい式(1)の最適解は下記のように表せます。
{\hat{\bf x}} = S_{\lambda}({\rm Re}({\bf W}^* {\bf y}))+ i S_{\lambda}({\rm Im}({\bf W}^* {\bf y})) \tag{6}

さて、ここで通常の離散フーリエ変換の係数は{\bf W}^* {\bf y}で得られることを思い出すと、式(6)は通常の離散フーリエ変換の結果の実部、虚部それぞれに対してソフトしきい値作用素を作用させていることになります。
これは式(5)により離散フーリエ係数の絶対値が\lambda/2未満であれば0にすることに対応します。l1ノルム正則化項の離散フーリエ変換なんて言っていますが、紐解いていくと離散フーリエの結果に対ししきい値未満ならば0にするという割りとよくやる操作とほぼ同じになるわけですね(しきい値を超えている要素には\pm \lambda/2を加える点は異なりますが)。
ということで、以下は上式によるスパース離散フーリエ変換を行う関数の実装です。

def sparse_dft_with_thresh(y, T, lam):
    Z = np.meshgrid(np.arange(T),np.zeros((T,1)))[0]
    u = np.arange(T)
    V = (Z.T * u) * (2 * np.pi /T)
    W= np.exp( 1j * V)/T
    return (soft_thresh((np.dot(W.T.conj(), y)).real, lam/2)
            +1j * soft_thresh((np.dot(W.T.conj(), y)).imag, lam/2))

def soft_thresh(b, lam):
    x_hat = np.zeros(b.shape[0])
    x_hat[b >= lam] = b[b >= lam] - lam
    x_hat[b <= -lam] = b[b <= -lam] + lam
    return x_hat

sparse_dft_with_threshがスパースな離散フーリエ変換を行う関数で、soft_threshがソフトしきい値作用素を施す関数です。soft_threshを解説と合わせるために、引数で渡す時にlam/2をしています。

複素数l1ノルム正則化

ここで1つ疑問が残ります。正則化として式(2)のノルムを使うとどうなるか?です。これまでの議論から、式(2)のノルムを使った場合、式(1)は次のように変形できます。
\displaystyle \sum_{i = 0}^{T - 1} \min_{x(i)}  \left\{ |(b(i) -  x(i))|^2+ \lambda| ( x(i)| \right\} \tag{7}
この式を見ると要素同士は独立して最適解を選ぶことができますが、実部と虚部は独立させることができません(|x(i)|  = \sqrt{{\bar x(i)} x(i)}なので)。
この最適解を得るために、まず以下のような式を考えます。
\displaystyle {\rm prox}_{\lambda}({\bf b})=  \min_{\bf x}  \frac{1}{2} \| {\bf b} -  {\bf x} \|_2^2+ \lambda \Omega({\bf x}) \tag{8}
この式(8)は近接作用素(proximal operator)と呼ばれるもので、近接勾配法などで重要な役割を持つ作用素です。今、式(7)について複素数を実部と虚部を並べた2次元ベクトルとして考えると(つまり、{\bf x} = \left[{\rm Re}(x) , {\rm Im}(x)\right]^T )、式(7)と式(8)の問題は(第1項の1/2を除いて)等価になります(\Omega({\bf x}) = \| x \|_2)。
ということで、式(8)を解くことを考えますが、実はこれはgroup lassoにおけるproximal operatorを求めることと等価です。group lassoのproximal operatorの導出やその基礎となる劣微分、凸解析、双対理論などは「ptimization with sparsity-inducing penalties」「Proximal Algorithms」にまとまっています。正直、まだ理解できていない部分が多いですが、これらの資料に記載されているgroup lassoにおけるproximal operatorを使うと式(7)の解は下記のようになります(たぶん)。
\begin{eqnarray} 
{\hat x} = S_{\lambda} (b) &=& \left\{
\begin{array}{ll}
 b -   \frac{\lambda b}{2|b|}& (|b| > \frac{\lambda}{2}) \\
0 & (|b| \leq  - \frac{\lambda}{2}) \\
\end{array}
 \right.
\end{eqnarray} \tag{9}

式(6)では、実部、虚部それぞれでしきい値未満だと0になりましたが、こちらは複素数の絶対値がしきい値未満だと0になります。つまり、実部と虚部は運命共同体で、しきい値未満となる要素は実部も虚部一緒に0になり、その逆もまたしかりです。ということで実装です。

def group_lasso_prox(b, lam):
    x_hat = np.zeros(b.shape[0], "complex")
    x_hat[np.abs(b) >= lam] = (b[np.abs(b) >= lam] -
                                lam *  b[np.abs(b) >= lam]/np.abs(b)[np.abs(b) >= lam])
    return x_hat

def sparse_dft_with_complex_l1_norm(y, T, lam):
    Z = np.meshgrid(np.arange(T),np.zeros((T,1)))[0]
    u = np.arange(T)
    V = (Z.T * u) * (2 * np.pi /T)
    W= np.exp( 1j * V)/T
    return (group_lasso_prox(np.dot(W.T.conj(), y), lam/2))

一応、実験した結果も載せておきます。実験には前回同様sin波とcos波とホワイトノイズを重ね合わせた疑似データを用います。

def main():
    T = 2**10
    y = (3 * np.cos(28* 2*np.pi * np.arange(T)/T) +2 * np.sin(28* 2*np.pi * np.arange(T)/T) +
            np.cos(400* 2*np.pi * np.arange(T)/T) + np.sin(125 * 2*np.pi * np.arange(T)/T) + 5 * randn(T))

    x1 = sparse_dft_with_thresh(y, T, 0.5)
    plt.figure()
    plt.xlim((0, T))
    plt.xticks(np.arange(0, T , 100))
    plt.ylim((-2, 2))
    plt.plot(x1.real, "r")
    plt.hold(True)
    plt.plot(x1.imag, "b")

    x2 = sparse_dft_with_complex_l1_norm(y, T, 0.5)
    plt.figure()
    plt.xlim((0, T))
    plt.xticks(np.arange(0, T , 100))
    plt.ylim((-2, 2))
    plt.plot(x2.real, "r")
    plt.hold(True)
    plt.plot(x2.imag, "b")

    plt.show()


実部と虚部を独立の変数として扱った正則化の結果
f:id:YamagenSakam:20180204165656p:plain


複素数l1ノルムによる正則化の結果
f:id:YamagenSakam:20180204165749p:plain


こうやって見ると、前者の方がスパース性という意味だと高いです。実部、虚部を運命共同体にする理由がない限り、前者のように独立して扱った方がよいように思えます。

最後に

ここまで見て感じたかと思いますが、結局のところ単なるしきい値での枝刈りに帰着するので、実装まで落としこんでしまうと正直あまり面白いものではありませんし、応用としてもあまり意味のあるものでもありません(これやるくらいなら、普通にFFTしてしきい値未満を0にする方が断然早いです)。
ただ、ここまで導出する過程で、ソフトしきい値作用素やproximal operator、劣微分、凸解析などなど最適化理論の色々なトピックスに触れることができました(まだ全然理解できていないですが)。
冒頭に書いた思いついたことはまず試してみるというスタンスは、試すために知らないトピックスを調べるきっかけになるというメリットがあります。やってみてダメなら何がダメか更に深追いするモチベーションにもなります。後はこういう試行錯誤は何より楽しいですね。ということで、今後は近接勾配法や凸解析の理論などの知識を深めていければと思います。

l1ノルム正則化でスパースな離散フーリエ変換

はじめに

離散フーリエ変換は離散化された時間軸の信号を周波数軸に写像する基底変換の一種で、周波数解析などによく使われます。この離散フーリエ変換の係数は、基底ベクトルの重み付き線形和で表したときの誤差の2乗和を最小にする重みと見なせます。つまり、離散フーリエ係数は一種の線形回帰問題の解として得られます。ということは、l1ノルム正則化項つけたLasso回帰として解けばスパースな離散フーリエ変換ができるのでは?ということでやってみました。

離散フーリエ変換の線形回帰的定式化

元の信号のデータ点数をT個、基底ベクトルの個数をF個とし、上述の通り回帰問題として離散フーリエ係数を求める問題を定式化すると下記のようになります。

\displaystyle \min_{\bf x} \| {\bf y} - {\bf W}{\bf x}\|^2_2 \tag{1}

ここで、{\bf y}は元の時間軸の信号で、{\bf y} = \left[ y(0), y(1), \cdots, y(T - 1) \right]^Tです。
{\bf W} = \left[{\bf w}_0, {\bf w}_1, \cdots, {\bf w}_{F - 1} \right]F個の基底ベクトルを劣ベクトルとして並べたT \times Fの行列です。
一般的に離散フーリエ変換では{\bf W}{\bf w}_k = \left[\exp(i\frac{2\pi k \times (0)}{T}) ,\exp(i\frac{2\pi k \times (1)}{T}), \cdots, \exp(i\frac{2\pi k \times (T - 1)}{T} ) \right]^Tという複素数の要素を持つ行列で、各列ベクトル同士は直交します({{\bf w}_i^{*} {\bf w}_j} = 0 \quad i \neq j)。
そして、\bf{x}が求めたい離散フーリエ係数であり、{\bf x} =  \left[ x(0), x(1), \cdots, x(F - 1) \right]^Tです。
通常、F = Tです。そのため、(1)の最適解は{\hat{\bf x}} = {\bf W}^{-1} {\bf y} =  {\bf W}^{*} {\bf y}で表され、最小2乗誤差は0になります。

l1ノルムによるスパース化

ということで、回帰問題としてみなせるならl1ノルムの正則化項をつける(つまりは、Lasso回帰を行う)とスパースな解が得られるだろうという考えです。(1)の目的関数に正則化項を付け加えます。

\displaystyle \min_{\bf x} \| {\bf y} - {\bf W}{\bf x}\|^2_2 + \lambda \|  {\bf x} \|_1 \tag{2}

こうなると最早この問題は解析的には解けず、正則化項は\|{\bf x} \|_1 \leq \gammaという制約条件を付けるのと等価なので、最適解による最小2乗誤差も一般的には0にはなりません。ということで、cvxpyというソルバを用いて上記最適化問題を解きます。
なお、この定式化は入口こそ離散フーリエ変換を回帰問題と見なしてLasso回帰を行うという考えですが、結局のところこれは、この資料等に記載されいるノイズあり圧縮センシングの一種です。 

複素数最適化問題から実数のみの最適化問題

離散フーリエ変換は上述のように、基底ベクトルが複素数で得られる解も一般的に複素数となります。ところが、このままだと実装上面倒くさい(というより、ソルバが動いてくれるかわからない)ので実数のみの式に変換します。といっても、\exp({i \theta}) = \cos\theta + i\sin\thetaより実部と虚部に分けてそれぞれ個別の変数として考えるだけです。
具体的には、{\bf W} {\bf W} = \left[ {\bf w}_0^c, {\bf w}_1^c, \cdots, {\bf w}_{T - 1}^c,  {\bf w}_0^s, {\bf w}_1^s \cdots, {\bf w}_{T - 1}^s  \right]というT \times 2Tの行列として再定義します。ここで、
{\bf w}_k^c = \left[\cos(\frac{2\pi k \times (0)}{T}) ,\cos(\frac{2\pi k \times (1)}{T}), \cdots, \cos(\frac{2\pi k \times (T - 1)}{T} )\right]^T
{\bf w}_k^s = \left[\sin(\frac{2\pi k \times (0)}{T})  ,\sin(\frac{2\pi k \times (1)}{T}), \cdots, \sin(\frac{2\pi k \times (T - 1)}{T} )\right]^T
です。あわせて説明変数{\bf x}を要素数2Tのベクトルとして式(2)を解けばよいです。この変換で基底ベクトルの数も説明変数の要素数も倍になったように見えますが、元々も複素数でそれぞれが実部と虚部を持っていたため本質的には変わりません。
また、離散フーリエ変換の場合、x(k) + {\bar x(k + \frac{T}{2}})という関係があるため、実質0\sim\frac{T}{2}-1番目の要素までを求めれば、残りの\frac{T}{2}\sim T-1番目は自動で決定されます。これは、l1ノルム正則化項をつけても変わりません(実際に試しました)。そのため、説明変数の要素数2T\rightarrow Tへとシュリンクすることも可能ですが、今回は理解を深めるためにもそのまま解くようにしました。

cvxpyを用いた実装

上述の通り、今回は最適化問題を解くためにcvxpyというツールを利用しました。このcvxpyの使い方はこちらの記事を参考にしました。cvxpyの使い方は他の解説記事に譲るとして、早速式(2)を解く関数のソースコードを示します。

from cvxpy import *
import numpy as np
from numpy.random import *
import matplotlib.pyplot as plt

#y:原信号、T:信号のデータ点数、lam:ハイパーパラメータλ
def sparse_dft(y, T, lam):
    #行列Wを生成
    V = np.meshgrid(np.arange(T),np.zeros((T,1)))[0]
    u = np.arange(T)
    V = (V.T * u) * (2 * np.pi /T)
    W = np.c_[np.cos(V), np.sin(V)]

    #変数xを定義
    x = Variable(T * 2)

    #ハイパーパラメータを定義
    lamda = Parameter(sign="positive")
    lamda.value =lam

    #目的関数を定義
    objective = Minimize(sum_squares((y - W * x))/T + lamda*norm(x, 1))
    p = Problem(objective)

    #最適化計算
    result = p.solve()

    return result, x

特に特別なことはしておらず、最初に{\bf W}を生成して、後はcvxpyの解説に従って変数、ハイパーパラメータ、目的関数を定義し、p.solve()で問題を解きます。

実験

以下のようなコードで疑似データを発生させて実験してみました。

    T = 2**10
    y = (3 * np.cos(28* 2*np.pi * np.arange(T)/T) +2 * np.sin(28* 2*np.pi * np.arange(T)/T) +
            np.cos(400* 2*np.pi * np.arange(T)/T) + np.sin(125 * 2*np.pi * np.arange(T)/T) + 5 * randn(T) )

いくつかの周波数が異なるsin波とcos波を重み付きで重ね合わせて、更に正規分布に従うノイズを乗せた1024点のデータです。疑似データの波形は以下の通り。
f:id:YamagenSakam:20180129000732p:plain

通常の離散フーリエ変換

まずは上記の疑似データにFFTにより離散フーリエ変換を行いました。

    Y = np.fft.fft(y)/T #スケーリングを合わすため要素数で割る
    plt.figure()
    plt.xlim((0, T))
    plt.xticks(np.arange(0, T , 100))
    plt.ylim((-1.5, 1.5))
    plt.plot(Y.real , "r")
    plt.hold(True)
    plt.plot(Y.imag , "b")

結果は下記の図です。
f:id:YamagenSakam:20180129000744p:plain
重ねたsin波、cos波に対応する周波数成分が大きな値になるのは当たり前ですが、ノイズが重畳しているせいで全周波数成分まんべんなく値が入ります。

Lassoによる離散フーリエ変換

次は今回記した手法での離散フーリエ変換\lambda = 0.5で実験)。上述のsparse_dft関数で計算します。

    (result, x) = sparse_dft(y, T, 0.5)

    plt.figure()
    plt.xlim((0, T))
    plt.xticks(np.arange(0, T , 100))
    plt.ylim((-1.5, 1.5))
    plt.plot(x.value[0:T,0], "r")
    plt.hold(True)
    plt.plot(x.value[T:T*2,0], "b")

f:id:YamagenSakam:20180129000758p:plain
通常の離散フーリエ変換とは異なりスパースな解が得られ、余分な周波数成分は0になっていることがわかります。

逆変換

得られたフーリエ係数から逆変換を行い元の信号を再構築してみました。こうすることで、余計なノイズ成分が除去された信号が得られるはずです。
なお、逆変換のため、得られた{\bf x}から\cosにかかる係数を実部、\sinにかかる係数を虚部として複素数に変換してから、ifftメソッドで逆変換を行っています。

    X=(x.value[0:T] + (x.value[T:T*2] * 1j))*T
    y_hat = np.fft.ifft(X.T)

    plt.figure()
    plt.xlim((0, T))
    plt.xticks(np.arange(0, T , 100))
    plt.ylim((-20, 20))
    plt.plot(y_hat.real.T, "r")

以下が原波形から今回の方法でノイズを除去したデータです。
f:id:YamagenSakam:20180129000809p:plain
また、以下が上述の疑似データ生成のコードからノイズ項を取り除いたデータの波形です(つまり、真に得たい波形)。
f:id:YamagenSakam:20180129002851p:plain
確かにノイズが除去された波形が得られています。実際のノイズ項を除いた信号と比べると若干スケーリングが小さくなっていますが、悪くない結果だと思います。

まとめ

離散フーリエ変換を回帰問題と見なし、目的関数にl1ノルムの正則化項をつけて誤差の最小化を行うことで、スパースな離散フーリエ係数が得られました。またそれを逆変換することで原信号からノイズを除去した信号を得ることもできました。また、最適化問題を解くためにcvxpyというツールを利用しました。
今回は何も考えずcvxpyに問題を突っ込んで解いただけでしたが、今後は、Lasso回帰の最適化等で使われる近接勾配法やcvxpyの詳しい使い方、凸解析など最適化の理論、実装両面で理解を深めていきたいと思います。

はてなダイアリーからはてなブログへ移行

このたび、本ブログをはてなダイアリーからはてなブログへ移行しました。旧はてなダイアリーの記事にアクセスすると、対応する本ブログの記事へのリダイレクトされるようなので、ブックマーク登録等の変更は不要です。
さて、この移行作業に当たって、元のダイアリーで書いていた数式がうまく表示されず修正に手間取った箇所がいくつかあったので、備忘録の意味でもここに記しておきます。(大体は元の書き方がよろしくなかっただけですが)

小中大括弧でエラー

  • 元々の書き方
[tex:\(x + y \)]
[tex:\{x + y \}]
[tex:\left{x + y \right}]
[tex:\[x + y \]]
[tex:\left[x + y \right]]

見え方
\(x + y \) :おかしい。
\{x + y \}:うまく表示される。
\left{x + y \right}:おかしい。
[x + y ]:うまく表示される。
\left[x + y \right]:おかしい。

小括弧は\エスケープだけではエラーになり、中大括弧は\left, \rightをつけるとエラーになります。そのため、小括弧には\left、\rightをつける、中大括弧は\left\{ \right}\}のように\left、\rightの後、更に\でエスケープする方法で対策しました。

  • 修正後の書き方:
[tex:\left(x + y \right)]
[tex:\{x + y \}]
[tex:\left\{x + y \right\}]
[tex:\[x + y \]]
[tex:\left\[x + y \right\]]

見え方
\left(x + y \right):\left, \right追加でOKになった。
\{x + y \}:変更なし。
\left\{x + y \right\}:\エスケープ追加でOKになった。
[x + y ]:変更なし。
\left[x + y \right]:\エスケープ追加でOKになった。

sum、min、max等の上置き,下置きが上付き、下付きになる

  • 元々の書き方
[tex:\sum_{i = 1}^{N}x_i]
[tex:\min_{x} f(x)]

見え方
\sum_{i = 1}^{N}x_i
\min_{x} f(x)

本当は\Sigma\minの真上、真下に添え字が来るようにしたいのに、右上、右下に来ています。これは\displaystyleにすることで解決しました。

  • 修正後の書き方
[tex:\displaystyle \sum_{i = 1}^{N}x_i]
[tex:\displaystyle \min_{x} f(x)]

見え方
\displaystyle \sum_{i = 1}^{N}x_i
\displaystyle \min_{x} f(x)

\bfにより意図しないところまでも太字化される

  • 元々の書き方
[tex:(\bf{x}_N)]

見え方
(\bf{x}_N):xの先にある添え字のNや)まで太字化される。

わかりにくいですが、xだけ太字にしたいのに、それ以降の文字まで太字になっています。これは元々の書き方がそもそも間違っており、中括弧の位置を\bf{x}→{\bf x}と変更すると修正できます。

  • 修正後の書き方
[tex:({\bf x}_N)]

見え方
({\bf x}_N):xのみが太字。


ということで、今回は以上。せっかくはてなブログに移行しましたが、次回更新は遠い未来になるかも知れません。