ProxQP API with examples

ProxQP solves convex quadratic programs, which minimize 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 infinite 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 be 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, 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 z projection 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 this criterion in general. 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.

Finally, note that ProxQP has a specific feature for handling primal infeasibility. More precisely, if the problem appears to be primal infeasible, it will solve the closest primal feasible problem in \(\ell_2\) sense, and (x,y,z) will satisfy.

$$\begin{equation}\label{eq:approx_closest_qp_sol_rel} \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}), \\ |A^\top(Ax-b)+C^\top[Cx-u]_+|_{\infty} \leq |A^\top 1 + C^\top 1|_{\infty}*\epsilon_{\text{abs}} +\epsilon_{\text{rel}}\max(|Ax|_{\infty},|b|_{\infty}), \\ \end{array} \right. \end{aligned} \end{equation}$$

You can find more details on these subjects (and how to activate this feature with ProxQP) in the subsections describing the Settings and Results classes. You can also find more technical references with this work.

ProxQP unified API for dense and sparse backends

ProxQP algorithm is implemented in two versions specialized for dense and sparse matrices. One simple and unified API has been designed for loading dense and sparse backends. Concretely, it contains three methods:

• init : for initializing the QP model, along with some parameters,
• solve : for solving the QP problem,
• update : for updating some parameters of the QP model.

In what follows, we will make several examples to illustrate how to use this API in C++ and Python. Some subtle differences exist, nevertheless, between the dense and sparse backends, and we will point them out when needed. We will also present all solver's possible settings and show where the results are stored. We will then give some recommendations about which backend to use, considering your needs.

The API structure

When creating a Qp object in C++ or Python, it automatically contains the following sub-classes:

• model: a class storing the QP problem model which we want to solve,
• results: a class storing the main solver's results,
• settings: a class with all solver's settings,
• work: a class (not exposed in Python), with auxiliary variables the solver uses for its subroutines.

For loading ProxQP with the dense backend, it is as simple as the following code below:

examples/cpp/loading_dense_qp.cpp	examples/python/loading_dense_qp.py
#include <iostream> #include "proxsuite/proxqp/dense/dense.hpp" using namespace proxsuite::proxqp; using T = double; int main() { dense::isize dim = 10; dense::isize n_eq(dim / 4); dense::isize n_in(dim / 4); dense::QP<T> qp(dim, n_eq, n_in); }	import proxsuite n = 10 n_eq = 2 n_in = 2 qp = proxsuite.proxqp.dense.QP(n, n_eq, n_in)

The dimensions of the problem (i.e., n is the dimension of primal variable x, n_eq the number of equality constraints, and n_in the number of inequality constraints) are used for allocating the space needed for the Qp object. The dense Qp object is templated by the floating precision of the QP model (in the example above in C++ a double precision). Note that for the model to be valid, the primal dimension (i.e., n) must be strictly positive. If it is not the case, an assertion will be raised precising this issue.

The dense backend also has a specific feature for efficiently handling box inequality constraints. To benefit from it, constructors are overloaded as follows:

examples/cpp/loading_dense_qp_with_box_constraints.cpp	examples/python/loading_dense_qp_with_box_constraints.py

Furthermore, the dense version of ProxQP has two different backends with different advantages:

• PrimalDualLDLT: it factorizes a regularized version of the KKT system and benefits from great accuracy and stability. Nevertheless, if the primal dimension (i.e., the one of x) is far smaller than the dimensions of the constraints, it will be slower than PrimalLDLT backend.
• PrimalLDLT: it factorizes at the beginning the matrix $$H+\rho I+\frac{1}{\mu_{eq}} A^\top A$$ and goes on then with rank one updates. It is less accurate than PrimalDualLDLT backend, but it will be far quicker if it happens that the primal dimension is much smaller than the ones of the constraints.

The QP constructor uses, by default, an automatic choice for deciding which backend suits a priori bests user's needs. It is based on a heuristic comparing a priori computational complexity of each backend. However, if you have more insights into your needs (e.g., accuracy specifications, primal dimension is known to be far larger than the one of the constraints etc.), we encourage you to specify directly in the constructor which backend to use. It is as simple as the following:

examples/cpp/loading_dense_qp_with_different_backend_choice.cpp

examples/python/loading_dense_qp_with_different_backend_choice.py

#include <iostream> #include "proxsuite/proxqp/dense/dense.hpp" using namespace proxsuite::proxqp; using T = double; int main() { dense::isize dim = 10; dense::isize n_eq(dim / 4); dense::isize n_in(dim / 4); // we load here a QP model with // n_eq equality constraints // n_in generic type of inequality constraints // and dim box inequality constraints // we specify PrimalDualLDLT backend dense::QP<T> qp( dim, n_eq, n_in, true, proxsuite::proxqp::DenseBackend::PrimalDualLDLT); // true specifies we take into accounts box constraints // n_in are any other type of inequality constraints // Other examples // we load here a QP model with // n_eq equality constraints // O generic type of inequality constraints // and dim box inequality constraints // we specify PrimalLDLT backend dense::QP<T> qp2( dim, n_eq, 0, true, proxsuite::proxqp::DenseBackend::PrimalLDLT); // true specifies we take into accounts box constraints // we don't need to precise n_in = dim, it is taken // into account internally // we let here the solver decide with Automatic // backend choice. dense::QP<T> qp3( dim, n_eq, 0, true, proxsuite::proxqp::DenseBackend::Automatic); // Note that it is the default choice }

import proxsuite n = 10 n_eq = 2 n_in = 2 # we load here a QP model with # n_eq equality constraints # n_in generic type of inequality constraints # and dim box inequality constraints # we specify PrimalDualLDLT backend qp = proxsuite.proxqp.dense.QP( n, n_eq, n_in, True, dense_backend=proxsuite.proxqp.dense.DenseBackend.PrimalDualLDLT, ) # true specifies we take into accounts box constraints # n_in are any other type of inequality constraints # Other examples # we load here a QP model with # n_eq equality constraints # O generic type of inequality constraints # and dim box inequality constraints # we specify PrimalLDLT backend qp2 = proxsuite.proxqp.dense.QP( n, n_eq, 0, True, dense_backend=proxsuite.proxqp.dense.DenseBackend.PrimalLDLT ) # true specifies we take into accounts box constraints # we don't need to precise n_in = dim, it is taken # into account internally # We let finally the solver decide qp3 = proxsuite.proxqp.dense.QP( n, n_eq, 0, True, dense_backend=proxsuite.proxqp.dense.DenseBackend.Automatic ) # Note that it is the default choice

For loading ProxQP with the sparse backend, they are two possibilities:

• one can use as before the dimensions of the QP problem (i.e., n, n_eq and n_in)
• or one can use the sparsity structure of the matrices defining the QP problem. More precisely, if H designs the quadratic cost of the model, A the equality constraint matrix, and C the inequality constraint matrix, one can pass in entry a boolean mask of these matrices (i.e., matrices with true value when one entry is non zero) for initializing the Qp object.

