本記事ではSDPを解くためのpythonライブラリであるPICOSとCVXPYのSDP実行時間を比較します．
SDP（半正定値計画問題）
**半正定値計画問題(SDP; SemiDefinite Programming)**は凸最適化問題の一種です．
数理計画やシステム制御の分野で注目を浴びており，現在盛んに研究されています．
SDPは例えば以下のような問題です．
正方行列$A$が与えられたとする．
$P\succ0$, $PA-A^TP\prec0$
を満たす行列Pを求めよ．
（$P \succ 0$は $P$ が正定行列であることを表す．）
SDPでは制約条件が**線形行列不等式（LMI;Linear Matrix Inequality）**で与えられます．
この問題に解Pが存在することとAの固有値の実部が全て負であることは同値です．
よって，この問題が解ければAの固有値が全て負であることを確認することができます．
この例だとAの固有値を直接計算した方が早いと言われてしまいますが，
複数のLMI制約が組み合わさったり，目的関数を最適化したい場合にはSDPとして解く必要があります．
システム制御の分野では制御仕様を複数のLMI制約で表現できる場合があり，
SDPを解けることは非常に有用です．
SDPを数値的に解くツールとしてMATLABやPythonが利用できます．
PythonのSDPライブラリ
PythonでSDPを解くためのライブラリは私の知る限りPICOS，CVXPY，PYLMI-SDPの3つです．
そのうちPYLMI-SDPは使用例やドキュメントを発見できなかったため扱いません．
PICOSやCVXPYは記述されたSDPをソルバーに渡すインターフェースの役割をしています．
ソルバーは実際に問題を解く作業を請け負っている部分です．
PICOSとCVXPYでサポートされているSDPソルバーを以下に示します．
表：PICOS, CVXPYが対応するSDPソルバー
| CVXOPT | MOSEK | SMCP | SCS | SDPA-P | |
|---|---|---|---|---|---|
PICOS |
○ | ○ | ○ | - | ○ |
CVXPY |
○ | ○ | - | ○ | - |
各ソルバーの説明
-
CVXOPT・・・
PICOS,CVXPYと一緒にインストールされる． - MOSEK・・・商用．アカデミックライセンスは無料．
- SMCP・・・この論文に基づいて作られた．
-
SCS・・・この論文に基づいて作られた．
CVXPYと一緒にインストールされる． - SDPA-P・・・多分この論文に基づいて作られた．
上記の通り，たくさんのライブラリおよびソルバーがあります．
多くの人が開発し充実しているのは良いのですが，どれを使うべきなのか分かりません．
そこで，これらを使ったSDPの実行時間を比較してみようと思います．
今回は上の表のうちSMCPとSDPA-Pを除いた以下の5つの組み合わせを比較します．
PICOS×CVXOPT, PICOS×MOSEK
CVXPY×CVXOPT, CVXPY×MOSEK, CVXPY×SCS
*他のソルバーも今後掲載予定
手順
- 定数行列 $A$ を
numpy.randomで生成． - 以下のSDPを解く．
$ ~~~~~~~~~~~~~~~~~ \rm{min}~~\rm{trace}(\it{P}) ~~~ \rm{s.t.} ~~ \it{P}\succ \rm{0},~\it{PA-A^TP}\prec \rm{0}.$ -
%timeitで平均実行時間を計測． - 以下のパラメータ毎に計測する．
・定数行列 $A$ および変数行列 $P$ の次数 $n$
$ ~~~~~~~~~~~~~~~~~ 1\leq n \leq20$
・許容誤差 $tol$
$ ~~~~~~~~~~~~~~~~~ tol={ 10^{-4},10^{-6},10^{-8},10^{-10},10^{-12} }$
実行コード
実行環境
| ライブラリ | バージョン |
|---|---|
| anaconda | 2019.03 |
| python | 3.7 |
| Jupyter Lab | |
| picos | 2.0.0 |
| cvxpy | 1.0.28 |
| cvxopt | 1.2.0 |
| mosek | 9.1.13 |
| scs | 2.1.1.2 |
結果
横軸は行列のサイズで，縦軸は実行時間です．
許容誤差は $tol=1e-08$ で固定しています．
一番速いのはPICOS×___MOSEK___です．
次点はPICOS×CVXOPTです．
それ以外はほぼ同じです．
横軸は許容誤差で，縦軸は実行時間です．
行列のサイズは $n=2$ で固定しています．
一番速いのはPICOS×___MOSEK___です．
次点はPICOS×CVXOPTです．
それ以外はほぼ同じです．
ただし，PICOS×MOSEKは $10^{-10}$ 以上の許容誤差では小さすぎる旨のエラーが出て実行できませんでした．
#まとめ
・変数や制約が大量にあったり，大量のSDPを解く場合には，PICOS×MOSEKがおすすめ．
・MOSEKの商用ライセンスまたは無料アカデミックライセンスを利用できない方はPICOS×CVXOPTを使いましょう．
※各ソルバーは異なるアルゴリズムを採用しているので，問題の性質によって結果が変わる可能性があります．
備忘録
・PICOSは最近ver2.0にアップデートされてバグあり(2020/3/22現在)．
→ライブラリを書き換えてしのぐ．
（searchTime の Zero Division Error ： 該当行をコメントアウトし, 適当な値(overheadTime=1)を入れる．）
・PICOS×MOSEKではobjectiveを'find'に設定するとエラーが出る．'min'で適当な関数を入れて解く．
・CVXPYで謎のverboseが発生するバグが発生，anacondaをアップデートしたら解決．
参考にした情報
・PICOSドキュメント
・CVXPYドキュメント
・PICOSメモ (半正定値計画問題)
・PICOSを使って錐計画問題を解く
・Semidefinite programming in Python
・LMIによるシステム制御 - ロバスト制御系設計のための体系的アプローチ
・行列不等式アプローチによる制御系設計 (システム制御工学シリーズ)