ProxQP solves convex quadratic programs, which consists in minimizing a convex quadratic cost under some linear constraints. It is mathematically described as:

$$\begin{equation}\label{eq:QP}\tag{QP} \begin{aligned} \min_{x\in\mathbb{R}^{d}} & \quad \frac{1}{2}x^{T}Hx+g^{T}x \\ \text{s.t.}&\left\{ \begin{array}{ll} Ax = b, \\ Cx \leq u. \\ \end{array} \right. \end{aligned} \end{equation}\\ \text{with } H\in\mathbb{R}^{d\times d}, A\in\mathbb{R}^{n_\text{eq}\times d}, C\in\mathbb{R}^{n_\text{in}\times d}, b\in\mathbb{R}^{n_\text{eq}}, u\in\mathbb{R}^{n_\text{in}}. $$ H is a real symmetric positive semi-definite matrix. d is the problem dimension (i.e., the number of primal variables), while n_eq and n_in are the numbers of equality and inequality constraints respectively.

For linearly constrained convex optimization problems such as \eqref{eq:QP}, strong duality holds and the associated KKT conditions are necessary and sufficient for ensuring a primal-dual point (x,y,z) to be optimal (see, e.g.,Section 5.2.3} and Section 2, page 5 for more details). For \eqref{eq:QP}, the KKT system is given by the set of equations:

$$\begin{equation}\label{qp:kkt}\tag{KKT} \begin{aligned} &\left\{ \begin{array}{ll} Hx+g+A^Ty+C^Tz = 0, \\ Ax-b = 0, \\ Cx \leq u, \\ z\odot[Cx-u] = 0,\\ \end{array} \right. \end{aligned} \end{equation}$$

where the last equation involves the Hadamard product (i.e., for two vectors u and v, the Hadamard product is the vector whose ith entry is u_i v_i).

In practice, we look for a triplet (x,y,z) satisfying these optimality conditions \eqref{qp:kkt} up to a certain level of absolute accuracy (dependent of the application), leading us to the following absolute stopping criterion on the primal and dual residuals:

$$\begin{equation}\label{eq:approx_qp_sol} \begin{aligned} &\left\{ \begin{array}{ll} |Hx+g+A^Ty+C^Tz|_{\infty} \leq \epsilon_{abs}, \\ |Ax-b|_{\infty} \leq \epsilon_{abs}, \\ |[Cx-u]_+|_{\infty}\leq \epsilon_{abs}. \\ \end{array} \right. \end{aligned} \end{equation}$$

The infite norm is preferred to the L2 norm as it is independent of the problem dimensions. It is also common to consider relative convergence criteria for early-stopping, as absolute targets might not bet reached due to numerical issues. ProxQP provides it in a similar way as OSQP (for more details see, e.g., OSQP's convergence criteria or section 3.4 in the corresponding paper). Hence more generally the following stopping criterion can be used:

$$\begin{equation}\label{eq:approx_qp_sol_relative_criterion} \begin{aligned} &\left\{ \begin{array}{ll} |Hx+g+A^Ty+C^Tz|_{\infty} \leq \epsilon_{\text{abs}} + \epsilon_{\text{rel}}\max(|Hx|_{\infty},|A^Ty|_{\infty},|C^Tz|_{\infty},|g|_{\infty}), \\ |Ax-b|_{\infty} \leq \epsilon_{\text{abs}} +\epsilon_{\text{rel}}\max(|Ax|_{\infty},|b|_{\infty}), \\ |[Cx-u]_+|_{\infty}\leq \epsilon_{\text{abs}} +\epsilon_{\text{rel}}\max(|Cx|_{\infty},|u|_{\infty}). \\ \end{array} \right. \end{aligned} \end{equation}$$

It is important to note that this stopping criterion on primal and dual residuals is not enough to guarantee that the returned solution satisfies all \eqref{qp:kkt} conditions. Indeed, as the problem has affine constraints and the objective is quadratic and convex, then as soon as the primal or the dual problem is feasible, then strong duality holds (see e.g., Theorem 2 from L. El Ghaoui's lesson) and to satisfy all optimality conditions we need to add a third criterion on the duality gap $r_g$:

$$\begin{equation}\label{eq:approx_dg_sol} \begin{aligned} r_g := | x^T H x + g^T x + b^T y + u^T [z]_+ + l^T [z]_- | \leq \epsilon^{\text{gap}}_{\text{abs}} + \epsilon^{\text{gap}}_{\text{rel}} \max(|x^T H x|, |g^T x|, |b^T y|, |u^T [z]_+|, |l^T [z]_-|), \ \end{aligned} \end{equation}$$

where $[z]_+$ and $[z]_-$ stand for the projection of z onto the positive and negative orthant. ProxQP provides the check_duality_gap option to include this duality gap in the stopping criterion. Note that it is disabled by default, as other solvers don't check in general this criterion. Enable this option if you want a stronger guarantee that your solution is optimal. ProxQP will then check the same termination condition as SCS (for more details see, e.g., SCS's optimality conditions checks as well as section 7.2 in the corresponding paper). The absolute and relative thresholds $\epsilon^{\text{gap}}_{\text{abs}}, \epsilon^{\text{gap}}_{\text{rel}}$ for the duality gap can differ from those $\epsilon_{\text{abs}}, \epsilon_{\text{rel}}$ for residuals because, contrary to residuals which result from an infinite norm, the duality gap scales with the square root of the problem dimension (thus it is numerically harder to achieve a given duality gap for larger problems). A recommended choice is $\epsilon^{\text{gap}}_{\text{abs}} = \epsilon_{\text{abs}} \sqrt{\max(n, n_{\text{eq}}, n_{\text{ineq}})}$. Note finally that meeting all residual and duality-gap criteria can be difficult for ill-conditioned problems.

A single solve function for dense and sparse backends

If if you don't want to pass through ProxQP API, it is also possible to use one single solve function. We will show how to do so with examples.

You just need to call a "solve" function with in entry the model of the convex QP you want to solve. We show you below examples in C++ and python for ProxQP sparse and dense backends. Note that the sparse and dense solvers take respectivaly entries in sparse and dense formats. Note finally that the dense backend benefits from a feature enabling it to handle more efficiently box inequality constraints. We provide an example below as well of how using it (you can find more details about it in ProxQP API).

examples/cpp/solve_without_api.cpp

examples/python/solve_without_api.py

#include <iostream>
#include "proxsuite/proxqp/sparse/sparse.hpp" // get the sparse backend of ProxQP
#include "proxsuite/proxqp/dense/dense.hpp"   // get the dense backend of ProxQP
#include "proxsuite/proxqp/utils/random_qp_problems.hpp" // used for generating a random convex qp
 
using namespace proxsuite::proxqp;
using T = double;
using Mat = dense::Mat<T>;
using Vec = dense::Vec<T>;
 
int
main()
{
  isize n = 10;
  isize n_eq(n / 4);
  isize n_in(n / 4);
  T p = 0.35;            // level of sparsity
  T conditioning = 10.0; // conditioning level for H
  auto H = utils::rand::sparse_positive_definite_rand(
    n, conditioning, p); // upper triangular matrix
  Mat H_dense = Mat(H);
  H_dense.template triangularView<Eigen::Lower>() = H_dense.transpose();
  Vec g = utils::rand::vector_rand<T>(n);
  auto A = utils::rand::sparse_matrix_rand<T>(n_eq, n, p);
  Mat A_dense = Mat(A);
  auto C = utils::rand::sparse_matrix_rand<T>(n_in, n, p);
  Mat C_dense = Mat(C);
  Vec x_sol = utils::rand::vector_rand<T>(n);
  Vec b = A * x_sol;
  Vec l = C * x_sol;
  Vec u = (l.array() + 10).matrix();
 
  // Solve the problem using the sparse backend
  Results<T> results_sparse_solver =
    sparse::solve<T, isize>(H, g, A, b, C, l, u);
  std::cout << "optimal x from sparse solver: " << results_sparse_solver.x
            << std::endl;
  std::cout << "optimal y from sparse solver: " << results_sparse_solver.y
            << std::endl;
  std::cout << "optimal z from sparse solver: " << results_sparse_solver.z
            << std::endl;
  // Solve the problem using the dense backend
  Results<T> results_dense_solver =
    dense::solve<T>(H_dense, g, A_dense, b, C_dense, l, u);
 
  // print an optimal solution x,y and z
  std::cout << "optimal x from dense solver: " << results_dense_solver.x
            << std::endl;
  std::cout << "optimal y from dense solver: " << results_dense_solver.y
            << std::endl;
  std::cout << "optimal z from dense solver: " << results_dense_solver.z
            << std::endl;
 
  // Solve the problem using the dense backend using its feature for handling
  // box constraints
 
  // some trivial boxes
  Vec u_box(n);
  Vec l_box(n);
  u_box.setZero();
  l_box.setZero();
  u_box.array() += 1.E10;
  l_box.array() -= 1.E10;
  // make sure to specify at least next 9 variables of the solve function after
  // u_box to make sure the overloading work and use the specific feature for
  // handling more efficiently box constraints
  Results<T> results_dense_solver_box = dense::solve<T>(H_dense,
                                                        g,
                                                        A_dense,
                                                        b,
                                                        C_dense,
                                                        l,
                                                        u,
                                                        l_box,
                                                        u_box,
                                                        proxsuite::nullopt,
                                                        proxsuite::nullopt,
                                                        proxsuite::nullopt,
                                                        proxsuite::nullopt,
                                                        proxsuite::nullopt,
                                                        proxsuite::nullopt,
                                                        proxsuite::nullopt,
                                                        proxsuite::nullopt,
                                                        proxsuite::nullopt);
  // print an optimal solution x,y and z
  std::cout << "optimal x from dense solver: " << results_dense_solver.x
            << std::endl;
  std::cout << "optimal y from dense solver: " << results_dense_solver.y
            << std::endl;
  std::cout << "optimal z from dense solver: " << results_dense_solver.z
            << std::endl;
  // Note that the last n elements of z corresponds to the multipliers
  // associated to the box constraints
}

import proxsuite
import numpy as np
import scipy.sparse as spa
from util import generate_mixed_qp
 
 
# load a qp object using qp problem dimensions
n = 10
n_eq = 2
n_in = 2
H, g, A, b, C, u, l = generate_mixed_qp(n, True)
 
# solve the problem using the sparse backend
results = proxsuite.proxqp.sparse.solve(H, g, A, b, C, l, u)
 
# solve the problem using the dense backend
 
results2 = proxsuite.proxqp.dense.solve(
    H.toarray(), g, A.toarray(), b, C.toarray(), l, u
)
# Note finally, that the matrices are in sparse format, when using the dense backend you
# should convert them in dense format
 
# print an optimal solution
print("optimal x: {}".format(results.x))
print("optimal y: {}".format(results.y))
print("optimal z: {}".format(results.z))
 
# solve the problem using the dense backend using its feature for handling box constraints
 
l_box = -np.ones(n) * 1.0e10
u_box = np.ones(n) * 1.0e10
 
# make sure to specify l_box=l_box, u_box=u_box in order to make work the
# overloading
results_dense_solver_box = proxsuite.proxqp.dense.solve(
    H.toarray(), g, A.toarray(), b, C.toarray(), l, u, l_box=l_box, u_box=u_box
)
# print an optimal solution
print("optimal x: {}".format(results_dense_solver_box.x))
print("optimal y: {}".format(results_dense_solver_box.y))
print("optimal z: {}".format(results_dense_solver_box.z))
# Note that the last n elements of z corresponds to the multipliers associated to the box
# constraints

The results of the solve call are available in the result object, which structure is described in ProxQP API with examples section (at "The results subclass" subsection).

Different options are available for the solve function. In the table below you have the list of the different parameters that can be specified in the solve function (without the model of the problem to solve). We precise from left to right their name, their default value, and a short description of it.

Option	Default value	Description
x	Value of the EQUALITY_CONSTRAINED_INITIAL_GUESS	Warm start value for the primal variable.
y	Value of the EQUALITY_CONSTRAINED_INITIAL_GUESS	Warm start value for the dual Lagrange multiplier for equality constraints.
z	0	Warm start value for the dual Lagrange multiplier for inequality constraints.
eps_abs	1.E-5	Asbolute stopping criterion of the solver.
eps_rel	0	Relative stopping criterion of the solver.
check_duality_gap	False	If set to true, include the duality gap in absolute and relative stopping criteria.
eps_duality_gap_abs	1.E-4	Asbolute duality-gap stopping criterion of the solver.
eps_duality_gap_rel	0	Relative duality-gap stopping criterion of the solver.
mu_eq	1.E-3	Proximal step size wrt equality constraints multiplier.
mu_in	1.E-1	Proximal step size wrt inequality constraints multiplier.
rho	1.E-6	Proximal step size wrt primal variable.
VERBOSE	False	If set to true, the solver prints information at each loop.
compute_preconditioner	True	If set to true, the preconditioner will be derived.
compute_timings	False	If set to true, timings in microseconds will be computed by the solver (setup time, solving time, and run time = setup time + solving time).
max_iter	10.000	Maximal number of authorized outer iterations.
initial_guess	EQUALITY_CONSTRAINED_INITIAL_GUESS	Sets the initial guess option for initilizing x, y and z.

All other settings are set to their default values detailed in ProxQP API with examples differents subsections of "The settings subclass" section.

Note that contrary to ProxQP API, the default value of the initial guess is the "EQUALITY_CONSTRAINED_INITIAL_GUESS" option. Indeed, there is no meaning in using the "WARM_START_WITH_PREVIOUS_RESULT" option as after the solve call, it is not possible to warm start any Qp object with previous results. Hence, the only meaningful initial guess options are:

EQUALITY_CONSTRAINED_INITIAL_GUESS,
NO_INITIAL_GUESS,
or a WARM_START explicitly provided by the user.

For the latter case, It is not necessary to specify the WARM_START initial guess in the solve function. It is sufficient to just add the warm start x, y and z in the solve function, and the solver will automatically make the setting change internally.

Finally, note that in C++, if you want to change one option, the order described above in the table matter. Any intermediary option not changed must be let to the std::nullopt value. In Python, you can just specify the name of the option you want to change in any order. We give an example below.

examples/cpp/solve_without_api_and_option.cpp

examples/python/solve_without_api_and_option.py

#include <iostream>
#include <proxsuite/proxqp/sparse/sparse.hpp> // get the sparse backend of ProxQP
#include <proxsuite/proxqp/dense/dense.hpp>   // get the dense backend of ProxQP
#include <proxsuite/proxqp/utils/random_qp_problems.hpp> // used for generating a random convex qp
 
using namespace proxsuite;
using namespace proxsuite::proxqp;
using T = double;
 
int
main()
{
  isize n = 10;
  isize n_eq(n / 4);
  isize n_in(n / 4);
  T sparsity_factor(0.15);
  T strong_convexity_factor(1.e-2);
  dense::Model<T> qp_random = utils::dense_strongly_convex_qp(
    n, n_eq, n_in, sparsity_factor, strong_convexity_factor);
  // Solve the problem using the dense backend
  // and suppose you want to change the accuracy to 1.E-9 and rho initial value
  // to 1.E-7
  proxsuite::proxqp::Results<T> results =
    proxsuite::proxqp::dense::solve<T>(qp_random.H,
                                       qp_random.g,
                                       qp_random.A,
                                       qp_random.b,
                                       qp_random.C,
                                       qp_random.l,
                                       qp_random.u,
                                       nullopt,
                                       nullopt,
                                       nullopt,
                                       T(1.E-9),
                                       nullopt,
                                       nullopt,
                                       nullopt,
                                       T(1.E-7));
  // print an optimal solution x,y and z
  std::cout << "optimal x: " << results.x << std::endl;
  std::cout << "optimal y: " << results.y << std::endl;
  std::cout << "optimal z: " << results.z << std::endl;
}

import proxsuite
import numpy as np
import scipy.sparse as spa
from util import generate_mixed_qp
 
 
# load a qp object using qp problem dimensions
n = 10
n_eq = 2
n_in = 2
H, g, A, b, C, u, l = generate_mixed_qp(n)
# solve the problem using the sparse backend
# and suppose you want to change the accuracy to 1.E-9 and rho initial value to 1.E-7
results = proxsuite.proxqp.dense.solve(
    H=H, g=g, A=A, b=b, C=C, l=l, u=u, rho=1.0e-7, eps_abs=1.0e-9
)
# Note that in python the order does not matter for rho and eps_abs
# print an optimal solution
print("optimal x: {}".format(results.x))
print("optimal y: {}".format(results.y))
print("optimal z: {}".format(results.z))

Note that if some elements of your QP model are not defined (for example a QP without linear cost or inequality constraints), you can either pass a None argument, or a matrix with zero shape for specifying it. We provide an example below in cpp and python (for the dense case, it is similar with sparse backend).

examples/cpp/initializing_with_none_without_api.cpp

examples/python/initializing_with_none_without_api.py

#include <iostream>
#include "proxsuite/proxqp/dense/dense.hpp"
#include <proxsuite/proxqp/utils/random_qp_problems.hpp> // used for generating a random convex qp
 
using namespace proxsuite::proxqp;
using T = double;
int
main()
{
  dense::isize dim = 10;
  dense::isize n_eq(0);
  dense::isize n_in(0);
  T strong_convexity_factor(0.1);
  T sparsity_factor(0.15);
  // we generate a qp, so the function used from helpers.hpp is
  // in proxqp namespace. The qp is in dense eigen format and
  // you can control its sparsity ratio and strong convexity factor.
  dense::Model<T> qp_random = utils::dense_strongly_convex_qp(
    dim, n_eq, n_in, sparsity_factor, strong_convexity_factor);
 
  Results<T> results =
    dense::solve<T>(qp_random.H,
                    qp_random.g,
                    qp_random.A,
                    qp_random.b,
                    qp_random.C,
                    qp_random.l,
                    qp_random.u); // initialization with zero shape matrices
  // it is equivalent to do dense::solve<T>(qp_random.H, qp_random.g,
  // nullopt,nullopt,nullopt,nullopt,nullopt);
  //  print an optimal solution x,y and z
  std::cout << "optimal x: " << results.x << std::endl;
  std::cout << "optimal y: " << results.y << std::endl;
  std::cout << "optimal z: " << results.z << std::endl;
}

import numpy as np
import proxsuite
 
H = np.array([[65.0, -22.0, -16.0], [-22.0, 14.0, 7.0], [-16.0, 7.0, 5.0]])
g = np.array([-13.0, 15.0, 7.0])
A = None
b = None
C = None
u = None
l = None
 
results = proxsuite.proxqp.dense.solve(
    H, g, A, b, C, l, u
)  # it is equivalent to do as well proxsuite.proxqp.dense.solve(H, g)
print("optimal x: {}".format(results.x))

Finally, note that you can also you ProxQP for solving QP with non convex quadratic. For doing so, you just need to provide to the solve function an estimate of the smallest eigenvalue of the quadratic cost H. The solver environment provides an independent function for estimating the minimal eigenvalue of a dense or sparse symmetric matrix. It is named "estimate_minimal_eigen_value_of_symmetric_matrix". You can find more details in ProxQP API with examples about the different other settings that can be used for setting other related parameters (e.g., for using a Power Iteration algorithm).