examples/cpp/loading_sparse_qp.cpp

examples/python/loading_sparse_qp.py

#include <iostream> #include <proxsuite/proxqp/sparse/sparse.hpp> // get the sparse API of ProxQP #include <proxsuite/proxqp/utils/random_qp_problems.hpp> // used for generating a random convex qp using namespace proxsuite::proxqp; using T = double; int main() { // design a qp object using QP problem dimensions isize n = 10; isize n_eq(n / 4); isize n_in(n / 4); sparse::QP<T, isize> qp(n, n_eq, n_in); // assume you generate these matrices H, A and C for your QP problem T p = 0.15; // level of sparsity T conditioning = 10.0; // conditioning level for H auto H = ::proxsuite::proxqp::utils::rand::sparse_positive_definite_rand( n, conditioning, p); auto A = ::proxsuite::proxqp::utils::rand::sparse_matrix_rand<T>(n_eq, n, p); auto C = ::proxsuite::proxqp::utils::rand::sparse_matrix_rand<T>(n_in, n, p); // design a qp2 object using sparsity masks of H, A and C proxsuite::proxqp::sparse::QP<T, isize> qp2( H.cast<bool>(), A.cast<bool>(), C.cast<bool>()); }

import proxsuite # load a qp object using qp problem dimensions n = 10 n_eq = 2 n_in = 2 qp = proxsuite.proxqp.sparse.QP(n, n_eq, n_in) import numpy as np import scipy.sparse as spa from util import generate_mixed_qp # load a qp2 object using matrix masks H, g, A, b, C, u, l = generate_mixed_qp(n, True) H_ = H != 0.0 A_ = A != 0.0 C_ = C != 0.0 qp2 = proxsuite.proxqp.sparse.QP(H_, A_, C_)

The sparse Qp object is templated by the floating precision of the QP model (in the example above, a double precision) and the integer precision used for the different types of non-zero indices used (for the associated sparse matrix representation used).

The init method

Once you have defined a Qp object, the init method enables you to set up the QP problem to be solved (the example is given for the dense backend, it is similar for the sparse backend).

examples/cpp/init_dense_qp.cpp

examples/python/init_dense_qp.py

#include <iostream> #include <proxsuite/proxqp/dense/dense.hpp> // load the dense solver backend #include <proxsuite/proxqp/utils/random_qp_problems.hpp> // used for generating a random convex qp using namespace proxsuite::proxqp; using T = double; int main() { isize dim = 10; isize n_eq(dim / 4); isize n_in(dim / 4); // generate a random qp T sparsity_factor(0.15); T strong_convexity_factor(1.e-2); dense::Model<T> qp_random = utils::dense_strongly_convex_qp( dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); dense::QP<T> qp(dim, n_eq, n_in); // create the QP object qp.init(qp_random.H, qp_random.g, qp_random.A, qp_random.b, qp_random.C, qp_random.l, qp_random.u); // initialize the model }

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 qp = proxsuite.proxqp.dense.QP(n, n_eq, n_in) # generate a random QP H, g, A, b, C, u, l = generate_mixed_qp(n) # initialize the model of the problem to solve qp.init(H, g, A, b, C, l, u)

If you use the specific feature of the dense backend for handling box constraints, the init method uses simply as follows:

examples/cpp/init_dense_qp_with_box.cpp

examples/python/init_dense_qp_with_box.py

#include <iostream> #include <proxsuite/proxqp/dense/dense.hpp> // load the dense solver backend #include <proxsuite/proxqp/utils/random_qp_problems.hpp> // used for generating a random convex qp using namespace proxsuite::proxqp; using T = double; int main() { isize dim = 10; isize n_eq(dim / 4); isize n_in(dim / 4); // generate a random qp T sparsity_factor(0.15); T strong_convexity_factor(1.e-2); dense::Model<T> qp_random = utils::dense_strongly_convex_qp( dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); // specify some trivial box constraints dense::Vec<T> u_box(dim); dense::Vec<T> l_box(dim); u_box.setZero(); l_box.setZero(); u_box.array() += 1.E10; l_box.array() -= 1.E10; dense::QP<T> qp(dim, n_eq, n_in, true); // create the QP object // and specify with true you take into account box constraints // in the model // you can check that there are constraints with method is_box_constrained std::cout << "the qp is box constrained : " << qp.is_box_constrained() << std::endl; // init the model qp.init(qp_random.H, qp_random.g, qp_random.A, qp_random.b, qp_random.C, qp_random.l, qp_random.u, l_box, u_box); // initialize the model }

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 qp = proxsuite.proxqp.dense.QP(n, n_eq, n_in, True) # you can check that there are constraints with method is_box_constrained print(f"the qp is box constrained: {qp.is_box_constrained()}") # generate a random QP H, g, A, b, C, u, l = generate_mixed_qp(n) l_box = -np.ones(n) * 1.0e10 u_box = np.ones(n) * 1.0e10 # initialize the model of the problem to solve qp.init(H, g, A, b, C, l, u, l_box, u_box)

Note that with its dense backend, ProxQP solver manipulates matrices in dense representations (in the same spirit, the solver with sparse backend manipulates entries in sparse format). In the example above, the matrices are originally in sparse format and eventually converted into dense format. 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 shapes 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.cpp

