OSQP: an operator splitting solver for quadratic programs

Abstract

We present a general-purpose solver for convex quadratic programs based on the alternating direction method of multipliers, employing a novel operator splitting technique that requires the solution of a quasi-definite linear system with the same coefficient matrix at almost every iteration. Our algorithm is very robust, placing no requirements on the problem data such as positive definiteness of the objective function or linear independence of the constraint functions. It can be configured to be division-free once an initial matrix factorization is carried out, making it suitable for real-time applications in embedded systems. In addition, our technique is the first operator splitting method for quadratic programs able to reliably detect primal and dual infeasible problems from the algorithm iterates. The method also supports factorization caching and warm starting, making it particularly efficient when solving parametrized problems arising in finance, control, and machine learning. Our open-source C implementation OSQP has a small footprint, is library-free, and has been extensively tested on many problem instances from a wide variety of application areas. It is typically ten times faster than competing interior-point methods, and sometimes much more when factorization caching or warm start is used. OSQP has already shown a large impact with tens of thousands of users both in academia and in large corporations.

Problem classes

In this section we describe the random problem classes used in the benchmarks and derive formulations with explicit linear equalities and inequalities that can be directly written in the form $Ax \in {\mathcal {C}}$ with ${\mathcal {C}} = [l, u]$.

1.1 Random QP

Consider the following QP

$$\begin{aligned} \begin{array}{ll} \text{ minimize } &{}\quad (1/2) x^TP x + q^Tx \\ \text{ subject } \text{ to } &{}\quad l \le A x \le u. \end{array} \end{aligned}$$

Problem instances The number of variables and constraints in our problem instances are n and $m=10n$. We generated random matrix $P= M M^T+ \alpha I$ where $M\in {\mathbf{R}}^{n \times n}$ and $15 \%$ nonzero elements $M_{ij}\sim {\mathcal {N}}(0,1)$. We add the regularization $\alpha I$ with $\alpha = 10^{-2}$ to ensure that the problem is not unbounded. We set the elements of $A\in {\mathbf{R}}^{m \times n}$ as $A_{ij} \sim {\mathcal {N}}(0,1)$ with only $15 \%$ being nonzero. The linear part of the cost is normally distributed, i.e., $q_i \sim {\mathcal {N}}(0,1)$. We generated the constraint bounds as $u_i \sim {\mathcal {U}}(0,1)$, $l_i \sim -{\mathcal {U}}(0,1)$.

1.2 Equality constrained QP

Consider the following equality constrained QP

$$\begin{aligned} \begin{array}{ll} \text{ minimize } &{}\quad (1/2) x^TP x + q^Tx \\ \text{ subject } \text{ to } &{}\quad A x = b. \end{array} \end{aligned}$$

This problem can be rewritten as (1) by setting $l = u = b$.

Problem instances The number of variables and constraints in our problem instances are n and $m=\lfloor n/2\rfloor$.

We generated random matrix $P= M M^T+ \alpha I$ where $M\in {\mathbf{R}}^{n \times n}$ and $15 \%$ nonzero elements $M_{ij}\sim {\mathcal {N}}(0,1)$. We add the regularization $\alpha I$ with $\alpha = 10^{-2}$ to ensure that the problem is not unbounded. We set the elements of $A\in {\mathbf{R}}^{m \times n}$ as $A_{ij} \sim {\mathcal {N}}(0,1)$ with only $15 \%$ being nonzero. The vectors are all normally distributed, i.e., $q_i, b_i \sim {\mathcal {N}}(0,1)$.

Iterative refinement interpretation Solution of the above problem can be found directly by solving the following linear system

$$\begin{aligned} \begin{bmatrix} P &{}\quad A^T\\ A &{}\quad 0\end{bmatrix} \begin{bmatrix}x \\ \nu \end{bmatrix} = \begin{bmatrix} -q \\ b \end{bmatrix}. \end{aligned}$$

(35)

If we apply the ADMM iterations (15)–(19) for solving the above problem, and by setting $\alpha =1$ and $y^0=b$, the algorithm boils down to the following iteration

