はじめに
前回、前々回とl1ノルム正則化項をつけた離散フーリエ変換により、スパースな解が得られること、そして基底ベクトルを並べた行列がユニタリ行列(直交行列)のため解析的に解けることを見てきました。それでは、が一般の行列の場合どうするか?今度こそ解析的に解けないので、反復法により最適解を求めます。しかし、目的関数にはl1ノルムという微分できない項があるため、最急降下法などのアルゴリズムは使えません。
このような微分できない項を含む最適化問題を効率的に解く手法の1つに近接勾配法(Proximal Gradient Method)があります。今回はこの近接勾配法と近接勾配法において重要な役割を果たすproximal operatorについて書いていきます。また、pythonで近接勾配法の実装を行いましたので、そのソースコードも載せておきます。
近接勾配法、proximal operatorとは
まずは、近接勾配法で解ける最適化問題の定義です。
ただし、は微分可能な凸関数、は微分可能とは限らない凸関数*1です。
l1ノルム正則化項付きの問題を式(1)に当てはめると、となります。
次にproximal operatorの説明です。とりあえず定義を見てみましょう。関数のproximal operatorは以下のように定義されます。
ここでは後述の更新式のステップ幅です。
このproximal operatorは何者か?というのは、として以下のような指示関数を考えると見えてきます。
ここでは閉凸集合です。この関数はが集合の要素なら0、の要素でないならに吹っ飛ばすという関数です。そして、この指示関数のproximal operatorは以下のようになります。
が集合の要素の場合、そのままとすればを最小化できます。一方で、が集合の要素ではない場合、集合の要素でとのユークリッド距離を最小とするものが選ばれます。すなわち、の集合への射影です(このような射影作用素をとして定義しています)。
よって、proximal operatorは射影作用素の関数への一般化と考えることができ、がによって定められる領域からはみ出していたら、その領域に入るように射影する作用素と見なせます。
proximal operatorの具体例を挙げると、がというl1ノルムであれば、前回述べたソフトしきい値作用素がproximal operatorになります。つまり、
です。そして、近接勾配法の更新式はこのproximal operatorを使って以下のようになります。
この更新式の収束性や正しさなどの説明は他の文献を参照いただくとして、直観的には次のような解釈ができると思います。
まず、は最急降下法の更新式です。また、上述のようにproximal operatorはによって定められる領域に射影する作用素です。従って、最急降下法の1ステップで得られる値がの領域からはみ出してたら、その領域へ射影するという操作を繰り返して最適解を探すアルゴリズムと見ることができます。
このように、関数の勾配と関数のproximal operatorがわかれば、近接勾配法により最適解が得られます。
なお、の決め方も収束のためには重要ですが、こちらについても別の文献を参照してください。
Moreau decompositionと双対ノルム
さて、ここまで見たように近接勾配法を使うためにはproximal operatorを求める必要があります。しかし、関数によっては、式(2)から直接的にproximal operatorを求めるのが難しいことがあります。例えば、とl2ノルムとした場合(2乗が付いていないことに注意)、式(2)を解こうとしてで微分しても、がかかる項が出てきてしまうので厄介です。
そこで便利なのが以下に示すMoreau decompositionが成り立つというproximal operatorの性質です。
ここで、はの共役関数のproximal operatorです。この共役関数とは下記の式で表される関数のことを言います。
よって、この共役関数のproximal operatorさえ分かれば、元の関数のproximal operatorも求められます。
では、というノルム関数の場合どうなるかを見ていきます。まずノルム関数の共役関数は以下の指示関数となります。
ここで、は双対ノルムと呼ばれるもので、以下の式で定義されます。
ということで、共役関数のproximal operatorは式(4)でみたように下記の射影作用素になります。
ただし、はで表される集合です。
つまり、をの範囲に射影する作用素がノルム関数の共役proximal operatorとなります。
おそらく式だけをザッと書いているので何故ノルムの共役関数が指示関数なのか?双対ノルムってそもそも何?どうやって求めるの?などの疑問が湧くかと思います。その辺りを少し補足しておきます。
双対ノルムについての補足
説明が前後しますが、後の説明のためにまずここで双対ノルムについて補足を書きます。そのための準備としてヘルダーの不等式という不等式を紹介します。
という2つの整数に、という関係があれば、
が成り立つというのがヘルダーの不等式です。実は、の双対ノルムはであり、その逆も成り立ちます。そのため、l2ノルムの双対ノルムはそのままl2ノルムですし、l1ノルムの双対ノルムはlノルムです。
このことはヘルダーの不等式から示すことができます。まずヘルダーの不等式より、任意のに対してが成り立ちます。lpノルムの双対ノルムはという条件下でのの最大値です。ここでをうまくとることにより、等式を成立させることができます。つまり、
が成立します。
何故、ノルム関数の共役関数が指示関数なのか?
まず、の時を考えます。式(8)の定義よりなわけですが、双対ノルムの定義式(10)と前述の前提より、かつとなるが存在することになります。これを満たすをとした場合、そのに定数をかけたも式(8)の最適解の候補として入ってきます。これはもちろん、とした場合も同様です。これらのことから、
となるため、のときはとなります。
一方で、の時は、式(12)のヘルダーの不等式よりなので、
です。前提より、なので、左辺(つまりは)が0以上になることはありません。ここで、という取り方もできるので、式(15)の左辺を0にすることができます。よって、のときはとなります。
結局のところproimal operatorはどうやって求めればよいか?
2つ方法があると思います。1つ目は式(2)から直接求める方法です。例えば、前回記したようの方法でl1ノルム正則化項のproximal operator式(5)を求めるというものです。
もう1つの方法としては、式(7)のMoreau decompositionの性質を使う方法です。特に正則化項がノルムなら、共役関数のproximal operatorはという超球への射影なので簡単に求められます。例えば、l2ノルムの双対ノルムはl2ノルムなので、
であり、求めたいproximal operatorは、
となります。
なお、式(2)にはステップ幅のが付いています。が付いていても基本的には上式のをに置き換えればよいだけですが、を含めたMoreau decompositionの一般形を書いておきましょう。
という形式になります。したがって、がノルム関数の場合は、
になります。ちなみに、多くの文献ではとしているので、として理論を展開している点に注意してください。
様々な正則化項のproximal operator
ここまでノルム関数が正則化項の時のproximal operator導出方法を見てきました。ここで、いくつかの正則化項に対するproximal operatorを記載しておきます。
l2ノルムの2乗正則化:
リッジ回帰等でよく使われる正則化項です。普通にproximal operatorの式を微分して0とおけば求められます。2乗しているので、以下のように式(16)とは異なる結果が得られます。
Elastic Net:
l1ノルム正則化とl2ノルム2乗正則化の混合です。l1ノルムは2つ変数が同程度の寄与度だったとしてもどちらかしか選ばれない(一方が0になってしまう)ケースが多いですが、Elastic Netではこの点を改善できるようです。このproximal operatorはl1ノルム正則化の時と同様、変数ごとに分離して求めることができます。
グループl1ノルム(グループに重複がない場合):
変数がいくつかのグループに分けられるとき、それぞれのグループのlpノルムを取って、更にその和を取るノルムです( ベクトルのうち、グループに所属する変数を並べたのが)。これを正則化項に用いるとグループごとにスパースになる、ならないが決まります。このグループl1ノルムを正則化項として使った回帰をgroup lassoと言います。
これはグループに重複がなく完全に分離ができるため、グループごとにlpノルムのproximal operatorを求めればよいです。p=2の時のグループのproximal operatorは下記の通りです。(式(16)とほぼ同じです。)
ちなみに、グループに重複がある場合も提案されていますが、ここでは割愛します(まだ追いきれてない)。
近接勾配法を実装
ということで式(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はにあたります。つまり、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回帰問題を解く関数
以下は、という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を引数としてとり、上記の目的関数を計算する関数として定義しています。
また、であり、(2)はこの計算を行うgrad_fをxに引数とる関数として定義しています。
(3)、(4)はを定めています。最大特異値の逆数より小さい値をとして定めると収束するようですが、その辺りはまだ理解ができておりません。とりあえず、文献の通りに実装します。
(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も変更する必要がありますが上記では省略しています)。
また、を別なものにしたかったらを別途定義して、grad_fとすればよいです。
行列関数の場合も同様で、grad_fとproximal operator (あと、objective)さえ定義してあげたら機能します。そのため、proximal operatorとしてトレースノルムを指定してあげると、低ランク解を得ることができます。
まとめ
今回は近接勾配法とproximal operatorについて見てきました。上記の実装については次の機会に応用事例で実験したいと思います(ちなみに、疑似データで実験したところそれっぽく動いていました)。
また、今回は取り上げませんでしたが、この界隈ではISTAの高速版FISTAや拡張ラグランジュ法、その改良であるAlternating Direction Method of Multipliers (ADMM) など様々なアルゴリズムが提案されており、非常に奥が深いです。
応用例(2018/6/15追記)
近接勾配法の応用例の記事をいくつか書いたのでリンクを貼っておきます。
参考文献
最後にこの記事を書くにあたって参考にした文献です(基本的に前回とほぼ同じです)。
- スパースモデリングのための凸最適化—近接勾配法による高速アルゴリズム:今回の件を調べるきっかけとなった解説論文
- Proximal Algorithms:共役関数、双対ノルム、Moreau decompositionについて非常に参考になった
- Optimization with Sparsity-Inducing Penalties:同上
- Convex Optimization:UCバークレーの講義資料?上2つの内容を補間するのに役にたった
また日本語の書籍で、proximal operator界隈が詳しく書かれているのは以下2冊かなと思います(他にもあったら教えてください)。
スパース性に基づく機械学習 (機械学習プロフェッショナルシリーズ)
- 作者: 冨岡亮太
- 出版社/メーカー: 講談社
- 発売日: 2016/12/23
- メディア: Kindle版
- この商品を含むブログを見る
機械学習のための連続最適化 (機械学習プロフェッショナルシリーズ)
- 作者: 金森敬文,鈴木大慈,竹内一郎,佐藤一誠
- 出版社/メーカー: 講談社
- 発売日: 2017/09/08
- メディア: Kindle版
- この商品を含むブログを見る
*1:もうひとつ、「単純」という条件も付きます。