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

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

凸最適化ソルバー

Pythonで使える有名なソルバー：

• CVXOPT
  • 特定の形式の最適化問題を解くことができる
  • 行列の形式に合わせる必要がある
    
    • https://cvxopt.org/userguide/coneprog.html
• CVXPY ←今回はこちらを紹介
  • 自由に目的関数・制約を書くことができる
  • matlabのCVXのpythonバージョン（両方Boydさんのソフトウェア）
    • BoydさんはConvex Optimizationの著者です

CVXPY

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

これから紹介するサンプルはgithubにアップしてあります。

例題

記述が単純なので例を見れば書き方がわかります.

線形計画問題（LP）

$$
\min_{x} c^\top x \mathrm{~s.t.~} Ax\le b
$$

 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

出力結果：

• 見た目通りに式を書くだけで双対ギャップを計算しながら解いてくれます
• 双対ギャップって何って方はラグランジュ双対問題について解説を参照

最小二乗法（二次計画問題：QP）

$$
\min_{x} \|Ax-b\|_2^2
$$

 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とも書ける

出力結果：

正定値射影（半正定値計画問題：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$

 # 行列用意
 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)

出力結果：

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

	X=\left[\begin{array}{cc}
		x&y\\
		y&z
		\end{array}\right]

として、この射影を図示するとこうなります。

LASSO

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

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)

スパースモデルなので、ちゃんと係数がゼロにつぶれていることが確認できます：

SVM

		\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}		

 # 適当にデータ作成
 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)

最適化結果をプロットすると、ちゃんと学習できていることが確認できます：

注意点

• @だけでなく*も行列の積として扱われる（要素ごとの積ではない）
  • 特にベクトル同士では内積になる
• 要素ごとの積は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を使うと数式の見た目通りに書ける
• 双対問題について知りたい方はラグランジュ双対問題について解説を参照