Solving quadratic programs to high precision using scaled iterative refinement

Quadratic optimization problems (QPs) are ubiquitous, and solution algorithms have matured to a reliable technology. However, the precision of solutions is usually limited due to the underlying floating-point operations. This may cause inconveniences when solutions are used for rigorous reasoning. We contribute on three levels to overcome this issue. First, we present a novel refinement algorithm to solve QPs to arbitrary precision. It iteratively solves refined QPs, assuming a floating-point QP solver oracle. We prove linear convergence of residuals and primal errors. Second, we provide an efficient implementation, based on SoPlex and qpOASES that is publicly available in source code. Third, we give precise reference solutions for the Maros and Mészáros benchmark library.


Introduction
Quadratic optimization problems (QPs) are optimization problems with a quadratic objective function and linear constraints. They are of interest directly, e.g., in portfolio optimization or support vector machines [1]. They also occur as subproblems in sequential quadratic programming, mixed-integer quadratic programming, and nonlinear model predictive control. Efficient algorithms are usually of active set, interior point, or parametric type. Examples of QP solvers are BQPD [5], CPLEX [2], Gurobi [12], qp solve [10], qpOASES [4], and QPOPT [7]. These QP solvers have matured to reliable tools and can solve convex problems with many thousands, sometimes millions of variables. However, they calculate and check the solution of a QP in floating-point arithmetic. Thus, the claimed precision may be violated and in extreme cases optimal solutions might not be found. This may cause inconveniences, especially when solutions are used for rigorous reasoning.
One possible approach is the application of interval arithmetic. It allows to include uncertainties as lower and upper bounds on the modeling level, see [13] for a survey for the case of linear optimization. As a drawback, all internal calculations have to be performed with interval arithmetic. Hence standard solvers can not be used any more, the computation times increase, and solutions may be very conservative.
We are only aware of one advanced algorithm that solves QPs exactly over the rational numbers. It is designed to tackle problems from computational geometry with a small number of constraints or variables, [6]. Based on the classical QP simplex method [16], it replaces critical calculations inside the QP solver by their rational counterparts. Heuristic decisions that do not affect the correctness of the algorithm are performed in fast floating-point arithmetic.
In this paper we propose a novel algorithm that can use efficient floating-point QP solvers as a black box. Our method is inspired by iterative refinement, a standard procedure to improve the accuracy of an approximate solution for a system of linear equalities, [15]: The residual of the approximate solution is calculated, the linear system is solved again with the residual as a right-hand side, and the new solution is used to refine the old solution, thus improving its accuracy. A generalization of this idea to the solution of optimization problems needs to address several difficulties: most importantly, the presence of inequality constraints; the handling of optimality conditions; and the numerical tolerances that floating-point solvers can return in practice.
For LPs this has first been developed in [9]. The approach refines primal-dual solutions of the Karush-Kuhn-Tucker (KKT) conditions and comprehends scaling and calculations in rational arithmetic. We generalize further and discuss the specific issues due to the presence of a quadratic objective function. The fact that the general approach carries over from LP to QP was remarked in [8]. Here we provide the details, provide a general lemma showing how the residuals bound the primal and dual iterates, and analyze the computational behavior of the algorithm based on an efficient implementation that is made publicly available in source code and can be used freely for research purposes.
The paper is organized as follows. In Section 2 we define and discuss QPs and their refined and scaled counterparts. We give one illustrating and motivating example for scaling and refinement. In Section 3 we formulate an algorithm and prove its convergence properties. In Section 4 we consider performance issues and describe how our implementation based on SoPlex and qpOASES can be used to calculate solutions for QPs with arbitrary precision. In Section 5 we discuss run times and provide solutions for the Maros and Mészáros benchmark library, [14]. We conclude in Section 6 with a discussion of the results and give directions for future research and applications of the algorithm.
In the following we will use · for the maximum norm · ∞ . The maximal entry of a vector max i {v i } is written as max{v}. Inequalities a ≤ b for a, b ∈ Q n hold componentwise.

