読者です 読者をやめる 読者になる 読者になる

MyEnigma

とあるエンジニアのブログです。#Robotics #Programing #C++ #Python #MATLAB #Vim #Mathematics #Book #Movie #Traveling #Mac #iPhone

等式制約下の凸二次計画法をソルバーを使わずにPythonで解いてみた

あたらしい数理最適化: Python言語とGurobiで解く

あたらしい数理最適化: Python言語とGurobiで解く

  • 作者: 久保幹雄,ペドロソジョア・ペドロ,村松正和,レイスアブドル
  • 出版社/メーカー: 近代科学社
  • 発売日: 2016/09/06
  • メディア: Kindle版
  • この商品を含むブログを見る

目次

はじめに

先日、凸最適化の概要を紹介しましたが、

myenigma.hatenablog.com

一般的にこの凸最適化を解くためには、

既存のソルバーと呼ばれるライブラリを使うことが一般的です。

 

Pythonでは以前紹介したcvxoptというソルバーが有名です。

myenigma.hatenablog.com

 

しかし、別の言語で最適化システムを実装したいときなど、

このソルバーを使わずに凸最適問題を解きたい時もあります。

 

そこで、まず手始めに、

等式制約条件のみの凸二次計画法

(不等式制約は無し)

の解法を数学的に導出し、

それを元にPythonでソルバーを使わずに

凸二次最適化を解くサンプルコードを

作成したので紹介したいと思います。

 

等式制約下の凸二次計画法の数学的解法

今回解きたい問題は、

下記のような等式制約を持つ、

凸二次計画問題であるとします。

f:id:meison_amsl:20170113152229p:plain

ここで凸性を担保するために、

Pは半正定(semi-difinite)

つまり, x'Pxは>=0であるとします。                                                            

まず、上記の問題を解くために

上記の問題をラグランジュの未定乗数法を使って

下記の式に変換します。

ラグランジュの未定乗数法に関しては下記を参照下さい。

myenigma.hatenablog.com

 

続いて、上記の式において、

KKT条件より、

上記の式における最適値x*は、

上記の式をそれぞれxとλで偏微分した値が0になります。

 

上記の二式より、

下記のような行列式を作ることができます。

f:id:meison_amsl:20170114071900p:plain

上記の式を変形すると、

f:id:meison_amsl:20170114072053p:plain

となるので、

xとλの最適値は下記の行列計算から計算できます。

f:id:meison_amsl:20170114072249p:plain

ここで、逆行列をとっている行列を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

 

上記のコードは下記でも公開しています。

github.com

 

参考資料

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

myenigma.hatenablog.com

あたらしい数理最適化: Python言語とGurobiで解く

あたらしい数理最適化: Python言語とGurobiで解く

  • 作者: 久保幹雄,ペドロソジョア・ペドロ,村松正和,レイスアブドル
  • 出版社/メーカー: 近代科学社
  • 発売日: 2016/09/06
  • メディア: Kindle版
  • この商品を含むブログを見る