examples/python/initializing_with_none.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); dense::QP<T> qp(dim, n_eq, n_in); 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); qp.init(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 qp.init(qp_random.H, qp_random.g, // nullopt,nullopt,nullopt,nullopt,nullopt); qp.solve(); // print an optimal solution x,y and z std::cout << "optimal x: " << qp.results.x << std::endl; std::cout << "optimal y: " << qp.results.y << std::endl; std::cout << "optimal z: " << qp.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 qp = proxsuite.proxqp.dense.QP(3, 0, 0) qp.init(H, g, A, b, C, l, u) # it is equivalent to do as well qp.init(H, g) qp.solve() print("optimal x: {}".format(qp.results.x))

With the init method, you can also setting-up on the same time some other parameters in the following order:

• compute_preconditioner: a boolean parameter for executing or not the preconditioner. The preconditioner is an algorithm used (for the moment, we use Ruiz equilibrator) for reducing the ill-conditioning of the QP problem, and hence speeding-up the solver and increasing its accuracy. It consists mostly of a heuristic involving linear scalings. Note that for very ill-conditioned QP problem, when one asks for a very accurate solution, the unscaling procedure can become less precise (we provide some remarks about this subject in section 6.D of the following paper). By default, its value is set to true.
• rho: the proximal step size wrt primal variable. Reducing its value speed-ups convergence wrt the primal variable (but increases as well the ill-conditioning of sub-problems to solve). The minimal value it can take is 1.e-7. By default, its value is set to 1.e-6.
• mu_eq: the proximal step size wrt equality constrained multiplier. Reducing its value speed-ups convergence wrt equality-constrained variables (but increases as well the ill-conditioning of sub-problems to solve). The minimal value it can take is 1.e-9. By default, its value is set to 1.e-3.
• mu_in: the proximal step size wrt inequality constrained multiplier. Reducing its value speed-ups convergence wrt inequality-constrained variable (but increases as well the ill-conditioning of sub-problems to solve). The minimal value it can take is 1.e-9. By default, its value is set to 1.e-1.

We provide below one example in C++ and Python.

examples/cpp/init_dense_qp_with_other_options.cpp

examples/python/init_dense_qp_with_other_options.py

#include <iostream> #include <proxsuite/proxqp/dense/dense.hpp> // load the dense solver backend #include <proxsuite/proxqp/utils/random_qp_problems.hpp> // used for generating a random convex qp using namespace proxsuite::proxqp; using T = double; int main() { isize dim = 10; isize n_eq(dim / 4); isize n_in(dim / 4); // generate a random qp T sparsity_factor(0.15); T strong_convexity_factor(1.e-2); dense::Model<T> qp_random = utils::dense_strongly_convex_qp( dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); dense::QP<T> qp( dim, n_eq, n_in); // create the QP // initialize the model, along with another rho parameter qp.init(qp_random.H, qp_random.g, qp_random.A, qp_random.b, qp_random.C, qp_random.l, qp_random.u, true, /*rho*/ 1.e-7); // in c++ you must follow the order speficied in the API for the parameters // if you don't want to change one parameter (here compute_preconditioner), // just let it be nullopt }

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 qp = proxsuite.proxqp.dense.QP(n, n_eq, n_in) # generate a random QP H, g, A, b, C, u, l = generate_mixed_qp(n) # initialize the model of the problem to solve with another rho parameter qp.init(H, g, A, b, C, l, u, rho=1.0e-7)

Furthermore, some settings must be defined before the init method takes effect. For example, if you want the solver to compute the runtime properly (the sum of the setup time and the solving time), you must set this option before the init method (which is part of the setup time). We provide below an example.

examples/cpp/init_dense_qp_with_timings.cpp

examples/python/init_dense_qp_with_timings.py

#include <iostream> #include <proxsuite/proxqp/dense/dense.hpp> // load the dense solver backend #include <proxsuite/proxqp/utils/random_qp_problems.hpp> // used for generating a random convex qp using namespace proxsuite::proxqp; using T = double; int main() { isize dim = 10; isize n_eq(dim / 4); isize n_in(dim / 4); // generate a random qp T sparsity_factor(0.15); T strong_convexity_factor(1.e-2); dense::Model<T> qp_random = utils::dense_strongly_convex_qp( dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); dense::QP<T> qp(dim, n_eq, n_in); // create the QP object qp.settings.compute_timings = true; // compute all timings qp.init(qp_random.H, qp_random.g, qp_random.A, qp_random.b, qp_random.C, qp_random.l, qp_random.u); // initialize the model }

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 qp = proxsuite.proxqp.dense.QP(n, n_eq, n_in) # generate a random QP H, g, A, b, C, u, l = generate_mixed_qp(n) # initialize the model of the problem to solve qp.settings.compute_timings # compute all timings qp.init(H, g, A, b, C, l, u)

The solve method

Once you have defined a Qp object and initialized it with a model, the solve method enables you solving the QP problem. The method is overloaded with two modes considering whether you provide or not a warm start to the method. We give below two examples (for the dense backend, with the sparse one it is similar).

examples/cpp/solve_dense_qp.cpp

examples/python/solve_dense_qp.py

#include <iostream> #include <proxsuite/proxqp/dense/dense.hpp> // load the dense solver backend #include <proxsuite/proxqp/utils/random_qp_problems.hpp> // used for generating a random convex qp using namespace proxsuite::proxqp; using T = double; int main() { isize dim = 10; isize n_eq(dim / 4); isize n_in(dim / 4); // generate a random qp T sparsity_factor(0.15); T strong_convexity_factor(1.e-2); dense::Model<T> qp_random = utils::dense_strongly_convex_qp( dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); dense::QP<T> qp(dim, n_eq, n_in); // create the QP object qp.init(qp_random.H, qp_random.g, qp_random.A, qp_random.b, qp_random.C, qp_random.l, qp_random.u); // initialize the model qp.solve(); // solve the problem without warm start auto x_wm = utils::rand::vector_rand<T>(dim); auto y_wm = utils::rand::vector_rand<T>(n_eq); auto z_wm = utils::rand::vector_rand<T>(n_in); qp.solve(x_wm, y_wm, z_wm); // if you have a warm start, put it here // print an optimal solution x,y and z std::cout << "optimal x: " << qp.results.x << std::endl; std::cout << "optimal y: " << qp.results.y << std::endl; std::cout << "optimal z: " << qp.results.z << std::endl; // Another example if you have box constraints (for the dense backend only for // the moment) dense::QP<T> qp2(dim, n_eq, n_in, true); // create the QP object // some trivial boxes dense::Vec<T> u_box(dim); dense::Vec<T> l_box(dim); u_box.setZero(); l_box.setZero(); u_box.array() += 1.E10; l_box.array() -= 1.E10; qp2.init(qp_random.H, qp_random.g, qp_random.A, qp_random.b, qp_random.C, qp_random.l, qp_random.u, l_box, u_box); // initialize the model with the boxes qp2.solve(); // solve the problem // An important note regarding the inequality multipliers // auto z_ineq = qp.results.z.head(n_in); contains the multiplier associated // to qp_random.C z_box = qp.results.z.tail(dim); the last dim elements // correspond to multiplier 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 qp = proxsuite.proxqp.dense.QP(n, n_eq, n_in) # generate a random QP H, g, A, b, C, u, l = generate_mixed_qp(n) # initialize the model of the problem to solve qp.init(H, g, A, b, C, l, u) # solve without warm start qp.solve() # solve with a warm start, for ex random one qp.solve(np.random.randn(n), np.random.randn(n_eq), np.random.randn(n_in)) # print an optimal solution print("optimal x: {}".format(qp.results.x)) print("optimal y: {}".format(qp.results.y)) print("optimal z: {}".format(qp.results.z)) # Another example if you have box constraints (for the dense backend only for the moment) qp2 = proxsuite.proxqp.dense.QP(n, n_eq, n_in, True) l_box = -np.ones(n) * 1.0e10 u_box = np.ones(n) * 1.0e10 qp2.init(H, g, A, b, C, l, u, l_box, u_box) qp2.solve() # An important note regarding the inequality multipliers z_ineq = qp.results.z[:n_in] # contains the multiplier associated to qp_random.C z_box = qp.results.z[ n_in: ] # the last dim elements correspond to multiplier associated to # the box constraints

Before ending this section, we will talk about how to activate some other settings before launching the solve method. To do so, you only need to define your desired settings (for example, the stopping criterion accuracy threshold eps_abs, or the verbose option) after initializing the Qp object. They will then be taken into account only if there are set before the solve method (otherwise, they will be taken into account when the next solve or update method is called). The full description of all the settings is provided in a dedicated section below. Here we just give an example to illustrate the mentioned notion above.

examples/cpp/solve_dense_qp_with_setting.cpp

examples/python/solve_dense_qp_with_setting.py

#include <iostream> #include <proxsuite/proxqp/dense/dense.hpp> // load the dense solver backend #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(dim / 4); dense::isize n_in(dim / 4); // generate a random qp T sparsity_factor(0.15); T strong_convexity_factor(1.e-2); dense::Model<T> qp_random = utils::dense_strongly_convex_qp( dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); dense::QP<T> qp(dim, n_eq, n_in); // create the QP object qp.init(qp_random.H, qp_random.g, qp_random.A, qp_random.b, qp_random.C, qp_random.l, qp_random.u); // initialize the model qp.settings.eps_abs = T(1.E-9); // set accuracy threshold to 1.e-9 qp.settings.verbose = true; // print some intermediary results qp.solve(); // solve the problem with previous settings // print an optimal solution x,y and z std::cout << "optimal x: " << qp.results.x << std::endl; std::cout << "optimal y: " << qp.results.y << std::endl; std::cout << "optimal z: " << qp.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 qp = proxsuite.proxqp.dense.QP(n, n_eq, n_in) # generate a random QP H, g, A, b, C, u, l = generate_mixed_qp(n) # initialize the model of the problem to solve qp.init(H, g, A, b, C, l, u) qp.settings.eps_abs = 1.0e-9 qp.settings.verbose = True # solve with previous settings qp.solve() # print an optimal solution print("optimal x: {}".format(qp.results.x)) print("optimal y: {}".format(qp.results.y)) print("optimal z: {}".format(qp.results.z))

The update method

The update method is used to update the model or a parameter of the problem, as for the init method. We provide below an example for the dense case.

examples/cpp/update_dense_qp.cpp

examples/python/update_dense_qp.py

#include <iostream> #include <proxsuite/proxqp/dense/dense.hpp> // load the dense solver backend #include <proxsuite/proxqp/utils/random_qp_problems.hpp> // used for generating a random convex qp using namespace proxsuite::proxqp; using T = double; int main() { isize dim = 10; isize n_eq(dim / 4); isize n_in(dim / 4); // generate a random qp T sparsity_factor(0.15); T strong_convexity_factor(1.e-2); dense::Model<T> qp_random = utils::dense_strongly_convex_qp( dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); dense::QP<T> qp(dim, n_eq, n_in); // create the QP object qp.init(qp_random.H, qp_random.g, qp_random.A, qp_random.b, qp_random.C, qp_random.l, qp_random.u); // initialize the model qp.solve(); // solve the problem // a new qp problem dense::Model<T> qp2 = utils::dense_strongly_convex_qp( dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); // re update the model qp.update(qp2.H, qp2.g, qp2.A, qp2.b, qp2.C, qp2.l, qp2.u); // solve it qp.solve(); // print an optimal solution x,y and z std::cout << "optimal x: " << qp.results.x << std::endl; std::cout << "optimal y: " << qp.results.y << std::endl; std::cout << "optimal z: " << qp.results.z << std::endl; // if you have boxes (dense backend only) you proceed the same way dense::QP<T> qp_box(dim, n_eq, n_in, true); // create the QP object dense::Vec<T> u_box(dim); dense::Vec<T> l_box(dim); u_box.setZero(); l_box.setZero(); u_box.array() += 1.E10; l_box.array() -= 1.E10; qp_box.init(qp_random.H, qp_random.g, qp_random.A, qp_random.b, qp_random.C, qp_random.l, qp_random.u, l_box, u_box); qp_box.solve(); u_box.array() += 1.E1; l_box.array() -= 1.E1; qp_box.update(qp2.H, qp2.g, qp2.A, qp2.b, qp2.C, qp2.l, qp2.u, l_box, u_box); qp_box.solve(); }

import proxsuite import numpy as np from util import generate_mixed_qp # load a qp object using qp problem dimensions n = 10 n_eq = 2 n_in = 2 qp = proxsuite.proxqp.dense.QP(n, n_eq, n_in) # generate a random QP H, g, A, b, C, u, l = generate_mixed_qp(n) # initialize the model of the problem to solve qp.init(H, g, A, b, C, l, u) # solve without warm start qp.solve() # create a new problem and update qp H_new, g_new, A_new, b_new, C_new, u_new, l_new = generate_mixed_qp(n, seed=2) qp.update(H_new, g_new, A_new, b_new, C_new, l_new, u_new) # solve it qp.solve() # print an optimal solution print("optimal x: {}".format(qp.results.x)) print("optimal y: {}".format(qp.results.y)) print("optimal z: {}".format(qp.results.z)) # if you have boxes (dense backend only) you proceed the same way qp2 = proxsuite.proxqp.dense.QP(n, n_eq, n_in, True) l_box = -np.ones(n) * 1.0e10 u_box = np.ones(n) * 1.0e10 qp2.init(H, g, A, b, C, l, u, l_box, u_box) qp2.solve() l_box += 1.0e1 u_box -= 1.0e1 qp2.update(H_new, g_new, A_new, b_new, C_new, l_new, u_new, l_box, u_box) qp2.solve()

Contrary to the init method, the compute_preconditioner boolean parameter becomes update_preconditioner, which enables you to keep the previous preconditioner (if set to false) to equilibrate the new updated problem, or to re-compute the preconditioner with the new values of the problem (if set to true). By default, the update_preconditioner parameter is set to false.

The major difference between the dense and sparse API is that in the sparse case only if you change the matrices of the model, the update will take effect only if the matrices have the same sparsity structure (i.e., the non-zero values are located at the same place). Hence, if the matrices have a different sparsity structure, you must create a new Qp object to solve the new problem. We provide an example below.

examples/cpp/update_sparse_qp.cpp

examples/python/update_sparse_qp.py

#include <iostream> #include <proxsuite/proxqp/sparse/sparse.hpp> // get the sparse API 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 p = 0.15; // level of sparsity T conditioning = 10.0; // conditioning level for H auto H = ::proxsuite::proxqp::utils::rand::sparse_positive_definite_rand( n, conditioning, p); auto g = ::proxsuite::proxqp::utils::rand::vector_rand<T>(n); auto A = ::proxsuite::proxqp::utils::rand::sparse_matrix_rand<T>(n_eq, n, p); auto C = ::proxsuite::proxqp::utils::rand::sparse_matrix_rand<T>(n_in, n, p); auto x_sol = ::proxsuite::proxqp::utils::rand::vector_rand<T>(n); auto b = A * x_sol; auto l = C * x_sol; auto u = (l.array() + 10).matrix().eval(); // design a qp2 object using sparsity masks of H, A and C proxsuite::proxqp::sparse::QP<T, int> qp( H.cast<bool>(), A.cast<bool>(), C.cast<bool>()); qp.init(H, g, A, b, C, l, u); qp.solve(); // update H auto H_new = 2 * H; // keep the same sparsity structure qp.update(H_new, nullopt, nullopt, nullopt, nullopt, nullopt, nullopt); // update H with H_new, it will work qp.solve(); // generate H2 with another sparsity structure auto H2 = ::proxsuite::proxqp::utils::rand::sparse_positive_definite_rand( n, conditioning, p); qp.update(H2, nullopt, nullopt, nullopt, nullopt, nullopt, nullopt); // nothing will happen // if only a vector changes, then the update takes effect auto g_new = ::proxsuite::proxqp::utils::rand::vector_rand<T>(n); qp.update(nullopt, g, nullopt, nullopt, nullopt, nullopt, nullopt); qp.solve(); // it solves the problem with another vector // to solve the problem with H2 matrix create a new qp object proxsuite::proxqp::sparse::QP<T, isize> qp2( H2.cast<bool>(), A.cast<bool>(), C.cast<bool>()); qp2.init(H2, g_new, A, b, C, l, u); qp2.solve(); // it will solve the new problem // print an optimal solution x,y and z std::cout << "optimal x: " << qp2.results.x << std::endl; std::cout << "optimal y: " << qp2.results.y << std::endl; std::cout << "optimal z: " << qp2.results.z << std::endl; }

import proxsuite import numpy as np from util import generate_mixed_qp # load a qp object using qp problem dimensions n = 10 n_eq = 2 n_in = 2 qp = proxsuite.proxqp.sparse.QP(n, n_eq, n_in) # generate a random QP H, g, A, b, C, u, l = generate_mixed_qp(n, True) # initialize the model of the problem to solve qp.init(H, g, A, b, C, l, u) qp.solve() H_new = 2 * H # keep the same sparsity structure qp.update(H_new) # update H with H_new, it will work qp.solve() # generate a QP with another sparsity structure # create a new problem and update qp H2, g_new, A_new, b_new, C_new, u_new, l_new = generate_mixed_qp(n, True) qp.update(H=H2) # nothing will happen qp.update(g=g_new) # if only a vector changes, then the update takes effect qp.solve() # it solves the problem with the QP H,g_new,A,b,C,u,l # to solve the problem with H2 matrix create a new qp object in the sparse case qp2 = proxsuite.proxqp.sparse.QP(n, n_eq, n_in) qp2.init(H2, g_new, A, b, C, l, u) qp2.solve() # it will solve the new problem # print an optimal solution print("optimal x: {}".format(qp.results.x)) print("optimal y: {}".format(qp.results.y)) print("optimal z: {}".format(qp.results.z))

Finally, if you want to change your initial guess option when updating the problem, you must change it in the setting before the update takes effect for the next solve (otherwise, it will keep the previous one set). It is important, especially for the WARM_START_WITH_PREVIOUS_RESULT initial guess option (set by default in the solver). Indeed, in this case, if no matrix is updated, the workspace keeps the previous factorization in the update method, which adds considerable speed up for the next solving. We provide below an example in the dense case.

examples/cpp/update_dense_qp_ws_previous_result.cpp

examples/python/update_dense_qp_ws_previous_result.py

#include <iostream> #include <proxsuite/proxqp/dense/dense.hpp> // load the dense solver backend #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() { dense::isize dim = 10; dense::isize n_eq(dim / 4); dense::isize n_in(dim / 4); // generate a random qp T sparsity_factor(0.15); T strong_convexity_factor(1.e-2); dense::Model<T> qp_random = utils::dense_strongly_convex_qp( dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); dense::QP<T> qp(dim, n_eq, n_in); // create the QP object qp.init(qp_random.H, qp_random.g, qp_random.A, qp_random.b, qp_random.C, qp_random.l, qp_random.u); // initialize the model qp.solve(); // solve the problem // re update the linear cost taking previous result qp.settings.initial_guess = InitialGuessStatus::WARM_START_WITH_PREVIOUS_RESULT; // it takes effect at the update because it is set before // (the workspace is not erased at the update method, hence // the previous factorization is kept) // a new linear cost slightly modified auto g = qp_random.g * 0.95; qp.update(nullopt, g, nullopt, nullopt, nullopt, nullopt, nullopt); qp.solve(); // print an optimal solution x,y and z std::cout << "optimal x: " << qp.results.x << std::endl; std::cout << "optimal y: " << qp.results.y << std::endl; std::cout << "optimal z: " << qp.results.z << std::endl; }

import proxsuite import numpy as np from util import generate_mixed_qp # load a qp object using qp problem dimensions n = 10 n_eq = 2 n_in = 2 qp = proxsuite.proxqp.dense.QP(n, n_eq, n_in) # generate a random QP H, g, A, b, C, u, l = generate_mixed_qp(n) # initialize the model of the problem to solve qp.init(H, g, A, b, C, l, u) # solve without warm start qp.solve() # create a new problem and update qp g_new = 0.95 * g # slightly different g_new qp.settings.initial_guess = ( proxsuite.proxqp.InitialGuess.WARM_START_WITH_PREVIOUS_RESULT ) qp.update(g=g_new) # solve it qp.solve() # print an optimal solution print("optimal x: {}".format(qp.results.x)) print("optimal y: {}".format(qp.results.y)) print("optimal z: {}".format(qp.results.z))

Note that one subsection is dedicated to the different initial guess options below.

The settings subclass

In this section, you will find first the solver's settings and then a subsection detailing the different initial guess options.

The solver's settings

In this table, you have the three columns from left to right: the name of the setting, its default value, and then a short description of it.

Setting	Default value	Description
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.
VERBOSE	False	If set to true, the solver prints information at each loop.
default_rho	1.E-6	Default rho parameter of result class (i.e., for each initial guess, except WARM_START_WITH_PREVIOUS_RESULT, after a new solve or update, the solver initializes rho to this value).
default_mu_eq	1.E-3	Default mu_eq parameter of result class (i.e., for each initial guess, except WARM_START_WITH_PREVIOUS_RESULT, after a new solve or update, the solver initializes mu_eq to this value).
default_mu_in	1.E-1	Default mu_in parameter of result class (i.e., for each initial guess, except WARM_START_WITH_PREVIOUS_RESULT, after a new solve or update, the solver initializes mu_in to this value).
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.
max_iter_in	1500	Maximal number of authorized inner iterations.
initial_guess	EQUALITY_CONSTRAINED_INITIAL_GUESS	Sets the initial guess option for initilizing x, y and z.
mu_min_eq	1.E-9	Minimal authorized value for mu_eq.
mu_min_in	1.E-8	Minimal authorized value for mu_in.
mu_update_factor	0.1	Factor used for updating mu_eq and mu_in.
eps_primal_inf	1.E-4	Threshold under which primal infeasibility is detected.
eps_dual_inf	1.E-4	Threshold under which dual infeasibility is detected.
update_preconditioner	False	If set to true, the preconditioner will be re-derived with the update method, otherwise it uses previous one.
compute_preconditioner	True	If set to true, the preconditioner will be derived with the init method.
alpha_bcl	0.1	alpha parameter of the BCL algorithm.
beta_bcl	0.9	beta parameter of the BCL algorithm.
refactor_dual_feasibility_threshold	1.E-2	Threshold above which refactorization is performed to change rho parameter.
refactor_rho_threshold	1.E-7	New rho parameter used if the refactor_dual_feasibility_threshold condition has been satisfied.
cold_reset_mu_eq	1./1.1	Value used for cold restarting mu_eq.
cold_reset_mu_in	1./1.1	Value used for cold restarting mu_in.
nb_iterative_refinement	10	Maximal number of iterative refinements.
eps_refact	1.E-6	Threshold value below which the Cholesky factorization is refactorized factorization in the iterative refinement loop.
safe_guard	1.E4	Safeguard parameter ensuring global convergence of the scheme. More precisely, if the total number of iterations is superior to safe_guard, the BCL scheme accepts always the multipliers (hence the scheme is a pure proximal point algorithm).
preconditioner_max_iter	10	Maximal number of authorized iterations for the preconditioner.
preconditioner_accuracy	1.E-3	Accuracy level of the preconditioner.
HessianType	Dense	Defines the type of problem solved (Dense, Zero, or Diagonal). In case the Zero or Diagonal option is used, the solver exploits the Hessian structure to evaluate the Cholesky factorization efficiently.
primal_infeasibility_solving	False	If set to true, it solves the closest primal feasible problem if primal infeasibility is detected.
nb_power_iteration	1000	Number of power iteration iteration used by default for estimating H lowest eigenvalue.
power_iteration_accuracy	1.E-6	If set to true, it solves the closest primal feasible problem if primal infeasibility is detected.
primal_infeasibility_solving	False	Accuracy target of the power iteration algorithm for estimating the lowest eigenvalue of H.
estimate_method_option	NoRegularization	Option for estimating the minimal eigen value of H and regularizing default_rho default_rho=rho_regularization_scaling*abs(default_H_eigenvalue_estimate). This option can be used for solving non convex QPs.
default_H_eigenvalue_estimate	0.	Default estimate of the minimal eigen value of H.
rho_regularization_scaling	1.5	Scaling for regularizing default_rho according to the minimal eigen value of H.

The different initial guesses

The solver has five different possible initial guesses for warm starting or not the initial iterate values:

• NO_INITIAL_GUESS,
• EQUALITY_CONSTRAINED_INITIAL_GUESS,
• WARM_START_WITH_PREVIOUS_RESULT,
• WARM_START,
• COLD_START_WITH_PREVIOUS_RESULT.

The different options will be commented below in the introduced order above.

No initial guess

If set to this option, the solver will start with no initial guess, which means that the initial values of x, y, and z are the 0 vector.

The different options for estimating H minimal Eigenvalue

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". In the sparse case, it uses a power iteration algorithm (with two options: the maximal number of iterations and the accuracy target for the estimate). In the dense case, we provide two options within the struct EigenValueEstimateMethodOption:

• PowerIteration: a power iteration algorithm will be used for estimating H minimal eigenvalue,
• ExactMethod: in this case, an exact method from EigenSolver is used to provide an estimate.

Estimating minimal eigenvalue is particularly usefull for solving QP with non convex quadratics. Indeed, if default_rho is set to a value strictly higher than the minimal eigenvalue of H, then ProxQP is guaranteed for find a local minimum to the problem since it relies on a Proximal Method of Multipliers (for more detail for example this work providing convergence proof of this property).

More precisely, ProxQP API enables the user to provide for the init or update methods estimate of the minimal eigenvalue of H (i.e., manual_minimal_H_eigenvalue). It the values are not empty, then the values of primal proximal step size rho will be updated according to: rho = rho + abs(manual_minimal_H_eigenvalue). It guarantees that the proximal step-size is larger than the minimal eigenvalue of H and hence to converging towards a local minimum of the QP. We provide below examples in C++ and python for using this feature appropriately with the dense backend (it is similar with the sparse one)

examples/cpp/estimate_nonconvex_eigenvalue.cpp

examples/python/estimate_nonconvex_eigenvalue.py

#include <iostream> #include <proxsuite/proxqp/dense/dense.hpp> // load the dense solver backend #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() { dense::isize dim = 10; dense::isize n_eq(dim / 4); dense::isize n_in(dim / 4); // generate a random qp T sparsity_factor(0.15); T strong_convexity_factor(1.e-2); dense::Model<T> qp_random = utils::dense_strongly_convex_qp( dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); // make the QP nonconvex dense::Vec<T> diag(dim); diag.setOnes(); qp_random.H.diagonal().array() -= 2. * diag.array(); // add some nonpositive values dense matrix Eigen::SelfAdjointEigenSolver<dense::Mat<T>> es(qp_random.H, Eigen::EigenvaluesOnly); T minimal_eigenvalue = T(es.eigenvalues().minCoeff()); // choose scaling for regularizing default_rho accordingly dense::QP<T> qp(dim, n_eq, n_in); // create the QP object // choose the option for estimating this eigenvalue T estimate_minimal_eigen_value = dense::estimate_minimal_eigen_value_of_symmetric_matrix( qp_random.H, EigenValueEstimateMethodOption::ExactMethod, 1.E-6, 10000); bool compute_preconditioner = false; // input the estimate for making rho appropriate qp.init(qp_random.H, qp_random.g, qp_random.A, qp_random.b, qp_random.C, qp_random.l, qp_random.u, compute_preconditioner, nullopt, nullopt, nullopt, estimate_minimal_eigen_value); // print the estimates std::cout << "ProxQP estimate " << qp.results.info.minimal_H_eigenvalue_estimate << std::endl; std::cout << "minimal_eigenvalue " << minimal_eigenvalue << std::endl; std::cout << "default_rho " << qp.settings.default_rho << 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 qp = proxsuite.proxqp.dense.QP(n, n_eq, n_in) # generate a random non convex QP H, g, A, b, C, u, l = generate_mixed_qp(n) # initialize the model of the problem to solve estimate_minimal_eigen_value = ( proxsuite.proxqp.dense.estimate_minimal_eigen_value_of_symmetric_matrix( H, proxsuite.proxqp.EigenValueEstimateMethodOption.ExactMethod ) ) qp.init(H, g, A, b, C, l, u, manual_minimal_H_eigenvalue=estimate_minimal_eigen_value) vals, _ = spa.linalg.eigs(H, which="SR") min_eigenvalue = float(np.min(vals)) # print the estimates print(f"{min_eigenvalue=}") print(f"{estimate_minimal_eigen_value=}") print(f"{qp.results.info.minimal_H_eigenvalue_estimate=}")

No initial guess

If set to this option, the solver will start with no initial guess, which means that the initial values of x, y, and z are the 0 vector.

Equality constrained initial guess

If set to this option, the solver will solve at the beginning the following system for warm starting x and y. $$\begin{bmatrix} H+\rho I & A^T \\ A & -\mu_{eq} I \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} -g \\ b \end{bmatrix}$$ z stays to 0. In general, this option warm starts well equality constrained QP.

Warm start with the previous result

If set to this option, the solver will warm start x, y, and z with the values of the previous problem solved and it will keep all the last parameters of the solver (i.e., proximal step sizes, for example, and the full workspace with the Cholesky factorization etc.). Hence, if the new problem to solve is the same as the previous one, the problem is warm-started at the solution (and zero iteration will be executed).

This option was initially thought to be used in optimal control-like problems when the next problem to be solved is close to the previous one. Indeed, if the problem changes only slightly, it is reasonable to warm start the new problem with the value of the previous one for speeding the whole runtime.

Note, however, that if your update involves new matrices or you decide to change parameters involved in the Cholesky factorization (i.e., the proximal step sizes), then for consistency, the solver will automatically refactorize the Cholesky with these updates (and taking into account the last values of x, y, and z for the active set).

Finally, note that this option is set by default in the solver. At the first solve, as there is no previous results, x, y and z are warm started with the 0 vector value.

Warm start

If set to this option, the solver expects then a warm start at the solve method.

Note, that it is unnecessary to set this option through the command below (for example, in C++) before the update or solve method call.

Qp.settings.initial_guess = proxsuite::qp::InitialGuessStatus::WARM_START;

It is sufficient to just add the warm start in the solve method, and the solver will automatically make the setting change internally.

Cold start with previous result

If set to this option, the solver will warm start x, y and z with the values of the previous problem solved. Contrary to the WARM_START_WITH_PREVIOUS_RESULT option, all other parameters of the solver (i.e., proximal step sizes for example, and the full workspace with the ldlt factorization etc.) are re-set to their default values (hence a factorization is reperformed taking into account of z warm start for the active set, but with default values of proximal step sizes).

This option has also been thought initially for being used in optimal control-like problems when the next problem to be solved is close to the previous one. Indeed, if the problem changes only slightly, it is reasonable to warm start the new problem with the value of the previous one for speeding the whole runtime.

Note finally that at the first solve, as there are no previous results, x, y, and z are warm started with the 0 vector value.

The results subclass

The result subclass is composed of the following:

• x: a primal solution,
• y: a Lagrange optimal multiplier for equality constraints,
• z: a Lagrange optimal multiplier for inequality constraints,
• se: the optimal shift in \(\ell_2\) with respect to equality constraints,
• si: the optimal shift in \(\ell_2\) with respect to inequality constraints,
• info: a subclass which containts some information about the solver's execution.

If the solver has solved the problem, the triplet (x,y,z) satisfies:

$$\begin{equation}\label{eq:approx_qp_sol_rel} \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}$$ accordingly with the parameters eps_abs and eps_rel chosen by the user.

