PaperlayerView exact revision on arXiv

Source-linked AI summary

OSQP: An Operator Splitting Solver for Quadratic Programs

Bartolomeo Stellato, Goran Banjac, Paul Goulart, Alberto Bemporad, Stephen Boyd

arXiv:1711.08013v4math.OC

TL;DR

The paper addresses the need for a reliable general-purpose solver for convex QPs that supports high accuracy and infeasibility detection despite limitations of existing first-order methods. It develops an ADMM-based operator-splitting method with efficient repeated-solve features, and reports broad benchmark improvements alongside extensive adoption of OSQP.

  • Problem

    Existing first-order methods typically cannot detect primal or dual infeasibility, while ADMM convergence depends strongly on problem data and step-size parameters.

  • Method

    The paper develops a general-purpose ADMM-based QP solver using a novel splitting and quasi-definite linear system, supplemented by polishing, warm starting, and factorization reuse.

  • Results

    Benchmarks across finance, control, machine learning, and difficult QP test sets showed great improvements over state-of-the-art academic and commercial QP solvers.

  • Takeaways & Limitations

    OSQP provides a robust solver that returns solutions for solvable problems or primal and dual infeasibility certificates for unsolvable ones, with efficient support for parametrized QPs.

  • Takeaways & Limitations

    Solution polishing can fail when the identified active constraints are redundant, because the associated linear system may be unsolvable in degenerate solutions.

Abstract

from arXiv · show

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.

1 Introduction

Convex quadratic programs arise across engineering, finance, control, machine learning, and other fields, but existing solution methods have important limitations. The paper introduces OSQP, an ADMM-based solver designed for robustness, high accuracy, infeasibility detection, and efficient repeated solves.

  • 1.1 The problem: Convex quadratic programs support applications including portfolio optimization, model predictive control, moving horizon estimation, support vector machines, Lasso, and Huber fitting.They also arise in signal processing and as subproblems in nonconvex optimization methods.
  • 1.2 Solution methods: Active-set methods can have worst-case complexity that grows exponentially with the number of constraints.They may need to investigate all possible active sets before reaching the optimal one.
  • 1.2 Solution methods: ADMM offers inexpensive iterations and practical convergence, but first-order methods typically cannot detect primal or dual infeasibility.ADMM parameter selection also strongly affects convergence, and reliable general-purpose operator-splitting QP solvers were lacking.
  • 1.3 Our approach: OSQP uses a novel ADMM splitting with a quasi-definite linear system that remains solvable for any convex problem data.The method avoids requirements such as strict objective convexity or linear independence of the constraints.
  • 1.3 Our approach: OSQP improves solution accuracy through polishing, which solves a smaller equality-constrained QP after identifying active constraints.When active constraints are identified correctly, the resulting accuracy can equal or exceed that of interior-point methods.
  • 1.3 Our approach: OSQP supports warm starts and factorization reuse, while benchmarks over diverse QP collections report up to an order-of-magnitude computational-time improvement.The open-source C implementation is library-free and intended for embedded and large-scale applications.

2 Optimality conditions

The paper reformulates a convex QP with an auxiliary variable and characterizes solutions through primal-dual optimality conditions. Strong alternatives then provide conditions under which vectors certify primal or dual infeasibility.

  • 2 Optimality conditions: Introducing z makes the QP equivalent to minimizing the quadratic objective subject to Ax = z and z ∈ C.This separates the linear mapping from membership in the convex constraint set.
  • 2 Optimality conditions: The optimality conditions use x and z as primal variables, y as the multiplier for Ax = z, and the normal cone N_C(z).Satisfying these conditions defines a primal-dual solution.
  • 2.1 Certificates of primal and dual infeasibility: Strong alternatives state that exactly one of the primal-feasibility and certificate sets is nonempty, under an applicable constraint qualification.The certificate formulation connects infeasibility detection to support-function conditions.
  • 2.1 Certificates of primal and dual infeasibility: For a general convex QP, a vector y satisfying A^T y = 0 and S_C(y) < 0 certifies primal infeasibility.Here S_C is the support function of C, and the certificate applies when the stated constraint qualification holds.
  • 2.1 Certificates of primal and dual infeasibility: For box-constrained QPs, primal and dual infeasibility certificates reduce to explicit vector conditions involving the lower and upper bounds.The paper gives the corresponding sign and zero conditions for finite and infinite bounds.