Refinement and Scaling of Quadratic Programs
In this section we collect some basic definitions and results that will be of use later on. We consider convex optimization problems of the following form.
Definition 1 (Quadratic optimization problem (QP)). Let a symmetric matrix Q ∈ Q n×n , a matrix A ∈ Q m×n , and vectors c ∈ Q n , b ∈ Q m , l ∈ Q n be given. We consider the quadratic optimization problem (QP) assuming that (P ) is feasible and bounded, and Q is positive semi-definite on the feasible set. The rounded version of this convex rational data QP will be denoted as (P ).
A point x * ∈ Q n is a global optimum of (P ) if and only if it satisfies the Karush-Kuhn-Tucker (KKT) conditions, [3], i.e., if multipliers y * ∈ Q m exist such that The pair (x * , y * ) is then called KKT pair of (P ). Primal feasibility is given by (1a-1b), dual feasibility by (1c), and complementary slackness by (1d). Refinement of this system of linear (in-)equalities is equivalent to the refinement of (P ).
Definition 2 (Refined QP). Let the QP (P ), scaling factors ∆ P , ∆ D ∈ Q + and vectors x * ∈ Q n , y * ∈ Q m be given. We define the refined QP as The rounded version of this refined and scaled rational data QP will be denoted as (P ∆ ).
The following theorem is the basis for our theoretical and algorithmic approaches. It is a generalization of iterative refinement for LP and was formulated and proven in [8,Theorem 5.2]. Again, primal feasibility refers to (1a-1b), dual feasibility to (1c), and complementary slackness to (1d).
Theorem 3 (QP Refinement). Let the QP (P ), scaling factors ∆ P , ∆ D ∈ Q + , vectors x * ∈ Q n , y * ∈ Q m , and the refined QP (P ∆ ) be given. Then for anyx ∈ R n ,ŷ ∈ R m and tolerances ǫ P , ǫ D , ǫ S ≥ 0: 1.x is primal feasible for (P ∆ ) within an absolute tolerance ǫ P if and only if x * +x ∆P is primal feasible for (P ) within ǫ P /∆ P .
2.ŷ is dual feasible for (P ∆ ) within an absolute tolerance ǫ D if and only if y * +ŷ ∆D is dual feasible for (P ) within ǫ D /∆ D .
For illustration, we investigate the following example.
Example 4 (QP Refinement). Consider the QP with two variables An approximate solution to a tolerance of 10 −6 is x * 1 = x * 2 = 0 with dual multiplier y * = 1. This solution is slightly primal and dual infeasible, but the solver can not recognize this on this scale. The situation is depicted in Figure 1 on the left.
The point x * seems to be the optimal solution satisfying the equality constraint and the brown circle representing the level curve of the objective function indicates the Figure 1: Illustration of nominal QP (left) and of refined QP (right) for Example 4. The scaled (and shifted) QP (P ∆ ) works like a zoom (and shift) for (P ), allowing to correct the solution x * (orange dot) from (0, 0) to (10 −6 , 0).

