Quadratic Program (QP) Solver that only depends on NumPy/SciPy?

Question

I would like students to solve a quadratic program in an assignment without them having to install extra software like cvxopt etc. Is there a python implementation available that only depends on NumPy/SciPy?

Halutz

If you could provide some links on what you mean by a quadratic program and maybe an example or two, it would allow more people to answer this question. Please update your question, because I am not too sure what you mean by QP and I might know how to write your program, although I don't know what it requires. Thank you! – Shingly

Sorry for not clarifying. QP is a special linear algebra problem, see Wikipedia (en.wikipedia.org/wiki/Quadratic_programming). – Halutz

I find it odd that a question asking for a python implemented QP solver that only depends on numpy/scipy and doesn't require additional software like cvxopt… has one answer that recommends cvxopt and another (the accepted answer) that recommends what's essentially unmaintained python bindings to another language (i.e. a non-python implementation). – Sickle

Answer 1

I ran across a good solution and wanted to get it out there. There is a python implementation of LOQO in the ELEFANT machine learning toolkit out of NICTA (http://elefant.forge.nicta.com.au as of this posting). Have a look at optimization.intpointsolver. This was coded by Alex Smola, and I've used a C-version of the same code with great success.

Richmound

I don't believe the project is active. The download link is broken, but this link works: elefant.forge.nicta.com.au/download/release/0.4/index.html There's a C++ fork of the project at users.cecs.anu.edu.au/~chteo/BMRM.html, but I don't believe it is active either. – Richmound

The links in this answer are broken and the suggested software is not pure Python (+Numpy/Scipy) – Isosteric

Answer 2

I'm not very familiar with quadratic programming, but I think you can solve this sort of problem just using scipy.optimize's constrained minimization algorithms. Here's an example:

import numpy as np
from scipy import optimize
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D

# minimize
#     F = x[1]^2 + 4x[2]^2 -32x[2] + 64

# subject to:
#      x[1] + x[2] <= 7
#     -x[1] + 2x[2] <= 4
#      x[1] >= 0
#      x[2] >= 0
#      x[2] <= 4

# in matrix notation:
#     F = (1/2)*x.T*H*x + c*x + c0

# subject to:
#     Ax <= b

# where:
#     H = [[2, 0],
#          [0, 8]]

#     c = [0, -32]

#     c0 = 64

#     A = [[ 1, 1],
#          [-1, 2],
#          [-1, 0],
#          [0, -1],
#          [0,  1]]

#     b = [7,4,0,0,4]

H = np.array([[2., 0.],
              [0., 8.]])

c = np.array([0, -32])

c0 = 64

A = np.array([[ 1., 1.],
              [-1., 2.],
              [-1., 0.],
              [0., -1.],
              [0.,  1.]])

b = np.array([7., 4., 0., 0., 4.])

x0 = np.random.randn(2)

def loss(x, sign=1.):
    return sign * (0.5 * np.dot(x.T, np.dot(H, x))+ np.dot(c, x) + c0)

def jac(x, sign=1.):
    return sign * (np.dot(x.T, H) + c)

cons = {'type':'ineq',
        'fun':lambda x: b - np.dot(A,x),
        'jac':lambda x: -A}

opt = {'disp':False}

def solve():

    res_cons = optimize.minimize(loss, x0, jac=jac,constraints=cons,
                                 method='SLSQP', options=opt)

    res_uncons = optimize.minimize(loss, x0, jac=jac, method='SLSQP',
                                   options=opt)

    print '\nConstrained:'
    print res_cons

    print '\nUnconstrained:'
    print res_uncons

    x1, x2 = res_cons['x']
    f = res_cons['fun']

    x1_unc, x2_unc = res_uncons['x']
    f_unc = res_uncons['fun']

    # plotting
    xgrid = np.mgrid[-2:4:0.1, 1.5:5.5:0.1]
    xvec = xgrid.reshape(2, -1).T
    F = np.vstack([loss(xi) for xi in xvec]).reshape(xgrid.shape[1:])

    ax = plt.axes(projection='3d')
    ax.hold(True)
    ax.plot_surface(xgrid[0], xgrid[1], F, rstride=1, cstride=1,
                    cmap=plt.cm.jet, shade=True, alpha=0.9, linewidth=0)
    ax.plot3D([x1], [x2], [f], 'og', mec='w', label='Constrained minimum')
    ax.plot3D([x1_unc], [x2_unc], [f_unc], 'oy', mec='w',
              label='Unconstrained minimum')
    ax.legend(fancybox=True, numpoints=1)
    ax.set_xlabel('x1')
    ax.set_ylabel('x2')
    ax.set_zlabel('F')

Output:

Constrained:
  status: 0
 success: True
    njev: 4
    nfev: 4
     fun: 7.9999999999997584
       x: array([ 2.,  3.])
 message: 'Optimization terminated successfully.'
     jac: array([ 4., -8.,  0.])
     nit: 4

Unconstrained:
  status: 0
 success: True
    njev: 3
    nfev: 5
     fun: 0.0
       x: array([ -2.66453526e-15,   4.00000000e+00])
 message: 'Optimization terminated successfully.'
     jac: array([ -5.32907052e-15,  -3.55271368e-15,   0.00000000e+00])
     nit: 3

Majormajordomo

I doubt that this is very efficient. I think an implementation of LOQO: An Interior Point Code for Quadratic Programming (citeseer.ist.psu.edu/viewdoc/summary?doi=10.1.1.39.2191) will be faster. – Halutz

How hard are the problems you need your students to solve? SLSQP solves my (admittedly rather simple) example in about 1.33msec. It can also handle any combination of bounds, inequality and equality constraints. If your heart is set upon using a particular solver that is optimised for QP then you will probably have to (A) have your students install extra dependencies, or (B) write it yourself. – Majormajordomo

Thanks for your follow up. The students should use it to solve an Support Vector Machine problem to compare it to a more efficient algorithm they should implement. It's a convex problem in about 100 variables. I might implement the LOQO, just thought I can't be the first. – Halutz

It's worth adding 'jac':(lambda x:-A) to the constraint definition, to make the solver more robust. – Knighthood

I was trying to implement some basic machine learning algorithms from scratch. SVM was on the todo list but I had no confident to pull it out. After reading your answer, I managed to write a svm of my own (github.com/Sacry/mla_sani/blob/master/mla_sani/supervised/…) and it works pretty as expected. I'm really really appreciated for your answer, thank you very much. – Uncertainty

Some useful Python documentation: scipy.optimize tutorial/docs (docs.scipy.org/doc/scipy/reference/tutorial/optimize.html, docs.scipy.org/doc/scipy/reference/optimize.html), scipy-lectures (scipy-lectures.org/advanced/mathematical_optimization/…). – Chicanery

Answer 3

This might be a late answer, but I found CVXOPT - http://cvxopt.org/ - as the commonly used free python library for Quadratic Programming. However, it is not easy to install, as it requires the installation of other dependencies.

Nedra

Well, as you described, it's not easy to install :-) Upvote as my thanks for the suggestion but I think I'll try another options first. – Taoism

@JimRaynor I have no problem installing cvxopt directly with pip install cvxopt in OS X. That's it. pip takes care of everything. And I have installed cvxopt in several machines already. Surely you need to have compilers installed, but that's also straightforward and if you are using scipy you most likely have them already. In case it helps, I use Anaconda as a Python distribution (which is fully free) and installing Anaconda is also straightforward. You don't need admin privileges and there isn't anything you need to config. Just download it, install it, and it's ready to go. – Boatyard

This library was one of the reasons I switched to Anaconda for the ease of managing the dependencies. I just couldn't install it with pip. If you already have Anaconda, use conda install -c https://conda.anaconda.org/omnia cvxopt and it's done. I'm on Windows 10 and Python 2.7. – Vulgus

Note that the question explicitly reqests a solver that does not require the installation of cvxopt – Isosteric

Answer 4

The qpsolvers package also seems to fit the bill. It only depends on NumPy and can be installed by pip install qpsolvers. Then, you can do:

from numpy import array, dot
from qpsolvers import solve_qp

M = array([[1., 2., 0.], [-8., 3., 2.], [0., 1., 1.]])
P = dot(M.T, M)  # quick way to build a symmetric matrix
q = dot(array([3., 2., 3.]), M).reshape((3,))
G = array([[1., 2., 1.], [2., 0., 1.], [-1., 2., -1.]])
h = array([3., 2., -2.]).reshape((3,))

# min. 1/2 x^T P x + q^T x with G x <= h
print "QP solution:", solve_qp(P, q, G, h)

You can also try different QP solvers (such as CVXOPT mentioned by Curious) by changing the solver keyword argument, for example solver='cvxopt' or solver='osqp'.

Nervous

qpsolvers is just a wrapper to other quadratic programming packages (like cvxopt), it requires a compiler for installation and is not in pure Python (+Numpy/Scypy) as requested – Isosteric

Answer 5

mystic provides a pure python implementation of nonlinear/non-convex optimization algorithms with advanced constraints functionality that typically is only found in QP solvers. mystic actually provides more robust constraints than most QP solvers. However, if you are looking for optimization algorithmic speed, then the following is not for you. mystic is not slow, but it's pure python as opposed to python bindings to C. If you are looking for flexibility and QP constraints functionality in a nonlinear solver, then you might be interested.

"""
Maximize: f = 2*x[0]*x[1] + 2*x[0] - x[0]**2 - 2*x[1]**2

Subject to: -2*x[0] + 2*x[1] <= -2
             2*x[0] - 4*x[1] <= 0
               x[0]**3 -x[1] == 0

where: 0 <= x[0] <= inf
       1 <= x[1] <= inf
"""
import numpy as np
import mystic.symbolic as ms
import mystic.solvers as my
import mystic.math as mm

# generate constraints and penalty for a nonlinear system of equations 
ieqn = '''
   -2*x0 + 2*x1 <= -2
    2*x0 - 4*x1 <= 0'''
eqn = '''
     x0**3 - x1 == 0'''
cons = ms.generate_constraint(ms.generate_solvers(ms.simplify(eqn,target='x1')))
pens = ms.generate_penalty(ms.generate_conditions(ieqn), k=1e3)
bounds = [(0., None), (1., None)]

# get the objective
def objective(x, sign=1):
  x = np.asarray(x)
  return sign * (2*x[0]*x[1] + 2*x[0] - x[0]**2 - 2*x[1]**2)

# solve    
x0 = np.random.rand(2)
sol = my.fmin_powell(objective, x0, constraint=cons, penalty=pens, disp=True,
                     bounds=bounds, gtol=3, ftol=1e-6, full_output=True,
                     args=(-1,))

print 'x* = %s; f(x*) = %s' % (sol[0], -sol[1])

Things to note is that mystic can generically apply LP, QP, and higher order equality and inequality constraints to any given optimizer, not just a special QP solver. Secondly, mystic can digest symbolic math, so the ease of defining/entering the constraints is a bit nicer than working with the matrices and derivatives of functions. mystic depends on numpy, and will use scipy if it is installed (however, scipy is not required). mystic utilizes sympy to handle symbolic constraints, but it's also not required for optimization in general.

Output:

Optimization terminated successfully.
         Current function value: -2.000000
         Iterations: 3
         Function evaluations: 103

x* = [ 2.  1.]; f(x*) = 2.0

Get mystic here: https://github.com/uqfoundation

Sickle

The suggested solution does not use a quadratic programming solver, but a nonlinear one. If a generic nonlinear solver is enough, a better answer is that by @Majormajordomo which only depends on Numpy/Scipy – Isosteric

@divenex: the OP didn't ask for a QP solver, they asked for something to solve a QP problem that only depended on numpy/scipy... and well, the solvers in mystic essentially only have a numpy dependency (note no scipy dependency!). You can see by the vote count that @Majormajordomo has a more direct solution (i.e. using scipy) -- and I knew that. My point is, mystic can solve QP problems as well as nonlinear problems... so, it's worth mentioning. – Sickle

My comment is meant to save time to other users coming to this page like me looking for a QP solver, as stated in the question title. mystic does not include such a QP solver. Actually, your generic optimizer mystic.fmin_powell is a copy of the Scipy code included in your package. I would suggest one calls Scipy directly. – Isosteric

@divenex: You are wrong, actually. While the fmin_powell code is indeed derived from scipy.optimize, it has been modified to accept arbitrary hard and soft constraints (and can thus it can execute in a space that effectively makes it a QP, LP, or MIP solver, if desired). If you create hard QP constraints, any of the mystic solvers that you apply the constraints to will only look for solutions in QP space. This was the entire point of my post, and one of the major features of mystic. Please feel free to contact me directly if you have questions on how mystic works. – Sickle

The fact that your modified Powell algorithm also optimizes a constrained quadratic or linear function does not make it "effectively" a Quadratic Programming (QP) solver. The advantage of a QP solver is that it exploits the quadratic form of the function for much faster and more robust convergence. The disadvantage is that a true QP solver, unlike the Powell algorithm, cannot solve other types of functions, like generic nonlinear ones. Here are links to true Python QP solvers, but none "only depends on NumPy/SciPy" as requested scaron.info/blog/quadratic-programming-in-python.html – Isosteric

Thanks for the link. It might help you to look into how constraints work in mystic, they are essentially coordinate transforms. – Sickle

Answer 6

can try python module quadprog for QP-tasks like OLS & try example

import numpy as np
import quadprog
from scipy import optimize
# e.g. least-squares problem
#print("    min. || G * x - a ||^2")
#print("    s.t. C * x <= b")
#print('\n')  

def solve_qp_scipy(G, a, C, b, meq):
    def f(x):
        return 0.5 * np.dot(x, G).dot(x) - np.dot(a, x)

    constraints = []
    if C is not None:
        constraints = [{
            'type': 'eq' if i < meq else 'ineq',
            'fun': lambda x, C=C, b=b, i=i: (np.dot(C.T, x) - b)[i]
        } for i in range(C.shape[1])]

    result = optimize.minimize(
        f, x0=np.zeros(len(G)), method='SLSQP', constraints=constraints,
        tol=1e-10, options={'maxiter': 2000})
    return result

# data: 
M = np.array([[1.0, 2.0, 0.0], [-8.0, 3.0, 2.0], [0.0, 1.0, 1.0]])
G = np.dot(M.T, M)  # this is a positive definite matrix
a = np.dot(np.array([3.0, 2.0, 3.0]), M)   # @ M
C = np.array([[1.0, 2.0, 1.0], [2.0, 0.0, 1.0], [-1.0, 2.0, -1.0]])
b = np.array([3.0, 2.0, -2.0])     #.reshape((3,))
meq = b.size

# for dense problem from https://github.com/qpsolvers/qpsolvers/blob/main/examples/quadratic_program.py
print("    min. 1/2 x^T P x + a^T x")
print("    s.t. G * x <= b")

# quadprog
xf, f, xu, iters, lagr, iact = quadprog.solve_qp(G, a, C, b, meq)
print(xf, '\n', f, '\n')

# SciPy
result = solve_qp_scipy(G, a, C, b, meq)
print(result.x, '\n', result.fun, '\n')

np.testing.assert_array_almost_equal(result.x, xf)
np.testing.assert_array_almost_equal(result.fun, f)

P.S. other qpsolvers, that I couldn't pip installwith msys32\mingw32 gcc 9.3.0, vs_BT2019 needed, I think

P.P.S. valuable references here - "How not to lie with statistics: the correct way to summarize benchmark results"

Sudanic

Answer 7

You can use cvxpy to solve a quadratic program (this package requires clarabel, ecos, numpy, osqp, pybind11, scipy, scs), and it has nice syntax that allows for equations to converted to code pretty intuitively. The authors are Steven Diamond, Eric Chu, and Stephen Boyd who wrote the well-known Convex Optimization textbook.

Just note that you may want to convert the quadratic optimization function (and affine constraints) into the standard form as shown in the documentation:

minimize       f_0(x) = x^T P x + q^T x
subject to     Gx ≤ h
               Ax = b

For example, if you want to solve the following QP:

minimize       x1^2 + x2^2 - x1x2 - x1
subject to     -x1 + x2 ≤ 3
               x1 - 4x2 ≤ 1
               2x1 + 5x2 ≤ 10

You could set up the objective and constraints, and extract optimal value p* with optimal primal variables x1*, x2* as follows:

import numpy as np
import cvxpy as cp

x = cp.Variable(2)

P = np.array([
    [1, -1/2],
    [-1/2, 1]
])
q = np.array([-1, 0])


A = np.array([
    [-1, 1],
    [1, -4],
    [2, 5]
])
b = np.array([3, 1, 10])

constraints = [
    A @ x <= b
]

obj = cp.Minimize(cp.quad_form(x, P) + q.T @ x)

# Form and solve problem.
prob = cp.Problem(obj, constraints)
prob.solve()
p_star = (cp.quad_form(x, P) + q.T @ x).value

x1, x2 = x.value

And you get optimal value p* = -1/3 with optimal primal variables x1*=2/3, x2*=1/3

Kristankriste