If the problem is primal infeasible and you have enabled the solver to solve the closest feasible problem, then (x,y,z) will satisfy. $$\begin{equation}\label{eq:approx_closest_qp_sol_rel_bis} \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}), \\ |A^\top(Ax-b)+C^\top[Cx-u]_+|_{\infty} \leq |A^\top 1 + C^\top 1|_{\infty}*\epsilon_{\text{abs}} +\epsilon_{\text{rel}}\max(|Ax|_{\infty},|b|_{\infty}), \\ \end{array} \right. \end{aligned} \end{equation}$$

(se, si) stands in this context for the optimal shifts in \(\ell_2\) sense which enables recovering a primal feasible problem. More precisely, they are derived such that

\begin{equation}\label{eq:QP_primal_feasible}\tag{QP_feas} \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+se, \\ Cx \leq u+si, \\ \end{array} \right. \end{aligned} \end{equation}\\ defines a primal feasible problem.

Note that if you use the dense backend and its specific feature for handling box inequality constraints, then the first \(n_{in}\) elements of z correspond to multipliers associated to the linear inequality formed with \(C\) matrix, whereas the last \(d\) elements correspond to multipliers associated to the box inequality constraints (see for example solve_dense_qp.cpp or solve_dense_qp.py).

The info subclass

In this table you have on the three columns from left to right: the name of the info subclass item, its default value, and then a short description of it.

