[PYTHON] Quadratic programming by the interior point method
Quadratic programming problem
Last time solved the linear programming problem using the interior point method (Karmarkar's algorithm). This time, I will try to solve the quadratic programming problem using the interior point method in the same way. The quadratic programming problem has the objective function in quadratic form and the constraints in linear form, and is expressed as a problem that minimizes the vector $ x $ as shown below.
\min \frac{1}{2} {\bf x}^T {\bf D} {\bf x} + {\bf c}^T {\bf x}\\
subject \ \ \ {\bf A}{\bf x} \ge {\bf b}
Here, if the dimension of $ {\ bf x} $ is $ n $ and the number of constraints is $ m $, then $ {\ bf x}, {\ bf c}, {\ bf b} $ are $ n $ dimension vectors. So, $ {\ bf A} $ is a matrix of $ n \ times m $, and $ {\ bf D} $ is a matrix of positive values of $ n \ times n $.
Interior-point method in quadratic programming
This time, it was introduced in Numerical Optimization as a method to solve the quadratic programming problem by the interior point method. I tried to implement the predictor-correction method. It is more complex than the linear programming problem.
#!/usr/bin/python
#coding:utf-8
import numpy as np
EPS_ZERO = 1.0e-9
def select_alpha(alphas):
return 1.0 if alphas.size == 0 or np.min(alphas) > 1.0 else np.min(alphas)
def inv_vec(vec):
inv = np.matrix(zeros((vec.shape[0], 1)))
cond = np.where(np.abs(vec) > EPS_ZERO)
inv[cond] = 1.0 / vec[cond]
return inv
def less_zero(vec):
return vec <= -EPS_ZERO
def quadratic_programming(x, gmat, c, amat, b,
tau=0.5, eps=1.0e-3, nloop=30):
"""
Interior point method(Predictor-Collector)Solving the quadratic programming problem by
object min z = 0.5 * x^T * G * x + cT * x
subject Ax >= b
"""
ndim = gmat.shape[0] #Number of dimensions
ncnst = amat.shape[0] #Number of constraints
x_idx, y_idx, l_idx = np.split(np.arange(0, ndim + 2 * ncnst),
[ndim, ndim + ncnst])
y_l_idx = np.r_[y_idx, l_idx]
zeros_d_c = np.zeros((ndim, ncnst))
zeros_c2 = np.zeros((ncnst, ncnst))
eye_c = np.identity(ncnst)
f = np.bmat([[-gmat * x - c],
[zeros((ncnst, 1))],
[amat * x - b]])
exmat = np.bmat([[gmat, zeros_d_c, -amat.T],
[zeros_d_c.T, zeros_c2, eye_c],
[-amat, eye_c, zeros_c2]])
delta_xyl0 = np.linalg.pinv(exmat) * f
abs_yl0 = np.abs(delta_xyl0[y_l_idx])
abs_yl0[abs_yl0 < 1.0] = 1.0
y, lmd = np.vsplit(abs_yl0, 2)
rd = gmat * x - amat.T * lmd + c
rp = amat * x - y - b
for n in range(nloop):
err = np.linalg.norm(rd)**2 + np.linalg.norm(rp)**2 + np.asscalar(y.T * lmd)
if err < eps:
break
# affine scaling step
iy = inv_vec(y)
f = np.bmat([[-rd],
[-lmd],
[rp]])
exmat = np.bmat([[gmat, zeros_d_c, -amat.T],
[zeros_d_c.T, np.diag(np.multiply(lmd, iy).A1), eye_c],
[-amat, eye_c, zeros_c2]])
d_xyl_aff = np.linalg.pinv(exmat) * f
#calculation of alpha
y_l_cnb = np.bmat([[y],
[lmd]])
cnd = less_zero(d_xyl_aff[y_l_idx])
alpha_aff = select_alpha(-y_l_cnb[cnd] / d_xyl_aff[y_l_idx][cnd])
#Central parameters
mu = np.asscalar(y.T * lmd / ncnst)
mu_aff = np.asscalar((y + alpha_aff * d_xyl_aff[y_idx]).T * (lmd + alpha_aff * d_xyl_aff[l_idx]) / ncnst)
sig = (mu_aff / mu)**3
# corrector step
dl_y_aff = np.multiply(np.multiply(d_xyl_aff[y_idx],
d_xyl_aff[l_idx]),
iy)
f = np.bmat([[-rd],
[-lmd - dl_y_aff + sig * mu * iy],
[rp]])
d_xyl = np.linalg.pinv(exmat) * f
# alpha_pri, alpha_dual calculation
cnd = less_zero(d_xyl)
alpha_pri = select_alpha(-tau * y[cnd[y_idx]] / d_xyl[y_idx][cnd[y_idx]])
alpha_dual = select_alpha(-tau * lmd[cnd[l_idx]] / d_xyl[l_idx][cnd[l_idx]])
alpha_hat = min([alpha_pri, alpha_dual])
# x, y,lmd update
x, y, lmd = np.vsplit(np.r_[x, y, lmd] + alpha_hat * d_xyl,
[ndim, ndim + ncnst])
rp = (1.0 - alpha_pri) * rp
rd = (1.0 - alpha_dual) * rd + (alpha_pri - alpha_dual) * gmat * d_xyl[x_idx]
return x, y, lmd
if __name__ == "__main__":
c = np.matrix([[-2.0],
[-4.0]])
dmat = np.matrix([[2.0, 0.0],
[0.0, 2.0]])
amat = np.matrix([[1.0, 1.0],
[1.0, -1.0],
[-3.0, 1.0]])
b = np.matrix([[-0.5],
[1.0],
[-3.0]])
#drawing
from pylab import *
ax = subplot(111, aspect='equal')
x = np.arange(-3.0, 3.01, 0.15)
y = np.arange(-3.0, 3.01, 0.15)
X,Y = meshgrid(x, y)
t = arange(-3.0, 3.01, 0.15)
func = lambda x, y : 0.5 * (dmat[0, 0] * x**2 + dmat[1, 1] * y**2) + c[0, 0] * x + c[1, 0] * y
const = [lambda x : -amat[i, 0] / amat[i, 1] * x + b[i, 0] / amat[i, 1] for i in range(amat.shape[0])]
Z = func(X, Y)
s = [const[i](t)foriinrange(amat.shape[0])]
pcolor(X, Y, Z)
for i in range(amat.shape[0]):
ax.plot(t, s[i], 'k')
x, y, lmd = quadratic_programming(np.matrix([[0.0], [0.0]]),
dmat, c, amat, b,
tau=0.5, eps=1.0e-3, nloop=25)
print x
ax.plot([x[0, 0]], [x[1, 0]],'yo')
axis([-3, 3, -3, 3])
show()
Like last time, I will try drawing with a two-dimensional problem. The area of the solution that satisfies the constraint is within the three lines, and the direction in which blue becomes darker is the optimum value. The yellow dots in the graph indicate the final optimal solution, and it can be seen that the appropriate solution can be obtained even if the initial value is given from outside the constraint.
Since the code created this time is Python, the calculation speed is slow, but I thought it would be a little easier to understand the algorithm. In this code, we used the pseudo inverse matrix when finding the steps, but it seems that it can be speeded up by using the conjugate gradient method or by using the sparseness of ʻexmat`.
Recommended Posts