7

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!

A.Torres
  • 413
  • 1
  • 6
  • 16
  • 2
    Did you look inside [`scipy.optimize`](https://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html)? 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. – abarnert Jul 03 '18 at 19:03
  • 1
    If there's nothing there, you'll need a third-party library. From a quick search, there's something called [`quadprog`](https://pypi.org/project/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. – abarnert Jul 03 '18 at 19:06
  • 1
    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). – abarnert Jul 03 '18 at 19:07
  • Did you look at cvxpy, a library that let you easily implement convex optimization (hence also quadratic programming)? – linello Aug 31 '21 at 09:07

4 Answers4

4

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]
4

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.

Lisa
  • 194
  • 7
2

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)
bstellato
  • 160
  • 1
  • 6
  • 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! – Prateek Jain Mar 03 '21 at 09:59
1

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)
Tastalian
  • 369
  • 4
  • 12