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

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

<img src="https://qiita-image-store.s3.amazonaws.com/0/202074/0c2267dc-c99f-df02-548f-f01a5c26027c.jpeg" width=25%>
[出版社ページ](http://www.coronasha.co.jp/np/isbn/9784339032222/)
[Amazonページ](https://www.amazon.co.jp/gp/product/4339032220/)

なぜMATLABなのかというと，制御工学や信号処理，通信などスパースモデリングの応用先の工学系分野では主流の言語がMATLABだからです．
<details>
著者がMATLAB以外の言語に詳しくない，というのが大きな理由かもしれません（笑）
</details>

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

# CVX について
スパースモデリングでは，<a href="https://ja.wikipedia.org/wiki/%E5%87%B8%E6%9C%80%E9%81%A9%E5%8C%96">凸最適化</a>をガンガン解きます．凸最適化を解くためのアルゴリズムを実装するのは，なかなかハードです．すぐにスパースモデリングを試してみたい方は，アルゴリズムが実装済みのツールボックスを使うことをおすすめします．

凸最適化問題を解くためのツールボックスはたくさんありますが，私は <a href="http://cvxr.com/cvx/">CVX</a> を強くお薦めします． CVX を使えば，凸最適化問題を非常に簡単にコーディングできます．

<font color="red">現時点でCVXは**MATLAB2017a以降には対応していません**．MATLAB2016b以前のバージョンをお使い下さい．</font>

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

<a href="http://cvxr.com/cvx/">http://cvxr.com/cvx/</a>

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

<a href="http://cvxr.com/cvx/download/">CVXダウンロードページ</a>

特別な目的がなければ**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を使えば，以下のようにして解くことが出来ます．

```matlab
cvx_begin
   variable x(1)
   minimize x^2 + x + 1
cvx_end
```
これを実行すると，以下のような結果が得られます．

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

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

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

```matlab
>> x

x =

   -0.5000
```

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

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

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

```math
\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を使えば，次のように直感的にコードが書けます．

```matlab
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
```

これを実行すると，

```matlab
Calling SDPT3 4.0: 11 variables, 5 equality constraints
（省略）
Status: Solved
Optimal value (cvx_optval): +0.416667
```
すぐに結果が得られます．

```matlab
>> x,y,z

x =

    0.3333


y =

    0.1667


z =

    0.5000
```

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

<img src="https://qiita-image-store.s3.amazonaws.com/0/202074/0c2267dc-c99f-df02-548f-f01a5c26027c.jpeg" width=25%>
[出版社ページ](http://www.coronasha.co.jp/np/isbn/9784339032222/)
[Amazonページ](https://www.amazon.co.jp/gp/product/4339032220/)