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.

ipqp.png

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

Quadratic programming by the interior point method
Saddle point search using the gradient method
The first Markov chain Monte Carlo method by PyStan
I tried increasing or decreasing the number by programming
The copy method of pandas.DataFrame is deep copy by default
[Pyro] Statistical modeling by the stochastic programming language Pyro ① ~ What is Pyro ~