はじめに
本記事で扱う内容：
本記事で扱わない内容：
- 数理最適化に関する理論的な内容
- 最適化問題に対する解法アルゴリズム
- CVXPY 以外の数理最適化モデリングツール(scipy.optimize，PuLP，Pyomo，etc.)との比較
- 「ソルバー」と「モデラー」の違い
- CVXPY で使える関数
数理最適化とは
初めに、数理最適化とは何かについてざっくりと述べておきたいと思う。ただし、このページに辿り着いたということは、恐らく実験や研究などで最適化問題を解く必要が生じた方が大半であると思うので、最小限の説明に留める。数理最適化とは、制約条件と呼ばれる関数で記述される領域において、目的関数と呼ばれる関数を最小化(もしくは最大化)する問題のことであり、一般に次のように記述される： $$ \begin{align} \mathrm{Minimize} \quad &f(x) \\ \mathrm{subject}\ \mathrm{to} \quad &g_{i}(x) \leqq 0 \quad (i=1,2,\ldots,m) \\ &h_{j}(x) = 0 \quad (j=1,2,\ldots,l) \\ \end{align} $$ ここに、各文字は以下を表す：
- $x$：$n$次元ベクトル(決定変数と呼ばれる)
- $f:\mathbb{R}^n \to \mathbb{R}$：目的関数
- $g_i:\mathbb{R}^n \to \mathbb{R}$：不等式制約関数
- $h_j:\mathbb{R}^n \to \mathbb{R}$：等式制約条件
最適化問題には様々な種類があり、それらは関数$f,\,g,\,h$の性質*1や決定変数$x\in\mathbb{R}^n$が取りうる値が連続値か離散値であるかなどによって分類される。例を挙げておこう。
| 最適化問題の名称 | 略称 | 決定変数 | 目的関数 | 不等式制約 | 等式制約 |
|---|---|---|---|---|---|
| 線形計画問題 | LP | 連続 | 線形 | 線形 | 線形 |
| 二次計画問題 | QP | 連続 | 二次関数 | 線形 | 線形 |
| 非線形最適化問題 | NLP*2 | 連続 | 非線形 | 非線形 | 非線形 |
| 凸最適化問題 | CO*3 | 連続 | 凸*4 | 凸 | 線形 |
| 整数線形計画問題 | ILP | 整数 | 線形 | 線形 | 線形 |
| 混合線形整数計画問題 | MLP | 連続 + 整数*5 | 線形 | 線形 | 線形 |
これらの他にも、制約条件に半正定値条件を含んだ半正定値計画問題(Semidefinite Programming Problem; SDP)や、制約条件に二次錐と呼ばれる構造を含んだ二次錐計画問題(Second Order Cone Programming Problem; SOCP)などが存在する。一般に、非線形最適化問題にはいくつもの局所的最適解*6が存在するため、大域的最適解を求めることは困難であり、多くのソルバーは局所的最適解を求解する*7。
CVXPY
上記で取り上げた数々の最適化問題に対する解法アルゴリズムの研究は盛んに行われており、近年は自分で実装を行う必要が少なくなってきている。CVXPYの公式ドキュメントによると、CVXPY は凸最適化問題を解くための Python 製モデリングツールであり、ソルバーとしてはオープンソースの ECOS，OSQP，SCS を使用しているようである。これら以外のソルバーにも対応しているが、別途インストールする必要がある。
導入
それでは、CVXPY をインストールしましょう。手順に関しては公式ドキュメントに書いてある。ここでは、Mac OS X の場合における例を示しておく。基本的には pip でインストールしてあげれば O.K. ですが、仮想環境にインストールしたい場合は pip ではなく、Anaconda を使ってインストールすることが推奨されているようだ。
$ pip install cvxpy
CVXPY をインストールすると、自動的に ECOS，OSQP，SCS もインストールされる。Numpy および Scipy がローカル環境にインストールされていない場合は、別途インストールする必要がある。また、各ライブラリの依存関係だが、現時点(2019年10月22日)では以下のようになっている：
Pyomo と比べると、ライブラリの導入だけでソルバーもインストールしてくれるので、手間は少ない。
数値実験
さて、ローカルに実行環境が構築されたので、実際にいくつかの最適化問題を CVXPY を使って解いていきたいと思う。今回は、以下の三種類の問題について取り扱う：
- 線形計画問題(Linear Programming Problem; LP)
- 最小二乗問題(Least-Squares Problem)
- 半正定値計画問題(Semidefinite Programming Problem; SDP)
また、実行環境は Jupyter を使用しており、numpyおよびcvxpyは以下のように実行済みであるとする。
$ import numpy as np$ import cvxpy
例1：線形計画問題(Linear Programming Problem; LP)
初めは、最も基本的な最適化問題である線形計画問題を解く。例として、次のような問題を考えよう。 $$ \begin{align} \mathrm{Minimize} \quad &2x_1 + 2x_2 \\ \mathrm{subject}\ \mathrm{to} \quad &3x_1 + 2x_2 \geqq 5 \\ &2x_1 + 6x_2 \geqq 8 \\ &x_1 \geqq 0,\, x_2 \geqq 0 \end{align} $$ この問題は、九州大学の神山先生のスライド*8で扱われていた問題です。行列とベクトルを用いて最適化問題を書き直すと次のようになる。 $$ \begin{align} \mathrm{Minimize} \quad &\begin{bmatrix}2 \\ 3\end{bmatrix}^\top \begin{bmatrix}x_1 \\ x_2\end{bmatrix} \\ \mathrm{subject}\ \mathrm{to} \quad &\begin{bmatrix} 3 & 2 \\ 2 & 6 \end{bmatrix} \begin{bmatrix}x_1 \\ x_2\end{bmatrix} \geqq \begin{bmatrix}5 \\ 8\end{bmatrix}, \quad \begin{bmatrix}x_1 \\ x_2\end{bmatrix} \geqq 0 \end{align} $$
Python によるサンプルコード
# Define the Data.
A = np.array([[3,2],[2,6]])
b = np.array([5,8])
c = np.array([2,3])
# Define and solve the CVXPY problem.
x = cvxpy.Variable(2) #決定変数の次元
objective = cvxpy.Minimize(c.T@x) #目的関数
constraints = [A@x >= b, x >= 0] #制約条件
prob = cvxpy.Problem(objective, constraints)
prob.solve()
# Print result.
print("\nThe optimal value is", prob.value) #最適値
print("A solution x is") #最適解
print(x.value)
print("A dual solution is") #双対問題の最適解
print(prob.constraints[0].dual_value)
出力は以下の通り：
The optimal value is 4.999999997381506
A solution x is
[1. 1.]
A dual solution is
[0.42857143 0.35714286]
例2：最小二乗問題(制約付き)
次に最小二乗問題について考える。公式のドキュメントでは、最もシンプルな制約のない最小二乗問題を扱っているが、ここでは制約条件を含んだ最小二乗問題を扱う。具体的には次のように記述される問題について考える。 $$ \begin{align} \mathrm{Minimize} \quad &\frac{1}{2} \| Cx - d \|_{2}^2 \\ \mathrm{subject\ to} \quad &Ax \leqq b \end{align} $$
係数$A,\,b,\,C,\,d$については、MathWorks のドキュメント*9で例として挙げられている数値を使う。
Python によるサンプルコード
# Define the data.
C = np.array([[0.9501, 0.7620, 0.6153, 0.4057],
[0.2311, 0.4564, 0.7919, 0.9354],
[0.6068, 0.0185, 0.9218, 0.9169],
[0.4859, 0.8214, 0.7382, 0.4102],
[0.8912, 0.4447, 0.1762, 0.8936]])
d = np.array([0.0578, 0.3528, 0.8131, 0.0098, 0.1388])
A = np.array([[0.2027, 0.2721, 0.7467, 0.4659],
[0.1987, 0.1988, 0.4450, 0.4186],
[0.6037, 0.0152, 0.9318, 0.8462]])
b = np.array([0.5251, 0.2026, 0.6721])
# Define and solve the CVXPY problem.
x = cvxpy.Variable(4) #決定変数の次元
objective = cvxpy.sum_squares(C*x - d) * 0.5 #目的関数
constraints = [A@x <= b] #制約条件
prob = cvxpy.Problem(cvxpy.Minimize(objective), constraints)
prob.solve()
# Print result.
print("\nThe optimal value is", prob.value)
print("The optimal x is")
print(x.value)
print("The norm of the residual is ", 0.5*cvxpy.norm(C*x - d, p=2).value)
出力は以下の通り：
The optimal value is 0.008792690396538318
The optimal x is
[ 0.12986199 -0.57569442 0.4251035 0.24384475]
The norm of the residual is 0.0663049447069704
例3：半正定値計画問題(Semidefinite Programming Problem; SDP)
最後に、半正定値計画問題について扱う。この問題の特筆すべき点は決定変数が行列であることであり、標準形は次のように記述される。
$$ \begin{align} \mathrm{Minimize} \quad & \left< C, X \right> \\ \mathrm{subject\ to} \quad &\left< A_j, X \right> = b_j \quad j=1,2,...,p \\ \quad &X \succcurlyeq O \end{align} $$
ここに、各文字は以下の通りである：
- $C$：$n\times n$ の定数行列
- $A_1, \, A_2, \ldots , A_p$：$n\times n$ の定数行列
- $b_1,\,b_2,\ldots,b_p$：定数
- $X$：$n\times n$の対称行列(決定変数)
- $X \succcurlyeq O$：正方行列$X$が半正定値行列*10であることを表す
また、$\left< C, X \right>$ は行列に対する内積を定義したものであり、以下で与えられる：
$$ \begin{align} \left< C, X \right> := \sum_{i=1}^{n} \sum_{j=1}^{n} C_{ij} X_{ij} \end{align} $$
Python によるサンプルコード
# Define the Data.
n = 3 #決定変数である行列Xのサイズ
p = 2 #制約条件の個数
C = np.array([[1, 2, 3],
[2, 9, 0],
[3, 0, 7]])
A = []
A.append(np.array([[1, 0, 1],[0, 3, 7],[1, 7, 5]]))
A.append(np.array([[0, 2, 8],[2, 6, 0],[8, 0, 4]]))
b = []
b.append(11)
b.append(19)
# Define and solve the CVXPY problem.
# Create a symmetric matrix variable.
X = cvxpy.Variable((n,n), symmetric=True)
# Define the objective and constraints
objective = cvxpy.Minimize(cvxpy.trace(C@X))
constraints = []
constraints = [X >> 0] # The operator >> denotes matrix inequality.
constraints += [cvxpy.trace(A[i]@X) == b[i] for i in range(p)]
prob = cvxpy.Problem(objective, constraints)
prob.solve()
# Print result.
print("The optimal value is", prob.value)
print("A solution X is")
print(X.value)
出力は次の通り：
The optimal value is 13.902227826690467
A solution X is
[[1.05592547 0.36918783 0.86829605]
[0.36918783 0.12908075 0.30358614]
[0.86829605 0.30358614 0.71400686]]
おわりに
最後に、今回の記事を作成するにあたって参考にした Web サイトを載せておく。
-
- CVXPY に関する先行記事。本記事では扱わなかった CVXPY で使うことのできる関数について扱っているので、本記事で分からないことがあったらこちらも参照されたい。
-
- 例1で取り上げた問題はこちらのスライドから流用させていただいた。関係ないが、このスライドは LaTeX Beamer で作成されていると思われるが、素晴らしくシンプルで美しい体裁をしていて感動した。
-
- 整数計画問題についてあまり勉強したことがないので、細かい分類を確認するために参照させていただいた。今回は全く取り上げなかったが、今後は離散凸解析にも明るくなりたい。
-
- 数理計画に関する用語が解説されている。かなり詳しくかかれているので、最適化に関する用語で分からないものがあったら、とりあえずこちらを参照すれば良いと思う。
追記情報
- 目次(Contents)を追加．(191103)
- コードのデザインを修正．(191103)
*1:一般にこれらの関数は実数値をとる微分可能な連続関数と仮定されることが多いが、応用上、微分不可能な関数を考えることもままある
*2:決して自然言語処理(Natural Language Processing; NLP)ではない
*3:“Convex Optimization”の略だが、あまり略されているところを見たことはない
*4:関数$f$が凸であるとは、エピグラフ(epigraph)が凸集合になるような関数である。これと同値な条件として、任意の二点$x,\,y$と任意の実数$\lambda \in (0,1)$に対して、$f(\lambda x + (1-\lambda)y) \leqq \lambda f(x) + (1-\lambda) f(y)$を満たすことを言う。
*5:連続変数と離散変数が混在することを表す
*6:最適化問題に対する「最適解」は「局所的最適解(local optimal solution)」と「大域的最適解(global optimal solution)」の二種類に大別される。最適化問題$\,\mathrm{min} \, f(x) \, \mathrm{s.t.} \, x \in \mathcal{F}\,$において、実行可能解$\,x^{\star} \in \mathcal{F}\,$が任意の近傍$\mathcal{B}(x^{\star}, \varepsilon) := \{ x \mid \|x - x^{\star}\| \leqq \varepsilon \}$に対して、$f(x^{\star}) \leqq f(x) \, (\forall x \in \mathcal{B}(x^{\star}, \varepsilon) )$を満たすとき、$x^{\star}\,$を局所的最適解という。一方、実行可能解$x^{\ast}$が考えている問題の任意の実行可能解に対して、$f(x^{\ast}) \leqq f(x) \, (\forall x \in \mathcal{F})$を満たすとき、$x^{\ast}$は大域的最適解という。
*10:正方行列$X\in\mathbb{R}^{n\times n}$が半正定値行列であるとは、任意の$n$次元ベクトル$d \in \mathbb{R}^n$に対して、二次形式が非負、すなわち$d^{\top} X d \geqq 0$ を満たすことを言う