Iterative Refinement for Quadratic Programming
To solve quadratic programs to arbitrary, a priori specified precisions, we apply the refinement idea from the previous section iteratively in Algorithm 1. Algorithm 1 expects QP data (Q, A, c, b, l) in rational precision, primal and dual termination tolerances (ǫ P , ǫ D ), complementary slack termination tolerance (ǫ S ), scaling limit α > 1 and iteration limit k max . First the rounded QP (P ) is solved with a floatingpoint QP solver oracle which returns optimal primal and dual solution vectors (Line 2 ). In Line 3 the main loop begins. The primal violations for constraints (b, Line 4 ) and for bounds (l, Line 5 ) are calculated. The maximal primal violation is saved as δ P in Line 6 . The reduced cost vectorĉ and its maximal violation δ D are calculated in Lines 7-8 . In Line 9 the scaling factor ∆ k is chosen as the maximum of α∆ k−1 and the inverses of the violations δ P and δ D . The complementary slack violation δ S is calculated in Line 10 . If the primal, dual and complementary slack violations are already below the specified tolerances the loop is stopped (Lines 11-12 ) and the optimal solution is returned (Line 17 ). Else (Line 13 ) the refined, scaled, and rounded QP (P ∆ k ) is solved with the floating-point QP oracle in Line 14 . We save the floating-point optimal primal and dual solution vectors (Line 15 ). We scale and add them to the current iterate (x k , y k ) to obtain (x k+1 , y k+1 ), Line 16 .
Note that all calculations except the expensive solves of the QPs are done in rational precision. Algorithm 1 uses only one scaling factor ∆ k for primal and dual infeasibility Algorithm 1 Iterative QP Refinement (IQPR) 1: Input: (P ) in rational precision, termination tolerances ǫ P and ǫ D and ǫ S , scaling limit α > 1, iteration limit k max 2: Initialization: ∆ 0 ← 1, solve (P ) approximately, save optimal (x 1 , y 1 ) 3: for k ← 1, 2, ..., k max do if δ P ≤ ǫ P and δ D ≤ ǫ D and δ S ≤ ǫ S then (x * , y * ) ← KKT pair returned as optimal 16: to avoid the scaling of the quadratic term of the objective. Keeping this matrix and the constraint matrix A fixed gives QP solvers the possibility to reuse the internal factorization of the basis system between refinements, as the transformation does not change the basis. This allows to use efficient hotstart techniques for all solves after the initial solve.
To investigate the performance of the algorithm we make, in analogy with the LP case [9, Assumption 2.9], the following assumption.
Assumption 5 (QP solver accuracy). We assume that there exists ǫ ∈ [0, 1) and σ ≥ 0 such that the QP solver oracle returns for all rounded QPs (P ∆ k ) solutions (x,ȳ) that satisfy Note that this ǫ corresponds to a termination tolerance passed to a floating-point solver, while the algorithm uses overall termination tolerances ǫ P and ǫ D and a scaling limit α > 1 per iteration. We denoteǫ = max{1/α, ǫ}.
Lemma 6 (Termination and residual convergence). Algorithm 1 applied to a primal and dual feasible QP (P ) and using a QP solver that satisfies Assumption 5 will terminate in at most k max = max { log(ǫ P )/ log(ǫ), log(ǫ D )/ log(ǫ), log(ǫ S /σ)/(2 log(ǫ)) + 1 } iterations. Furthermore, after each iteration k = 1, 2, ... the primal-dual iterate (x k , y k ) and the scaling factor ∆ k satisfy Proof. We prove (3) by induction over k, starting with k = 1. Asǫ ≥ ǫ, the claims (3b-3e) follow directly from Assumption 5. Using Lines 6, 4-5 , and Assumption 5 we obtain and with Lines 8,7 and Assumption 5 Thus from Line 9 we have and hence claim (3a) for the first iteration. Assuming (3) holds for k we know that δ P,k , δ D,k ≤ǫ k and ∆ k ≥ 1/ǫ k . With the scaling factor ∆ k using x * = x k and y * = y k we scale the QP (P ) as in Theorem 3 and hand it to the QP solver. By Theorem 3 this scaled QP is still primal and dual feasible and by Assumption 5 the solver hands back a solution (x,ŷ) with tolerance ǫ ≤ǫ. Therefore using Theorem 3 again the next refined iterate (x k+1 , y k+1 ) has a tolerance in QP (P ) ofǫ/∆ k ≤ǫ k+1 , which proves (3b-3d).
The results show that even though we did not use the violation of the complementary slackness to choose the scaling factor in Algorithm 1, the complementary slackness violation is bounded by the square ofǫ.

Remark 7 (Nonconvex QPs). Algorithm 1 can also be used to calculate high precision KKT pairs of nonconvex QPs. If the black box QP solver hands back local solutions of the quality specified in Assumption 5 Lemma 6 holds as well for nonconvex QPs and Algorithm 1 returns a high precision local solution.
However, assuming strict convexity, an even stronger result holds. Inspired by the result for the equality-constrained QP [3, Proposition 2.12] we investigate how this righthand side convergence of the KKT conditions is related to the primal-dual solution.
Lemma 8 (Primal and dual solution accuracy). Let QP (P ) be given and be strictly convex, the minimal and maximal eigenvalues of Q be λ min (Q) and λ max (Q), respectively, and the minimal nonzero singular value of A be σ min (A). Let the KKT conditions (1) hold for (x * , y * , z * ), i.e., and the disturbed KKT conditions for disturbances e ∈ Q m , g, f, h ∈ Q n , and i ∈ Q hold for (x, y, z), i.e., Then and Proof. By (4a) and (5a) we have that A(x − x * ) = e and taking the Moore-Penrose pseudoinverse A + of A we define δ = A + e with Aδ = e and δ 2 ≤ σ min (A) −1 e 2 . Using this we can start to derive the dual bound by taking the difference of (4b) and (5b) Multiplying from the left with Q −1 (A T (y − y * ) + (z − z * )) transposed gives The second term of (9) can be expressed as