Info item	Default value	Description
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.
iter	0	Total number of iterations.
iter_ext	0	Total number of outer iterations.
mu_updates	0	Total number of mu updates.
rho_updates	0	Total number of rho updates.
status	PROXQP_NOT_RUN	Status of the solver.
setup_time	0	Setup time (takes into account the equilibration procedure).
solve_time	0	Solve time (takes into account the first factorization).
run_time	0	the sum of the setup time and the solve time.
objValue	0	The objective value to minimize.
pri_res	0	The primal residual.
dua_res	0	The dual residual.

Note finally that when initializing a QP object, by default, the proximal step sizes (i.e., rho, mu_eq, and mu_in) are set up by the default values defined in the Setting class. Hence, when doing multiple solves, if not specified, their values are re-set respectively to default_rho, default_mu_eq, and default_mu_in. A small example is given below in C++ and Python.

examples/cpp/init_with_default_options.cpp

examples/python/init_with_default_options.py

#include <iostream> #include <proxsuite/proxqp/dense/dense.hpp> // load the dense solver backend #include <proxsuite/proxqp/utils/random_qp_problems.hpp> // used for generating a random convex qp using namespace proxsuite::proxqp; using T = double; int main() { isize dim = 10; isize n_eq(dim / 4); isize n_in(dim / 4); // generate a random qp T sparsity_factor(0.15); T strong_convexity_factor(1.e-2); dense::Model<T> qp_random = utils::dense_strongly_convex_qp( dim, n_eq, n_in, sparsity_factor, strong_convexity_factor); dense::QP<T> qp( dim, n_eq, n_in); // create the QP // initialize the model, along with another rho parameter qp.settings.initial_guess = InitialGuessStatus::NO_INITIAL_GUESS; qp.init(qp_random.H, qp_random.g, qp_random.A, qp_random.b, qp_random.C, qp_random.l, qp_random.u, true, /*rho*/ 1.e-7, /*mu_eq*/ 1.e-4); // Initializing rho sets in practive qp.settings.default_rho value, // hence, after each solve or update method, the qp.results.info.rho value // will be reset to qp.settings.default_rho value. qp.solve(); // So if we redo a solve, qp.settings.default_rho value = 1.e-7, hence // qp.results.info.rho restarts at 1.e-7 The same occurs for mu_eq. qp.solve(); // There might be a different result with WARM_START_WITH_PREVIOUS_RESULT // initial guess option, as by construction, it reuses the last proximal step // sizes of the last solving method. }

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 qp = proxsuite.proxqp.dense.QP(n, n_eq, n_in) # generate a random QP H, g, A, b, C, u, l = generate_mixed_qp(n) # initialize the model of the problem to solve qp.init(H, g, A, b, C, l, u, rho=1.0e-7, mu_eq=1.0e-4) qp.solve() # If we redo a solve, qp.settings.default_rho value = 1.e-7, hence qp.results.info.rho restarts at 1.e-7 # The same occurs for mu_eq. qp.solve() # There might be a different result with WARM_START_WITH_PREVIOUS_RESULT initial guess option, as # by construction, it reuses the last proximal step sizes of the last solving method.

