Regularized Interior Point Methods for Constrained Optimization and Control

Mathematical optimization plays an important role in model-based and data-driven control systems, forming the basis for advanced techniques such as optimal control, nonlinear model predictive control (MPC) and parameters estimation. Significant research effort on computationally efficient real-time optimization algorithms contributed to the success of MPC over the years and yet the demand for fast and reliable methods for a broad spectrum of applications is growing; see [29, 27] and references therein. In order to tackle these challenges, it is desirable to have an algorithm that benefits from warm-starting information, can cope with infeasibility, is robust to problem scaling, and exploits the structure of the problem. In order to reduce computations and increase robustness, a common approach is to relax the requirements on the solutions, in terms of optimality, constraint violation, or both [14, 27]. In this work, we propose to address such features by combining proximal regularization and interior point techniques, for developing a stabilized, efficient and robust numerical method. We advocate for this strategy by bringing together and combining a variety of ideas from the nonlinear programming literature.

Let us consider the constrained nonconvex problem

where functions $f\colon\mathbb{R}^{n}\to\mathbb{R}$ and $c\colon\mathbb{R}^{n}\to\mathbb{R}^{m}$ are (at least) continuously differentiable. Problems with general inequality constraints on $x$ and $c(x)$ can be reformulated as (P) by introducing auxiliary variables. Nonlinear programming (NLP) problems such as (P) have been extensively studied and there exist several approaches for their numerical solution. Interior point (IP) [31, 32], penalty and augmented Lagrangian (AL) [8, 2, 6] and sequential programming [16] schemes are predominant ideas and have been re-combined in many ways [9, 5, 3].

Starting from linear programming, IP methods had a significant impact on the field of mathematical optimization [18]. By solving a sequence of barrier subproblems, they can efficiently handle inequality constraints and scale well with the problem size. The state-of-the-art solver Ipopt, described by [32], is an emblem of this remarkable success. However, relying on Newton-type schemes for approximately solving the subproblems, IP algorithms may suffer degeneracy and lack of constraint qualifications if suitable counter-mechanisms are not implemented. On the contrary, proximal techniques naturally cope with these scenarios thanks to their inherent regularizing action. Widely investigated in the convex setting [26], their favorable properties have been exploited to design (primal-dual) stabilized methods building on the proximal point algorithm [17, 20, 11]. The analysis of their close connection with the AL framework [25] led to the development of techniques applicable to more general problems [22, 10, 24].

The combination of IP and proximal strategies has been successfully leveraged in the context of convex quadratic programming [1, 7] and for linear [13] and nonlinear [28] least-squares problems. With this work we address general NLPs and devise a method for their numerical solution, which can be seen as an extension of a regularized Lagrange–Newton method to handle bound constraints via a barrier function [10], or as a proximally stabilized IP algorithm, generalizing the ideas put forward by [7].¹¹1Although beyond the scope of this paper, topics such as infeasibility detection [3], ill-posed subproblems [2], and feasibility restoration [32, §3.3] are of great relevance in this context. Appropriate mechanisms should be incorporated in practical implementations, as they can significantly improve the performances.

In this section we introduce some optimality concepts, following [6, Ch. 3] and [12, §2].

Notice that (1c) provides a condition modeling the approximate satisfaction of the (elementwise) complementarity condition $\min\{x,-z\}=0$ within some tolerance $\varepsilon\geq 0$. The IP algorithm discussed in Section 3 satisfies a stronger version of these conditions, since the iterates it generates meet the constraints $x\geq 0$ and $z\leq 0$ by construction. Furthermore, we point out that the condition $\min\{x_{i},-z_{i}\}\leq\varepsilon$ is analogous to $-x_{i}z_{i}\leq\varepsilon$, more typical for interior point methods, but does not depend on a specific barrier function, e.g., the logarithmic barrier in [32, §2.1].

We shall consider the limiting behavior of approximate KKT stationary points when the tolerance $\varepsilon$ vanishes. In fact, having $x^{k}\to x^{\ast}$ with $x^{k}$ $\varepsilon_{k}$-KKT stationary for (P) and $\varepsilon_{k}\searrow 0$ does not guarantee KKT stationarity of a limit point $x^{\ast}$ of $\{x^{k}\}$. This issue raises the need for defining KKT stationarity in an asymptotic sense [6, Def. 3.1].

Any local minimizer $x^{\ast}$ for (P) is AKKT stationary, independently of constraint qualifications [6, Thm 3.1].

