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

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

近接勾配法と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:もうひとつ、「単純」という条件も付きます。