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$ のベクトルとして式（２）を解けばよいです。この変換で基底ベクトルの数も説明変数の要素数も倍になったように見えますが、元々も複素数でそれぞれが実部と虚部を持っていたため本質的には変わりません。
また、離散フーリエ変換の場合、 $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の使い方は他の解説記事に譲るとして、早速式（２）を解く関数のソースコードを示します。

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の詳しい使い方、凸解析など最適化の理論、実装両面で理解を深めていきたいと思います。