The methodology presented in this section builds upon the AL framework, interpreted as a proximal point scheme in the nonconvex regime, and IP methods. The basic idea is to construct a sequence of proximally regularized subproblems and to approximately solve each of them as a single barrier subproblem, effectively merging the AL and IP outer loops. Reduced computational cost can be achieved with an effective warm-starting of the IP iterations and with the tight entanglement of barrier and proximal penalty strategies, by monitoring and updating the parameters’ values alongside with the inner tolerance.

A classical approach is to consider a sequence of bound-constrained Lagrangian (BCL) subproblems [8, 6]

where $\rho_{k}>0$ and $\hat{y}^{k}\in\mathbb{R}^{m}$ are some given penalty parameter and dual estimate, respectively. The nonlinearly-constrained Lagrangian (NCL) scheme [22] considers equality-constrained subproblems by introducing an auxiliary variable $s\in\mathbb{R}^{m}$ and the constraint $c(x)=s$. Analogously, a proximal point perspective yields the equivalent reformulation

recovering the dual regularization term obtained, e.g., by [24, 10, 11]. By construction, these regularized subproblems are always feasible and satisfy a strong constraint qualification, namely the LICQ, at all points.

The regularized subproblems (3)–(4) can be numerically solved via IP algorithms. Let us consider a barrier parameter $\mu_{k}>0$ and barrier functions $b_{i}\colon\mathbb{R}\to\overline{\mathbb{R}}$, $i=1,\ldots,n$, each with domain $\operatorname{dom}b_{i}=(0,\infty)$, and such that $b_{i}(t)\to\infty$ as $t\to 0^{+}$ and $b_{i}^{\prime}\leq 0$. Exemplarily, the logarithmic function $x\mapsto-\ln(x)$ is one of such barrier functions. Other choices can be considered as well, e.g., to handle bilateral constraints [4]. We collect these barrier functions to define $b\colon\mathbb{R}^{n}\to\overline{\mathbb{R}}$, $b\colon x\mapsto\sum_{i=1}^{n}b_{i}(x_{i})$, whose domain is $\operatorname{dom}b=(0,\infty)^{n}$. Thus, analogously to [3], a barrier counterpart for the BCL subproblem (3) reads

whereas for the constrained subproblem (4) this leads to

which is a regularized version of [32, Eq. 3] and reminiscent of [5, Eq. 2]. It should be stressed that, in stark contrast with classical AL and IP schemes, we intend to find an (approximate) solution to the regularized subproblem (4) by (approximately) solving only one barrier subproblem (6). Inspired by [9, 7], our rationale is to drive $\rho_{k},\mu_{k}$ and the inner tolerance $\epsilon_{k}$ concurrently toward zero, effectively knitting together proximal and barrier strategies.

It should be noted that a primal (Tikhonov-like) regularization term is not explicitly included in (3)–(6). In fact, the original objective $f$ could be replaced by a (proximal) model of the form $x\mapsto f(x)+\frac{\sigma_{k}}{2}\|x-\hat{x}^{k}\|^{2}$, with some given primal regularization parameter $\sigma_{k}\geq 0$ and reference point $\hat{x}^{k}\in\mathbb{R}^{n}$. However, as this term can be interpreted as an inertia correction, we prefer the subsolver to account for its contribution; cf. [32, §3.1]. In this way, the subsolver can surgically tune the primal correction term as needed, possibly improving the convergence speed, and surpassing the issue that suitable values for $\sigma_{k}$ are unknown a priori.

The overall procedure is detailed in Algorithm 3.1. At every outer iteration, indexed by $k$, Algorithm 3.1 requires to compute an approximate stationary point, with the associated Lagrange multiplier, for the regularized barrier subproblem (6). As the dual estimate $\hat{y}^{k}$ is selected from some bounded set $Y\subset\mathbb{R}^{m}$ at Algorithm 3.1, the AL scheme is safeguarded and has stronger global convergence properties [6, Ch. 4]. The assignment of $z^{k}$ at Algorithm 3.1 follows from comparing and matching the stationarity conditions for (P) and (6). After checking termination, we monitor progress in constraint violation and complementarity, based on (1), and update parameters $\rho_{k}$ and $\mu_{k}$ accordingly, as well as the inner tolerance $\varepsilon_{k}$. At Algorithms 3.1 and 3.1 we consider relaxed conditions for satisfactory feasibility and complementarity as it is preferable to have the sequences $\{\rho_{k}\}$, $\{\mu_{k}\}$, and $\{\epsilon_{k}\}$ bounded away from zero, in order to avoid unnecessary ill-conditioning and tight tolerances. Sufficient conditions to guarantee boundedness of the penalty parameter $\{\rho_{k}\}$ away from zero are given, e.g., by [2, §5]. Remarkably, as established by 4.2 in Section 4, there is no need for the barrier parameter $\mu_{k}$ to vanish in order to achieve $\varepsilon$-complementarity in the sense of (1c), for $\varepsilon>0$.

