OpenOpt, FuncDeisignerを使って2次計画法

2次計画法を扱う機会があったのでPythonで書かれているフレームワークOpenOptとモジュールFuncDesignerを使ってやってみました。

なお、OpenOptとFuncDesigner自体は

$pip install openopt funcdesigner

で簡単にインストールできます。（要NumPy）

以下、使い方まとめ、というかメモ。

2次計画法とは

2次計画法(Quadratic Problems, Quadratic Programming)とは「線形制約に従って、2次式で表される目的関数を最小化するベクトルxを見つける問題」と定義されるそうです。具体的には
\(
\begin{align*}
\min ~ &\frac{1}{2}\mathbf{x}^TH\mathbf{x} + \mathbf{f}^T x \\
\text{subject to} ~ ~
&\mathbf{lb} \le \mathbf{x} \le \mathbf{ub} \\
&A\mathbf{x} \le \mathbf{b} \\
&A_{eq}\mathbf{x} = \mathbf{b}_{eq}
\end{align*}
\)
の式を満たすようなxを見つければいいわけです。（上の式はOpenOptのページにある表記を用いています。ややこしいのですが"ub"と"lb"は、"u","l"と"b"の積を表しているわけではなく、"ub","lb"でひとつのベクトルです。）

今回はこのサイトの問題
\(
\begin{align*}
\min ~ &2x_1^2 + x_2^2 +x_1x_2 + x_1 + x_2 \\
\text{subject to} ~ ~
&x_1 \ge 0 \\
&x_2 \ge 0 \\
&x_1 + x_2 = 1 \\
\end{align*}
\)
を例題として考えます。もちろん上のサイトで紹介されてる通り、CVXOPTというパッケージだけを使ってもいいのですが、今回は様々なソルバーを使うことが出来るフレームワーク、OpenOptを用いてみます。

OpenOpt

OpenOptはいくつかのソルバーを含む数値計算最適化パッケージですが、その一番の特徴はその他のソルバーと接続することができることです。
早速ですが、このOpenOptを用いて上記の問題を解くコードを書いてみます。
問題の式を、上記で示した行列を用いた形に変形した上で、その行列やベクトルを変数として代入していきます。

…とても短いですね。なお実行結果は以下の通りです。

$ python qp_notfd.py

------------------------- OpenOpt 0.5310 ------------------------- solver: qlcp problem: unnamed type: QP iter objFunVal log10(maxResidual) 0 0.000e+00 0.00 1 0.000e+00 -100.00 istop: 1000 Solver: Time Elapsed = 0.02 CPU Time Elapsed = 0.0 objFunValue: 1.875 (feasible, MaxResidual = 0) x = [ 0.25 0.75]

今回はソルバーはOpenOptに標準で付属している"qlcp"を用いました。
（他にも2次計画法で使えるソルバーがこのページにまとめてあります。）
OpenOptのサイトの例では行列をNumPyの関数を用いて表していますが、ネストしたリストでも問題なく動いてくれます。なお、コード中の"r.xf"の型はnumpy.ndarrayです。

FuncDesigner

さて、今回のような問題をOpenOptを用いて解く際に面倒なのは、式を行列を用いた形に変形することだと思います。このような煩わしさを解決するためにFuncDesignerという便利なモジュールが存在します。これを用いればOpenOptをより直観的に使えるようになります。
FuncDesignerを用いてコードを書き換えたのが以下です。

変数を宣言して、目的関数、制約条件を書くだけ。とても簡単です。
気をつけないといけない点としては、初期点を指定してあげないといけないことぐらいです。
たぶんこのレベルの2次計画法なら初期点はなんでもOKだと思います。（もしも{x1:0, x2:0}でダメだったら{x1:1, x2:1}にするかー、ぐらいの感覚で使ってます…。ちょっと初期点についてはよく分かりません。）

さて、場合によっては変数がとても多くなる場合もあります。
そんな時は変数がベクトルで宣言できた方が便利ですし、やっぱり行列を用いた式で書きたくなると思います。上のコードと同様にFuncDesignerを用いながら、変数をベクトルにして書き換えます。

ここでは行列・ベクトルの表現にNumPyの関数を使いましたが、たぶん使わななくても大丈夫です。
また、fd.dot(x,y)は行列x,yの積を表し、fd.sum(x)は行列xの要素の総和を表します。共にFuncDesignerの関数です。FuncDesignerには他にも様々な数式操作のための関数があるので、必要に応じて使いたいものです。

ちなみに、FuncDesignerではn次元ベクトルの変数はfd.oovar(size=n)で宣言するのが一般的ですが、残念ながら2次計画法ではfd.oovars(n)で宣言しないと動かないようです。ただの線形計画法ならfd.oovar(size=n)でもOKなのですが…。

CVXOPT

問題をCVXOPTのページから取っておきながら、言及しないのはどうかと思うので一応説明すると、CVXOPTは凸最適化パッケージです。このパッケージだけで最適化問題を解くことも可能で、ソルバーも含んでいます。もちろんこのソルバーはOpenOptから呼び出すことも可能で、2次計画法なら今まで紹介したコード中のp.solve("qlcp")をp.solve("cvxopt_qp")に変えれば使うことができます。

参考サイト

• QP - OpenOpt
• FuncDesignerDoc - OpenOpt
• Solving a quadratic program — CVXOPT
• satomacoto: Pythonで線形計画法を解いてみる
• OpenOpt の線形計画で圧縮センシング
• ソフトマージンSVMを書いてみた - 唯物是真 @Scaled_Wurm

以上、OpenOptを使った2次計画法についてでした。
なお、今回始めてGistとMathJaxを使ってみました。うまく表示されてるといいのですが。