With this and (9) we bound from above the term
Taking the norm on the right and reordering terms gives

This is a quadratic expression in
It has two roots, but only one is greater than zero and bounds This can be expressed as where a and d are defined as above. This proves (6). To prove the primal bound we multiply equation (8) Taking norms gives the inequality Combining the dual bound (11) and (12) we get the final primal bound which proves (7).
Note that λ max (Q)λ min (Q) is the condition number of Q. The above assumption and lemmas can be summarized to a statement about the convergence of the algorithm for a strictly convex QP.
Theorem 9 (Rate of convergence). Algorithm 1 with corresponding input and using a QP solver satisfying Assumption 5 solving the QP (P ) that is also strictly convex has a linear rate of convergence with a factor ofǫ 1/2 for the primal iterates, i.e.
with x * being the unique solution of (P ).
Proof. By Assumption 5 and Lemma 6 we know that the right-hand side errors of the KKT conditions are bounded by Here we set the violations h of the inequality KKT multipliers z to zero and count them as additional dual violations g for simplicity. Also note that in Lemma 8 the bound is just depending on the norm of the right-hand side violation vectors, two different violation vectors with the same norm give the same bound. Therefore we just consider the norms. Combining the above with Lemma 8 we get for the primal iterate in iteration k with constants Looking at the quotient and seeing that This theoretical investigation shows us two things. First, we have linear residual convergence with a rate ofǫ. In contrast to usual convergence results our algorithm achieves this rate in practice by the use of rational computations if the floating-point solver delivers solutions of the quality specified in Assumption 5. This is also checked by the rational residual calculation in our algorithm in every iteration. Second, this residual convergence implies primal iterate convergence with a linear rate ofǫ 1/2 for strictly convex QPs.

Implementation
Following previous work [9] on the LP case we implemented Algorithm 1 in the same framework within the SoPlex solver [17], version 2.2.1.2, using the GNU multiple precision library (GMP) [11] for rational computations, version 6.1.0. As underlying QP solver we use the active-set solver qpOASES [4] version 3.2. This version of qpOASES was originally designed for small to medium QPs (up to 1,000 variables and constraints). Furthermore, we implemented an interface to a pre-release version of qpOASES 4.0, which can handle larger, sparse QPs of a size up to 40,000 variables and constraints. Compared to the matured qpOASES 3.2, this version is not yet capable of hotstarts and in some cases less robust. Nevertheless, it allows us to study the viability of iterative refinement on larger QPs. The source code of our implementation is available for download in a public repository. 1 In order to treat general QPs with inequalities, our implementation recovers the form (P ) by adding one slack variable per inequality constraint. Note that not only lower, but also upper bounds on the variables need to be considered. However, this is a straightforward modification to our algorithm and realized in the implementation.
One advantage of using the active-set QP solver qpOASES is the returned basis information. We use the basis in three aspects: first, to calculate dual and complementary slack violations; second, to explicitly set nonbasic variables to their lower bounds after the refinement step in Line 16 of Algorithm 1; and third, to compute a rational solution primal tolerance (ǫ P ) 1e-100 1e-100 1e-100 1e-100 1e-10 dual tolerance (ǫ D ) 1e-100 1e-100 1e-100 1e-100 1e-10 maxscaleincrement (α) 1e12 1e12 1e12 1e12 1e12 sparse no no yes yes no max num backstepping (l max ) 10 10 10 10 1 refinement limit (k max ) 300 50 50 50 10 ratfac minstalls 2 0 0 51 30 defined by the corresponding system of linear equations. This is solved by a standard LU factorization in rational arithmetic. If the resulting primal-dual solution is verified to be primal and dual feasible, the algorithm can terminate early with an exact optimal basic solution.
Since the LU factorization can be computationally expensive, we only perform this step if we believe the basis to be optimal. When the QP solver returns the same basis as "optimal" for several iterations this can be used as a heuristic indicator that the basis might be truly optimal, even if the iteratively corrected numerical solution is not yet exact. Hence, the number of consecutive iterations with the same basis is used to trigger a rational basis system solve. This can be controlled by a threshold parameter called "ratfac minstalls", see Table 1.
If the floating-point solver fails to compute an approximately optimal solution, we decrease the scaling factor by two orders of magnitude and try to solve the resulting QP again. The scaling factor is reduced either until the maximum number of backstepping rounds is reached or until the next backstepping round would result in a scaling factor lower than in the last refinement iteration (k − 1).
Currently, the implementation has no detection mechanism for infeasible or unbounded QPs. Because the testset at hand contains only convex and feasible QPs, this does not affect the numerical experiments. The default parameter set (s1) of our implementation is given in Table 1. The other four parameter sets (s2-s5) are used for our numerical experiments to derive either exact or inexact solutions.
We exploit the different features of the two qpOASES versions. Version 3.2 has hotstart capabilities that allow to reuse the internal basis system factorization of the preceding optimal basis. Therefore we start in the old optimal basis and build on the progress made in the previous iterations instead of solving the QP from scratch at every iteration. Additionally we increase the termination tolerance and relax other parameters that ensure a reliable solve. This speeds up the solving process and is possible because the inaccuracies, introduced by this in the floating-point solution, are detected anyway and handed back to the QP solver in the next iteration for correction. If the QP solver fails we simply change to reliable settings and resolve the same QP from the same starting basis before downscaling. Hence, in Algorithm 1 each 'solve' statement means: try fast settings first and if this fails switch to slow and reliable settings of qpOASES 3.2. These two sets of options are given in Table 2. For the pre-release version 4.0 we use default settings and no resolves. We either factorize after each iteration or not at all (see Table 1). We perform two different experiments. The goal of the first experiment is to solve as many QPs from the testset as precisely as possible in order to analyze the iterative refinement procedure computationally and to provide exact reference solutions for future research on QP solvers. In the second experiment we want to compare qpOASES (version 3.2, no QP refinement, one solve, default settings) to low accuracy refinement (low tolerance of 1e-10 in Algorithm 1, using also qpOASES 3.2). This allows us to investigate whether refinement could also be beneficial in cases that do not require extremely high accuracy, but a strictly guaranteed solution tolerance in shortest possible runtime.

