The solver

Problem statement

OSQP solves convex quadratic programs (QPs) of the form

\[\begin{split}\begin{array}{ll} \mbox{minimize} & \frac{1}{2} x^T P x + q^T x \\ \mbox{subject to} & l \leq A x \leq u \end{array}\end{split}\]

where \(x\in\mathbf{R}^{n}\) is the optimization variable. The objective function is defined by a positive semidefinite matrix \(P \in \mathbf{S}^{n}_{+}\) and vector \(q\in \mathbf{R}^{n}\). The linear constraints are defined by matrix \(A\in\mathbf{R}^{m \times n}\) and vectors \(l\) and \(u\) so that \(l_i \in \mathbf{R} \cup \{-\infty\}\) and \(u_i \in \mathbf{R} \cup \{+\infty\}\) for all \(i \in \{1,\ldots,m\}\).

Algorithm

The solver runs the following ADMM algorithm (for more details see the related papers at the Citing OSQP section):

\[\begin{split}(x^{k+1}, \nu^{k+1}) & \gets \text{solve linear system}\\ \tilde{z}^{k+1} & \gets z^{k} + \rho^{-1}(\nu^{k+1} - y^{k})\\ z^{k+1} &\gets \Pi(\tilde{z}^{k} + \rho^{-1}y^{k})\\ y^{k+1} &\gets y^{k} + \rho (\tilde{z}^{k+1} - z^{k+1})\end{split}\]

where \(\Pi\) is the projection onto the hyperbox \([l,u]\). \(\rho\) is the ADMM step-size.

Linear system solution

The linear system solution is the core part of the algorithm. It can be done using a direct or indirect method.

With a direct linear system solver we solve the following linear system with a quasi-definite matrix

\[\begin{split}\begin{bmatrix} P + \sigma I & A^T \\ A & -\rho^{-1}I \end{bmatrix} \begin{bmatrix} x^{k+1} \\ \nu^{k+1} \end{bmatrix}= \begin{bmatrix}\sigma x^{k} - q \\ z^{k} - \rho^{-1} y^k \end{bmatrix}.\end{split}\]

With an indirect linear system solver we solve the following linear system with a positive definite matrix

\[\left(P + \sigma I + \rho A^T A \right)x^{k+1} = \sigma x^{k} - q + A^T (\rho z^{k} - y^{k}).\]

OSQP core is designed to support different linear system solvers. For their installation see this section. To specify your favorite linear system solver see this section.

Convergence

At each iteration \(k\) OSQP produces a tuple \((x^{k}, z^{k}, y^{k})\) with \(x^{k} \in \mathbf{R}^{n}\) and \(z^{k}, y^{k} \in \mathbf{R}^{m}\).

The primal and and dual residuals associated to \((x^{k}, z^{k}, y^{k})\) are

\[\begin{split}\begin{aligned} r_{\rm prim}^{k} &= Ax^{k} - z^{k}\\ r_{\rm dual}^{k} &= Px^{k} + q + A^{T} y^{k}. \end{aligned}\end{split}\]

Complementary slackness is satisfied by construction at machine precision. If the problem is feasible, the residuals converge to zero as \(k\to\infty\). The algorithm stops when the norms of \(r_{\rm prim}^{k}\) and \(r_{\rm dual}^{k}\) are within the specified tolerance levels \(\epsilon_{\rm prim}>0\) and \(\epsilon_{\rm dual}>0\)

\[ \| r_{\rm prim}^{k} \|_{\infty} \le \epsilon_{\rm prim}, \quad \| r_{\rm dual}^{k} \|_{\infty} \le \epsilon_{\rm dual}.\]

We set the tolerance levels as

\[\begin{split}\epsilon_{\rm prim} &= \epsilon_{\rm abs} + \epsilon_{\rm rel} \max\lbrace \|Ax^{k}\|_{\infty}, \| z^{k} \|_{\infty} \rbrace \\ \epsilon_{\rm dual} &= \epsilon_{\rm abs} + \epsilon_{\rm rel} \max\lbrace \| P x^{k} \|_{\infty}, \| A^T y^{k} \|_{\infty}, \| q \|_{\infty} \rbrace.\end{split}\]

