Embedded Code Generation with CVXPY

Convex optimization is used in many domains, including signal and image processing [1, 2], control [3, 4], and finance [5, 6], to mention just a few. A (parametrized) convex optimization problem can be written as

where $x\in\mathbf{R}^{n}$ is the optimization variable, $f_{0}$ is the objective function to be minimized, $f_{1},\ldots,f_{p}$ are the inequality constraint functions, and $g_{1}\ldots,g_{r}$ are the equality constraint functions. We require that $f_{0},\ldots,f_{p}$ are convex functions, and $g_{1},\ldots,g_{r}$ are affine functions [7]. The parameter $\theta\in\mathbf{R}^{d}$ specifies data that can change, but is constant and given when we solve an instance of the problem. We refer to the parametrized problem (1) as a problem family; when we specify a fixed value of $\theta$, we refer to it as a problem instance. We let $x^{\star}$ denote an optimal point for the problem (1), assuming it exists.

The problem family can be specified using a domain-specific language (DSL) for convex optimization. Such systems allow the user to specify the functions $f_{i}$ and $g_{j}$ in a simple format that closely follows the mathematical description of the problem. Examples include YALMIP [8] and CVX [9] (in Matlab), CVXPY [10] (in Python), Convex.jl [11] and JuMP [12] (in Julia), and CVXR [13] (in R). We focus on CVXPY, which also supports the declaration of parameters, enabling it to specify problem families, not just problem instances.

DSLs parse the problem description and translate (or canonicalize) it to an equivalent problem that is suitable for a solver that handles some generic class of problems, such as linear programs (LPs), quadratic programs (QPs), second-order cone programs (SOCPs), semidefinite programs (SDPs), and others such as exponential cone programs [7]. We focus on solvers that are suitable for embedded applications, i.e., are single-threaded, can be statically compiled, and do not make system calls: OSQP [14] handles QPs, SCS [15] and ECOS [16] handle cone programs that include SOCPs. After the canonicalized problem is solved, a solution of the original problem is retrieved from a solution of the canonicalized problem.

It is useful to think of the whole process as a function that maps $\theta$, the parameter that specifies the problem instance, into $x^{\star}$, an optimal value of the variable. With a DSL, this process consists of three steps. First the original problem description is canonicalized to a problem in some standard (or canonical) form; then the canonicalized problem is solved using a solver; and finally, a solution of the original problem is retrieved from a solution of the canonicalized problem.

Most of these DSLs are organized as parser-solvers, which carry out the canonicalization each time the problem is solved (with different parameter values). This simple setting is illustrated in figure 1(a). We are interested in applications where we solve many instances of the problem, possibly in an embedded application with hard real-time constraints. For such applications, a code generator makes more sense. A code generator takes as input a description of a problem family, and generates specialized solver source code for that specific family. That source code is then compiled, and we have an efficient solver for the specific family.

This workflow is illustrated in figure 1(b). The compiled solver has a number of advantages over parser-solvers. First, the compiled solver can be deployed in embedded systems, fulfilling rules for safety-critical code [17]. Second, by caching canonicalization and exploiting the problem structure, the compiled solver is faster.

A well known code generator is CVXGEN [18]. It handles problems that can be transformed to QPs and includes a custom interior-point solver. CVXGEN is used in many applications, including autonomous driving, dynamic energy management, and real-time trading in finance. SpaceX uses CVXGEN to generate flight code for high-speed onboard convex optimization for precision landing of space rockets [19].

CVXGEN was designed for use in real-time control systems, where the problems solved are not too big, either in terms of the number of variables or number of parameters. CVXGEN unrolls for-loops in its generated source code files to increase the solving speed, but this can also result in large compiled code size. Due to the flat and explicit code generated, CVXGEN only handles problems with up to around a few thousand parameters. (More accurately, CVXGEN is limited to 4000 nonzero entries in the linear system of equations solved in each iteration.)

CVXPYgen is based on the open-source Python-embedded DSL CVXPY. CVXPY handles many types of conic programs and certain types of nonconvex problems, whereas we focus on LPs, QPs, and SOCPs for code generation. CVXPY provides modeling instructions that follow the mathematical description for convex optimization problems. It ensures that the modeled problems are convex, using disciplined convex programming (DCP). DCP is the process of constructing convex functions by assembling given base functions in mathematical expressions using a simple set of rules [24]. DCP ensures that the resulting problem is convex, and also, readily canonicalized to a standard form.

In DCP, parameters are treated as constants, optionally with specified sign, and there are no restrictions about how these constants appear in the expressions defining the problem family. The recently developed concept of disciplined parametrized programming (DPP) puts additional restrictions on how parameters can enter a problem description. If a problem family description is DPP-compliant, then canonicalization and retrieval can be represented as affine mappings [25]. Thus DPP-compliant problems are reducible to ASA-form, which stands for Affine-Solve-Affine [25]. This is the key property we exploit in CVXPYgen. More about the DCP and DPP rules can be found in the aforementioned papers, or at https://www.cvxpy.org.

After CVXPY has reduced the DPP-compliant problem to ASA-form, CVXPYgen extracts a sparse matrix $C$ that canonicalizes the user-defined parameters $\theta$ to the parameters $\tilde{\theta}$ appearing in the standard form solver:

CVXPYgen analyzes $C$ to determine the user-defined parameters (i.e., components of $\theta$) that every standardized form parameter depends on. This information is used when generating the custom solver, where only slices of the above mapping are computed if not all user-defined parameters are updated between solves. In addition, it is very useful to know the set of updated canonical parameters when using the OSQP solver or SCS, as detailed below.

In a similar way the retrieval of the solution $x^{\star}$ for the original problem from a solution $\tilde{x}^{\star}$ of the canonicalized problem is an affine mapping,