The solver's status

The solver has five status:

• PROXQP_SOLVED: the problem is solved.
• PROXQP_MAX_ITER_REACHED: the maximum number of iterations has been reached.
• PROXQP_PRIMAL_INFEASIBLE: the problem is primal infeasible.
• PROXQP_SOLVED_CLOSEST_PRIMAL_FEASIBLE: the closest feasible problem in L2 sense is solved.
• PROXQP_DUAL_INFEASIBLE: the problem is dual infeasible.
• PROXQP_NOT_RUN: the solver has not been run yet.

Infeasibility is detected using the necessary conditions exposed in section 3.4. More precisely, primal infeasibility is assumed if the following conditions are matched for some non zeros dy and dz (according to the eps_prim_inf variable set by the user):

$$\begin{equation}\label{eq:approx_qp_sol_prim_inf}\tag{PrimalInfeas} \begin{aligned} &\left\{ \begin{array}{ll} |A^Tdy|_{\infty} \leq \epsilon_{\text{primal_inf}} |dy|_{\infty} , \\ b^T dy \leq -\epsilon_{\text{primal_inf}} |dy|_{\infty}, \\ |C^Tdz|_{\infty} \leq \epsilon_{\text{primal_inf}} |dz|_{\infty}, \\ u^T [dz]_+ - l^T[-dz]_+\leq -\epsilon_{\text{primal_inf}} |dz|_{\infty}. \\ \end{array} \right. \end{aligned} \end{equation}$$

Dual infeasibility is assumed if the following two conditions are matched for some non-zero dx (according to the eps_dual_inf variable set by the user):

$$\begin{equation}\label{eq:approx_qp_sol_dual_inf}\tag{DualInfeas} \begin{aligned} &\left\{ \begin{array}{ll} |Hdx|_{\infty} \leq \epsilon_{\text{dual_inf}} |dx|_{\infty} , \\ g^T dx \leq -\epsilon_{\text{dual_inf}} |dx|_{\infty}, \\ |Adx|_{\infty}\leq \epsilon_{\text{dual_inf}} |dx|_{\infty}, \\ (Cdz)_i \geq \epsilon_{\text{dual_inf}} |dx|_{\infty}, \mbox{ if } u_i = +\infty, \\ (Cdz)_i \leq \epsilon_{\text{dual_inf}} |dx|_{\infty}, \mbox{ otherwise }. \end{array} \right. \end{aligned} \end{equation}$$

If the problem turns out to be primal or dual infeasible, then x, y, and z stored in the results class will be the certificate of primal or dual infeasibility. More precisely:

• if the problem is dual infeasible, Qp.results.x will be the certificate dx of dual infeasibility satisfying \eqref{eq:approx_qp_sol_dual_inf} at precisifion Qp.settings.eps_dual_inf specified by the user,
• if the problem is primal infeasible, Qp.results.y and Qp.results.z will be, respectively, the certificates dy and dz of primal infeasibility satisfying \eqref{eq:approx_qp_sol_prim_inf} at precision Qp.settings.eps_primal_inf specified by the user.

Which backend to use?

We have the following generic advice for choosing between the sparse and dense backend. If your problem is not:

• too large (less than some thousand variables),
• and too sparse (a sparsity ratio of your matrices greater than 0.1),

then we recommend using the solver with the dense backend.

The sparsity ratio of matrix A is defined as:

$$ \text{sparsity}(A) = \frac{\text{nnz}(A)}{\text{number}_{\text{row}}(A) * \text{number}_{\text{col}}(A)}, $$

which accounts for the percentage of non-zero elements in matrix A.

Some important remarks when computing timings

We first provide some details about what is measured in the setup and solve time of ProxQP, which is of some importance when doing benchmarks with other solvers, as they can measure different things in a feature with a similar name.

Then we conclude this documentation section with some compilation options for ProxQP which can considerably speed up the solver, considering your OS architecture.

What do the timings take into account?

An important remark about quadratic programming solver is that they all rely at some point on a factorization matrix algorithm, which constitutes the time bottleneck of the solver (as the factorization has a cubic order of complexity wrt dimension of the matrix to factorize).

Available solvers often have a similar API as the one we propose, with first an "init" method for initializing the model, and then a "solve" method for solving the QP problem. For not biasing the benchmarks, it is important to know where is done the first matrix factorization, as it constitutes a considerable cost. In our API, we have decided to make the first factorization in the solve method, as we consider that factorizing the problem is part of the solving part. Hence in terms of timing:

• results.info.setup_time: measures only the initialization of the model and also the preconditioning procedure (if it has been activated),
• results.info.solve_time: measures everything else (including the first factorization),
• results.info.run_time = results.info.setup_time + results.info.solve_time.

It is important to notice that some other solvers API have made different choices. For example, OSQP measures in the setup time the first factorization of the system (at the time ProxQP algorithm was published). Hence we recommend that for benchmarking ProxQP against other solvers, you should compare ProxQP runtime against the other solvers' runtime (i.e., everything from what constitutes their setup to their solve method). Otherwise, the benchmarks won't take into comparable account timings.

Architecture options when compiling ProxSuite

We highly encourage you to enable the vectorization of the underlying linear algebra for the best performance. You just need to activate the cmake option BUILD_WITH_SIMD_SUPPORT=ON, like:

mkdir build && cd build
cmake .. -DCMAKE_BUILD_TYPE=Release -DBUILD_TESTING=OFF -DBUILD_WITH_SIMD_SUPPORT=ON
make
make install

ProxQP can be compiled more precisely with two SIMD instructions options for x86 instruction set architectures: AVX-2 and AVX-512. They can considerably enhance the speed of ProxQP, and we encourage you to use them if your OS is compatible with them. For your information, our benchmarks for ProxQP algorithm have been realised with AVX-2 compilation option.