20 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
解
説
スパースモデリングのための凸最適化
—
近接勾配法による高速アルゴリズム
永原 正章*
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||1.
はじめに
近年，情報学の分野でスパースモデリングが注目を集 めている[1]．これは，大量の高次元データからデータ フィッティングと説明変数の選択を自動的に行う方法で， ビッグデータなどの大規模なデータであっても本質的に は少数の説明変数しか存在しないというスパース性に着 目する． スパースモデリングにおいては，スパース性を測る指 標としてℓ0ノルム，すなわち高次元ベクトルの非ゼロ 要素の総数が用いられる．あるデータbが，スパースな （すなわち小さいℓ0ノルムをもつ）高次元ベクトルxの 線形変換b = Axで与えられたとする．その線形変換は 圧縮的，すなわちbの次元はxの次元に比べて小さいと する．言い換えれば，行列Aは横長であるとする．この とき，連立方程式b = Axは無数に解をもち，データbだ けからもとのベクトルxは復元できない．しかし，その 無数の解の中から最もスパースなものを選べば，それが もとの信号xとなる．これがスパースモデリングのアイ デアであり，圧縮センシングともよばれる[2–4]． スパースモデリングの技術が盛んに研究され，応用さ れている背景に最適化手法の発展がある．「最もスパー スなものを選ぶ」という最適化は，本質的に組合せ最適 化であり，ディジタル画像のような数百万次元の高次元 データに対して，そのまま総当たり的に最適解を求める ことは（現在のコンピュータでは）不可能である．スパー スモデリングにおける最も重要なアイデアは，スパース 性を測るℓ0ノルムを凸関数であるℓ1ノルムで近似する ことである．このアイデアは1994年ごろからStanford 大学のDonohoらによって提唱され，その有効性が計算 機シミュレーションにより検証された[5,6]．当時はなぜ ℓ1ノルムによる近似がうまくいくのかよくわかっていなかったが，CandesとTao[7]やDonoho[8]など2000 年代の研究によって，その正当性が（部分的に）証明さ れた． このℓ1ノルムを用いた最適化は凸最適化となり，内 点法など標準的な数値最適化手法[9]を用いれば，総当 たり法よりはるかに早く解が求まる．しかし，画像処理 ∗ 北九州市立大学 環境技術研究所
Key Words: sparse modeling, convex optimization,
proxi-mal splitting, image processing.
などデータが大規模になると，内点法でも解をみつける のが難しくなる．このような背景のもと，近年のスパー スモデリングの発展は凸最適化理論の発展と足並みをそ ろえている．具体的には，ℓ1ノルムなど，凸ではある が微分ができない関数を含む凸最適化問題を極めて効率 的に解くために近接勾配法というアルゴリズムが提案さ れ[10,11]，さらにその加速法が示された[12]．この高速 アルゴリズムは，内点法に比べてはるかに高速である． 本稿では，スパースモデリングのための凸最適化の基 礎を理解していただくために，最小2乗法と正則化から 話をはじめて，多項式曲線フィッティングを題材にℓ1正 則化（LASSOともよばれる）によるスパースモデリン グの定式化と，それを解くための近接勾配法にもとづく 高速アルゴリズムを説明する．
2.
最小 2 乗法と正則化
本節では，凸最適化のうち最も基本的な最小2乗法と 正則化について復習する．2.1
劣決定系と最小ノルム解 まず，つぎの連立方程式を考えよう． x1+ x2+ x3= 3 x1−x3= 0 (1) 変数の数は三つであるが，方程式の数は二つであるので， この連立方程式は「解けない」．もう少し詳しくいうと， この方程式は無数の解をもち，それはt∈ Rをパラメー タとして， x1= t, x2=−2t+3, x3= t (2) と書ける．このような連立方程式を劣決定系 (underde-termined system)とよぶ． この状況は，探偵のたとえでいえば，証拠や証言が 足りず，犯人を一人に絞り込めないという状況に似てい る．ここで，捜査により「犯人（解）は容疑者（解の候 補）の中で一番小さい」という情報が新たに得られた とする．この情報を使い，解の候補の中から，ベクトル x = (x1,x2,x3)のℓ2ノルム（ユークリッドノルム）が一 番小さいものを選ぶ．すなわち，∥x∥2 2= x21+ x22+ x23 = t2+ (−2t+3)2+ t2 = 6(t−1)2+ 3 を最小化するtはt = 1であるので，(2)式より解は (x1,x2,x3) = (1,1,1)と一意に求まる． 以上のアイデアを一般化してみよう．行列とベクトル を使い，連立方程式を Ax = b (3) と書く．(1)の例では， A = [ 1 1 1 1 0−1 ] , x = x1 x2 x3 , b = [ 3 0 ] である．行列Aのサイズはm× nで，劣決定系を考え， m < nと仮定する．また，行列Aの横ベクトルはすべて 一次独立，すなわち行フルランクであると仮定する．行 フルランクでない場合は，連立方程式の中に冗長な方程 式が含まれているということであり，このような冗長性 はあらかじめ排除されていると仮定する．この仮定のも とで，連立方程式(3)には解が無数に存在する．そのう ち，最もℓ2ノルムが小さい解を求めよう．これはつぎの 最適化問題として定式化される． min x∈Rn∥x∥ 2 2 subject to Ax = b (4) この最適化問題の解を最小ノルム解 (minimum-norm solution)とよぶ．Lagrangeの未定乗数法を用いれば最 小ノルム解x∗は以下のように閉形式で求まる． x∗= A⊤(AA⊤)−1b (5)
2.2
回帰問題と最小2
乗法 つぎに，2次元データ D ={(t1,y1),(t2,y2),...,(tm,ym)} (6) が与えられたとき，このデータを表現するn−1次多項式 y = f (t) = a0+ a1t + a2t2+···+an−1tn−1 (7) を求める問題を考える．このような問題を回帰問題 (re-gression)または多項式曲線フィッティング(polynomial curve fitting)とよぶ． まずはじめに，与えられたデータ点の上を確実に通 る補間多項式(interpolating polynomial)を求めてみよ う．多項式(7) が2次元データ(6)のすべての点の上を 通ることより，係数a0,a1,...,an−1に関するつぎの連立 方程式が得られる． 2 4 6 8 10 12 14 -60 -40 -20 0 20 interpolating polynomial 第1図 データにノイズが加わったときの14次補間多項式 a0+ a1t1+ a2t12+···+an−1tn−11 = y1 a0+ a1t2+ a2t22+···+an−1tn−12 = y2 .. . a0+ a1tm+ a2tm2+···+an−1tnm−1= ym (8) ここで， A = 1 t1 ... tn1−1 1 t2 ... tn−12 .. . ... . .. ... 1 tm ... tn−1m , x = a0 a1 .. . an−1 , b = y1 y2 .. . ym とおくと，(3)式と同じ方程式Ax = bが得られる．行列 AはVandermonde行列とよばれる．いま，m = nとす れば，行列Aは正方行列となり，その行列式は det(A) = ∏ 1≤i<j≤m (tj−ti) で与えられる．これより，もしti̸= tj (i̸= j)ならば，行 列Aは逆行列をもち，連立方程式(8)の解はA−1bとな る．すなわち，データのサイズmに対してm−1次多項 式を選べば，完全にデータ点上を通る多項式（の係数） が得られることがわかる． 具体的な例題で補間多項式を求めてみよう．いま，つ ぎのデータが与えられているとする． t 1 2 3 . . . 14 15 y 2 4 6 . . . 28 30 このデータに対して上の方法で補間多項式を求めると， a1= 2, ai= 0 (i̸= 1) という答えが得られ，補間多項式 はy = 2tとなる．これは確かにすべてのデータの上を通 る多項式である． データにはノイズがつきものである．そこで，上の データyに平均0，標準偏差0.5の正規分布に従うGauss ノイズを加え，このデータをもとに補間多項式を求めて みる．得られた曲線を第1図に示す．曲線はデータ点の 上を確かに通ってはいるものの，ノイズが影響して振動 し，線形の関係y = 2tは失われている．このような現象 を過学習または過適合(over fitting) とよぶ．22 2 4 6 8 10 12 14 0 5 10 15 20 25
30 least squares solution
第2図 最小2乗法による近似直線（実線）と14次補間多 項式（破線） 過学習は，多項式の次数がデータ数に比べて大きすぎ ることが原因である．データの背景に線形性があるとい うことがあらかじめわかっていれば，多項式を1次と仮 定し，y = a0+ a1tとしたうえで，データに最も近くな る係数a0とa1を求めるべきであろう．この場合，すべ てのデータ点の真上を多項式が通ることは不可能となる． しかし，データにノイズが乗っていることを考えると， 多項式がすべてのデータ点の真上を通る必要はまったく ない． 多項式とデータとの距離をℓ2ノルム（ユークリッドノ ルム）で測るとすると，問題は以下で定式化される． min x∈Rn∥Ax−b∥ 2 2 (9) いま，n < m，すなわち多項式の次数をm− 2以下と仮 定すると，一般に連立方程式Ax = bの解は存在しなく なる．しかし，このような状況でも，ti̸= tj (i̸= j) が 満たされていれば，最適化問題(9)の解は一意に定まる． 実際，Vandermonde行列の性質より，上記が満たされ れば行列Aは列フルランクとなり，その一意解は x∗= (A⊤A)−1A⊤b (10)
で与えられる．これを最小2乗解(least squares solution)
とよぶ．最小ノルム解(5)と同様に最小2乗解も閉形式 で得られる．第2図に最小2乗法で得られた直線を示す． すべてのデータ点の上を通ってはいないが，過学習は起 きていない．
2.3
正則化法 前節でみたように，多項式曲線フィッティングでは，多 項式の次数をデータ数よりもはるかに小さくし，最小2 乗法を用いることで過学習を避けることができる．しか し，どの程度次数を小さくすべきかがわからない場合は どうすればよいであろうか．これを解決する方法をここ で説明する．まず，例題から考えてみよう． 2次元のデータを正弦波から生成する．具体的には， ti= i−1, yi= sin(ti) + ϵ, i = 1,2,...,11でデータ{ti,yi} 0 2 4 6 8 10 -1.5 -1 -0.5 0 0.5 11.5 least square solution (deg=10)
第3図 正弦波データの補間多項式（10次），破線はもとの 正弦波 0 2 4 6 8 10 -1.5 -1 -0.5 0 0.5 1
1.5 least square solution (deg=6)
第4図 正弦波データの最小2乗多項式（6次），破線はも との正弦波 を生成する．ただし，ϵは平均0，標準偏差0.2の正規分 布に従う確率変数(Gaussノイズ)とする．この11点の データに対する補間多項式を第3図に，また6次の多項 式による最小2乗解を第4図にそれぞれ示す． 第3図の 補間多項式では過学習がみられるが，6次の多項式はう まくフィッティングできている．しかし，6次多項式が よいというのは，オリジナルの正弦波の曲線がみえてい るからわかることであり，実際には何次の多項式がよい かはデータだけからはわからない．これを解決するため に，10次多項式と6次多項式の係数ベクトルc10とc6を みてみよう． c10= −0.0248 0.9783 −4.7147 8.5862 −6.6640 2.7528 −0.6727 0.1007 −0.0091 0.0005 −0.0000 , c6= −0.0406 −0.6730 2.0942 −1.2085 0.2650 −0.0249 0.0008
0 2 4 6 8 10 -1.5 -1 -0.5 0 0.5 1
1.5 regularized least square solution (deg=10)
第5図 正則化最小2乗法で求まった多項式（10次） ここで，各ベクトルにおいて絶対値の大きい要素を順に 三つ選び，太文字で示している．過学習を行う10次補 間多項式の係数は，6次多項式の係数に比べて大きいこ とがわかる．この観察より，誤差を小さくするとともに 係数そのものも小さくすればよいのではないかと思いつ く．実際，その方法はうまくいく．すなわち，つぎの最 適化問題を解けばよい． min x∈Rn∥Ax−b∥ 2 2+ λ∥x∥22 (11)
これを正則化最小2乗法(regularized least squares)，ま たはリッジ回帰(ridge regularization)とよび，加えられ たλ∥x∥2 2の項を正則化項 (regularization term)とよぶ． 係数λは正則化パラメータ(regularization parameter) とよばれ，誤差と解の大きさとのバランスを調整する重 要なパラメータである． 最適化問題(11)に対しては，最小2 乗解と同じく Lagrangeの未定乗数法により，下記のように閉形式で 解が得られる． x∗= (λI + A⊤A)−1A⊤b (12) 上の多項式補間の例題で，次数を10とし，正則化パ ラメータをλ = 1として，正則化最小2乗法により求め た多項式を第5図に示す．10次多項式であるが，6次の 最小2乗解とほぼ同等のフィッティングが達成されてい ることがわかる．
3.
スパースモデリングと ℓ
1正則化
回帰問題のもう一つの例を考えてみよう．80次多項式 y =−t80+ t (13) から，ti= 0.1i (i = 0,1,...,10) として11点のデータ D ={(ti,yi)}i=0,1,...,10を生成する．なお，多項式は80 次以下の多項式ということはわかっているものとする． 与えられたデータDと「80次以下」という条件から，も との多項式(13)を復元することは可能であろうか？ 明 らかに解は無限個存在し，一意に多項式を定めることは できない．実際，連立方程式Ax = bで書けば，Aは横 長の行列となっており，解は無限個存在する． もとの 80次多項式 (13)は，その係数のほとんど がゼロである．すなわち，係数を要素とするベクトル x = (a0,a1,...,an−1)はスパースである．この情報も使 えると仮定しよう．すなわち，ここでは， もとの多項式は係数がスパースである という情報を加味したうえで，11個のデータだけからも との80次多項式を復元することを考える． この目標のために，つぎのℓ0ノルムとよばれる量を定 義する． ∥x∥0≜ xの非零要素の数 (14) このℓ0ノルムを用いて，つぎの最適化問題を考える． min x∈Rn 1 2∥Ax−b∥ 2 2+ λ∥x∥0 (15) これをℓ0正則化 (ℓ0 regularization)とよぶ．このアイ デアは，前節で述べた正則化最小2乗法と同じである． 正則化最小2乗法では，解のℓ2ノルムを正則化項に用い たが，今回はスパース性を使うためℓ0ノルムを正則化 項に用いている．この最適化問題は有限回の最小2乗法 を解くことにより厳密解を求めることが可能である（い わゆる総当たり法）．しかし，その計算回数はnの指数 オーダで増大し，規模の大きな問題に対して，この方法 ではまったく歯が立たない．その原因は，ℓ0正則化項に ある．これを強引にℓ1ノルムに変更したつぎの最適化問 題を考えよう． min x∈Rn 1 2∥Ax−b∥ 2 2+ λ∥x∥1 (16) ただし，ベクトルxのℓ1ノルムを ∥x∥1≜ n ∑ i=1 |xi| で定義する．(16)式の最適化問題をℓ1正則化(ℓ1 regu-larization) またはLASSO（「ラッソ」と読む）とよぶ． なお，LASSOは1996年にTibshirani[13]によって提案 されたとされているが，それに先駆けて1980年代末に 同様のアイデアを九工大の石川真澄名誉教授1が提案し ていたようである[2]． 正則化項λ∥x∥1は凸関数であり，最適化問題(16)は 凸最適化問題となる．したがって，数値最適化法（たと えば内点法）を用いれば，容易に解が求まる．また，非 常に興味深いことに，もとの最適化問題(15)の解とそ の凸緩和である(16)式の解が一致する条件（十分条件） が知られている[14]．多くの応用でℓ1正則化が有効であ 1石川先生は，筆者も発足メンバーとして関わっている ひびきのAI社会実装研究会の現会長である．24 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 第6図 11次の補間多項式 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 第7図 スパースモデリング：ℓ1正則化による曲線y = t80−t の復元 る例が示されており，スパースモデリングといえばℓ1ノ ルムを用いた最適化を指すことがあるほど，ポピュラー な手法になっている1． 実際に11点のデータから多項式 (13)を復元してみ よう．まず，11点のデータから10次の補間多項式を求 めたのが第6図である．ノイズがない状況を考えている ので過学習特有の振動現象はみられないが，t = 0.9から t = 1のあたりでもとの曲線と大きくずれている．いっぽ う，正則化パラメータをλ = 10−7としたときのℓ1正則 化による復元（スパースモデリング）の結果を第7図に 示す．もとの曲線y = t80 −tが正確に復元できているこ とがわかる． 1もちろん，スパースモデリングには，ℓ1正則化やその ほかℓ1ノルムを用いた凸最適化に帰着する方法以外 に，たとえば欲張り法を用いたものやℓpノルム（ただ し，0 < p < 1）を用いたものなど，さまざまな方法が ある．文献[15,4]などを参照されたい．
0
λ
−
λ
S
λ(b)
b
第8図 ソフトしきい値作用素Sλ(b)4.
高速アルゴリズム
さて，(16)式は凸最適化問題であり，内点法などの標 準的な数値最適化法により効率的に解が求まる．しかし， 規模が非常に大きい場合，たとえば画像処理のようにサ イズnが数百万となった場合，標準的な内点法では計算 に膨大な時間がかかってしまう．また，ℓ1正則化を制御 ループの中で使うような場合（たとえば[16,17]など）， 最適解を極めて高速に求める必要がある．これらの場合 に有効な高速アルゴリズムを本節では解説する．4.1
簡単なケース まず，最も簡単なケースでℓ1正則化を考えてみよう． 次元をm = n = 1としたスカラーのℓ1正則化 min x∈Rf (x), f (x)≜ 1 2(x−b) 2+ λ |x|, x ∈ R (17) を考える．関数f (x)をx≥ 0とx < 0で場合分けし，さ らにb≥ λ, −λ < b < λ, b < −λの三つの場合に分けて計 算すれば，関数(17)を最小化するx∗∈ Rは x∗= Sλ(b) = b−λ, b ≥ λ 0, −λ < b < λ b + λ, b≤ −λ (18) となる．関数 Sλ(b) をソフトしきい値作用素 (soft-thresholding operator) とよぶ．第8図にソフトしき い値作用素を示す．以上より，スカラーの場合，ℓ1正則 化(16)の解はソフトしきい値作用素を用いて，閉形式 (18)で書けることがわかった． つぎに簡単なケースは，行列Aが直交行列Qのとき， すなわち，A = Q, QQ⊤= Q⊤Q = Iが成り立つ場合で ある．このとき，任意のベクトルx∈ Rnに対して， ∥x∥2 2= x⊤x = x⊤Q⊤Qx =∥Qx∥22=∥Q⊤x∥22 が成り立つ．これより， ∥Qx−b∥2 2=∥Q⊤(Qx−b)∥22 =∥Q⊤Qx−Q⊤b∥22 =∥x−Q⊤b∥2 2となり，β = Q⊤bとおくと，(16)式の評価関数は f (x) =1 2∥Qx−b∥ 2 2+ λ∥x∥1 =1 2∥x−β∥ 2 2+ λ∥x∥1 = n ∑ i=1 {1 2(xi−βi) 2+ λ |xi| } = n ∑ i=1 fi(xi), fi(xi) = 1 2(xi−βi) 2+ λ |xi| (19) と書き換えることができる．ただし，βi, xiはそれぞれ ベクトルβおよびxの第i要素を表す．これより，最小 化問題は min x∈Rnf (x) = minx1,...,xn n ∑ i=1 fi(xi) と書ける．各xiは独立に選べるので，スカラーのℓ1正 則化(17)に帰着し， min x1,...,xn n ∑ i=1 fi(xi) = n ∑ i=1 min xi fi(xi) = n ∑ i=1 fi(x∗i) となる．ただし，x∗ i= Sλ(βi)である．ここで，ベクト ルβ∈ Rmに対して，ソフトしきい値作用素Sλ(β)をベ クトルβの各要素に(18)式の関数が作用するように定義 しなおす．すると，行列Aが直交行列Qの場合のℓ1正 則化(16)の解はつぎの閉形式で得られることがわかる． x∗= Sλ(β) = Sλ(Q⊤b) (20) つぎに行列Aが二つのn次直交行列PとQにより， A =[P Q]∈ Rn×2n (21) と書ける場合を考える．行列Aの分割に合わせて， x = [ v w ] , v = x1 .. . xn , w = xn+1 .. . x2n と表すと，以下の等式が成り立つ． Ax =[P Q] [ v w ] = P v + Qw ∥x∥1= n ∑ i=1 |xi|+ 2n ∑ i=n+1 |xi| = ∥v∥1+∥w∥1 これより， f (x) =1 2∥Ax−b∥ 2 2+ λ∥x∥1 =1 2∥P v +Qw −b∥ 2 2+ λ∥v∥1+ λ∥w∥1 と書ける．ここで， g(v,w)≜1 2∥P v +Qw −b∥ 2 2+ λ∥v∥1+ λ∥w∥1 とおき，f (x)を最小化する代わりに，vとwを別々の変 数とみて，g(v,w)を最小化することを考えてみる．すな わち，初期値v[0]とw[0]を適当に与えて， v[k + 1] = arg min z g(z,w[k]) w[k + 1] = arg min z g(v[k],z), k = 0,1,2,... を繰り返すことを考える．ここで， arg min z g(z,w[k]) = arg min z {1 2∥P z −(b−Qw[k])∥ 2 2+ λ∥z∥1 } であり，行列 Pは直交行列であるから，この解は(20) 式よりソフトしきい値作用素を用いて arg min z g(z,w[k]) = Sλ ( P⊤(b−Qw[k])) と書けることがわかる．同様にして， arg min z g(v[k],z) = Sλ ( Q⊤(b−P v[k])) が成り立つ．以上をまとめると， [ v[k + 1] w[k + 1] ] = Sλ ([ P⊤(b−Qw[k]) Q⊤(b−P v[k]) ]) = Sλ ([ P⊤(b−Ax[k]+P v[k]) Q⊤(b−Ax[k]+Qw[k]) ]) = Sλ ([ P⊤(b−Ax[k])+v[k]) Q⊤(b−Ax[k])+w[k]) ]) k = 0,1,2,... すなわち， x[k + 1] = Sλ(x[k]−A⊤(Ax[k]−b)) k = 0,1,2,... (22) と簡略化して書くことができる．行列Aが(21)式を満 たす場合のこの繰り返し解法をブロック座標緩和(block coordinate relaxation, BCR)とよび，任意の初期値x[0] に対して，ℓ1正則化の最適値に収束することが知られて いる[18,19]． 繰り返しのアルゴリズム(22)は，一般の行列Aに対 して計算することが可能である．一般の場合でも，この アルゴリズムでℓ1正則化の最適解を得ることができる のではないか，と考えるのは自然な発想である．実際， 上記の繰り返しアルゴリズムをほんの少し修正するだけ で，一般的なℓ1正則化を解く高速アルゴリズムが得られ る．以下，その仕組みをみていこう．
4.2
近接勾配法による高速アルゴリズム まずはじめに，凸最適化に必要な基礎概念を復習して おこう．26 関数 f :Rn→ R ∪ {+∞}をプロパーな閉凸関数とす る．すなわち，関数fのエピグラフ(epigraph) epi(f )≜ {(x,t) ∈ Rn×R : f(x) ≤ t} が空でない閉凸集合であるとする．関数fの実効定義域 (effective domain) を dom(f )≜ {x ∈ Rn: f (x) < +∞} で定義する．また，関数fの近接作用素 (proximal op-erator)を
proxf(v)≜ arg min x { f (x) +1 2∥x−v∥ 2 2 } (23) で定義する．この近接作用素がℓ1正則化の高速アルゴリ ズムに重要な役割を果たす． 近接作用素の意味を考えるため，閉凸集合C ⊂ Rnに 対する指示関数(indicator function) iC(x)≜ 0, if x∈ C +∞, if x̸∈ C (24) を考えよう．指示関数は，Cが閉凸集合なら，そのエピ グラフを描いてみればわかる通り，プロパーな閉凸関数 となる．このとき，指示関数iCの近接作用素は
proxiC(v) = arg min x∈Rn { iC(x) +1 2∥x−v∥ 2 2 } = arg min x∈C ∥x−v∥ 2 2 ≜ ΠC(v) (25) と計算できる．ここで，ΠC(v)は点vの集合Cへの射影 である．(25)式より，近接作用素(23)は射影作用素の一 般化（すなわち，iCを一般のfに置き換えたもの）と考 えることができる．なお，関数fをγ > 0でスケーリン グした以下の近接作用素も重要である．
proxγf(v)≜ arg min x { f (x) + 1 2γ∥x−v∥ 2 2 } (26) 以上の準備のもとで，つぎの最適化問題を考える． min x∈Rn{g(x)+h(x)} (27) ここで，関数gはdom(g) =Rnを満たす微分可能な凸 関数，関数hはプロパーな閉凸関数であるとする．関数 hは微分可能であるとは限らず，また(24)式のような指 示関数であってもよい．この最適化問題に対する近接勾 配法(proximal gradient method)のアルゴリズムは以 下で与えられる． x[k + 1] = proxγh ( x[k]−γ∇g(x[k])), k = 0,1,2,... (28) ここで，γ > 0はこのアルゴリズムのステップサイズで あり，∇g(x)は関数gのx∈ Rnにおける勾配を表す．こ のアルゴリズムの意味を考えてみよう．(28)式の右辺の 関数を ϕ(x)≜ proxγh ( x−γ∇g(x)) とおく．このとき，近接作用素の定義(23)，または(26) 式より ϕ(x) = arg min z { h(z) + 1 2γ∥z −x+γ∇g(x)∥ 2 2 } = arg min z {˜g(z,x)+h(z)} ˜ g(z,x)≜ g(x)+∇g(x)⊤(z−x)+ 1 2γ∥z −x∥ 2 2 と書ける．なお上の式変形で，g(x)はϕ(x)の中の最小 化に関しては定数であることに注意する．関数 g(z,x)˜ は，点x∈ Rnの近傍における関数g(z)の2次近似であ る．すなわち，各ステップにおいてx[k]の近傍でg(z) を2次近似して関数˜g(·,x[k])+hを逐次最小化している ことになる．さらに，もし勾配∇gがLipschitz連続で あれば，(27)式は少なくとも一つ解をもち，その解をx∗ とおくと， x∗= ϕ(x∗) = proxγh ( x−γ∇g(x)) が成り立つ（たとえば[20]の4.2節をみよ）．すなわち， (27)式の最適解はϕの不動点である．したがって，反復 アルゴリズム(28)が収束すれば，その収束先は(27)式 の最適解（の一つ）であることがわかる．実際，近接勾 配法のアルゴリズム(28)の収束に関して，つぎの定理が 成り立つ[21]． 【定理 1】 関数gの勾配∇gはRn上でLipschitz連 続とし，LをそのLipschitz定数とする．すなわち，Lは ∥∇g(x)−∇g(y)∥2≤ L∥x−y∥2, ∀x,y ∈ Rn
を満たすもののうち最小のものとする．ステップサイ ズを γ≤1 L (29) とする．このとき，近接勾配法により得られるベクトル 列{x[k]}は最適化問題(27)の解x∗に収束し， ∥x[k +1]−x∗∥ 2≤ ∥x[k]−x∗∥2, k = 0,1,2,... が成り立つ．またf (x) = g(x) + h(x)とおいたとき f (x[k])−f(x∗)≤L∥x[0]−x∗∥ 2 2 2k , k = 0,1,2,... (30) が成り立つ． さて，では具体的にℓ1正則化(16)に対して，近接勾 配法のアルゴリズムを導出してみよう．いま， g(x) =1 2∥Ax−b∥ 2 2, h(x) = λ∥x∥1 であり，g(x)の勾配は
∇g(x) = A⊤(Ax−b) (31)
またγh(x) = γλ∥x∥1の近接作用素は，
proxγh(v) = arg min x ( λ∥x∥1+ 1 2γ∥x−v∥ 2 2 ) = arg min x ( γλ∥x∥1+1 2∥x−v∥ 2 2 ) = Sγλ(v) となる（Q = Iの場合の(19)の最小化を考えよ）．これ より，ℓ1正則化(16)を解く近接勾配法のアルゴリズムは x[k + 1] = Sγλ(x[k]−γA⊤(Ax[k]−b)) k = 0,1,2,... (32) となる．このアルゴリズムを反復縮小しきい値アル ゴリズム(iterative shrinkage thresholding algorithm, ISTA)とよぶ．いま，(31)式より∇gのLipschitz定数 はL = λmax(A⊤A)となる．ただし，λmaxは最大固有値 を表す．したがって定理1の(29)式より， γ≤λ 1 max(A⊤A) となるγを選べば，(32)式の反復によりℓ1正則化(16) の最適解を得ることができる．行列Aが(21)式と表さ れる場合の反復法(22)との類似性に注意してほしい． 定理1の(30)式より，ISTAのアルゴリズムはステップ 数kに対して誤差がO(1/k)のオーダで減少する．ISTA のアルゴリズムは第kステップでx[k + 1]を計算するた めに現在の推定値x[k]のみを使用するが，一つ前の推定 値x[k−1]も使えば，さらに速くなる．実際，つぎのア
ルゴリズムはFISTA (Fast ISTA)とよばれ，O(1/k2)
で収束することが知られている[21,19]．
x[k] = Sγλ(y[k]−γA⊤(Ay[k]−b)) t[k + 1] =1 + √ 1 + 4t[k]2 2 y[k + 1] = x[k] +t[k]−1 t[k + 1](x[k]−x[k −1]) k = 0,1,2,... (33) なお，収束オーダO(1/k2)は最適であり，これ以上速く することは（gの2階微分などの情報を使わない限り）理 論的には不可能である．詳しくは，[22]を参照されたい．
5.
おわりに
本稿では，スパースモデリングで用いる凸最適化につ いて，その基礎を述べた．まず，正則化法とよばれる方 法を最小2乗法の枠組みで説明し，その後，スパースモ デリングの概念とℓ1ノルムを使った定式化を多項式補間 の文脈で紹介した．そして，スパースモデリングでよく 用いられるℓ1正則化（LASSOともよばれる）の最適解 を求めるための近接勾配法を用いた高速アルゴリズムを その導出も含め詳しく述べた．スパースモデリングで使 われる最適化はℓ1正則化に限らず，また高速アルゴリズ ムも近接勾配法に限らない．これらの展開は，日本語の 参考書[15,22,23]や解説記事[1,2,24]などを参考にされ たい． 謝 辞 本研究の一部はJSPS科研費15K14006, 15H02668, 16H01546の助成を受けたものです． (2016年9月8日受付) 参 考 文 献 [1] 山内: 特集 スパースモデリングの発展—原理から応用 まで—;電子情報通信学会誌, Vol. 99, No. 5 (2016) [2] 田中: 圧縮センシングの数理; IEICE FundamentalsReview, Vol. 4, No. 1, pp. 39–47 (2010)
[3] Y. C. Eldar and G. Kutyniok: Compressed
Sens-ing: Theory and Applications, Cambridge University
Press (2012)
[4] K. Hayashi, M. Nagahara and T. Tanaka: A user’s guide to compressed sensing for communications sys-tems; IEICE Trans. on Communications, Vol. E96-B, No. 3, pp. 685–712 (2013)
[5] S. Chen and D. Donoho: Basis pursuit; Signals,
Systems and Computers; Conference Record of the Twenty-Eighth Asilomar Conference on, Vol. 1, pp.
41–44 (1994)
[6] S. S. Chen, D. L. Donoho and M. A. Saun-ders: Atomic decomposition by basis pursuit; SIAM
J. Sci. Comput., Vol. 20, No. 1, pp. 33–61 (1998)
[7] E. J. Candes and T. Tao: Near-optimal signal re-covery from random projections: Universal encoding strategies?; IEEE Trans. Inf. Theory, Vol. 52, No. 12, pp. 5406–5425 (2006)
[8] D. L. Donoho: Compressed sensing; IEEE Trans. Inf.
Theory, Vol. 52, No. 4, pp. 1289–1306 (2006)
[9] S. Boyd and L. Vandenberghe: Convex Optimization, Cambridge University Press (2004)
[10] P. L. Combettes and V. R. Wajs: Signal recovery by proximal forward-backward splitting; SIAM J. on
Multiscale Modeling and Simulation, Vol. 4, No. 4,
pp. 1168–1200 (2005)
[11] P. L. Combettes and J.-C. Pesquet: Proximal
Split-ting Methods in Signal Processing, Springer New
York, pp. 185–212 (2011)
[12] A. Beck and M. Teboulle: A fast iterative shrinkage-thresholding algorithm for linear inverse problems;
SIAM J. Imaging Sci., Vol. 2, No. 1, pp. 183–202
(2009)
[13] R. Tibshirani: Regression shrinkage and selection via the LASSO; J. R. Statist. Soc. Ser. B, Vol. 58, No. 1, pp. 267–288 (1996)
28
by ℓ1 minimization; Ann. Statist., Vol. 37, No. 5A, pp. 2145–2177 (2009)
[15] M. Elad（玉木訳）: スパースモデリング, 共立出版 (2016)
[16] M. Nagahara, D. Quevedo and J. Østergaard: Sparse packetized predictive control for networked control over erasure channels; IEEE Trans. Autom. Control, Vol. 59, No. 7, pp. 1899–1905 (2014)
[17] M. Nagahara, D. E. Quevedo and D. Neˇsi´c: Maxi-mum hands-off control: a paradigm of control effort minimization; IEEE Trans. Autom. Control, Vol. 61, No. 3, pp. 735–747 (2016)
[18] S. Sardy, A. G. Bruce and P. Tseng: Block co-ordinate relaxation methods for nonparametric sig-nal denoising with wavelet dictionaries; J.
Com-put. Graph. Statist., Vol. 9, No. 2, pp. 361–379 (2000)
[19] M. Zibulevsky and M. Elad: L1-L2 optimization in signal and image processing; IEEE Signal Process.
Mag., Vol. 27, pp. 76–88 (2010)
[20] N. Parikh and S. Boyd: Proximal algorithms;
Foun-dations and Trends in Optimization, Vol. 1, No. 3,
pp. 123–231 (2013)
[21] A. Beck and M. Teboulle: Gradient-based algorithms with applications to signal-recovery problems;
Con-vex Optimization, Cambridge University Press (2010)
[22] 鈴木: 確率的最適化,講談社(2015) [23] 富岡: スパース性に基づく機械学習,講談社(2015) [24] 小野: 近接分離による分散凸最適化—交互方向乗数法に 基づくアプローチを中心として—;計測と制御, Vol. 55, No. 11 (2016) 著 者 略 歴 なが 永 はら原 まさ正 あき章 (正会員) 愛媛県生まれ．2003年京都大学大学院 情報学研究科博士課程修了．京都大学助 手，助教，講師を経て，2016年より北九 州市立大学環境技術研究所教授．同年より インド工科大学ムンバイ校(IIT Bombay) の客員教授を兼任．専門は自動制御と人工 知能．博士（情報学）．2012年，IEEE制御システム部門よ り国際賞であるTransition to Practice Awardを受賞．同賞 の受賞は日本人初である．そのほか，計測自動制御学会や電
子情報通信学会の論文賞など，受賞多数．IEEEの上級会員
(Senior Member)．著書に「マルチエージェントシステムの 制御」（共著，2016年度SICE著述賞受賞）などがある．ひ
