Convex Optimization in Python: 3 Essential Packages
Python has taken a big role in scientific computation and engineering. Packages, such as NumPy and SciPy, are heavily relied on in various computational applications and day-to-day research. How about computational optimization? We will have a look at the 3 most popular Python packages for computational optimization. More specifically, we will have a look at the essential convex optimization tools for Python.
Convex optimizations is a well-established field within computational optimization. Convex problems have a particular property, where an optimal solution can be guaranteed to be found under fairly mild restrictions. Various efficient algorithms like the interior-point method have been used for decades to solve complex problems. Around these algorithms, there are a lot of solvers that can handle a wide range of convex problems. Thus, solving convex problems is not really problematic. The problematic and time-consuming part is to write the optimization problems in a way that can be fed to the solvers.
A case could be made that Matlab package CVX was a groundbreaking framework, which brought high-level convex optimization modeling to anyone's reach. The researchers were able to apply convex optimization easily to complex problems, without having to spend tons of time bringing the problems into a standard form that could be fed to the solver libraries.
In this post, we will go through 3 popular Python packages that make convex optimization more accessible to researchers and developers familiar with Python. These packages are
Example problem
We will consider a simple example problem as a test case for each of the packages. The problem is a simple quadratic optimization problem with linear constraints
\[ \begin{array}{rl} \underset{x, y}{\text{min.}} & \displaystyle x^2 + y^2 \\ \text{s.t.} & \displaystyle x + y = 1, \\ & \displaystyle x \geq 0, y \geq 0 \end{array} \]
This problem has a trivial solution \(x = 0.5, y = 0.5\), which each of the packages should be able to find.
CVXPY
CVXPY is a modeling framework for convex problems. It closely resembles the famous CVX Matlab package, which is not too surprising as it originates from the same research groups at Stanford University. We have also used this one before in our solution for orthogonal channel sum-rate maximization.
Similar to CVX, CVXPY is a high-level modeling tool that translates complex convex problems into various different kinds of standard forms that are then fed to a solver for processing. It follows the disciplined convex programming approach, where you can freely use a set of functions and sets like exponential function, logarithm, or conic sets to formulate a problem. CVXPY can then translate any combination of these functions and sets into a standard form that is supported by a convex solver. It has also partial support for mixed-integer and quasiconvex problems.
It supports a wide range of solvers: OSQP, ECOS, GLOP, MOSEK, CBC, CVXOPT, NAG, PDLP, GUROBI, and SCS. Not all solvers are equal, and sometimes you have to try different ones before you find the one that can actually handle your problem. This is particularly the case when the complexity and numerical accuracy requirements grow.
Implementing our example problem in CVXPY is really easy. It follows more or less the same notation as in the mathematical model. We can directly write the implementation just by following the mathematical formulation of our problem
import cvxpy as cp
x = cp.Variable(1) # Real variable x
y = cp.Variable(1) # Real variable y
objective = cp.Minimize(x**2 + y**2) # Objective
constraints = []
constraints.append(x >= 0)
constraints.append(y >= 0)
constraints.append(x + y == 1)
problem = cp.Problem(objective, constraints)
opt = problem.solve()
print("Optimal solution: %f (x = %f, y = %f)" % (opt, x.value, y.value))Solution
Optimal solution: 0.500000 (x = 0.500000, y = 0.500000)PICOS
PICOS stands for Python Interface to Conic Optimization Solvers. Even though the name specifies conic problems, the framework is very powerful and flexible. It can be used to easily model a wide variety of problems. In addition to continuous convex problems, it has support for mixed-integer optimization.
In many ways, PICOS is very similar to CXVPY as they both support similar high-level approaches for modeling convex problems. I'd say that for most cases, it is pretty much just up to personal preference, which one to choose.
PICOS supports all the standard solvers, namely, CPLEX, CVXOPT, ECOS, GLPK, Gurobi, MOSEK, QSQP, SCIP, and SMCP.
How would we model our example problem with PICOS? This is actually quite similar to the CVXPY implementation. Since PICOS is a high-level modeling framework, we can pretty much write our problem in the same form that was used in the mathematical presentation.
import picos as pc
x = pc.RealVariable("x", 1) # Real variable x
y = pc.RealVariable("y", 1) # Real variable y
problem = pc.Problem()
problem.minimize = x**2 + y**2 # Objective
problem += x >= 0
problem += y >= 0
problem += x + y == 1
problem.solve(solver="cvxopt")
opt = problem.value
print("Optimal solution: %f (x = %f, y = %f)" % (opt, x, y))Solution
Optimal solution: 0.500000 (x = 0.500000, y = 0.500000)CVXOPT
CVXOPT takes the more traditional approaches, where you need to manually transform the optimization into a standard for the can be fed to the solver. It is actually used under the hood by CVXPY and PICOS.
If you are willing to your problem into the standard form that is accepted by solvers, you can greatly optimize the run time of your code by using CVXOPT directly.
However, translating a complex problem into a standard can be error-prone. Thus, it is highly encouraged to use the aforementioned higher-level modeling frameworks for experimentation.
How would we feed our example problem into the solver? As mentioned, it is a quadratic problem and CVXOPT accepts quadratic problems in the form
\[ \begin{array}{rl} \underset{\mathbf{x}}{\text{min.}} & \displaystyle {1 \over 2}\mathbf{x}^T \mathbf{Q} \mathbf{x} + \mathbf{p}^T \mathbf{x} \\ \text{s.t.} & \displaystyle \mathbf{G}\mathbf{x} \preceq \mathbf{h}, \\ & \displaystyle \mathbf{A}\mathbf{x} = b \end{array} \]
In our case, this would mean that \(\mathbf{x} = \begin{bmatrix} x\\ y \end{bmatrix}\), \(\mathbf{Q} = \begin{bmatrix} 2 & 0\\ 0 & 2\end{bmatrix}\), \(\mathbf{p} = \begin{bmatrix} 0\\ 0 \end{bmatrix}\), \(\mathbf{G} = \begin{bmatrix} -1 & 0\\ 0 & -1\end{bmatrix}\), \(\mathbf{h} = \begin{bmatrix} 0\\ 0 \end{bmatrix}\), \(\mathbf{A} = [1, 1]\) and \(b = 1\).
Once the correct form has been achieved, the implementation becomes quite straightforward. Just plug in the values and pass them on to the solver
from cvxopt import matrix, solvers
Q = matrix([ [2, 0.0], [0.0, 2] ])
p = matrix([0.0, 0.0])
G = matrix([[-1.0,0.0],[0.0,-1.0]])
h = matrix([0.0,0.0])
A = matrix([1.0, 1.0], (1,2))
b = matrix(1.0)
sol = solvers.qp(Q, p, G, h, A, b)
print("Optimal solution: %f (x = %f, y = %f)" % (sol['primal objective'], sol['x'][0], sol['x'][1]))Solution
Optimal solution: 0.500000 (x = 0.500000, y = 0.500000)Which one to choose?
Between CVXPY and PICOS, it is a pretty close call and really just up to your personal preference. Both are great high-level tools that allow you to experiment with optimization problems.
If you need to go after the performance, the high-level modeling does add an extra interpretation layer, which will have an impact on the performance. How much impact? It really depends on the problem's complexity. If performance is a key criterion go to the lower levels and choose CVXOPT. It requires some manual work to translate your problem into one of the standard forms accepted by a solver. But if you are ready to do this, you will be paid off by a faster run time.