$$\begin{aligned} \begin{bmatrix} x^{k+1} \\ \nu ^{k+1} \end{bmatrix} = \begin{bmatrix} x^k \\ \nu ^k \end{bmatrix} + \begin{bmatrix} P + \sigma I &{}\quad A^T\\ A &{}\quad -\rho ^{-1} I \end{bmatrix}^{-1}\left( \begin{bmatrix} -q \\ b \end{bmatrix} - \begin{bmatrix} P &{}\quad A^T\\ A &{}\quad 0\end{bmatrix} \begin{bmatrix} x^k \\ \nu ^k \end{bmatrix} \right) , \end{aligned}$$

which is equivalent to (31) with $g = (-q, b)$ and ${\hat{t}}^k = (x^k,\nu ^k)$. This means that Algorithm 1 applied to solve an equality constrained QP is equivalent to applying iterative refinement [32, 92] to solve the KKT system (35). Note that the perturbation matrix in this case is

$$\begin{aligned} \varDelta K = \begin{bmatrix} \sigma I \\ &{}\quad -\rho ^{-1} I \end{bmatrix}, \end{aligned}$$

which justifies using a low value of $\sigma$ and a high value of $\rho$ for equality constraints.

1.3 Optimal control

We consider the problem of controlling a constrained linear time-invariant dynamical system. To achieve this, we formulate the following optimization problem [12]

$$\begin{aligned} \begin{array}{ll} \text{ minimize } &{}\quad x_T^TQ_T x_T + \sum \limits _{t=0}^{T-1} x_t^TQ x_t + u_t^TR u_t \\ \text{ subject } \text{ to } &{}\quad x_{t+1} = Ax_t + Bu_t \\ &{} \quad x_t \in {\mathcal {X}}, u_t \in {\mathcal {U}} \\ &{} \quad x_0 = x_{\mathrm{init}}. \end{array} \end{aligned}$$

(36)

The states $x_t\in {\mathbf{R}}^{n_x}$ and the inputs $u_k\in {\mathbf{R}}^{n_u}$ are subject to polyhedral constraints defined by the sets ${\mathcal {X}}$ and ${\mathcal {U}}$. The horizon length is T and the initial state is $x_{\mathrm{init}}\in {\mathbf{R}}^{n_x}$. Matrices $Q\in {\mathbf{S}}^{n_x}_+$ and $R\in {\mathbf{S}}^{n_u}_{++}$ define the state and input costs at each stage of the horizon, and $Q_T\in {\mathbf{S}}^{n_x}_+$ defines the final stage cost.

By defining the new variable $z = (x_0, \ldots , x_{T}, u_0, \ldots , u_{T-1})$, problem (36) can be written as a sparse QP of the form (2) with a total of $n_x(T+1) + n_u T$ variables.

Problem instances We defined the linear systems with $n = n_x$ states and $n_u = 0.5n_x$ inputs. We set the horizon length to $T=10$. We generated the dynamics as $A = I + \varDelta$ with $\varDelta _{ij} \sim {\mathcal {N}}(0, 0.01)$. We chose only stable dynamics by enforcing the norm of the eigenvalues of A to be less than 1. The input action is modeled as B with $B_{ij} \sim {\mathcal {N}}(0, 1)$.

The state cost is defined as $Q = {\mathbf{diag}}(q)$ where $q_i \sim {\mathcal {U}}(0, 10)$ and $70\%$ nonzero elements in q. We chose the input cost as $R = 0.1I$. The terminal cost $Q_T$ is chosen as the optimal cost for the linear quadratic regulator (LQR) applied to A, B, Q, R by solving a discrete algebraic Riccati equation (DARE) [12]. We generated input and state constraints as

$$\begin{aligned} {\mathcal {X}} = \{x_t \in {\mathbf{R}}^{n_x} \mid -{\overline{x}} \le x_t \le {\overline{x}} \}, \qquad {\mathcal {U}}=\{u_t \in {\mathbf{R}}^{n_u} \mid -{\overline{u}} \le u_t \le {\overline{u}}\}, \end{aligned}$$

