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$ がかかっていないため、最適解も $\lambda$ に $1/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、劣微分、凸解析などなど最適化理論の色々なトピックスに触れることができました（まだ全然理解できていないですが）。
冒頭に書いた思いついたことはまず試してみるというスタンスは、試すために知らないトピックスを調べるきっかけになるというメリットがあります。やってみてダメなら何がダメか更に深追いするモチベーションにもなります。後はこういう試行錯誤は何より楽しいですね。ということで、今後は近接勾配法や凸解析の理論などの知識を深めていければと思います。