目次
- 目次
- はじめに
- 二次計画法(Quadratic Programming)とは？
- 二次計画法のPythonソルバーcvxopt
- インストール方法
- CVXOPTで簡単な二次計画問題を解いてみる。
- cvxoptにおける等式制約の注意点
- 制約の緩和のためのスラック変数の導入
- 大規模な二次計画問題のベンチマークテスト
- 参考資料
- MyEnigma Supporters
はじめに
先日、Googleが開発した
非線形最適化ソルバーCeres Solverの使い方を説明しましたが、
同じ非線形最適化でも、ある特別な形をしている最適化問題の場合、
より効率的に解ける場合があります。
その中でも有名なのが、二次計画法(Quadratic Programming)です。
今回はこの二次計画法の概要の説明と、
二次計画法を解くための
Pythonソルバーcvoptの使い方を説明したいと思います。
二次計画法(Quadratic Programming)とは？
二次計画法(Quadratic Programming:QP)は、
制約付き非線形最適化問題の代表的な問題の一つです。
二次計画法の問題は、複数個のパラメータベクトルに対して、
最適化(最小化)するコスト関数が下記のような形で、
表せるものを指します。
上記の式において、二次形式のコスト行列Pは、
正定である必要があります。(凸関数である条件のため)
下記が二次計画法を二次元平面上に図示したものです。
Pは制約条件による、状態ベクトルの範囲になります。
点線はコスト関数値の等高線です。
ちなみに、上記の二次形式のコスト行列Pが正定でない場合、
下記のように、コスト関数が楕円ではなくなり、
凸関数ではなくなってしまいます。
上記の数式の内、
一つ目の式が、最小化するコスト関数、
二つ目が制約条件を表す式になります。
コスト関数の部分が、行列の二次形式になっているため、
二次計画法と呼ばれます。
コスト関数において、
xが求めたいパラメータのベクトル(n次元)、
Qはn x nの実数値対称行列,
cはn x 1の実数ベクトルになります。
制約条件の式においては、
制約条件の数とmとすると、
Aはm x nの行列、
bはm次元のベクトルです。
制約条件を上手く定式化し、制約条件の行列式に入れることで、
制約条件付きの非線形最適化を実現することができます。
二次計画法のPythonソルバーcvxopt
二次計画法のソルバーとして一番有名なのは、
MATLABのツールボックスのOptimization Toolboxですが、
Python製の二次計画法ソルバーcvxoptというのも有名です。
CVXOPTは凸最適化用のフリーソフトウェアです。
Python製なので、インタラクティブシェルだけでなく、
Pythonスクリプトの中でも使えますし、
Pythonの拡張機能を使うことで他の言語のシステムに連結することもできます。
インストール方法
下記のリンクからコードの圧縮ファイルをダウンロードし、
解凍したディレクトリの直下で、
$sudo python setup.py install
とすればOKです。
なぜかpipだとエラーが出て、インストールできませんでした。。
CVXOPTで簡単な二次計画問題を解いてみる。
実際にこのCVXOPTを使って
簡単な二次計画問題を解いてみたいと思います。
まず、CVXOPTは下記のような二次計画の式を想定しています。
基本的には、前述した二次計画法の式と同じですが、
制約条件の不等号制約と等号制約をそれぞれ別々の
行列として考える必要があります。
CVXOPTは上記の式のP,q,G,h,A,bを与えることで、
上記の二次計画問題を解いてくれます。
Pとqは必ず与える必要がありますが、
それ以外のパラメータはオプションです。
では実際に下記の問題を解いてみたいと思います。
上記の式を、先程の二次計画法の式に合わせるように
行列化すると、下記のようになります。
ここで求めるパラメータを[x,y]としていることと、
複数の不等式を一つの行列式として合わせるために、
両辺の符号を変更していることなどに注意して下さい。
また、今回は等式制約は無いので、A,bは不要です。
では実際にCVXOPTで解いてみます。
まず初めに先程の行列式を元に、
各行列を作成します。
ちなみにCVXOPTは独自のmatrixクラスをもっており、
行列はそちらのクラスで設定する必要があります。
そのまま、matrixクラスに入力することもできますが、
今回はnumpyの行列ライブラリを使って、
行列を指定してみたいと思います。
import numpy import cvxopt from cvxopt import matrix P=matrix(numpy.diag([1.0,0.0])) q=matrix(numpy.array([3.0,4.0])) G=matrix(numpy.array([[-1.0,0.0],[0,-1.0],[-1.0,-3.0],[2.0,5.0],[3.0,4.0]])) h=matrix(numpy.array([0.0,0.0,-15.0,100.0,80.0]))
numpyのdiag関数などを使えば、
大規模な対角行列などが作りやすくなります。
あとは、実際に解くだけです。
二次計画法を解くメソッドは、
cvxopt.solvers.qpです。
先程の行列を入れるだけです。
Aとbがある場合は、Aとbをそれぞれhの後ろに入れるだけです。
sol=cvxopt.solvers.qp(P,q,G,h) print(sol) print(sol["x"]) print(sol["primal objective"])
返り値は結果の情報がpythonの辞書形式で格納されています。
xは最終的なパラメータベクトル
primal objectiveは最終的なコスト関数になります。
上記のコードを解くと、
下記のような実行結果が得られます。
8回の繰り返し計算の結果、
x=0, y=5の時に最小であることがわかりました。
その他のパラメータの説明は下記のドキュメントを参照下さい。
上記のコードはこちらでも公開しています。
cvxoptにおける等式制約の注意点
等式制約を入れる時に、
Ax=bで表されるように等式制約を入れますが、
なんとなくAをn x nの行列にしたくなります。
そして等式制約に関係の無い部分は0にしたくなりますが、
そのように設定して、qp関数を呼ぶと、
ValueError: Rank(A) < p or Rank([G; A]) < n
というエラーが出てしまいます。
しかし、cvxoptの場合は、制約の数をmとした時に、
Aをm x nの行列にし、bをm x 1のベクトルにする必要があります。
つまり、Aの中にすべての要素が0になる
行はいれてはいけないということです。
制約の緩和のためのスラック変数の導入
下記の記事で述べた通り、
二次計画法を解く際に制約条件が厳しすぎると
最適化計算ができない場合があります。
そこで、スラック変数という疑似変数を追加することで、
できるだけ制約条件を満たしつつ、
どうしても満たせない場合は、多少制約を無視する形で、
コスト関数を最小化するパラメータを得ることができます。
(このような制約をソフト制約といいます。)
∞ノルム法によるスラック変数
このスラック変数の追加の方法には２つ方法があります。
一つ目は一つのスラック変数sを追加する方法です.
(∞ノルム法とも呼ばれます)
前述の二次計画法の標準形に下記のようにスラック変数を追加します。
上記の式において、sはスラック変数(s>0の条件付きです),
ρはスラック変数のコストパラメータ,
1_vはすべて1のベクトルです。
上記の式を見ると、
ρを大きな値にしておけば、
制約条件を満たす最小値が得られる場合は、
s->0となり、通常のハード制約として機能します。
しかし、
どうしてもその制約条件の中で最適値が得られない場合、
コストが上がってしまいますが、ｓが大きくなり、
制約条件が緩和する方向になります。
これにより制約条件をできるだけ満たす形での
最適値を計算することができるのです。
1ノルム法によるスラック変数
しかし、上記の計算方法ですと、
一つの変数の制約条件が緩和されると、
他の変数の制約条件も緩和されてしまうので、
できるだけ制約条件を守りたい場合は、
もう一つの方法であるスラック変数を複数持つ方法があります。
(1ノルム法と呼ばれます)
下記のように、スラック変数をベクトルにし、
コスト関数をベクトルの和にすることで、
それぞれの個々の制約条件を緩和することができます。
この1ノルム法だと、
それぞれの制約条件を個別に緩和することができるのですが、
追加するスラック変数が、制約条件の分だけ増えてしまうので、
計算コストが増大してしまうという問題があります。
大規模な二次計画問題のベンチマークテスト
前述のように、二次計画問題は、
大規模な最適問題でも高速に大域的な最適値を計算できます。
そこで比較的大規模な二次計画問題をどれほどの計算時間で解くことができるかを、
ベンチマークテストしてみました。
使用するソルバーはC製のECOSを利用し、
解く最適問題としては、
ホライズンの長いMPC問題を使いました。
従って、比較的スパースな二次計画問題ですので、
ソルバーにとっては、高速に解きやすい問題であるといえます。
また、計算のプラットフォームとしては、
下記の少し古いMacを使っています。
CPU: 1.7 GHz Intel Core i7
メモリ: 8 GB 1600 MHz DDR3
下記が複数二次計画問題における
ベンチマークテスト結果です。
| Benchmark index | State Size | Equality constraints size | inequality constraints size | Time[s] |
|---|---|---|---|---|
| 1 | 1500 | 1000 | 3000 | 0.093 |
| 2 | 3000 | 2000 | 6000 | 0.227 |
| 3 | 9000 | 6000 | 18000 | 0.973 |
上記の結果を見ると分かる通り、
1万近い状態数で、かつ、等式制約が約5千,
不等式制約が約2万の問題でも、
1秒以下で解くことができていることがわかります。
すごいスピードですね。
参考資料
MyEnigma Supporters
もしこの記事が参考になり、
ブログをサポートしたいと思われた方は、
こちらからよろしくお願いします。