where ${\overline{x}}_i\sim {\mathcal {U}}(1, 2)$ and ${\overline{u}}_i \sim {\mathcal {U}}(0, 0.1)$. The initial state is uniformly distributed with $x_{\mathrm{init}} \sim {\mathcal {U}}(-0.5{\overline{x}}, 0.5{\overline{x}})$.

1.4 Portfolio optimization

Portfolio optimization is a problem arising in finance that seeks to allocate assets in a way that maximizes the risk adjusted return [13, 15, 67, 17, §4.4.1],

$$\begin{aligned} \begin{array}{ll} \text{ maximize } &{}\quad \mu ^Tx - \gamma (x^T\varSigma x) \\ \text{ subject } \text{ to } &{}\quad {\mathbf {1}}^Tx = 1 \\ &{} \quad x \ge 0, \end{array} \end{aligned}$$

where the variable $x \in {\mathbf{R}}^{n}$ represents the portfolio, $\mu \in {\mathbf{R}}^{n}$ the vector of expected returns, $\gamma > 0$ the risk aversion parameter, and $\varSigma \in {\mathbf{S}}_+^n$ the risk model covariance matrix. The risk model is usually assumed to be the sum of a diagonal and a rank $k < n$ matrix

$$\begin{aligned} \varSigma = FF^T+ D, \end{aligned}$$

where $F\in {\mathbf{R}}^{n\times k}$ is the factor loading matrix and $D\in {\mathbf{R}}^{n\times n}$ is a diagonal matrix describing the asset-specific risk.

We introduce a new variable $y=F^Tx$ and solve the resulting problem in variables x and y

$$\begin{aligned} \begin{array}{ll} \text{ minimize } &{}\quad x^TD x + y^Ty - \gamma ^{-1} \mu ^Tx \\ \text{ subject } \text{ to } &{}\quad y=F^Tx \\ &{}\quad {\mathbf {1}}^Tx = 1 \\ &{}\quad x \ge 0, \end{array} \end{aligned}$$

(37)

Note that the Hessian of the objective in (37) is a diagonal matrix. Also, observe that $FF^T$ does not appear in problem (37).

Problem instances We generated portfolio problems for increasing number of factors k and number of assets $n=100k$. The elements of matrix F were chosen as $F_{ij} \sim {\mathcal {N}}(0,1)$ with $50 \%$ nonzero elements. The diagonal matrix D is chosen as $D_{ii} \sim {\mathcal {U}}[0, \sqrt{k}]$. The mean return was generated as $\mu _i \sim {\mathcal {N}}(0,1)$. We set $\gamma = 1$.

1.5 Lasso

The least absolute shrinkage and selection operator (Lasso) is a well known linear regression technique obtained by adding an $\ell _1$ regularization term in the objective [19, 90]. It can be formulated as

$$\begin{aligned} \begin{array}{ll} \text{ minimize }&\quad \Vert Ax - b\Vert _2^2 + \lambda \Vert x\Vert _1, \end{array} \end{aligned}$$

where $x\in {\mathbf{R}}^{n}$ is the vector of parameters and $A\in {\mathbf{R}}^{m \times n}$ is the data matrix and $\lambda$ is the weighting parameter.

We convert this problem to the following QP

$$\begin{aligned} \begin{array}{ll} \text{ minimize } &{}\quad y^Ty + \lambda {\mathbf {1}}^Tt\\ \text{ subject } \text{ to } &{}\quad y = Ax - b\\ &{}\quad -t \le x \le t, \end{array} \end{aligned}$$

where $y\in {\mathbf{R}}^{m}$ and $t\in {\mathbf{R}}^{n}$ are two newly introduced variables.

Problem instances The elements of matrix A are generated as $A_{ij} \sim {\mathcal {N}}(0,1)$ with $15 \%$ nonzero elements. To construct the vector b, we generated the true sparse vector $v\in {\mathbf{R}}^{n}$ to be learned

