凸最適化ソルバーCVXPYの紹介

この記事はslideshareにアップした[凸最適化 〜双対定理とソルバーCVXPYの紹介〜](https://www.slideshare.net/ssuser2a48ea/cvxpy)から抜粋したものです。
slideshareで全画面にしないと日本語が消えてしまったのでQiitaに移植しました。

# 凸最適化ソルバー
Pythonで使える有名なソルバー：

- CVXOPT
    - 特定の形式の最適化問題を解くことができる
    - 行列の形式に合わせる必要がある
![cvxopt.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/294661/b0c6ed47-0023-48bf-f293-d3d92670b8f8.png)
        - https://cvxopt.org/userguide/coneprog.html
- **[CVXPY](https://www.cvxpy.org/index.html) ←今回はこちらを紹介**
    - **自由に目的関数・制約を書くことができる**
    - matlabのCVXのpythonバージョン（両方Boydさんのソフトウェア）
        - Boydさんは[Convex Optimization](https://web.stanford.edu/~boyd/cvxbook/bv_cvxbook.pdf)の著者です

# CVXPY
- インストール
    - pipの場合: `pip install cvxpy`
    - condaの場合: `conda install -c conda-forge cvxpy`
- インポート
    - import cvxpy as cvx

これから紹介するサンプルは[github](https://github.com/birdwatcherYT/cvxpy_example)にアップしてあります。
## 例題
記述が単純なので例を見れば書き方がわかります. 

### 線形計画問題（LP）
$$
	\min_{x} c^\top x \mathrm{~s.t.~} Ax\le b
$$

```python
 np.random.seed(1)
 # 次元数
 m = 3
 n = 2
 # 各種定数・変数
 A = np.random.randn(m, n)
 b = np.random.randn(m)
 c = np.random.randn(n)
 x = cvx.Variable(n) # 変数定義
 # 問題設定
 obj = cvx.Minimize(c.T @ x) # 最小化
 constraint = [A @ x <= b] # 制約
 prob = cvx.Problem(obj, constraint)# 問題
 prob.solve(verbose=True) # 解く
 # 表示
 print("obj: ", prob.value)
 print("x: ", x.value)
```
- `cvx.Variable(n)`: `n`次元の変数を作成
- `cvx.Minimize(・)`: 最小化したい目的関数を設定
- `cvx.Problem(・, ・)`: 目的関数と制約を設定
    - 制約はリストでわたす
- `.solve()`: 設定した問題を解く
- `@`で行列の積（Python3.5以降のnumpyと同等）
- `seed`によっては解なしになります: 
    - obj:  -inf
    - x:  None

出力結果：
![スクリーンショット 2021-03-13 11.33.48.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/294661/51d10985-73c3-d37b-ab9e-68a44e95a2a6.png)

- 見た目通りに式を書くだけで双対ギャップを計算しながら解いてくれます
- 双対ギャップって何って方は[ラグランジュ双対問題について解説](https://qiita.com/birdwatcher/items/b23209f06177373c6df4)を参照

### 最小二乗法（二次計画問題：QP）
$$
		\min_{x} \\|Ax-b\\|_2^2
$$

```python
 np.random.seed(1)
 # 次元数 
 m = 10
 n = 4
 # 定数・変数
 A = np.random.randn(m, n)
 b = np.random.randn(m)
 x = cvx.Variable(n)
 # 問題定義
 obj = cvx.Minimize(cvx.sum_squares(A @ x - b))
 prob = cvx.Problem(obj)
 prob.solve(verbose=True)
 # 結果表示
 print("obj: ", prob.value)
 print("x: ", x.value)
```
- `cvx.sum_squares(A @ x - b)`は`cvx.norm(A @ x - b, p=2)**2`とも書ける

出力結果：
![スクリーンショット 2021-03-13 11.38.05.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/294661/53124c42-019a-3c38-0d81-e1b62c247040.png)


### 正定値射影（半正定値計画問題：SDP）
$$
		\min_{X\in \mathbb{R}^{n\times n}}\\|A-X\\|_F^2 \mathrm{~s.t.~} X\succeq O
$$
ただし、$X$は対称行列

- $X\succeq O ~\Leftrightarrow~ a^\top X a\ge 0, \forall a ~\Leftrightarrow~ \lambda_i \ge 0$

```python
 # 行列用意
 n = 2 
 A = np.random.randn(n, n)
 A = A + A.T
 X = cvx.Variable((n,n), symmetric=True)
 # 問題定義
 constraints = [X >> 0] # 半正定値制約
 obj = cvx.Minimize(cvx.norm(A-X, "fro")**2)
 prob = cvx.Problem(obj, constraints)
 prob.solve(verbose=True)
 # 結果表示
 print("obj: ", prob.value)
 print("X: ")
 print(X.value)
```

出力結果：
![スクリーンショット 2021-03-13 11.40.56.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/294661/26d40df8-50f4-43f9-22fb-e4d458328a98.png)

- 自動的に適切なソルバーが選ばれていることが確認できます

```math
	X=\left[\begin{array}{cc}
		x&y\\
		y&z
		\end{array}\right]
```
として、この射影を図示するとこうなります。
![スクリーンショット 2021-03-13 11.42.32.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/294661/8d1b1d7a-de08-226d-5256-7519ff434307.png)

### LASSO
$$
	\min_{w,b} \sum_{i}(w^\top x_i + b - y_i)^2+\lambda\\|w\\|_1
$$

```python
data_boston = sklearn.datasets.load_boston()
X = data_boston.data
y = data_boston.target
(n,d) = X.shape
lam = 3000
# 変数定義
w = cvx.Variable(d)
b = cvx.Variable()
# 問題定義
obj = 0
for i in range(n):
    obj += (X[i]@w+b-y[i])**2
obj += lam * cvx.norm(w, p=1)    
obj = cvx.Minimize(obj)
prob = cvx.Problem(obj)
prob.solve(verbose=True)
# 結果表示
print("obj: ", prob.value)
print("w: ")
print(w.value)
print("b: ", b.value)
```

スパースモデルなので、ちゃんと係数がゼロにつぶれていることが確認できます：
![スクリーンショット 2021-03-13 11.45.59.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/294661/cfecdda6-32a0-59ee-2458-20eca2e16ee1.png)



### SVM
```math
		\begin{split}
			\min_{w,b,\xi} ~&~ \frac{1}{2}\|w\|_2^2+C\sum_{i}\xi_i\\
			\mathrm{s.t.} ~&~ \xi_i \ge 1 - y_i(w^\top x_i + b), \\
						   ~&~ \xi_i\ge0
		\end{split}		
```

```python
 # 適当にデータ作成
 n = 100
 X, y = make_classification(random_state=8,n_samples=n,n_features=2, n_redundant=0, n_informative=2,n_clusters_per_class=1,n_classes=2)
 y[y == 0] = -1
 C = 1
 # 変数
 w = cvx.Variable(2)
 b = cvx.Variable()
 xi = cvx.Variable(n)
 # 問題定義
 obj = cvx.Minimize( cvx.sum_squares(w) / 2 + C * cvx.sum(xi) )
 constraint = [xi >= 0] + [(xi[i] >= 1 - y[i]*(w@X[i]+b)) for i in range(n)]
 prob = cvx.Problem(obj, constraint)
 prob.solve(verbose=True)
```

最適化結果をプロットすると、ちゃんと学習できていることが確認できます：
![スクリーンショット 2021-03-13 11.49.28.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/294661/1c58ce60-9414-64c1-ee7a-341176fd9912.png)

## 注意点
- `@`だけでなく`*`も行列の積として扱われる（要素ごとの積ではない）
    - 特にベクトル同士では内積になる
- **要素ごとの積はcvx.multiply**を使う
    - https://www.cvxpy.org/tutorial/functions/index.html?highlight=elementwise

SVMの例

- $\xi_i \ge 1 - y_i(w^\top x_i + b)$
    - [OK] `constraint = [(xi[i] >= 1 - y[i]*(w@X[i]+b)) for i in range(n)]`
- $\xi \ge 1 - y\circ (X w + b)$ （$\circ$: 要素ごとの積）
    - [OK] `constraint = [xi >= 1 - cvx.multiply(y, X@w+b)]`
    - [NG] `constraint = [xi >= 1 - y*(X@w+b)]`


# まとめ
- CVXPYを使うと数式の見た目通りに書ける
- 双対問題について知りたい方は[ラグランジュ双対問題について解説](https://qiita.com/birdwatcher/items/b23209f06177373c6df4)を参照