\(\rho\) step-size

To ensure quick convergence of the algorithm we adapt \(\rho\) by balancing the residuals. In default mode, the inteval (i.e., number of iterations) at which we update \(\rho\) is defined by a time measurement. When the iterations time becomes greater than a certain fraction of the setup time, i.e. adaptive_rho_fraction, we set the current number of iterations as the interval to update \(\rho\). The update happens as follows

\[\rho^{k+1} \gets \rho^{k} \sqrt{\frac{\|r_{\rm prim}\|_{\infty}}{\|r_{\rm dual}\|_{\infty}}}.\]

Note that \(\rho\) is updated only if it is sufficiently different than the current one. In particular if it is adaptive_rho_tolerance times larger or smaller than the current one.

Infeasible problems

OSQP is able to detect if the problem is primal or dual infeasible.

Primal infeasibility

When the problem is primal infeasible, the algorithm generates a certificate of infeasibility, i.e., a vector \(v\in\mathbf{R}^{m}\) such that

\[A^T v = 0, \quad u^T v_{+} + l^T v_{-} < 0,\]

where \(v_+=\max(v,0)\) and \(v_-=\min(v,0)\).

The primal infeasibility check is then

\[\left\|A^T v \right\|_{\infty} \le \epsilon_{\rm prim\_inf}, \quad u^T v_{+} + l^T v_{-} < 0.\]

Dual infeasibility

When the problem is dual infeasible, OSQP generates a vector \(s\in\mathbf{R}^{n}\) being a certificate of dual infeasibility, i.e.,

\[\begin{split}P s = 0, \quad q^T s < 0, \quad (As)_i = \begin{cases} 0 & l_i \in \mathbf{R}, u_i\in\mathbf{R} \\ \ge 0 & l_i\in\mathbf{R}, u_i=+\infty \\ \le 0 & u_i\in\mathbf{R}, l_i=-\infty. \end{cases}\end{split}\]

The dual infeasibility check is then

\[\begin{split}\| P s \|_{\infty} \le \epsilon_{\rm dual\_inf} , \quad q^T s < 0, \\ (A s)_i \begin{cases} \in \left[-\epsilon_{\rm dual\_inf}, \epsilon_{\rm dual\_inf}\right] & u_i, l_i \in \mathbf{R}\\ \ge -\epsilon_{\rm dual\_inf} &u_i = +\infty\\ \le \epsilon_{\rm dual\_inf} &l_i = -\infty.\end{cases}\end{split}\]

Polishing

Polishing is an additional algorithm step where OSQP tries to compute a high-accuracy solution. We can enable it by turning the setting polish to 1.

Polishing works by guessing the active constraints at the optimum and solving an additional linear system. If the guess is correct, OSQP returns a high accuracy solution. Otherwise OSQP returns the ADMM solution. The status of the polishing phase appears in the information status_polish.

Note that polishing requires the solution of an additional linear system and thereby, an additional factorization if the linear system solver is direct. However, the linear system is usually much smaller than the one solved during the ADMM iterations.

The chances to have a successful polishing increase if the tolerances eps_abs and eps_rel are small. However, low tolerances might require a very large number of iterations.

Implementations

The OSQP library provides a C language implementation of the above ADMM algorithm, and also interfaces to several high level languages to enable those languages to solve QPs using OSQP.

There are also several community-developed implementations of this ADMM algorithm.

Implementation	Language	Link
jaxopt.OSQP	JAX	jaxopt.github.io/stable/quadratic_programming.html#osqp
jOSQP	Java	github.com/quantego/josqp

Note that these implementations may not have the same API or features as the OSQP library, and any questions or support requests should be directed to the authors of the implementations.