Quadratic Programming in Python using Numpy?

Question

I am in the process of translating some MATLAB code into Python. There is one line that is giving me a bit of trouble:

[q,f_dummy,exitflag, output] = quadprog(H,f,-A,zeros(p*N,1),E,qm,[],[],q0,options);

I looked up the documentation in MATLAB to find that the quadprog function is used for optimization (particularly minimization).

I attempted to find a similar function in Python (using numpy) and there does not seem to be any.

Is there a better way to translate this line of code into Python? Or are there other packages that can be used? Do I need to make a new function that accomplishes the same task?

Thanks for your time and help!

Ladder

Did you look inside scipy.optimize? NumPy only has the basic stuff for working with arrays of numbers; SciPy builds a bunch of other math and science stuff on top of it, including a variety of optimization algorithms. – Elihu

If there's nothing there, you'll need a third-party library. From a quick search, there's something called quadprog, and there's Python bindings for something called QuadProg++, and I'm sure there are others as well. (Some of them may not be built on top of NumPy, however.) But at that point, this becomes a library-recommendation question, and Stack Overflow is unfortunately not the place to find those. – Elihu

One last thing: If you do need to go searching, you may first want to decide whether you want an implementation of the same quadratic programming algorithm MATLAB uses (whichever one that is; I'm sure it's documented, and then you can search for Python implementations of that algorithm), or just anything that works for your data (in which case you could try lots of them and see). – Elihu

Did you look at cvxpy, a library that let you easily implement convex optimization (hence also quadratic programming)? – Infatuation

Answer 1

There is a library called CVXOPT that has quadratic programming in it.

def quadprog_solve_qp(P, q, G=None, h=None, A=None, b=None):
    qp_G = .5 * (P + P.T)   # make sure P is symmetric
    qp_a = -q
    if A is not None:
        qp_C = -numpy.vstack([A, G]).T
        qp_b = -numpy.hstack([b, h])
        meq = A.shape[0]
    else:  # no equality constraint
        qp_C = -G.T
        qp_b = -h
        meq = 0
    return quadprog.solve_qp(qp_G, qp_a, qp_C, qp_b, meq)[0] 

Kolb

As Lisa pointed out, this snippet is for quadprog, not CVXOPT. – Organometallic

Answer 2

I will start by mentioning that quadratic programming problems are a subset of convex optimization problems which are a subset of optimization problems.

There are multiple python packages which solve quadratic programming problems, notably

1. cvxopt -- which solves all kinds of convex optimization problems (including quadratic programming problems). This is a python version of the previous cvx MATLAB package.

2. quadprog -- this is exclusively for quadratic programming problems but doesn't seem to have much documentation.

3. scipy.optimize.minimize -- this is a very general minimizer which can solve quadratic programming problems, as well as other optimization problems (convex and non-convex).

You might also benefit from looking at the answers to this stackoverflow post which has more details and references.

Note: The code snippet in user1911226' answer appears to come from this blog post: https://scaron.info/blog/quadratic-programming-in-python.html which compares some of these quadratic programming packages. I can't comment on their answer, but they claim to be mentioning the cvxopt solution, but the code is actually for the quadprog solution.

Paine

Answer 3

OSQP is a specialized free QP solver based on ADMM. I have adapted the OSQP documentation demo and the OSQP call in the qpsolvers repository for your problem.

Note that matrices H and G are supposed to be sparse in CSC format. Here is the script

import numpy as np
import scipy.sparse as spa
import osqp


def quadprog(P, q, G=None, h=None, A=None, b=None,
             initvals=None, verbose=True):
    l = -np.inf * np.ones(len(h))
    if A is not None:
        qp_A = spa.vstack([G, A]).tocsc()
        qp_l = np.hstack([l, b])
        qp_u = np.hstack([h, b])
    else:  # no equality constraint
        qp_A = G
        qp_l = l
        qp_u = h
    model = osqp.OSQP()
    model.setup(P=P, q=q,
                A=qp_A, l=qp_l, u=qp_u, verbose=verbose)
    if initvals is not None:
        model.warm_start(x=initvals)
    results = model.solve()
    return results.x, results.info.status


# Generate problem data
n = 2   # Variables
H = spa.csc_matrix([[4, 1], [1, 2]])
f = np.array([1, 1])
G = spa.csc_matrix([[1, 0], [0, 1]])
h = np.array([0.7, 0.7])
A = spa.csc_matrix([[1, 1]])
b = np.array([1.])

# Initial point
q0 = np.ones(n)

x, status = quadprog(H, f, G, h, A, b, initvals=q0, verbose=True)

Drift

hi @bstellao, I am getting the following error while executing it with Inequality constraints. ERROR in LDL_factor: Error in KKT matrix LDL factorization when computing the nonzero elements. The problem seems to be non-convex ERROR in osqp_setup: KKT matrix factorization. The problem seems to be non-convex. ERROR : Workspace allocation error! – Sal

Answer 4

You could use the solve_qp function from qpsolvers. It can be installed, along with a starter kit of open source solvers, by pip install qpsolvers[open_source_solvers]. Then you can replace your line with:

from qpsolvers import solve_qp

solver = "proxqp"  # or "osqp", "quadprog", "cvxopt", ...
x = solve_qp(H, f, G, h, A, b, initvals=q_0, solver=solver, **options)

There are many solvers available in Python, each coming with their pros and cons. Make sure you try different values for the solver keyword argument to find the one that fits your problem best.

Here is a standalone example based on your question and the other comments:

import numpy as np

from qpsolvers import solve_qp

H = np.array([[4.0, 1.0], [1.0, 2.0]])
f = np.array([1.0, 1])
G = np.array([[1.0, 0.0], [0.0, 1.0]])
h = np.array([0.7, 0.7])
A = np.array([[1.0, 1.0]])
b = np.array([1.0])
q_0 = np.array([1.0, 1.0])

solver = "cvxopt"  # or "osqp", "proxqp", "quadprog", ...
options = {"verbose": True}
x = solve_qp(H, f, G, h, A, b, initvals=q_0, solver=solver, **options)

Organometallic