posted with カエレバ

北村充森北出版 2015-02-27

Amazonで探す

楽天市場で探す

Yahooショッピングで探す

はじめに

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

myenigma.hatenablog.com

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

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

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

myenigma.hatenablog.com

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

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

そこで、まず手始めに、

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

(不等式制約は無し)

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

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

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

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

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

今回解きたい問題は、

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

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

f:id:meison_amsl:20170113152229p:plain

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

Pは半正定(semi-difinite)

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

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

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

下記の式に変換します。

$L(x,\lam)=\frac{1}{2}x^TPx+q^Tx-\lam(Ax+b)$

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

myenigma.hatenablog.com

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

KKT条件より、

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

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

$\frac{\partial L}{\partial x}=x*-A^T\lam*+q^T=0$

$\frac{\partial L}{\partial \lam}=-A\lam*+b=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