この記事は、学生時代にお世話になった数理最適化論の世界でしょっちゅう目にするであろう、2次錐計画問題(Second Order Cone Programming、以後SOCPと記述)の簡単な紹介と、それをPythonを使って解いてみる実装例を述べたいと思います。理工系の大学1年で習得する線形代数の知識があれば、数式を含めて理解できる記事の構成にしました。
なお、私が大学院から離れてから現在までに仕事上全く数理最適化論に関わっておらず、当時の知識を思い出しながら記事を書きました。そのため、誤った記述がありましたらコメントで指摘して頂けると嬉しいです。
#SOCPの解説
そもそもSOCPって何？
1文で説明すると
$n$個の変数$\mathbf{x}:= ( x_1 , \ldots , x_n )$において、
$i = 1 , \ldots , m$に対し、ある与えられた定数$\mathbf{A}_i \in \mathbb{R}^{n \times n_i} , \mathbf{b}_i \in \mathbb{R}^{n_i} , \mathbf{c}_i \in \mathbb{R}^n , d_i \in \mathbb{R} , \mathbf{e} \in \mathbb{R}^n$に関して、
$m$個の2次錐制約条件$|| \mathbf{A}_i \mathbf{x} + \mathbf{b}_i || \leq \mathbf{c}_i^\top \mathbf{x} + d_i$を満たすとき、
線形の目的関数$\mathbf{e}^\top \mathbf{x}$の値を最小化しなさい
…みたいな数学の問題です。数学的には、これを以下のように記述できます。
\begin{aligned}
& \text{minimize} & & \mathbf{e}^\top \mathbf{x} \\
& \text{subject to} & & \left\| \mathbf{A}_i \mathbf{x} + \mathbf{b}_i \right\| \leq \mathbf{c}_i^\top \mathbf{x} + d_i \ \ \ (i = 1 , \ldots , m)
\end{aligned}
なお、この記事ではユークリッドノルムの定義として、ノルム$|| \cdot ||$を記述します。上記よりさらに一般的にしたSOCPの記述もありますが、この記事では有限次元上の実数空間に絞って記述していきます。
##ようわからん
数理最適化論に触れたことのないほとんどの人は、上記の説明だけでは一体なんのことやら…と思っているはずです。ここで、高校数学の教科書に書いてある内容を思い出しましょう。例えば、以下のような線形計画問題(Linear Programming、以後LPと記述)を思い出してみてください。
変数$x_1 , x_2$において、
$x_1 \geq 0 , x_2 \geq 0 , 2 x_1 + x_2 \leq 6 , x_1 + 3 x_2 \leq 6$をすべて満たすとき、
$5 x_1 + 6 x_2$の最大値を求めなさい
このような問題に対し、解き方を忘れた人向けに解説しましょう。まず、以下のように$x_1 \geq 0 , x_2 \geq 0 , 2 x_1 + x_2 \leq 6 , x_1 + 3 x_2 \leq 6$を満たす領域を図示します。
次に、最大にしたい$5 x_1 + 6 x_2$に対し、$k = 5 x_1 + 6 x_2$を満たすように定数$k$を用意します。この式を変形すると$x_2 = - \frac{5}{6} x_1 + \frac{k}{6}$となるので、先程の領域に含まれるように傾き$- \frac{5}{6}$の直線を平行移動させ、その$x_2$切片である$\frac{k}{6}$が最も大きくなる場合を考えます。
上記の大まかな説明を図で表わしてみました。水色の領域に含まれるように傾き$- \frac{5}{6}$の赤色の破線を徐々に上へ平行移動させ、Max Valueの矢印が示されるように、最も$x_2$切片の大きくなる場合が図での赤色の太い直線となっています。その場合での$x_1$と$x_2$の値は、水色の領域と赤色の太い直線との交点、すなわち$2 x_1 + x_2 = 6$と$x_1 + 3 x_2 = 6$との交点である$x_1 = \frac{12}{5} , x_2 = \frac{6}{5}$となり、これがそのまま$5 x_1 + 6 x_2$の最大値を与える点となります。よって、最大値は$5 \times \frac{12}{5} + 6 \times \frac{6}{5} = \frac{96}{5}$と求められます。
##LPはすべてSOCPに書き換えられる！
以上の問題は、
\begin{aligned}
& \text{maximize} & & 5 x_1 + 6 x_2 \\
& \text{subject to} & & x_1 \geq 0 \\
& & & x_2 \geq 0 \\
& & & 2 x_1 + x_2 \leq 6 \\
& & & x_1 + 3 x_2 \leq 6
\end{aligned}
と数学的に記述できます。$5 x_1 + 6 x_2$を最大化することは、$- 5 x_1 - 6 x_2$を最小化することと明らかに同値なので、最初の一行目は以下の表記と全く同じ意味です。
\text{miniimize} \ \ \ \ \ \ - 5 x_1 - 6 x_2
ここで、もう一度SOCPの数学的な記述を見てみましょう。
\begin{aligned}
& \text{minimize} & & \mathbf{e}^\top \mathbf{x} \\
& \text{subject to} & & \left\| \mathbf{A}_i \mathbf{x} + \mathbf{b}_i \right\| \leq \mathbf{c}_i^\top \mathbf{x} + d_i \ \ \ (i = 1 , \ldots , m)
\end{aligned}
1行目の数式を見ると$\mathbf{e}^\top \mathbf{x}$となっており、この2つのベクトルを$\mathbf{e} = (-5, -6)$と$\mathbf{x} = ( x_1, x_2 )$とおけば、内積を計算すること先程のLPの問題に帰着できることがわかります。また、2行目の数式は、$|| \mathbf{A}_i \mathbf{x} + \mathbf{b}_i || \leq \mathbf{c}_i^\top \mathbf{x} + d_i$という不等式の制約条件が$m$個ありますよ、という意味なので、これもLPの問題で出てきた$( m = ) 4$つの不等式に帰着できそうです。実際に$i = 1$番目の不等式である$x_1 \geq 0$は、
\mathbf{A}_1 =
\begin{pmatrix}
0 & 0
\end{pmatrix}
, \mathbf{b}_1 =
\begin{pmatrix}
0
\end{pmatrix}
, \mathbf{c}_1 =
\begin{pmatrix}
1 \\
0
\end{pmatrix}
, d_1 = 0
とおいて、$|| \mathbf{A}_1 \mathbf{x} + \mathbf{b}_1 || \leq \mathbf{c}_1^\top \mathbf{x} + d_1$に代入すると一致しますし、$\mathbf{c}_1$の0と1を逆にすることで$x_2 \geq 0$も表わすことができます。$i = 3$番目の不等式である$2 x_1 + x_2 \leq 6$は、例えば
\mathbf{A}_3 =
\begin{pmatrix}
0 & 0
\end{pmatrix}
, \mathbf{b}_3 =
\begin{pmatrix}
0
\end{pmatrix}
, \mathbf{c}_3 =
\begin{pmatrix}
-2 \\
-1
\end{pmatrix}
, d_3 = 6
とおいて代入すると一致しますし(これ以外のおき方もあります)、最後の$x_1 + 3 x_2 \leq 6$も$\mathbf{c}_3$の値を変えれば一致させることができます。以上より、$|| \mathbf{A}_i \mathbf{x} + \mathbf{b}_i || \leq \mathbf{c}_i^\top \mathbf{x} + d_i$の不等式を特別に$ \mathbf{A}_i \mathbf{x} = \mathbf{0}$とすればLPの任意の不等式制約を表わせられるので、すべてLPはSOCPとして定式化することが可能であると言えます(本当は$ \mathbf{A}_i \mathbf{x} + \mathbf{b}_i = \mathbf{0}$でも大丈夫ですが、$ \mathbf{A}_i \mathbf{x} = \mathbf{0}$とするだけで不等式制約の左辺は定数項となりますので、$d_i$と同様に考えることができます)。
##んじゃ、LPに帰着できないSOCPって何よ？
すべてのLPはSOCPに帰着できると例を挙げて説明しましたが、その逆は成り立ちません。つまり、反例としてLPに帰着できないSOCPの問題が存在することを意味します。まぁ、$ \mathbf{A}_i \mathbf{x} \neq \mathbf{0}$の場合とか明らかに怪しいですもんね。以下の例を考えてみましょう。
\mathbf{A}_i =
\begin{pmatrix}
1 & 0 \\
0 & 0
\end{pmatrix}
, \mathbf{b}_i =
\begin{pmatrix}
0 \\
1
\end{pmatrix}
, \mathbf{c}_i =
\begin{pmatrix}
0 \\
1/4
\end{pmatrix}
, d_i = 0
これを不等式制約$|| \mathbf{A}_i \mathbf{x} + \mathbf{b}_i || \leq \mathbf{c}_i^\top \mathbf{x} + d_i$に代入すると、$\sqrt{ x_1^2 + 1 } \leq \frac{1}{4} x_2$と計算できます。この不等式は両辺を2乗して移項すると$\frac{1}{16} x_2^2 - x_1^2 \geq 1$となりますので、以下のように双曲線で区切られた上側の領域としてLPに帰着できない不等式制約を図示することができます。
なお、双曲線の下側の薄く描かれた領域は$- \frac{1}{4} x_2 \leq \sqrt{ x_1^2 + 1 }$の領域なので間違いです($\frac{1}{16} x_2^2 - x_1^2 \geq 1 \Longleftrightarrow - \frac{1}{4} x_2 \leq \sqrt{ x_1^2 + 1 } \leq \frac{1}{4} x_2$を満たしますが、$\frac{1}{16} x_2^2 - x_1^2 \geq 1 \Longleftarrow \sqrt{ x_1^2 + 1 } \leq \frac{1}{4} x_2$を満たすことに注意です)。
##変数が2個だけとかつまんねぇーーーーー
もう1つだけ不等式制約の例を挙げてみましょう。しかし、ここから変数の数を$(n=)3$個にしてみましょう。
\mathbf{A}_i =
\begin{pmatrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 0
\end{pmatrix}
, \mathbf{b}_i =
\begin{pmatrix}
0 \\
0 \\
0
\end{pmatrix}
, \mathbf{c}_i =
\begin{pmatrix}
0 \\
0 \\
1
\end{pmatrix}
, d_i = 0
先程と同様に代入すると、$\sqrt{x_1^2 + x_2^2} \leq x_3$となります。この不等式によって表される領域は、以下の図で示される青い面の内側にあたります。
まるでアイスクリームのコーンみたいな形ですね。このように、コーン型などのような2次形式の不等式制約で表わされる領域を$m$個用意し、それらすべて重なった共通境域上で線形目的関数の値を最小化する問題がSOCPなんだと思って下さい。
##ちょっと待って、等式制約があるSOCPの一般形があったんだけど
記事を投稿した当時のWikipediaでは、以下のようにSOCPの一般形を記述されていました。
\begin{aligned}
& \text{minimize} & & \mathbf{e}^\top \mathbf{x} \\
& \text{subject to} & & \left\| \mathbf{A}_i \mathbf{x} + \mathbf{b}_i \right\| \leq \mathbf{c}_i^\top \mathbf{x} + d_i \ \ \ (i = 1 , \ldots , m) \\
& & & \mathbf{Fx} = \mathbf{g}
\end{aligned}
ここで、$\mathbf{F} \in \mathbb{R}^{l \times n} , \mathbf{g} \in \mathbb{R}^l$は定数の行列(ベクトル)です。先程のSOCPの一般形に$\mathbf{Fx} = \mathbf{g}$が追加されていますが、これは$2 x_1 - x_2 + 4 x_3 = 3$のような線形方程式で表わされる制約が計$l$個存在していることを意味します。その$l$個の線形方程式を、行列とベクトルを用いて$\mathbf{Fx} = \mathbf{g}$のように簡潔に表現しています。
さて、どっちがSOCPの一般形なのかと言いますと、どっちもSOCPの一般形として正しいです。なぜなら、$\mathbf{Fx} = \mathbf{g}$は、実質$\mathbf{Fx} \leq \mathbf{g}$かつ$\mathbf{Fx} \geq \mathbf{g}$であるからです。そしてLPをSOCPに帰着させるように、これら2つの線形不等式は$|| \mathbf{A}_i \mathbf{x} + \mathbf{b}_i || \leq \mathbf{c}_i^\top \mathbf{x} + d_i$における$ \mathbf{A}_i \mathbf{x} = \mathbf{0}$のパターンに帰着できます。無理矢理かもしれませんが、$\mathbf{Fx} = \mathbf{g}$は入れても、入れなくても、どちらでも構わないのです。
#SOCPの求解方法
##高校数学の知識だけではSOCPに太刀打ちできん！
上記のような変数が2個のLPに対する解き方は高校数学の教科書レベルで対応できますが、より一般的な変数が$n$個のSOCPの場合になると無理になる問題の方が圧倒的に多くなります。よって高校数学とは異なったアプローチを使ってSOCPを解くことになるのですが、現在SOCPへのアプローチとして最も広く知られている主双対内点法(Primal-Dual Interior Point Method)とよばれるアルゴリズム群を簡単に紹介していこうと思います。
主双対内点法とは、不等式制約を満たすある初期点を用意し、そこから反復して点を更新していき、主問題と双対問題を同時に解くアルゴリズムとなります。主双対内点法の最も重要な特徴は、定数として与えられる$\mathbf{A}_i , \mathbf{b}_i , \mathbf{c}_i , d_i , \mathbf{e}$のサイズである、$n , n_i , m$の多項式で表わされる反復回数によって最適解に収束する点です。この特徴はLPを解く際に使われるシンプレックス法や、メタヒューリスティクスなアルゴリズムである遺伝的アルゴリズムを凌駕する性質です。
なお、これまでに出てきた主問題、双対問題のようなキーワードを知らなくても、とりあえずSOCPを解いてみる程度であれば深追いする必要はないので、ここでは触れていません。この近辺の導入知識は、例えばこれなら分かる最適化数学―基礎原理から計算手法までの本などに詳しく書かれています。また、SOCPを解く主双対内点法の具体的な仕組みもここでは触れないようにします。この知識は、理工系の大学1年で習得する線形代数などの知識をベースとした高度な数学的知識を要求されるからです。興味のある人は、内点法(経営科学のニューフロンティア)などが参考になるのではないかと思います。
##んじゃ、どうしろと？
手っ取り早い話、今まで説明してきた定数の$\mathbf{A}_i , \mathbf{b}_i , \mathbf{c}_i , d_i , \mathbf{e}$を入力したら、そこから先は自動的にSOCPを解いて、最適解などを出力してくれるライブラリがあれば良いのです。機械学習やディープラーニングの世界でもそうだと思いますが、数学的な理論を理解し、一からアルゴリズムをちまちまと実装するのは学習上は良いのでしょうが、限られた時間の上では骨折りの作業と化します。なので、事前に用意した定数を放り込めば、あとは勝手にSOCPを解いてくれるようなライブラリを使えるようになれることが一番大事です。
そんなライブラリとして、例えばPython言語をサポートしているものの一部を、以下に簡単にまとめました。
| ライブラリ名 | 特徴 |
|---|---|
| Gurobi | Gurobi Optimization社による商用ライブラリ。Python、C++、MATLABなどの言語をサポート。 |
| CPLEX | IBM社による商用ライブラリ(無料版も存在する)。Python、Java、C#などの言語をサポート。 |
| SCIP | 混合整数計画問題にも強いOSS。Python、Java、MATLABなどの言語をサポート。 |
| CVXOPT | インストールが容易で、凸計画問題の求解に特化したOSS。 |
| OpenOpt | 非線形計画問題などにも対応できるOSS。 |
本記事では、インストールが容易なOSSであるCVXOPTに絞ってSOCPを解いていきます。インストール方法は、こちらからどうぞ。
##CVXOPTの行列構造はNumPyのやつと違うぞ！
CVXOPXで独自に定義された行列のインスタンスを用いてSOCPを扱うのですが、その行列のデータ構造に気をつけなければなりません。NumPyなどの行列データは行志向でデータを格納しますが、CVXOPTでは列志向でデータを格納します。以下の例のように、NumPyのnumpy.array、CVXOPTのcvxopt.matrixクラスにより、2重にネストされたリストなどを用いて行列を表せられますが、その構造は完全に逆であることに注意が必要です。
>>> import numpy as np
>>> import cvxopt as co
>>> A = np.array([[1, 2], [3, 4], [5, 6]])
>>> print(A)
[[1 2]
[3 4]
[5 6]]
>>> B = co.matrix([[1, 2], [3, 4], [5, 6]])
>>> print(B)
[ 1 3 5]
[ 2 4 6]
*による行列同士の算術演算子にも注意が必要です。NumPyでは2つの行列の各要素ごとに積を計算します。
>>> import numpy as np
>>> A = np.array([[0, 1], [2, 3]])
>>> B = np.array([[3, 2], [1, 0]])
>>> print(A)
[[0 1]
[2 3]]
>>> print(B)
[[3 2]
[1 0]]
>>> print(A*B)
[[0 2]
[2 0]]
ところが、CVXOPTでは2つの行列そのものの積を計算します。
>>> import cvxopt as co
>>> C = co.matrix([[0, 2], [1, 3]])
>>> D = co.matrix([[3, 1], [2, 0]])
>>> print(C)
[ 0 1]
[ 2 3]
>>> print(D)
[ 3 2]
[ 1 0]
>>> print(C*D)
[ 1 0]
[ 9 4]
なお、NumPyでprint(C*D)の出力結果のような行列を出力したい場合はnp.dot(A,B)とすれば良いんでしたね。
また、以下の例のようにNumPyで定義した行列構造をCVXOPTで定義した行列構造に変換することができます。
>>> import numpy as np
>>> import cvxopt as co
>>> A = np.array([[1, 2], [3, 4], [5, 6]])
>>> print(A)
[[1 2]
[3 4]
[5 6]]
>>> B = co.matrix(A)
>>> print(B)
[ 1 2]
[ 3 4]
[ 5 6]
逆の変換も可能です。
>>> import numpy as np
>>> import cvxopt as co
>>> C = co.matrix([[1, 2], [3, 4], [5, 6]])
>>> print(C)
[ 1 3 5]
[ 2 4 6]
>>> D = np.array(C)
>>> print(D)
[[1 3 5]
[2 4 6]]
このとき、NumPyの行志向の行列構造とCVXOPTの列志向の行列構造が入れ替わらず、そのまま行数と列数が維持されることに注意が必要です。その他諸々の使い方はこちらを参考にして下さい。
##そんなことより、早く実装例を教えてくれ
ぶっちゃけ、ここにすべて書いてあります…が、一番最初に説明したSOCPの一般形とは結構異なる形でSOCPを表わしているように見えるので、正直分かりにくいです。こういう時は、例題を見ながら1つずつ値を突っ込んでみるのが常套手段でしょう。3変数$x_1 , x_2 , x_3$に関する以下の例題を考えてみましょう。
\begin{aligned}
& \text{minimize} & & - x_1 + 3 x_2 + 2 x_3\\
& \text{subject to} & & \sqrt{ 4 x_1^2 + 3 x_2^2 + x_3^2 - 2 \sqrt{3} x_2 x_3 - 16 x_1 - 6 x_2 + 2 \sqrt{3} x_3 + 19 } \leq x_2 + \sqrt{3} x_3 - 1 \\
& & & x_1 - x_2 \geq 0
\end{aligned}
上記の1つ目の不等式制約は、先程説明したアイスクリームのコーンのような形をした2次形式の不等式$\sqrt{x_1^2 + x_2^2} \leq x_3$において、$x_2 x_3$平面で$\pi / 6$だけ回転させ、さらに$x_1$の方向へ2、$x_2$の方向へ1だけ平行移動させたものになります。興味のある人は調べてみてください。
さて、SOCPの一般形に持ち込む必要があるので、1つ目の不等式制約の左辺の平方根の中身を以下のように2乗和になるように変形させます。
\begin{aligned}
& 4 x_1^2 + 3 x_2^2 + x_3^2 - 2 \sqrt{3} x_2 x_3 - 16 x_1 - 6 x_2 + 2 \sqrt{3} x_3 + 19 \\
= & 4 ( x_1^2 - 4 x_1 ) + ( 3 x_2^2 - 2 \sqrt{3} x_2 x_3 + x_3^2 ) - 6 x_2 + 2 \sqrt{3} x_3 + 19 \\
= & 4 \{ (x_1 - 2)^2 - 4 \} + ( \sqrt{3} x_2 - x_3 )^2 - 6 x_2 + 2 \sqrt{3} x_3 + 19 \\
= & 4 (x_1 - 2)^2 + ( \sqrt{3} x_2 - x_3 )^2 - 2 ( \sqrt{3} x_2 - x_3 ) \sqrt{3} + 3 \\
= & ( 2 x_1 - 4 )^2 + ( \sqrt{3} x_2 - x_3 - \sqrt{3} )^2
\end{aligned}
これより1つ目の不等式制約は、以下のように変形できます。
\begin{aligned}
\sqrt{ 4 x_1^2 + 3 x_2^2 + x_3^2 - 2 \sqrt{3} x_2 x_3 - 16 x_1 - 6 x_2 + 2 \sqrt{3} x_3 + 19 } & \leq x_2 + \sqrt{3} x_3 - 1 \\
\sqrt{ ( 2 x_1 - 4 )^2 + ( \sqrt{3} x_2 - x_3 - \sqrt{3} )^2 } & \leq x_2 + \sqrt{3} x_3 - 1 \\
\left\|
\begin{pmatrix}
2 x_1 - 4 \\
\sqrt{3} x_2 - x_3 - \sqrt{3}
\end{pmatrix}
\right\| & \leq x_2 + \sqrt{3} x_3 - 1 \\
\left\|
\begin{pmatrix}
2 & 0 & 0 \\
0 & \sqrt{3} & -1
\end{pmatrix}
\begin{pmatrix}
x_1 \\
x_2 \\
x_3
\end{pmatrix}
+
\begin{pmatrix}
-4 \\
- \sqrt{3}
\end{pmatrix}
\right\| & \leq
\begin{pmatrix}
0 & 1 & \sqrt{3}
\end{pmatrix}
\begin{pmatrix}
x_1 \\
x_2 \\
x_3
\end{pmatrix}
+ (-1)
\end{aligned}
2つ目の不等式制約$x_1 - x_2 \geq 0$は以下のように考えると、$\mathbf{A}_2 , \mathbf{b}_2$が零行列(ベクトル)になるので簡単です。
\left\|
\begin{pmatrix}
0 & 0 & 0
\end{pmatrix}
\begin{pmatrix}
x_1 \\
x_2 \\
x_3
\end{pmatrix}
+
\begin{pmatrix}
0
\end{pmatrix}
\right\| \leq
\begin{pmatrix}
1 & -1 & 0
\end{pmatrix}
\begin{pmatrix}
x_1 \\
x_2 \\
x_3
\end{pmatrix}
+ 0
目的関数は内積の形に変形します。
- x_1 + 3 x_2 + 2 x_3 =
\begin{pmatrix}
-1 & 3 & 2
\end{pmatrix}
\begin{pmatrix}
x_1 \\
x_2 \\
x_3
\end{pmatrix}
よって、以下のように考えることができます。
\begin{aligned}
\mathbf{A}_1 =
\begin{pmatrix}
2 & 0 & 0 \\
0 & \sqrt{3} & -1
\end{pmatrix}
, \mathbf{b}_1 =
\begin{pmatrix}
-4 \\
- \sqrt{3}
\end{pmatrix}
, \mathbf{c}_1 =
\begin{pmatrix}
0 \\
1 \\
\sqrt{3}
\end{pmatrix}
, d_1 = -1
, \\
\mathbf{A}_2 =
\begin{pmatrix}
0 & 0 & 0
\end{pmatrix}
, \mathbf{b}_2 =
\begin{pmatrix}
0
\end{pmatrix}
, \mathbf{c}_2 =
\begin{pmatrix}
1 \\
-1 \\
0
\end{pmatrix}
, d_2 = 0
, \mathbf{e} =
\begin{pmatrix}
-1 \\
3 \\
2
\end{pmatrix}
\end{aligned}
では、CVXOPTを用いてSOCPを解いてみましょう。
>>> import cvxopt as co
>>> import numpy as np
>>> from math import sqrt
>>>
>>> e = np.array([-1.0, 3.0, 2.0])
>>> c = co.matrix(e)
>>>
>>> g1 = np.array([[0.0, -1.0, -sqrt(3)], [-2.0, 0.0, 0.0], [0.0, -sqrt(3), 1.0]])
>>> g2 = np.array([[-1.0, 1.0, 0.0], [0.0, 0.0, 0.0]])
>>> g = [co.matrix(g1), co.matrix(g2)]
>>>
>>> h1 = np.array([-1.0, -4.0, -sqrt(3)])
>>> h2 = np.array([0.0, 0.0])
>>> h = [co.matrix(h1), co.matrix(h2)]
>>>
>>> sol = co.solvers.socp(c, Gq=g, hq=h)
pcost dcost gap pres dres k/t
0: 1.6667e+00 3.3333e+00 7e+00 4e-01 9e-01 1e+00
1: 1.0249e+00 1.1074e+00 2e-01 2e-02 4e-02 6e-02
2: 1.0001e+00 1.0011e+00 1e-03 2e-04 4e-04 7e-04
3: 1.0000e+00 1.0000e+00 1e-05 2e-06 4e-06 7e-06
4: 1.0000e+00 1.0000e+00 1e-07 2e-08 4e-08 7e-08
Optimal solution found.
>>> print(sol['x'])
[ 2.00e+00]
[ 1.00e+00]
[ 6.00e-08]
最初の変数であるeで目的関数の係数ベクトル$\mathbf{e}$の定義を行っており、eからcvxopt.matrixの行列の形式に変換したものをcに代入しています。値は整数値の係数でも浮動小数点型として入れるようにして下さい。
注意するのは次以降の変数のg1とg2です。これも浮動小数点型として入れるのですが、2つの不等式制約の係数行列(ベクトル)$\mathbf{A}_1 , \mathbf{A}_2 , \mathbf{c}_1 , \mathbf{c}_2$に対し、以下の式で表される行列をそれぞれg1とg2に代入しています。
\begin{pmatrix}
- \mathbf{c}_1^\top \\
- \mathbf{A}_1
\end{pmatrix}
=
\begin{pmatrix}
0 & -1 & - \sqrt{3} \\
-2 & 0 & 0 \\
0 & - \sqrt{3} & 1
\end{pmatrix}
,
\begin{pmatrix}
- \mathbf{c}_2^\top \\
- \mathbf{A}_2
\end{pmatrix}
=
\begin{pmatrix}
-1 & 1 & 0 \\
0 & 0 & 0
\end{pmatrix}
なぜ、このような行列になるのかはここの数式を見れば把握できますが、ここでは触れないようにします。
これより、上記のg1とg2をリストを用いて変数gに格納しています。gをco.matrixのリストにする理由は複数の不等式制約に対応するためであり、リストの0番目に1番目の不等式制約の係数の一部の情報であるg1を、リストの1番目にg2を入れています。
gの次の変数はh1とh2であり、残りの係数ベクトルである$\mathbf{b}_1 , \mathbf{b}_2 , d_1 , d_2$を用いて、以下のようにそれぞれ代入しています。
\begin{pmatrix}
d_1 \\
\mathbf{b}_1
\end{pmatrix}
=
\begin{pmatrix}
-1 \\
-4 \\
- \sqrt{3}
\end{pmatrix}
,
\begin{pmatrix}
d_2 \\
\mathbf{b}_2
\end{pmatrix}
=
\begin{pmatrix}
0 \\
0
\end{pmatrix}
これも、上記の行列になる理由は触れません。gと同様、h1とh2もリストを使ってhに格納します。
こうして定義した3つの変数c、g、hをcvxopt.solvers.socp()の引数に指定して、SOCPを解いてくれる流れになっています。解いた結果を変数solに返します。solは辞書となっており、$x_1 , x_2 , x_3$の最適解はsol['x']で分かります。sol['x']の出力結果より、$x_1$が2.00e+00、$x_2$が1.00e+00、$x_3$が6.00e-08に対応します。しかし、主双対内点法は反復的に最適解を求めるので、6.00e-08は数値誤差であり、0であるとみなします。
以上より、最適解は$x_1 = 2 , x_2 = 1 , x_3 = 0$となり、最小値は$-1 \times 2 + 3 \times 1 + 2 \times 0 = -1$となります。
##んじゃ、適当な値を突っ込んでみたら何でも解けるんだな！！！
残念ながら、そんなわけがありません。2変数$x_1 , x_2$に関する以下の例題を見て下さい(式変形を以後省略し、SOCPの一般形で表すようにしています)。
\begin{aligned}
& \text{minimize} & &
\begin{pmatrix}
1 & 1
\end{pmatrix}
\begin{pmatrix}
x_1 \\
x_2
\end{pmatrix}
\\
& \text{subject to} & & \left\|
\begin{pmatrix}
1 & 0 \\
0 & \frac{1}{2}
\end{pmatrix}
\begin{pmatrix}
x_1 \\
x_2
\end{pmatrix}
+
\begin{pmatrix}
0 \\
0
\end{pmatrix}
\right\| \leq
\begin{pmatrix}
0 & 0
\end{pmatrix}
\begin{pmatrix}
x_1 \\
x_2
\end{pmatrix}
+ 1 \\
& & & \left\|
\begin{pmatrix}
0 & \frac{1}{\sqrt{3}} \\
0 & 0
\end{pmatrix}
\begin{pmatrix}
x_1 \\
x_2
\end{pmatrix}
+
\begin{pmatrix}
0 \\
1
\end{pmatrix}
\right\| \leq
\begin{pmatrix}
\frac{1}{2} & 0
\end{pmatrix}
\begin{pmatrix}
x_1 \\
x_2
\end{pmatrix}
+ 0
\end{aligned}
このSOCPをCVXOPTで求めてみると、以下のようになります。
>>> import cvxopt as co
>>> import numpy as np
>>> from math import sqrt
>>>
>>> e = np.array([1.0, 1.0])
>>> c = co.matrix(e)
>>>
>>> g1 = np.array([[0.0, 0.0], [-1.0, 0.0], [0.0, -0.5]])
>>> g2 = np.array([[-0.5, 0.0], [0.0, -1.0 / sqrt(3)], [0.0, 0.0]])
>>> g = [co.matrix(g1), co.matrix(g2)]
>>>
>>> h1 = np.array([1.0, 0.0, 0.0])
>>> h2 = np.array([0.0, 0.0, 1.0])
>>> h = [co.matrix(h1), co.matrix(h2)]
>>>
>>> sol = co.solvers.socp(c, Gq=g, hq=h)
pcost dcost gap pres dres k/t
0: 0.0000e+00 -2.1725e+00 1e+01 2e+00 8e-01 1e+00
1: 1.2143e+00 1.1223e+00 3e+00 4e-01 2e-01 6e-01
2: -1.2088e+01 2.9060e+01 1e+04 2e+01 7e+00 7e+01
3: 2.8664e-01 1.6374e+02 3e+03 1e+00 5e-01 2e+02
4: 2.8957e-01 1.6456e+04 3e+05 1e+00 5e-01 2e+04
5: 2.8958e-01 1.6457e+06 3e+07 1e+00 5e-01 2e+06
6: 2.8958e-01 1.6457e+08 3e+09 1e+00 5e-01 2e+08
Certificate of primal infeasibility found.
>>> print(sol['x'])
None
co.solvers.socp(c, Gq=g, hq=h)を実行すると、最後にCertificate of primal infeasibility found.と出力されており、sol['x']にはNoneが格納されている事が分かります。詳細な説明は省きますが、これは主問題の不等式制約をすべて満たす領域である実行可能領域が存在しないことを意味します。実際、このSOCPの不等式制約を変形しますと$\sqrt{x_1^2 + \frac{x_2^2}{4}} \leq 1$と$\sqrt{\frac{1}{3} x_2^2 + 1} \leq \frac{1}{2} x_1$となり、それぞれを図示すると以下の青色と赤色の領域になります。
この図を見ると分かりますが、2つの領域の共通領域がありません。要するに2つの不等式制約を同時に満たす$X_1 , x_2$の組は1つも存在しない(=実行不能、Infeasibleである)ことを意味します。これだと、当然最適解なんて求められませんよね。ですので、sol[x]もNoneを格納しているわけです。
では、Certificate of primal infeasibility found.が出力されなければ大丈夫…というわけにもいきません。今度は別のSOCPの例題を見てみましょう。
\begin{aligned}
& \text{minimize} & &
\begin{pmatrix}
-1 & -1 & -1
\end{pmatrix}
\begin{pmatrix}
x_1 \\
x_2 \\
x_3
\end{pmatrix}
\\
& \text{subject to} & & \left\|
\begin{pmatrix}
1 & 0 & 0 \\
0 & 1 & 0
\end{pmatrix}
\begin{pmatrix}
x_1 \\
x_2 \\
x_3
\end{pmatrix}
+
\begin{pmatrix}
0 \\
0
\end{pmatrix}
\right\| \leq
\begin{pmatrix}
0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
x_1 \\
x_2 \\
x_3
\end{pmatrix}
+ 0
\end{aligned}
このSOCPも、以下のようにCVXOPTで求めることができません。
>>> import cvxopt as co
>>> import numpy as np
>>>
>>> e = np.array([-1.0, -1.0, -1.0])
>>> c = co.matrix(e)
>>>
>>> g1 = np.array([[0.0, 0.0, -1.0], [-1.0, 0.0, 0.0], [0.0, -1.0, 0.0]])
>>> g = [co.matrix(g1)]
>>>
>>> h1 = np.array([0.0, 0.0, 0.0])
>>> h = [co.matrix(h1)]
>>>
>>> sol = co.solvers.socp(c, Gq=g, hq=h)
pcost dcost gap pres dres k/t
0: 0.0000e+00 -0.0000e+00 2e+00 1e+00 2e+00 1e+00
1: -9.9000e+01 -0.0000e+00 2e+02 1e+00 2e+00 1e+02
2: -9.9990e+03 -0.0000e+00 2e+04 1e+00 2e+00 1e+04
3: -1.0000e+06 -0.0000e+00 2e+06 1e+00 2e+00 1e+06
4: -1.0000e+08 -0.0000e+00 2e+08 1e+00 2e+00 1e+08
Certificate of dual infeasibility found.
>>> print(sol['x'])
[-3.18e-16]
[-5.89e-16]
[ 1.00e+00]
しかし、co.solvers.socp(c, Gq=g, hq=h)の実行した標準出力が前回のCertificate of primal infeasibility found.とは1単語だけ違っています。またsol['x']の中身も違っており、数値誤差を考えて$( x_1 , x_2 , x_3 ) = (0, 0, 1)$と求められているかのように見えますが、これは途中で反復解を求めることを途中でうち切っただけです。実は、このただ1つだけの不等式制約は係数は違えどもアイスクリームのコーンの形であり、目的関数値$- x_1 - x_2 - x_3$を最小化することは$x_1 + x_2 + x_3$を最大化することと同値であることを考えればすぐ分かります。アイスクリームのコーンの形の中に含まれる任意の点$( x_1 , x_2 , x_3 )$において、$x_1$と$x_2$を固定して$x_3$をどこまでも大きくしてもアイスクリームのコーンの形の中に入っています。なので、$x_1 + x_2 + x_3$の値を$x_3$のせいでどこまでも大きくすることができ、目的関数値$- x_1 - x_2 - x_3$もまた、どこまでも小さくすることができてしまいます。つまり、目的関数値が発散してしまう(=非有界、Unboundedである)わけです。
SOCPを解くときは、実行不能であったり非有界であったりすることを常に意識する必要があります。これ以外にも反復解を計算する際に、これ以上反復できずにストップするケースもあります。最適解が求められたかどうかはsol['status']で判別できます。最適解が見つかった場合は'optimal'が格納されています。
##せっかくなので、CVXOPTでSOCPの性能測定をやってみる
最後に、$n$個の変数$\mathbf{x} = ( x_1 , \ldots , x_n )$に対し、$m = 20$個のSOCPの不等式制約がある場合、$n = 5 , 10, \ldots , 35$によって計算時間がどのぐらい変化するのか数値実験してみました。数値実験に用いたマシンのスペックを以下に簡単にまとめました。
| スペック名 | 情報 |
|---|---|
| OS | Debian Buster(開発版) |
| プロセッサ | Intel(R) Core(TM) i5-5200U CPU 2.20GHz |
| Python | 3.5.3 |
| NumPy | 1.13.0 |
| CVXOPT | 1.1.9 |
数値実験は$m = 20$、$n = 5 , 10, \ldots , 35$において、異なる$\mathbf{A}_i , \mathbf{b}_i , \mathbf{c}_i , d_i , \mathbf{e}$の値を乱数で100回生成し、それぞれ実行時間を計測した平均値を求めてみます。しかし、ただ闇雲に乱数を発生しても実行不能や非有界のケースが多くなってしまいます。そこで、以下のような乱数の発生条件を設けました。
- $\mathbf{A}_i$の対角成分以外はすべて0
- $\mathbf{b}_i$は零ベクトル
- 上記以外は25%の確率で0、残75%の確率で$[-10 , 10]$の一様分布から乱数を生成
実際に使用したPythonのコードは以下の通りです。
import cvxopt as co
import numpy as np
from random import random
from time import time
import matplotlib.pyplot as plt
m = 20
n_min = 5
n_max = 35
dn = 5
limit_times = 100
all_result = []
co.solvers.options['show_progress'] = False
for n in range(n_min, n_max + dn, dn):
print()
print("------ n = " + str(n) + ", m = " + str(m) + " ------")
times = 0
result = []
while times < limit_times:
threshold = np.random.rand(n)
c = co.matrix(np.array(
[0.0 if threshold[i] > 0.75 else 20.0 * random() - 10.0 for i in range(n)]
))
threshold = [[[1.0 if i > 0 and i - 1 != j else random() for j in range(n)] for i in range(n + 1)] for _ in range(m)]
g = [co.matrix(np.array(
[[0.0 if threshold[i][j][k] > 0.75 else 20.0 * random() - 10.0 for k in range(n)] for j in range(n + 1)]
)) for i in range(m)]
threshold = np.array([[1.0 if i > 0 else random() for i in range(n + 1)] for _ in range(m)])
h = [co.matrix(np.array(
[0.0 if threshold[i][j] > 0.75 else 20.0 * random() - 10.0 for j in range(n + 1)]
)) for i in range(m)]
try:
start = time()
sol = co.solvers.socp(c, Gq=g, hq=h)
goal = time()
if sol['status'] == 'optimal':
result.append((goal - start) * 1000)
times += 1
print(str(times) + "/" + str(limit_times) + " : " + str((goal - start) * 1000) + "(ms)")
except ValueError:
print(" -> Retry : Rank(A) is lesser than the number of A's rows or Rank([g[i]; A]) is lesser than n")
all_result.append(sum(result) / len(result))
print(" -> Finish : Average = " + str(sum(result) / len(result)) + "(ms)")
fig = plt.figure(figsize=(14, 9))
plt.bar(range(n_min, n_max + dn, dn), all_result, align='center')
plt.xticks(range(n_min, n_max + dn, dn), range(n_min, n_max + dn, dn))
plt.title("$m = " + str(m) + "$", fontsize=30)
plt.xlabel("$n$", fontsize=30)
plt.ylabel("CPU Time (ms)", fontsize=30)
plt.tick_params(labelsize=25)
fig.savefig("SOCP_n" + str(n_min) + "-" + str(n_max) + "_m" + str(m) + ".png", dpi=60)
co.solvers.options['show_progress'] = Falseとは、co.solvers.socp()を実行した時にコンソール上で出力される途中の計算結果を表示させない設定を行っています。また、co.solvers.socp()によって実行される主双対内点法の中で連立1次方程式を解く過程があるのですが、$\mathbf{A}_i , \mathbf{b}_i , \mathbf{c}_i , d_i , \mathbf{e}$の乱数値によっては連立1次方程式の解が存在しないようなケースがあり得ます。この場合、co.solvers.socp()でValueErrorが発生するので、そのような乱数値は無視するようtry-exceptで囲んでいます。
結果は以下のようになりました。
$n$が大きくなるにつれ、計算時間の増加量が大きくなっているような傾向が読み取れますね。
最後に
ここまでSOCPがどういうものなのか、SOCPを解くためには何を使ってどう注意すれば良いのかの基本を書きました。どこかの参考になって頂ければ嬉しいです。
なお、本記事で使ったPythonのコードはこちらでもまとめています。