An infeasible-start framework for convex quadratic optimization, with application to constraint-reduced interior-point and other methods

A framework is proposed for solving general convex quadratic programs (CQPs) from an infeasible starting point by invoking an existing feasible-start algorithm tailored for inequality-constrained CQPs. The central tool is an exact penalty function scheme equipped with a penalty-parameter updating rule. The feasible-start algorithm merely has to satisfy certain general requirements, and so is the updating rule. Under mild assumptions, the framework is proved to converge on CQPs with both inequality and equality constraints and, at a negligible additional cost per iteration, produces an infeasibility certificate, together with a feasible point for an (approximately) \(\ell _1\)-least relaxed feasible problem, when the given problem does not have a feasible solution. The framework is applied to a feasible-start constraint-reduced interior-point algorithm previously proved to be highly performant on problems with many more inequality constraints than variables (“imbalanced”). Numerical comparison with popular codes (OSQP, qpOASES, MOSEK) is reported on both randomly generated problems and support-vector machine classifier training problems. The results show that the former typically outperforms the latter on imbalanced problems. Finally, application of the proposed infeasible-start framework to other feasible-start algorithms is briefly considered, and is tested on a simplex iteration.

1. Note however that, in the context of feasible-direction methods for general nonlinear optimization problems, methods that craftily combine the two phases have been known for decades; see [34, 36].

2. One option, not considered in the present paper, may be to resort to homogeneous self-dual (HSD) embedding methods, first proposed in [44] for the solution of linear optimization problems and later extended in [43] to that of linear complementarity problems (which include convex quadratic programming as a special case). We are unaware of studies combining HSD with a simplex algorithm or a constraint-reduction method. Development, analysis, and numerical testing of such combinations may be of much interest.

3. Earlier algorithms that include such features include, in the context of sequential quadratic programming (in nonlinear optimization), that of [10].

4. Only Proposition 1 invokes boundedness of the primal optimal solution set. It is not clear at this point whether such assumption is necessary indeed. In any case, if it is, it of course can be achieved by imposing large bounds to the components of \(\mathbf{x}\).

5. Note that, while, for primal linear optimization problems in standard form, dual strict feasibility implies boundedness of the primal optimal solution set (see, e.g., [42, Theorem 2.3]), it does not for problem (P).

6. An \(\ell _\infty \) penalty function can be substituted for this \(\ell _1\) penalty function with minor adjustments: see [22, 23] for details.

7. Iterations that require a strictly primal-feasible \(x^\ell \) are automatically catered to, subject to providing a strictly feasible initial \(({\mathbb {x}},{\mathbb {z}})\) in the Initialization step of Meta-Algorithm IS. Due to the relaxation, such strictly feasible points for (P\(_\varphi \)) are readily available.

8. A constraint-reduced version of Mehrotra’s Predictor Corrector.

9. While the authors of [28] (who are also the authors of the present paper) were not able to do away with such linear-independence assumption in proving that Theorem 1 in [28] holds, intuition and extensive numerical testing suggest that Assumption 3 can be dropped.