3 Solution with ADMM

The ADMM splitting reduces each iteration to a uniquely solvable quasi-definite linear system, while convergence theory covers both solvable and infeasible quadratic programs. Cached factorization, warm starts, and simple parallelizable updates improve repeated-solve efficiency.

  • ADMM formulation: ADMM solves the split quadratic program through auxiliary variables and projection-based iterations.The splitting introduces auxiliary variables and indicator functions, with updates using vector operations and projection onto the constraint set.
  • Linear-system solution: The KKT system has a unique solution for arbitrary P and A, without requiring strong convexity or linearly independent constraints.Positive parameters in the splitting make the KKT matrix full rank for all problem data.
  • Linear-system solution: A single factorization can be reused across ADMM iterations because the KKT matrix remains unchanged when ρ is fixed.The cached factors support efficient forward and backward substitution at subsequent iterations.
  • Linear-system solution: The quasi-definite KKT matrix admits an LDLT factorization and can support division-free solution steps after storing D^-1.Its positive-definite and negative-definite diagonal blocks ensure nonzero diagonal factors.
  • Convergence and infeasibility detection: The algorithm converges to optimality conditions for solvable problems and provides infeasibility certificates otherwise.For infeasible problems, the iterates need not converge, but their difference sequence converges and yields primal or dual infeasibility certificates.
  • Termination criteria: Termination criteria stop iterations when a primal-dual solution or infeasibility certificate is obtained within prescribed tolerances.Optimality uses primal and dual residual tolerances, while separate criteria detect primal and dual infeasibility.

4 Solution polishing

Solution polishing uses an approximate primal-dual solution to infer active constraints, then solves a smaller system to recover a high-accuracy candidate. Regularization makes this system solvable under degeneracy, while iterative refinement reduces the resulting error.

  • High-accuracy QP solutions can be recovered from final ADMM iterates by solving one additional system of equations.
  • Lower- and upper-active constraints are identified from the signs of the dual variables.L contains indices with negative dual variables, while U contains indices with positive dual variables.
  • The polishing system enforces active constraints at their corresponding lower or upper bounds and uses the associated rows of A.
  • The candidate is accepted as a primal-dual solution only if it satisfies the optimality conditions.This procedure is called solution polishing, and its system is usually smaller than the full KKT system because at most n constraints are active in non-degenerate QPs.
  • Regularization is needed because the polishing system may be unsolvable when correctly identified active constraints contain redundant constraints.The regularized system is quasi-definite and therefore always solvable, but it perturbs the solution slightly.
  • Iterative refinement compensates for regularization error using repeated forward and backward solves without another matrix factorization.The iterates converge to the true solution when it exists, and the procedure is run for a fixed number of passes.

5 Preconditioning and parameter selection

The solver addresses ill-conditioning and scaling through heuristic equilibration, problem rescaling, and adaptive parameter selection. Defaults and residual-based updates are chosen from extensive numerical testing, while adaptation trades extra factorizations for faster convergence.

  • Ill-conditioned or badly scaled data can impair first-order methods, motivating preconditioning and parameter selection to speed convergence.
  • Preconditioning: Matrix equilibration heuristically rescales problem data to reduce the condition number of a symmetric matrix representing the QP.
  • Preconditioning: Cost scaling multiplies the cost function by a positive scalar to prevent excessively large dual variables.
  • Preconditioning: The scaled problem is transformed so primal and dual residuals and termination tolerances can still be evaluated in unscaled terms.
  • Preconditioning: Modified Ruiz equilibration normalizes matrix columns and applies cost scaling using the largest relevant norm.
  • Parameter selection: Extensive testing on millions of problem instances informed the algorithm’s default parameter choices.
  • Parameter selection: σ = 10^-6 regularizes the linear system, while α = 1.6 lies within the empirically favorable range [1.5, 1.8].
  • Parameter selection: Constraint-specific ρ values improve performance, and adaptive updates adjust their common scale using primal and dual residuals.

6 Parametric programs

