スパースモデリング (2): CVXの使い方

MATLABについて

人工知能・機械学習の分野ではPythonがプログラミング言語としては主流ですが，書籍「スパースモデリング- 基礎から動的システムへの応用 -」ではMATLABを使用したスパースモデリングのコードを紹介しています．

[出版社ページ](http://www.coronasha.co.jp/np/isbn/9784339032222/) [Amazonページ](https://www.amazon.co.jp/gp/product/4339032220/)

なぜMATLABなのかというと，制御工学や信号処理，通信などスパースモデリングの応用先の工学系分野では主流の言語がMATLABだからです．

著者がMATLAB以外の言語に詳しくない，というのが大きな理由かもしれません（笑）

なお，MATLABと言っても特殊なことはしておらず，例えばPythonでもほぼ同じことができると思います．

CVX について

スパースモデリングでは，凸最適化をガンガン解きます．凸最適化を解くためのアルゴリズムを実装するのは，なかなかハードです．すぐにスパースモデリングを試してみたい方は，アルゴリズムが実装済みのツールボックスを使うことをおすすめします．

凸最適化問題を解くためのツールボックスはたくさんありますが，私は CVX を強くお薦めします． CVX を使えば，凸最適化問題を非常に簡単にコーディングできます．

現時点でCVXはMATLAB2017a以降には対応していません．MATLAB2016b以前のバージョンをお使い下さい．

CVX のインストール

CVX のホームページに行って，CVXをダウンロードしましょう．

http://cvxr.com/cvx/

このホームページにアクセスすると，上にDownloadというメニューがありますので，クリックしてください．または，直接ダウンロードページにアクセスしてください．

CVXダウンロードページ

特別な目的がなければStandard bundlesをダウンロードしてください（ご自分のPCのOSに合わせて，ダウンロードするファイルを選んでください）．

例えば，Mac OSなら，

cvx-maci64.zip

をダウンロードして，解凍します．すると，cvxというフォルダができますので，MATLABを実行するフォルダに移動させてください．Mac OSであれば自分の書類フォルダ (Documentsフォルダ）の中にMATLABフォルダがありますので，その下に解凍したcvxフォルダを移動させてください．

つぎに，MATLABを起ち上げて，cvxフォルダに移動し，cvx_startup.mを実行します．すなわち，MATLABのコマンドラインで

cvx_startup

を実行します．すると，自動的にパスが設定されて，CVXツールボックスが使えるようになります．うまくインストールできたかどうかは，コマンドラインで

cvx_setup

と打ってみてください．

Testing with a simple model...done!

と出れば，インストールはうまく行っています．

CVXで簡単な最適化問題を解いてみる

ここでは，CVXを使って，簡単な最適化問題を解いてみましょう．

次の問題を考えます．

$$
\mathrm{minimize} ~~x^2+x+1
$$

1次元の最適化問題です．CVXを使えば，以下のようにして解くことが出来ます．

cvx_begin
   variable x(1)
   minimize x^2 + x + 1
cvx_end

これを実行すると，以下のような結果が得られます．

Calling SDPT3 4.0: 3 variables, 1 equality constraints
------------------------------------------------------------

 num. of constraints =  1
 dim. of sdp    var  =  2,   num. of sdp  blk  =  1
*******************************************************************
   SDPT3: Infeasible path-following algorithms
*******************************************************************
 version  predcorr  gam  expon  scale_data
   HKM      1      0.000   1        0    
it pstep dstep pinfeas dinfeas  gap      prim-obj      dual-obj    cputime
-------------------------------------------------------------------
 0|0.000|0.000|4.5e+00|6.1e+00|2.0e+02| 1.000000e+01  0.000000e+00| 0:0:00| chol  1  1 
 1|1.000|1.000|1.8e-15|6.4e-02|1.6e+01| 7.458648e+00 -8.032648e+00| 0:0:00| chol  1  1 
 2|0.996|0.989|0.0e+00|7.0e-03|1.7e-01|-1.659194e-01 -3.225147e-01| 0:0:00| chol  1  1 
 3|0.989|0.989|0.0e+00|7.0e-04|1.9e-03|-2.490721e-01 -2.495622e-01| 0:0:00| chol  1  1 
 4|0.988|0.990|0.0e+00|7.0e-05|2.2e-05|-2.499890e-01 -2.498726e-01| 0:0:00| chol  1  1 
 5|0.983|0.990|2.2e-16|7.2e-07|3.2e-07|-2.499998e-01 -2.499987e-01| 0:0:00| chol  1  1 
 6|1.000|1.000|8.9e-16|1.0e-12|1.9e-08|-2.500000e-01 -2.500000e-01| 0:0:00|
  stop: max(relative gap, infeasibilities) < 1.49e-08
-------------------------------------------------------------------
 number of iterations   =  6
 primal objective value = -2.49999991e-01
 dual   objective value = -2.50000009e-01
 gap := trace(XZ)       = 1.85e-08
 relative gap           = 1.24e-08
 actual relative gap    = 1.24e-08
 rel. primal infeas (scaled problem)   = 8.88e-16
 rel. dual     "        "       "      = 1.00e-12
 rel. primal infeas (unscaled problem) = 0.00e+00
 rel. dual     "        "       "      = 0.00e+00
 norm(X), norm(y), norm(Z) = 1.3e+00, 2.5e-01, 1.3e+00
 norm(A), norm(b), norm(C) = 2.0e+00, 2.0e+00, 2.2e+00
 Total CPU time (secs)  = 0.31  
 CPU time per iteration = 0.05  
 termination code       =  0
 DIMACS: 8.9e-16  0.0e+00  1.1e-12  0.0e+00  1.2e-08  1.2e-08
-------------------------------------------------------------------
 
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +0.75

この内容は，凸最適化理論を良く知っている方なら「なるほど，なるほど」と思いながら読めると思いますが，素人だと，ほとんど意味がわからないと思います．でもそれで結構です．

この結果を解読してみましょう． **SDPT3**というソルバを使って，繰り返し計算により最小点を見つけています．ここで使われているアルゴリズムは**Infeasible path-following algorithms**というものです．6回の繰り返しで計算が終了しています．**primal**と**dual**はそれぞれ，主問題および双対問題のことで，このギャップ (**gap**) が小さくなれば，制度の良い解が得られていることになります．今，**relative gap**は1.853e-08なので，十分小さく，良い精度の解が得られました．この繰り返し計算に0.31秒かかっています．最適値 (**Optimal value**) は0.75です．

コマンドプロンプトで，xと入力してください．計算結果がわかります．

>> x

x =

   -0.5000

最適化問題の答えは$x=-0.5$です．

$x^2+x+1$の最小値は以下のように手計算で求められます． $$ x^2+x+1 = (x+0.5)^2+0.75 $$ これより，$x=-0.5$のとき最小値$0.75$をとることがわかります．

もう少し複雑な最適化問題も解いてみましょう．

\begin{align*}
 &\mathrm{minimize~~}  x^2 + 2y^2 +z^2\\
 &\mathrm{subject~to~~}\\
 &\qquad   x+y+z\geq 1\\
 &\qquad   x+y\leq z
\end{align*}

CVXを使えば，次のように直感的にコードが書けます．

cvx_begin
    variables x y z
    minimize( x^2 + 2*y^2 + z^2 )
    subject to
      x + y + z >= 1
      x + y <= z
cvx_end

これを実行すると，

Calling SDPT3 4.0: 11 variables, 5 equality constraints
（省略）
Status: Solved
Optimal value (cvx_optval): +0.416667

すぐに結果が得られます．

>> x,y,z

x =

    0.3333


y =

    0.1667


z =

    0.5000

CVXを使えば，書籍「スパースモデリング- 基礎から動的システムへの応用 -」に記載されている様々な最適化問題が解けます．ぜひインストールして使ってみてください．

[出版社ページ](http://www.coronasha.co.jp/np/isbn/9784339032222/) [Amazonページ](https://www.amazon.co.jp/gp/product/4339032220/)