$$\begin{aligned} v_i \sim {\left\{ \begin{array}{ll} 0 &{}\quad \text{ with } \text{ probability } p=0.5\\ {\mathcal {N}}(0, 1/n) &{} \quad \text{ otherwise }. \end{array}\right. } \end{aligned}$$

Then we let $b=Av + \varepsilon$ where $\varepsilon$ is the noise generated as $\varepsilon _i \sim {\mathcal {N}}(0, 1)$. We generated the instances with varying n features and $m = 100n$ data points. The parameter $\lambda$ is chosen as $(1/5)\Vert A^Tb\Vert _{\infty }$ since $\Vert A^Tb\Vert _{\infty }$ is the critical value above which the solution of the problem is $x=0$.

1.6 Huber fitting

Huber fitting or the robust least-squares problem performs linear regression under the assumption that there are outliers in the data [55, 56]. The fitting problem is written as

$$\begin{aligned} \begin{array}{ll} \text{ minimize }&\quad \sum \nolimits _{i=1}^m \phi _{\mathrm{hub}}(a_i^T x - b_i), \end{array} \end{aligned}$$

(38)

with the Huber penalty function $\phi _{\mathrm{hub}}:{\mathbf{R}}\rightarrow {\mathbf{R}}$ defined as

$$\begin{aligned} \phi _{\mathrm{hub}}(u) = {\left\{ \begin{array}{ll} u^2 &{}\quad |u| \le M \\ M(2|u|-M) &{}\quad |u| > M. \end{array}\right. } \end{aligned}$$

Problem (38) is equivalent to the following QP [66, Eq. (24)]

$$\begin{aligned} \begin{array}{ll} \text{ minimize } &{}\quad u^Tu + 2M {\mathbf {1}}^T(r + s) \\ \text{ subject } \text{ to } &{}\quad Ax - b - u = r - s \\ &{} \quad r \ge 0 \\ &{} \quad s \ge 0. \end{array} \end{aligned}$$

Problem instances We generate the elements of A as $A_{ij} \sim {\mathcal {N}}(0,1)$ with $15 \%$ nonzero elements. To construct $b\in {\mathbf{R}}^m$ we first generate a vector $v\in {\mathbf{R}}^n$ as $v_i \sim {\mathcal {N}}(0,1/n)$ and a noise vector $\varepsilon \in {\mathbf{R}}^m$ with elements

$$\begin{aligned} \varepsilon _i \sim {\left\{ \begin{array}{ll} {\mathcal {N}}(0,1/4) &{}\quad \hbox { with probability}\ p=0.95 \\ {\mathcal {U}}[0,10] &{} \quad \text {otherwise}. \end{array}\right. } \end{aligned}$$

We then set $b = Av + \varepsilon$. For each instance we choose $m=100n$ and $M=1$.

1.7 Support vector machine

Support vector machine problem seeks an affine function that approximately classifies the two sets of points [21]. The problem can be stated as

$$\begin{aligned} \begin{array}{ll} \text{ minimize }&\quad x^Tx + \lambda \sum \nolimits _{i=1}^m \max (0, b_i a_i^T x + 1), \end{array} \end{aligned}$$

where $b_i \in \{ -1, +1 \}$ is a set label, and $a_i$ is a vector of features for the ith point. The problem can be equivalently represented as the following QP

$$\begin{aligned} \begin{array}{ll} \text{ minimize } &{}\quad x^Tx + \lambda {\mathbf {1}}^Tt \\ \text{ subject } \text{ to } &{}\quad t \ge {\mathbf{diag}}(b)Ax + {\mathbf {1}}\\ &{}\quad t \ge 0, \end{array} \end{aligned}$$

where ${\mathbf{diag}}(b)$ denotes the diagonal matrix with elements of b on its diagonal.

Problem instances We choose the vector b so that

$$\begin{aligned} b_i = {\left\{ \begin{array}{ll} +1 &{}\quad i \le m/2 \\ -1 &{}\quad \text {otherwise}, \end{array}\right. } \end{aligned}$$

and the elements of A as

$$\begin{aligned} A_{ij} \sim {\left\{ \begin{array}{ll} {\mathcal {N}}(+1/n,1/n) &{} \quad i \le m/2 \\ {\mathcal {N}}(-1/n,1/n) &{} \quad \text {otherwise}, \end{array}\right. } \end{aligned}$$

with $15\%$ nonzeros per case.