Experiment 1
We use the three different parameter sets (s2-s4) given in Table 1 to calculate exact solutions. The first set (s2) contains a primal and dual termination tolerance of 1e-100, enables rational factorization in every iteration, and allows for 50 refinements and 10 backsteppings using a dense QP formulation with qpOASES version 3.2. In contrast the other two sets (s3, s4) with qpOASES version 4.0 use a sparse QP formulation, either with factorization in every iteration or without factorization. Table 3 states for each setting the number of instances which it solved exactly, for which tolerance 1e-100 was reached, and for which it failed to produce a high-precision solution. In total these three strategies could solve 91 out of the 138 QPs in the testset exactly and 39 instances within tolerance 1e-100. For eight instances no high-precision solution was computed. These "virtual best" results stated in the fifth column consider for each QP the result of the individual parameter sets that resulted in the smallest violation. It should be emphasized that for each of the three parameter sets there exists at least one instance for which it produced the most accurate solution.
The last column reports the average number of nonzeros of the QPs in the three "virtual best" categories. This suggests that for problems with less nonzeros a higher accuracy was reached. Furthermore, one sees that the entries in the second column (set s2) do not sum to 138 because with the dense QP formulation some of the problems do not return from the cluster due to memory limitations. Hence the set s2 fails in total on 32 instances. For the parameter set s4 without rational factorization we see that one QP is solved by chance exactly while for all others the algorithm terminates with violations greater zero.
In order to solve the 197 (=33+45+118) QPs to high precision the algorithm needed on average 8.84 refinements. This confirms the linear convergence because we bounded the increase of the scaling factor in each iteration by α = 10 12 and terminate after reaching a tolerance of 10 −100 . If qpOASES would consistently return solutions with an accuracy of 10 −12 we would expect the algorithm to need 9 iterations (100/12 ≈ 8.33 . . . rounded up). We see that qpOASES usually delivers solutions of a tolerance below 10 −12 .
Detailed results can be found in the Appendix in Tables 5, 6, and 7. If an exact solution is found qpOASES usually returns the optimal basis in the first three iterations (refinements). Subsequently, the corresponding basis system is solved exactly by a rational LU factorization. For six problems we found that the objective values given in [14] differ from our results by more than 1e-7: GOULDQP2, HS268, S268, HUESTIS, HUES-MOD, and LISWET8. This might be due to the use of a floating-point QP solver with termination tolerance about 1e-7 when originally computing the values reported. The precise objective value can be found in the online material associated with this paper. Experiment 2 In the following the iterative refinement algorithm is set to a termination tolerance 10 −10 and the rational factorization of the basis system is disabled. The refinement limit is set to 10 and the backstepping limit is set to one (parameter set s5). We compare this implementation to qpOASES 3.2 with the three predefined qpOASES settings (MPC, Default, Reliable) that include termination tolerances of 2.2210e-7, 1.1105e-9, and 1.1105e-9, respectively. For these fast solves we select only part of the testset, including the 73 problems that have no more than 1,000 variables and constraints. This corresponds to the problem sizes for which qpOASES 3.2 was originally designed. In order to allow for a meaningful comparison of runtimes, the evaluation only considers QPs which were solved by all three qpOASES 3.2 settings and by refinement to "optimality", where optimality was certified by qpOASES 3.2 (with its internal floating-point checks) or rational checks in our algorithm, respectively.
An overview over the performance results is given in Table 4. We report runtime, QP solver iterations, and the final tolerance reached, each time as arithmetic and shifted geometric mean. To facilitate a more detailed analysis, we consider the series of subsets "> t" of instances, for which at least one algorithm took more than t seconds. Equivalently, we exclude the QPs for which all settings took at most 0.01 s, 0.1 s, 1 s, and 10 s seconds. Defining the exclusion by all instead of one method only avoids a biased definition of these sets of increasing difficulty.
The results show that in no case the mean runtime of the refinement algorithm is larger than the runtime of qpOASES with reliable setting. At the same time, the accuracy reached is always significantly higher. Compared to qpOASES Default, which results in an even lower level of precision, refinement is faster in arithmetic and slightly slower in shifted geometric mean. The QP solver iterations of the refinement are comparable to the MPC setting. When looking at the different subsets we see that for QPs with larger runtime the refinement approach performs relatively better (smaller runtime, iterations and lower tolerance) than the three qpOASES 3.2 standard settings. The refinement guarantees the tolerance of 1e-10 if it does not fail. To achieve this tolerance, for 9 QPs two refinements were necessary, for 21 QPs only one refinement was necessary, and for 35 instances no refinement was necessary at all. The rational computation overhead stated in brackets after the runtime and is well below 2%. The details are shown in Table 8 in the Appendix. Also note that due to exclusion of fails (which mainly occur with the qpOASES MPC settings) the summarized results have a slight bias towards qpOASES.