OSQP accelerates repeated solutions of parametric programs by reusing factorizations and warm-starting from prior primal-dual solutions. The reuse strategy depends on whether only vectors or also matrix values change.

  • Parametric programs arise when QPs are solved repeatedly with varying data in control, statistics, finance, and sequential quadratic programming.
  • The reuse framework assumes fixed problem dimensions and fixed sparsity patterns for P and A.
  • When only q, l, and u vary, OSQP stores one KKT factorization and reuses it across subsequent solves.Because factorization is the most expensive computational step, this substantially reduces computation.
  • When nonzero values in P and A change without changing their locations, OSQP reuses the symbolic factorization but recomputes the numerical factorization.This yields a modest reduction in computation time.
  • Warm starting initializes primal and dual iterates from the previous OSQP solution and is particularly effective when successive QP solutions change little.The ρ estimation can also be warm-started using the previous primal-to-dual residual ratio.

7 OSQP

OSQP is an open-source C implementation designed to solve convex QPs with sparse linear algebra and configurable solver backends. Its reporting separates preprocessing, factorization, and ADMM iteration costs, including reuse cases.

  • OSQP is an open-source C package that solves convex QPs without assumptions beyond convexity.
  • The solver is accessible from multiple programming languages and modeling interfaces, including C++, Fortran, Python, Matlab, R, Julia, Ruby, Rust, CVXPY, JuMP, and YALMIP.
  • OSQP uses compressed-sparse-column matrices and provides an interface for switching among linear-system solution algorithms.
  • The package ships with the open-source QDLDL direct solver and supports dynamically loaded advanced solvers such as MKL Pardiso.
  • OSQP uses automatic ρ selection by default and sets the fixed number of iterative-refinement steps to 3.
  • Reported computation time distinguishes preprocessing, matrix factorization, and ADMM iterations, and reused factorizations report only ADMM solve time.

8 Numerical examples

OSQP was benchmarked against open-source and commercial interior-point and active-set solvers on broad, application-oriented QP test sets. It achieved strong performance across several benchmarks, while setup and polishing costs varied with problem structure and accuracy requirements.

  • Benchmark problems: OSQP was compared with ECOS, qpOASES, GUROBI, and MOSEK at low accuracy εabs = εrel = 10^-3 and high accuracy εabs = εrel = 10^-5.The benchmark used solver-specific feasibility tolerances aligned with these accuracy settings.
  • Benchmark problems: 1400 QP instances spanning 7 problem classes were generated from real or non-trivial random data, covering control, portfolio optimization, and machine learning applications.The instances ranged from n=101–10^4, m=10^2–10^5, and 10^2–10^8 nonzeros.
  • Benchmark problems: OSQP showed the best performance across the broad benchmarks, with MOSEK better at lower accuracy and GUROBI better at higher accuracy in the reported comparisons.OSQP was competitive with or faster than GUROBI for several problem classes, while ECOS was generally slower than the other interior-point solvers and qpOASES had issues with many constraints.
  • Solver statistics: 22–32% of total solution time was spent on polishing for SuiteSparse problems, while polishing increased solution time by a median of 10–20% in another benchmark class.On hard Maros-Mészáros problems, polishing rarely succeeded because the problems were badly scaled and degenerate, but its median time increase stayed below 10%.
  • SuiteSparse matrix problems: OSQP showed the best performance on SuiteSparse least-squares-derived QPs, with GUROBI slightly slower and MOSEK third.GUROBI and MOSEK had higher failure rates because their reported solutions did not satisfy the original problem’s optimality conditions.
  • Maros-Mészáros problems: GUROBI was fastest on the hard Maros-Mészáros test set, while OSQP remained competitive at both accuracy levels and MOSEK was slowest.More than 95% of cases required fewer than 5 ρ updates, despite some worst-case issues from bad scaling.
  • Warm start and factorization caching: Warm starting and factorization caching produced more than 10-fold time improvement and 8–11-fold fewer iterations for Lasso, 2.6–4-fold time improvement for MPC, and 5.8–7-fold time reduction for portfolio problems.MPC iteration counts changed little because cold-start iteration counts were already low; portfolio iterations fell by 2.9–3.6 times.

9 Conclusions

