二次計画法を久しぶりにやる - 驚異のアニヲタ社会復帰の予備

二次計画法をRでやる必要が出てきたので久しぶりにやってみる。
二次計画法はRはでは

ただし、行列 $D$ が正定値行列でないと

matrix D in quadratic function is not positive definite!

と言われてしまうので、逆行列を持つように微妙に変化させた行列を無理やり作る。

R のときの最適化は $min(-d^T b + \frac{1}{2} b^T D b)$ with the constraints $A^T b \geq b_0$ を行うので、それに合わせて引数となる行列とベクトルを書き換える。
meq の引数は上からmeq 個の条件式が等号の場合であり、すべてが不等式ならmeq=0 とする。

例題　http://courses.csail.mit.edu/6.867/wiki/images/a/a7/Qp-cvxopt.pdf

$\min(\frac{1}{2}x^2+3x+4y$ は
$\min\frac{1}{2}\begin{bmatrix}x\\y\end{bmatrix}^T \begin{bmatrix}1&0\\0&0\end{bmatrix}\begin{bmatrix}x\\y\end{bmatrix}+\begin{bmatrix}3\\4\end{bmatrix}^T \begin{bmatrix}x\\y\end{bmatrix}$
と書けるのでこれに合うように書く。

library(quadprog)
library(Matrix)
min(-d^T b + 1/2 b^T D b) with the constraints A^T b >= b_0

Amat <- rbind(c(1, 0, 1, -2, -3), c(0, 1, 3, -5, -4))
bvec <- c(0, 0, 15, -100, -80)
dvec <- -c(3, 4)
Dmat <- matrix(c(1, 0, 0, 0), 2)
Dmat <- nearPD(Dmat)$mat
solve.QP(Dmat, dvec, Amat, bvec=bvec, meq=0)$solution

$solution
[1] 0 5

from cvxopt import matrix
from cvxopt import solvers
P = matrix([[1.0,0.0],[0.0,0.0]])
q = matrix([3.0,4.0])
G = matrix([[-1.0,0.0,-1.0,2.0,3.0],[0.0,-1.0,-3.0,5.0,4.0]])
h = matrix([0.0,0.0,-15.0,100.0,80.0])
sol = solvers.qp(P,q,G,h)
list(sol["x"])

[7.131816408856328e-07, 5.000001008391809]