where $R$ is a sparse matrix. Typically $R$ is a selector matrix, with only one nonzero entry in each row, equal to one, in which case this step can be handled via simple pointers in C.

CVXPYgen generates allocation-, library-, and division-free C code for the canonicalization and retrieval steps, which in essence are nothing more than sparse matrix-vector multiplication, with some logic that exploits pointers or partial updates. Sparse matrices are stored in compressed sparse column format [26] and dense matrices are stored as vectors via column-major flattening.

Any solver can be used to solve the canonicalized problem, which provides the final link:

where $\mathcal{S}$ denotes the mapping from the canonicalized parameters to a solution of the canonicalized problem. (We assume here that the problem instance is feasible, and that when there are multiple solutions, we simply pick one.) If available, CVXPYgen uses the canonical solver’s code generation method to produce C code for canonical solving. As of now, only OSQP provides this functionality [27]. Otherwise, the solver’s C code is simply copied, possibly modified for use in embedded applications.

OSQP and SCS provide a set of C functions for updating their parameters. This way, when only canonical vector parameters are updated, the factorization of the linear system involved in the OSQP or SCS algorithms can be cached and re-used, which can lead to substantial speed up and division-free code. In the same way as only the parts of $\tilde{\theta}$ are re-canonicalized that depend on the updated parts of $\theta$, only the OSQP or SCS update functions associated with these parameters are called before the canonicalized problem is solved.

The code and the full documentation for CVXPYgen with its generated solvers are available at

https://pypi.org/project/cvxpygen.

We consider the nonnegative least squares problem

where $x\in\mathbf{R}^{n}$ is the variable and $G\in\mathbf{R}^{m\times n}$, $h\in\mathbf{R}^{m}$ are parameters, so $\theta=(G,h)$. We will canonicalize this to the standard form accepted by OSQP,

where $\tilde{x}\in\mathbf{R}^{\tilde{n}}$ is the canonical variable and all other symbols are canonical parameters, i.e., $\tilde{\theta}=(P,q,A,l,u)$. (In this form, entries of $l$ can be $-\infty$, and entries of $u$ can be $+\infty$.)

The naïve canonicalization of (2) to (3) takes $\tilde{x}=x$ and

In this canonicalization, $\tilde{\theta}$ is not an affine function of $\theta$, since some entries of $\tilde{\theta}$ are products of entries of $\theta$.

The canonicalization that uses DPP first expresses problem (2) as

with variable $\tilde{x}=(\tilde{x}_{1},\tilde{x}_{2})$, where $\tilde{x}_{1}=x$ and $\tilde{x}_{2}\in\mathbf{R}^{m}$. We can express this as (3) with parameters

where the second part of $u$ has $\infty$ in every entry, i.e., there is no upper bound on the second part of $A\tilde{x}$. In this canonicalization, $\tilde{\theta}$ is indeed an affine function of $\theta$. The retrieval map has the simple (linear) form $x^{\star}=[I~{}0]\tilde{x}^{\star}$.

We generate code for this problem as shown in figure 2. The problem is modeled with CVXPY in lines 5–9. The actual code generation is done in line 11. The variable and parameters are named in lines 5–7 via their name attributes. These names are used for C variable and function naming.

We compare CVXPYgen to CVXGEN for a model predictive control (MPC) problem family, described in Appendix A. In particular, we compare solve times of the C interface and executable sizes. The MPC problems are parametrized by their horizon length $H\in\left\{6,12,18,30,60\right\}$; the number of variables is around $10H$. Figure 3 shows the resulting solve times averaged over 100 simulation steps. CVXGEN is not able to generate code for $H>18$. Together with the automatically chosen OSQP solver, CVXPYgen outperforms CVXGEN for all problem sizes.

The bottom of figure 3 presents the example executable sizes for CVXGEN and CVXPYgen, respectively. For all values of $H$, the executables corresponding to CVXPYgen are considerably smaller.

The execution times cited above are on a MacBook Pro 2.3GHz Intel i5. We have also used these generated solvers to control the position of a custom-built 14-by-14 cm quadcopter. The generated code was compiled in a robot operating system (ROS) node, and run on the drone’s Intel Atom x5-Z8350 processor, at 30 Hz. We provide a video of the quadcopter following a circle trajectory at

https://polybox.ethz.ch/index.php/s/MARR9CGaLqmQaJ0.

Here we compare CVXPYgen to CVXPY, for a (financial) portfolio optimization problem family, described in Appendix B. This family of problems is parametrized by the number of assets in the portfolio $N\in\{10,20,40,60,100\}$. The number of variables in these problems is around $2N$. Figure 4 gives the results for 500 solves, a two year back-test using historical data. We see that the average solver speed is about 6 times faster with CVXPYgen, for $N=10$, with the ratio dropping to 2.5 for $N=100$. The execution times are measured on the same MacBook described in the previous section.

An interesting metric is the break-even point, which is the number of instances that need to be solved before CVXPYgen is faster than CVXPY, when we include the code generation and compilation time. This number is around 5000, and not too dependent on $N$. A typical back-test might involve daily trading, with around 250 trading days in each year, over 4 years, with hundreds of different hyper-parameter values, which gives on the order of 100,000 solves, well above this break-even point.

We have described CVXPYgen, a tool for generating custom C code that solves instances of a family of convex optimization problems specified within CVXPY. This gives a seemless path from prototyping an application using Python and CVXPY, to a final embedded implementation in C. In addition to CVXPYgen supporting a wider variety of problems (such as SOCPs) than the state-of-the-art code generator CVXGEN, numerical experiments show that it outperforms CVXGEN in terms of allowable problem size, compiled code size, and solve times. For applications running on general purpose machines, we obtain a significant speedup over CVXPY when many problem instances are to be solved.