We shall mention that considering equivalent yet different subproblem formulations may affect the practical performance of the subsolver. It is enlightening to pinpoint the effect of the dual regularization in (6) and to appreciate its interactions with the linear algebra routines used to solve the linear systems arising in Newton-type methods. Although (6) has more (possibly many more) variables than (5), a simple reordering yields matrices with the same structure [10, 24]. Let us have a closer look. Defining the Lagrangian function $\mathcal{L}_{k}(x,y):=f(x)+\mu_{k}b(x)+\langle y,c(x)\rangle$, the stationarity condition for (5) reads $0=\nabla_{x}\mathcal{L}_{k}(x,y_{k}(x))$, where $y_{k}(x):=\hat{y}^{k}+\rho_{k}^{-1}c(x)$, and the corresponding Newton system is

where $H_{k}(x,y)\in\mathbb{R}^{n\times n}$ denotes the Hessian matrix $\nabla_{xx}^{2}\mathcal{L}_{k}(x,y)$ or a symmetric approximation thereof. A linear transformation yields the equivalent linear system

Analogous Newton systems for (6) read

and formally solving for $\delta\lambda$ gives the condensed system

The resemblances between these linear systems are apparent, as well as the differences. The AL relaxation in (5) introduces a dual regularization for both the linear algebra and nonlinear solver, whose hidden constraint $c(x)+\rho_{k}(\hat{y}^{k}-y)=0$ holds pointwise due to the identity $y=y_{k}(x)$. We remark that, entering the (2,2)-block, the dual regularization prevents issues due to linear dependence. Furthermore, the primal regularization is left to the inertia correction strategy of the subsolver, affecting the (1,1)-block as in [32, §3.1]. If the approximation $H_{k}(x,y)$ is positive definite, e.g., by adopting suitable quasi-Newton techniques, the matrix in (10) is symmetric quasi-definite and can be efficiently factorized with tailored linear algebra routines [30].

In this section we analyze the asymptotic properties of the iterates generated by Algorithm 3.1 under the following blanket assumptions:

1. (a1)

   Functions $f\colon\mathbb{R}^{n}\to\mathbb{R}$ and $c\colon\mathbb{R}^{n}\to\mathbb{R}^{m}$ in (P) are continuously differentiable.

2. (a2)

   Subproblems (6) are well-posed for all parameters’ values, namely for any $\mu_{k}\leq\mu_{0}$, $\rho_{k}\leq\rho_{0}$, and $\hat{y}^{k}\in Y$.

First, we characterize the iterates in terms of stationarity.

Patterning [12, Thm 4.2(ii)], we establish asymptotic complementarity.

Like all penalty-type methods in the nonconvex setting, Algorithm 3.1 may generate limit points that are infeasible for (P). Patterning standard arguments, the following result gives sufficient conditions for the feasibility of limit points; cf. [6, Ex. 4.12].

These conditions are generally difficult to check a priori. Nevertheless, in the situation where each iterate $x^{k}$ is actually a (possibly inexact) global minimizer of (6), then limit points generated by Algorithm 3.1 have minimum constraint violation and tend to minimize the objective function subject to minimal infeasibility [6, Thm 5.1, Thm 5.3]. In particular, limit points are indeed feasible if (P) admits feasible points. However, these properties cannot be expected by solving the subproblems only up to stationarity. Nonetheless, even in the case where a limit point is not necessarily feasible, the next result shows that it is at least a stationary point for a feasibility problem associated to (P).

Finally, we qualify the output of Algorithm 3.1 in the case of feasible limit points. In particular, it is shown that any feasible limit point is AKKT stationary for (P) in the sense of 2.3. Under some additional boundedness conditions, feasible limit points are KKT stationary, according to 2.2.

Provided that the iterates admit a feasible limit point, finite termination of Algorithm 3.1 with an $\varepsilon$-KKT stationary point can be established as a direct consequence of 4.5.