10. Numerical experimentation, including tests with small-size problems that do not satisfy Assumption 4, suggests that mere boundedness of the optimal solution set of (P\('\)), as implied by Assumption 1\('\), is sufficient for the results to hold. A proof of this is elusive at this time though. In any case, boundedness of the feasible set of course can be achieved by imposing appropriately large bounds to the components of \(\mathbf{x}\).

11. Alternatively, \(\mathbf{x}^{\hat{k}}\) is also feasible for the adjusted (rather than relaxed) problem obtained by still replacing \(\mathbf{b}\) with \(\mathbf{b}'\) but then including instead the equality constraints \(\mathbf{C}\mathbf{x}=\mathbf{d}'\), with \(\mathbf{d}'=\mathbf{C}\mathbf{x}^{\hat{k}}\).

12. Approximate primal feasibility \((\mathbf{s},\mathbf{t}_+,\mathbf{t}_-)\ge \mathbf{0}\) (or \(\approx \mathbf{0}\)) is implicitly taken into account in the last three terms in the numerator.

13. A possibility would be to freeze \(\mathbf{z}\) at zero when \(\mathbf{x}^0\) is primal feasible and indeed, when a component \(z_i^k\) of \(\mathbf{z}^k\) reaches zero at some iteration k, freeze that component to zero thereafter. This was not done in the tests reported here.

14. Here we denote the CR-MPC algorithm with constraint reduction turned off as MPC\(^{*}\) (rather than MPC) to avoid confusion with the original MPC algorithm in [31].

15. To avoid biases due to different stopping criteria, we also experimented with tolerances set to \(10^{-6}\) (while keeping \(\texttt {tol} = 10^{-8}\) for all MPC versions) and observed results with slightly fewer iterations and nearly identical computation time.

16. Here the iteration counts of OSQP and qpOASES are significantly higher than other solver, which is as expected. However, OSQP failed to converge within 10 000 iterations in several problem instances. We observed that OSQP converges on all problem instances when the tolerance is relaxed from \(10^{-8}\) to \(10^{-4}\), but the computation time is still higher than other compared solvers.

17. We used an implementation due to Luke Winternitz, who kindly made it available to us.

18. Two-phase RPS is not included in Fig. 1b because when \(p\not =0\), the dual of (P) is not in standard dual form, a form it requires.

19. Alternatively, following the suggestion made at the very end of Sect. 5, an \((n+1)\)-variable problem with feasible start could be solved.

20. A Matlab-formatted version of these data sets was kindly made available to us by Jin Jung, first author of [26].

21. Typical base iterations accept, some require, problems to be expressed in standard form.

22. Slightly modified versions of penalty_parameter_update(), of the Reqs PU’ of Sect. 4 and, accordingly, of the penalty-parameter updating rule of Sect. 6.1 are to be used here. Specifically, all instances of \(\mathbf{z}\) and the associated dual variable \(\varvec{\xi }\) must be removed.

23. This unbounded feasible primal sequence can be constructed by taking \(\{{\mathbb {x}}^{\ell }\}\) such that \({\mathbb {x}}_B^{\ell }:= \hat{x}_B + t^\ell \Delta \hat{x}_B\) and \({\mathbb {x}}_{B^c}^{\ell }:= t^\ell e_j\), where \(\hat{x}\), \(\Delta \hat{x}_B\), and B are respectively the primal BFS, direction, and basis at which unboundedness is detected in Iteration RPS (Pseudocode 6), \(e_j\) is a standard unit vector that takes value 1 at the entering column index j and 0 elsewhere, and \(t^\ell \rightarrow \infty \) as \(\ell \rightarrow \infty \). Here the direction \(\Delta \hat{x}_B\) and the entering column index j would be appended to the output of RPS_iteration() as part of \(\texttt {info}^{\mathrm{out}}\) (together with \(\varvec{\lambda }_1^{\mathrm{out}}\) and \(\varvec{\lambda }_2^{\mathrm{out}}\)) when unboundedness is detected. In Meta-Algorithm IS, this unbounded sequence would trigger penalty_parameter_update() to increase the penalty parameter \(\varphi \). After the increment of \(\varphi \), the primal variable then needs to be reinitialized to the BFS \(\hat{x}\) (at which unboundedness was detected) so that RPS_iteration() can be applied to solve the updated, possibly bounded problem.

Authors and Affiliations

1. Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN, 37831, USA

   M. Paul Laiu

2. Department of Electrical and Computer Engineering and Institute for Systems Research, University of Maryland, College Park, MD, 20742, USA

   André L. Tits

Corresponding author

Correspondence to M. Paul Laiu.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This manuscript has been authored, in part, by UT-Battelle, LLC, under Contract No. DE-AC0500OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for the United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).

M. Paul Laiu’s research was sponsored by the Office of Advanced Scientific Computing Research and performed at the Oak Ridge National Laboratory, which is managed by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725.

Appendix

Two algorithms considered in this paper are presented here for the sake of completeness. Some notation used in the appendix is “local”, i.e., only applies to this appendix; indeed, some of the same symbols are used in the main body of the paper for unrelated quantities.

A Feasible-start algorithm CR-MPC

In this appendix, we present a brief description of Algorithm CR-MPC in [28], as applied to a generic inequality-constrained CQP of the form (6). We refer the reader to [28] for the complete version of Algorithm CR-MPC in full detail.

For ease of reference, we restate the generic inequality-constrained CQP (6) here:

Algorithm CR-MPC is a feasible-start, descent CQP variant of S. Mehrotra’s predictor–corrector algorithm [31]; computational efficiency is enhanced on imbalanced problems by solving a mere approximate Newton-KKT system for search directions (“constraint reduction”). Specifically, the Newton-KKT system is approximated by including only a selected subset of constraints, which aims to estimate the active constraint set at a solution to (58). Also, an adaptive positive definite regularization W of the original Hessian H is substituted in that system.

Let Q denote the index set of the selected set of constraints, then the affine-scaling search direction is given by the solution of the approximate Newton-KKT system for (58)

and the corrector direction is computed by solving

where \(\mu _{(Q)}:={\left\{ \begin{array}{ll} s_Q^T\lambda _Q/|Q|\,, &{} \quad \text {if } Q\ne \emptyset \\ 0\,, &{} \quad \text {otherwise} \end{array}\right. }\) and \(\sigma =(1-\alpha ^{\mathrm{a}})^3\) with \(\alpha ^{\mathrm{a}}\) the maximum feasible step for the affine-scaling direction. A pseudocode for an iteration in Algorithm CR-MPC is given in Pseudocode 5, which can serve as base_iteration() in Meta-Algorithm IS with minimal modifications that tailor the code for penalized relaxed problems of the form (P\(_\varphi \)).

B Revised primal simplex (RPS) algorithm

In this appendix, we give a brief description of the revised primal simplex method, which is considered in Sect. 6.3.1 in the comparison to the CR-MPC method discussed in Appendix A, as well as in Sect. 7 as a base iteration in the IS framework. We consider the RPS method that handles linear optimization problems in primal standard form

The RPS method considered here starts from an initial primal BFS for (61) and, at each iteration, it reduces the objective function value by moving to an adjacent primal BFS with lower value. Specifically, each RPS iteration forms a descent search direction by adding an entering column to the basis associated to the BFS, and takes a step towards the descent direction until one of the BFS entries becomes zero, the index of which is then identified as the leaving column and is removed from the basis. Pseudocode 6 presents the main steps of an RPS iteration, which can be used as iteration_RPS() in the wrapped RPS base iteration in Pseudocode 4 by incorporating the penalty parameter \(\varphi \) into the input.

Reprints and permissions