Conclusion
We presented a novel refinement algorithm and proved linear convergence of residuals and errors. Notably, this theoretical convergence result also carries over to our implementation due to the use of exact rational calculations. We provided high-precision solutions for most of the QPs in the Maros and Mészáros testset, correcting inaccuracies in optimal solution values reported in the literature. This is beneficial for future research on QP solvers that are evaluated on this testset. In a second experiment we saw that iterative refinement provides proven tolerance solutions with smaller or equal computation times compared to qpOASES with "Reliable" solver settings. It can therefore be used as tool to increase the reliability and speed of standard floating-point QP solvers. If optimal solutions are needed for rigorous reasoning or to make decisions in the real world the algorithm presented is useful because it is able to fully ensure a specified tolerance. This tolerance then can be adapted to the necessity of the application at hand. At the same time this comes with little overhead in rational computation time, which is important for practical applications.
Regarding algorithmic research and solver development, our framework also provides the possibility to compare different floating-point QP solvers by looking at the number of refinements needed with each solver to detect optimal bases or solutions of a specified tolerance as a measure for solver accuracy. Solver robustness can be checked precisely because violations are computed in rational precision. In the future, the implementation should be extended for the detection of unbounded or infeasible QPs. Also one could try more general variable transformations, e.g. having a different scaling factor for each variable. As a concluding remark, we hope that the idea of checking numerical results of floating-point algorithms in exact or safe arithmetic will become a future trend when applying or analyzing numerical algorithms.