目次
はじめに
先日、凸最適化の概要を紹介しましたが、
一般的にこの凸最適化を解くためには、
既存のソルバーと呼ばれるライブラリを使うことが一般的です。
Pythonでは以前紹介したcvxoptというソルバーが有名です。
しかし、別の言語で最適化システムを実装したいときなど、
このソルバーを使わずに凸最適問題を解きたい時もあります。
そこで、まず手始めに、
等式制約条件のみの凸二次計画法
(不等式制約は無し)
の解法を数学的に導出し、
それを元にPythonでソルバーを使わずに
凸二次最適化を解くサンプルコードを
作成したので紹介したいと思います。
等式制約下の凸二次計画法の数学的解法
今回解きたい問題は、
下記のような等式制約を持つ、
凸二次計画問題であるとします。
ここで凸性を担保するために、
Pは半正定(semi-difinite)
つまり, x'Pxは>=0であるとします。
まず、上記の問題を解くために
上記の問題をラグランジュの未定乗数法を使って
下記の式に変換します。
ラグランジュの未定乗数法に関しては下記を参照下さい。
続いて、上記の式において、
KKT条件より、
上記の式における最適値x*は、
上記の式をそれぞれxとλで偏微分した値が0になります。
上記の二式より、
下記のような行列式を作ることができます。
上記の式を変形すると、
となるので、
xとλの最適値は下記の行列計算から計算できます。
ここで、逆行列をとっている行列をKKT行列と呼ぶようです。
Pythonサンプルコード
上記の式から、
下記のようなPythonコードで
等式制約下の凸二次計画法を解くことができます。
下記のコードのsolve_qp_with_ep_const関数は、
上記の凸二次計画法のパラメータ行列であるP,q,A,bを与えると、
最適値のxベクトルを計算してくれます。
def solve_qp_with_ep_const(P, q, A, b): """ solve quadratic programming with only equality constraints min 0.5*x*P*x + q.T*x s.t Ax = b """ # input check if not isinstance(P, np.matrix): raise TypeError("'P' must be a np.matrix") if not isinstance(q, np.matrix): raise TypeError("'q' must be a np.matrix") if not isinstance(A, np.matrix): raise TypeError("'A' must be a np.matrix") if not isinstance(b, np.matrix): raise TypeError("'b' must be a np.matrix") if P.shape[0] != P.shape[1]: raise ValueError("'P' must be a square matrix") if P.shape[1] != q.shape[1]: raise ValueError("'P' or 'q' is invalid matrix size") if A.shape[0] != b.shape[1]: raise ValueError("'A' or 'b' is invalid matrix size") K1 = np.concatenate((P, A.T), axis=1) K2 = np.concatenate((A, np.zeros((A.shape[0], A.shape[0]))), axis=1) K = np.concatenate((K1, K2), axis=0) d = np.concatenate((-q.T, b), axis=0) star = np.linalg.solve(K, d) x_star = star[0:A.shape[1], :] return x_star def test_solve_qp_with_ep_const(): P = np.matrix(np.diag([1.0, 0.0])) q = np.matrix(np.array([3.0, 4.0])) A = np.matrix([1.0, 1.0]) b = np.matrix(1.0) x = solve_qp_with_ep_const(P, q, A, b) assert x[0] - 1.0 < 0.0001 assert x[1] - 0.0 < 0.0001
上記のコードは下記でも公開しています。
参考資料
MyEnigma Supporters
もしこの記事が参考になり、
ブログをサポートしたいと思われた方は、
こちらからよろしくお願いします。