In this section we test an instance of the proposed regularized interior point approach, denoted RegIP, on the CUTEst benchmark problems [19]. RegIP is compared in terms of robustness against the IP solver Ipopt [32] and the AL solver Percival [15], which is based on a BCL method [8] coupled with a trust-region matrix-free solver [21] for the subproblems. We do not report runtimes nor iteration counts since a fair comparison would require close inspection of heuristics and fallbacks [32, §3].

We implemented RegIP in Julia and set up the numerical experiments adopting the JSO software infrastructure by [23]. The IP solver Ipopt acts as subsolver to execute Algorithm 3.1, warm-started at the current primal $(x^{k-1},y^{k-1})$ and dual $(y^{k-1}$, $z^{k-1})$ estimates. We use its parameter tol to set the (inner) tolerance $\epsilon_{k}$, disabling other termination conditions, and let Ipopt control the barrier parameter as needed to approximately solve the regularized subproblem.²²2Solving a sequence of barrier subproblems may hinder the computational efficiency of RegIP compared to the approach behind Algorithm 3.1, but does not degrade its reliability. Ongoing research focuses on solving (6) and letting the IP subsolver update the barrier parameter after warm-starting at the current primal-dual estimate, in the spirit of [7, Alg. 3]. We let the safeguarding set be $Y:=\{v\in\mathbb{R}^{m}\,|\,\|v\|_{\infty}\leq 10^{20}\}$ and choose $\hat{y}^{k}$ by projecting the current estimate $y^{k-1}$ onto $Y$. We set the initial penalty parameter to $\rho_{0}=10^{-6}$, the inner tolerance $\epsilon_{0}=\sqrt[3]{\varepsilon}$, and parameters $\theta_{\rho}=0.5$, $\kappa_{\rho}=0.5$, and $\kappa_{\epsilon}=0.5$. RegIP declares success, and returns a $\varepsilon$-KKT stationary point, as soon as $\epsilon_{k}\leq\varepsilon$ and $C^{k}\leq\varepsilon$.³³3The condition $\epsilon_{k}\leq\varepsilon$ implies both $\varepsilon$-stationarity (1a) and $\varepsilon$-complementarity (1c) required for $\varepsilon$-KKT stationarity. This follows from the observation that, in RegIP, the subsolver approximately solves (4) at Algorithm 3.1, not (6); see Footnote 2. Instead, if $\epsilon_{k}\leq\varepsilon$, $C^{k}>\varepsilon$ and $\rho_{k}\leq\rho_{\min}:=10^{-20}$, RegIP stops declaring (local) infeasibility. For Ipopt, we set the tolerance tol to $\varepsilon$, remove the other desired thresholds, and disable termination based on acceptable iterates. For Percival, we set absolute and relative tolerances atol, rtol, and ctol to $\varepsilon$. We select all CUTEst problems with at most 1000 variables and constraints, obtaining a test set with 609 problems. All solvers are provided with the primal-dual initial point available in CUTEst, a time limit of $60$ seconds, the maximum number of iterations set to $10^{9}$, and a tolerance $\varepsilon\in\{10^{-3},10^{-5}\}$. A solver is deemed to solve a problem if it returns with a successful status; it fails otherwise. The source codes for the numerical experiments have been archived on Zenodo at

doi: 10.5281/zenodo.7109904.

Table 1 summarizes the results, stratified by solver, termination tolerance ($\varepsilon$) and problem size ($n$, $m$). For each combination, we indicate the number of times RegIP solves a problem that the other solver fails (“W”) or solves (“T+”) and the number of times RegIP fails on a problem that the other one fails (“T-”) or solves (“L”). The results show that RegIP succeeds on more problems than the other solvers, consistently for both low and high accuracy, indicating that the underlying regularized IP approach can form the basis for reliable and scalable solvers.⁴⁴4We agree with [5, §5] on the fact that “strong statements concerning the relative efficiency or robustness [$\ldots$] are not possible in nonlinear optimization.”

This paper has presented a regularized interior point approach to solving constrained nonlinear optimization problems. Operating as an outer regularization layer, a quadratic proximal penalty provides robustness whilst consuming minimal computation effort once embedded into existent interior point codes as a principled inertia correction strategy. Furthermore, regularizing the equality constraints allows to safely adopt more efficient linear algebra routines, while waiving the need for an infeasibility detection mechanism within the subsolver. Preliminary numerical results indicate that a close integration of proximal regularization within interior point schemes is key to provide efficient and robust solvers. We encourage further research in this direction.