cvxpyを使った凸最適化

Pythonで凸最適化を行うための便利なライブラリcvxpyを使う機会がありましたので、使い方を整理しておきます。

凸最適化

凸最適化 (convex optimization) は、制約条件がある中で目的関数の最大化 (または最小化) を行う最適化問題の1つですが、特に以下の特徴を持ちます。

目的関数が2次式 (線形を仮定していない)
局所的な最大値 (最小値) が大域的な最大値 (最小値) と一致する

変数の集合 $\vec{x}$ とすると、最大化 (または最小化) したい目的関数 $f(\vec{x})$ 、m個の不等式制約 $g_i(\vec{x})$ 、n個の等式制約 $h_j(\vec{x})$ によって一般化されます。

max. or min.
- s.t. $g_1(\vec{x}) \leq 0, ..., g_m(\vec{x}) \leq 0$
- s.t. $h_1(\vec{x}) = 0, ..., h_n(\vec{x}) = 0$

凸最適化は、ラグランジュ関数L ( $\lambda_i, \mu_j$ はラグランジュ乗数) を立式し、

$L(\vec{x}, \vec{\lambda}, \vec{\mu}) = f(\vec{x}) + \sum_{i=0}^{m} \lambda_i g_i(\vec{x}) + \sum_{j=0}^{n} \mu_j h_j(\vec{x})$

以下のカルーシュ・キューン・タッカー条件 (KKT条件) を使って方程式を立てることで、解析的に導出できます。

$\frac{\partial L}{\partial \hat{\vec{x}}} = 0$
$\frac{\partial L}{\partial \hat{\vec{\mu}}} = 0$
$\lambda_i \geq 0$
$g_i(\hat{\vec{x}}) \leq 0$
$\lambda_i g_i(\hat{\vec{x}}) = 0$

cvxpy

cvxpyは上の凸最適化を簡単に計算してくれるPythonライブラリです。

ドキュメント: https://www.cvxpy.org/
GitHub: https://github.com/cvxgrp/cvxpy

pipで簡単にインストールできます。

$ pip install cvxpy

cvxpyを使った実装

それではcvxpyを使って凸最適化を行ってみましょう。ここでは以下の問題で最適化を行います。

min. $x_1^2 + 2x_2^2$
s.t. $x_2 \geq - x_1 + 1$ , $x_2 \geq \frac{1}{2} (x_1 - 1)$

この問題を解く実装が↓です。

変数はVariableで宣言
- 第1引数shapeで次元数を渡せるので、ベクトルや行列でもOKです
目的関数は、最小化ならMinimize、最大化ならMaximizeで、定義
- Variableとの計算で使える演算子は、こちらを参照
Problemで目的関数・制約条件を渡してインスタンスを作り、solveメソッドで最適化
- 制約条件は複数あるためリストにまとめて渡す
最適化後の変数は、valueプロパティで取得 ( となってます)
- 制約条件から、dual_valueプロパティでラグランジュ乗数も取得できます ( $\lambda_2 = 0$ ですので2番目の制約条件は使われなかったことになります)

import numpy as np
import cvxpy as cp  # pip install cvxpy

# 変数
x_1 = cp.Variable()
x_2 = cp.Variable()

# 最小化したい関数 (目的関数)
objective = cp.Minimize(
    cp.power(x_1, 2) + 2 * cp.power(x_2, 2)
)

# 制約条件 (不等式制約)
constraints = [
    x_2 >= - x_1 + 1,
    x_2 >= (x_1 - 1) / 2.0
]

# 問題を定義
problem = cp.Problem(objective, constraints)

# 最適化 (戻り値は最適化後に得られた値=最小値)
result = problem.solve()
print(result)  # 0.6666666666666667

# 最小値の時のx_1とx_2の値
print(x_1.value, x_2.value)  # 0.6666666666666667 0.33333333333333337

# 最小値の時のラグランジュ定数
print(constraints[0].dual_value, constraints[1].dual_value)  # 1.3333333333333337 0.0

上では1番目の制約 ( $x_2 \geq - x_1 + 1$ ) 上に最適点がありましたが、制約式上に無い、つまり制約がなかった場合の目的関数の最小値と一致する場合を見ておきます。

目的関数を以下のように書き換えます。この場合、 $(\hat{x_1}, \hat{x_2}) = (3, 3)$ となり、制約条件は実質不要となります。

# 最小化したい関数 (目的関数)
objective = cp.Minimize(
    cp.power(x_1 - 3, 2) + 2 * cp.power(x_2 - 3, 2)
)

これで最適化すると、予想通り $(\hat{x_1}, \hat{x_2}) = (3, 3)$ が得られました。さらにラグランジュ乗数が全て0となっており、制約式上に無いことがわかります。

# 最適化
result = problem.solve()
print(result)  # 8.467107550921043e-48

# 最小値の時のx_1とx_2の値
print(x_1.value, x_2.value)  # 3.0 3.0

# 最小値の時のラグランジュ定数
print(constraints[0].dual_value, constraints[1].dual_value)  # 0.0 0.0

もうちょっと見てみます。最小化式を最大化式に無理やり書き換えてみます。この場合、目的関数は無限まで発散します。

# 最小化から最大化する式に変える
objective = cp.Maximize(
    cp.power(x_1, 2) + 2 * cp.power(x_2, 2)
)

実行すると、以下のエラーとなりました。

DCPError: Problem does not follow DCP rules.

まとめ

今回は凸最適化を行うPythonライブラリcvxpyについて、簡単に紹介しました。