The proposed ADMM-based QP solver combines broad robustness and infeasibility detection with polishing, warm starts, cached factorizations, and an open-source implementation. Benchmarks across multiple domains showed substantial timing and failure-rate improvements over competing solvers, while OSQP gained tens of thousands of users.

  • The solver uses a novel ADMM splitting whose quasi-definite linear system is always solvable independently of problem data.The method imposes no assumptions beyond convexity and is designed as a general-purpose, robust algorithm.
  • It provides primal and dual infeasibility certificates without homogeneous self-dual embedding or additional iteration complexity.
  • Solution polishing can produce accuracy equal to or higher than interior-point methods when the active constraints are identified correctly.Polishing solves an additional small equality-constrained QP after guessing the active constraints.
  • Warm starts reduce iterations, while cached factorizations greatly improve computation time when problem matrices remain unchanged, especially for parametric QPs.
  • The open-source C implementation uses sparse linear algebra and can be configured as library-free and division-free after initial factorization for embedded applications.
  • Benchmarks across finance, control, machine learning, and test sets showed great timing and failure-rate improvements over academic and commercial QP solvers.OSQP also has tens of thousands of users in academia and large corporations.

A Problem classes

The benchmark section uses random problem classes with explicit linear equalities and inequalities represented in the standard QP constraint form. These formulations support direct benchmarking of the solver.

  • The benchmarks use random problem classes to evaluate the solver.
  • The formulations derive explicit linear equalities and inequalities for direct representation as Ax ∈ C.
  • The constraint set is represented as C = [l, u] in the form Ax ∈ C.

A.1 Random QP

The random QP instances use sparse, normally distributed problem data with m = 10n constraints. A regularization term αI with α = 10^-2 is added to P to ensure boundedness.

  • Each random QP instance has n variables and m = 10n constraints.
  • The quadratic matrix is generated as P = MM^T + αI with 15% nonzero entries in M and α = 10^-2.The regularization is added to ensure that the problem is not unbounded.
  • The matrices A and M are sparse with normally distributed nonzero entries, while q is normally distributed and bounds satisfy ui ∼ U(0, 1), li ∼ −U(0, 1).

A.2 Equality constrained QP

Equality-constrained QPs are expressed as the standard bounded-constraint form by setting l = u = b. The benchmark instances use m = floor(n/2) and normally distributed sparse data.

  • An equality-constrained QP minimizes (1/2)x^TPx + q^Tx subject to Ax = b.
  • The equality constraint is rewritten in the standard QP form by setting l = u = b.
  • The equality-constrained instances use n variables and m = floor(n/2) constraints.
  • The instances use sparse normally distributed P, A, q, and b, with P = MM^T + αI and α = 10^-2.
  • For equality-constrained QPs, the ADMM algorithm is equivalent to iterative refinement applied to the KKT system.
  • The iterative-refinement interpretation justifies using a low σ and a high ρ for equality constraints.

A.3 Optimal control

The paper formulates control, portfolio, regression, robust-fitting, and classification applications as quadratic programs and benchmarks OSQP across problem classes. These formulations use variable transformations and sparse or structured representations to support efficient solution.

  • Optimal control: Control of constrained linear time-invariant systems is formulated as a sparse QP with nx(T + 1) + nuT variables.The formulation includes polyhedral state and input constraints over horizon T.
  • Portfolio optimization: Portfolio optimization is reformulated with y = F^T x, yielding a QP whose objective Hessian is diagonal and omits FF^T.The reformulation introduces factor-exposure variables y alongside portfolio variables x.
  • Lasso: Lasso regression is converted to a QP by introducing residual variables y and auxiliary variables t for the absolute-value constraints.The resulting formulation minimizes y^T y + λ1^Tt subject to y = Ax − b and −t ≤ x ≤ t.
  • Huber fitting: Huber fitting is represented as a QP with residual decomposition variables r and s constrained to be nonnegative.The formulation minimizes u^T u + 2M1^T(r + s) subject to Ax − b − u = r − s.
  • Support vector machine: Support vector machine classification is equivalently formulated as a QP using hinge-slack variables t and the signed feature matrix diag(b)A.The constraints require t ≥ diag(b)Ax + 1 and t ≥ 0.
  • Benchmark results: OSQP is competitive or faster than GUROBI for several problem classes, while benchmark performance profiles rank OSQP best overall across the reported comparisons.MOSEK performs better at lower accuracy and GUROBI at higher accuracy in the cited benchmark results.
Loading 1711.08013v4…