Adv Syst Sci Appl 2022; 02:73–84 Published online at https://ijassa.ipu.ru. Accelerating Sequential Quadratic Programming for Inequality-Constrained Optimization near Critical Lagrange Multipliers Alexey F. Izmailov1*, Ivan S. Rodin1 1Lomonosov Moscow State University, VMK Faculty, IO Department, Moscow, Russia Abstract: We consider a sequential quadratic programming algorithm for optimization problems with equality and inequality constraints, equipped with the standard Armijo linesearch procedure for a nonsmooth exact penalty function, intended for globalization of convergence. We are interested in the case when the standard assumptions for local superlinear convergence of the method may not hold. Specifically, we allow for violation of standard constraint qualifications and second-order sufficient optimality conditions, in which case attraction to so-called critical Lagrange multipliers is known to have a negative impact on convergence rate. In these circumstances, some known acceleration techniques can be expected to take effect only provided the true Hessian and the full SQP step are asymptotically accepted, and these are the main issues addressed in this work. The presented constructions extend some previously known ones to the case when inequality constraints are involved. Keywords: constrained optimization, Karush–Kuhn–Tucker optimality system, critical Lagrange multiplier, 2-regularity, Newton-type methods, sequential quadratic programming, linesearch globalization of convergence, merit function, nonsmooth exact penalty function, true Hessian, unit stepsize, extrapolation 1. INTRODUCTION In this work we consider the equality- and inequality-constrained optimization problem minimize f(x) subject to h(x) = 0, g(x) ≤ 0, (1.1) where the objective function f : Rn → R and the constraint mappings h : Rn → Rl and g : Rn → Rm are sufficiently smooth, and we are concerned with some crucial peculiarities of performance of Newton-type methods near solutions or stationary points of problem (1.1) that are in a sense singular. Let L : Rn × Rl × Rm → R be the Lagrangian of problem (1.1), i.e., L(x, λ) = f(x) + 〈λ, h(x)〉+ 〈µ, g(x)〉, where 〈·, ·〉 stands for the Euclidian inner product. Then the stationary points and associated Lagrange multipliers of problem (1.1) are characterized by the Karush–Kuhn–Tucker (KKT) optimality system ∂L ∂x (x, λ, µ) = 0, h(x) = 0, µ ≥ 0, g(x) ≤ 0, 〈µ, g(x)〉 = 0. (1.2) ∗Corresponding author: izmaf@cs.msu.ru 74 A. F. IZMAILOV, I. S. RODIN For a feasible point x̄ of problem (1.1), the linear independence constraint qualification (LICQ) consists of saying that the gradients h′j(x̄), j ∈ {1, . . . , l}, g′i(x̄), i ∈ A(x̄), are linearly independent, where A(x̄) = {i ∈ {1, . . . , m} | gi(x̄) = 0} is the set of indices of inequality constraints active at x̄. If x̄ is a local solution of (1.1), satisfying LICQ, then there exists the unique pair (λ̄, µ̄) ∈ Rl × Rm such that the triple (x̄, λ̄, µ̄) solves (1.2). In this work, we are mostly interested in the case when LICQ does not hold, but x̄ is stationary for (1.1) with some (possibly nonunique) associated Lagrange multiplier. In this case, the set of Lagrange multipliers may naturally contain some special instances, called critical multipliers (to be defined formally in Section 2 below), that are known to strongly attract dual sequences of various primal-dual optimization algorithms, and this phenomenon has a strong negative impact on convergence rate; see [14–16] and the summaries of this research in [18, Section 7.1] and [19]. The effect of attraction to critical multipliers can be locally avoided by using some dual stabilization mechanisms, like the one of the stabilized sequential quadratic programming method; see the original proposals in [22, 24], as well as a more recent treatment in [17] and [18, Section 7.2.2]. However, even such modified algorithms typically have large domains of convergence to critical multipliers when they exist [11], and this becomes even more of an issue when globalization of convergence is concerned. Generally, there can be at least two different approaches to the specified criticality issue. One is indeed to try to avoid convergence to a critical multiplier, and some further tools for this were developed very recently in [5]. The essence of an alternative approach, adopted in this work, can be expressed as follows: since it is difficult to avoid convergence to critical multipliers, this convergence can be accelerated by employing the knowledge of its rather special character. This point will be explained in Section 2, first for generic systems of nonlinear equations and for optimization problems with equality constraints only, using the results from [6–9] and [10], respectively. One of the outcomes of this analysis is a very simple extrapolation strategy for accelerating convergence of the basic sequential quadratic programming (SQP) method to critical Lagrange multipliers, that can be easily incorporated into globalization schemes based on linesearch. The main goal of this paper is to extend this approach to optimization problems involving inequality constraints as well, and this is accomplished in Section 3. Section 4 presents some numerical results confirming the use of the approach, and Section 5 contains some concluding remarks and discussion of directions for further research. Some words about our notation and terminology. The Euclidian (l2), l1, and l∞ norms will be denoted by ‖ · ‖, ‖ · ‖1, and ‖ · ‖∞, respectively. For a function ϕ : Rn → R, by ϕ′(x; ξ) we denote the standard directional derivative of ϕ at x ∈ Rn in a direction ξ ∈ Rn. The null space and the range space of a linear operator Λ : Rn → Rl are denoted by ker Λ and im Λ, respectively. For a given z ∈ Rm and an index set I ⊂ {1, . . . , m}, the symbol zI stands for the subvector of z with components corresponding to i ∈ I . A set U ⊂ Rp is called starlike with resect to ū ∈ Rp if for every u ∈ U and t ∈ (0, 1], it holds that tu+ (1− t)ū ∈ U . Any v ∈ Rp is called an excluded direction for a set U starlike with respect to ū if ū+ tv 6∈ U for all t > 0. A set which is starlike with respect to a given point is called asymptotically dense if the complement of the corresponding set of excluded directions is open and dense. 2. NONLINEAR EQUATIONS AND EQUALITY-CONSTRAINED OPTIMIZATION Consider a system of nonlinear equations Φ(u) = 0 (2.3) with a smooth mapping Φ : Rp → Rp. A solution ū of (2.3) is called (non)singular if the Jacobian Φ′(ū) is a (non)singular matrix. For a given current iterate uk ∈ Rp, the basic Copyright © 2022 ASSA. Adv Syst Sci Appl (2022) ACCELERATING SEQUENTIAL QUADRATIC PROGRAMMING FOR INEQUALITY-CONSTRAINED 75 Newton method (NM) defines the next iterate as uk+1 = uk + vk, where vk is a solution of the linearized (at uk) equation Φ(uk) + Φ′(uk)v = 0. (2.4) The classical theory says that when u0 ∈ Rp is taken close enough to a nonsingular solution of (2.3), then the NM defined this way generates the unique sequence of iterates {uk}, and this sequence converges to ū superlinearly; see, e.g., [18, Section 2.1.1]. Here we are mainly interested in the cases of when the solution in question can be singular, perhaps even nonisolated. To that end, we make use of the following second-order regularity concept assuming that Φ is twice differentiable at ū: the mapping Φ is said to be 2-regular at ū in a direction v ∈ Rp if im Φ′(ū) + Φ′′(ū)[v, ker Φ′(ū)] = Rp. If Φ′(ū) is nonsingular, then Φ is 2-regular at ū in every direction, including v = 0. At the same time, 2-regularity may hold at singular (and even nonisolated) solutions in nonzero directions, including those from ker Φ′(ū). The following theorem summarizes the results on local convergence of the basic NM, obtained in [8, Theorem 6.1] and [9, Theorem 2.1]. Theorem 2.1: Let Φ : Rp → Rp be twice differentiable near ū ∈ Rp, with its second derivative Lipschitz- continuous with respect to ū, that is, Φ′′(u)− Φ′′(ū) = O(‖u− ū‖) as u→ ū. Let ū be a singular solution of equation (2.3), and assume that there exists v̄ ∈ ker Φ′(ū) such that Φ is 2-regular at ū in this direction. Then there exists a set U ⊂ Rp starlike with respect to ū, asymptotically dense at ū, and such that for every starting point u0 ∈ U , there exists the unique sequence {uk} ⊂ Rp such that for all k, the step vk = uk+1 − uk is the solution of (2.4), and uk 6= ū, {uk} converges to ū, lim k→∞ ‖uk+1 − ū‖ ‖uk − ū‖ = 1 2 , (2.5) and the sequence {(uk − ū)/‖uk − ū‖} converges to some v ∈ ker Φ′(ū). The convergence pattern specified by (2.5) suggests that acceleration of convergence can be achieved by doubling the basic NM step. This, however, should be done with due care, as simply doubling each step may force the resulting sequence to leave the domain of convergence, which is asymptotically dense but need not contain a full neighborhood of the solution in question. To that end, somehow more involved acceleration techniques have been developed in [7,9], such as extrapolation and overrelaxation. The simplest extrapolation technique consists of generating, along with the basic NM sequence {uk}, an auxiliary sequence {ûk} obtained by doubling the NM step: for each k compute ûk+1 = uk + 2vk. (2.6) According to [9, Theorem 4.1], under the assumptions of Theorem 2.1 above, this auxiliary sequence {ûk} converges linearly to ū, with the asymptotic ratio of 1/4 (instead of 1/2 for {uk}, given by (2.5)), from all points in the domain of convergence of NM sequences {uk}. We now turn our attention to globalization issues, and recall that convergence of the NM can be globalized in standard ways, e.g., by the Armijo linesearch procedure for the residual ‖Φ(·)‖ used as a merit function; see [18, Section 5.1.1]. The following is a prototype algorithm of this kind. Copyright © 2022 ASSA. Adv Syst Sci Appl (2022) 76 A. F. IZMAILOV, I. S. RODIN Algorithm 2.1: Fix the parameters σ ∈ (0, 1) and θ ∈ (0, 1). Choose u0 ∈ Rp, and set k = 0. 1. Compute vk ∈ Rp solving the linear system (2.4). If (2.4) cannot be solved, stop with failure. 2. Set α = 1. If the inequality ‖Φ(uk + αvk)‖ ≤ (1− σα)‖Φ(uk)‖ (2.7) holds, set αk = α and go to Step 3. Otherwise, keep replacing α by θα until (2.7) is satisfied. 3. Set uk+1 = uk + αkv k, increase k by 1, and go to Step 1. According to [18, Theorem 5.3], this algorithm converges globally in a sense that any accumulation point ū of any sequence {uk} generated by this algorithm satisfies (Φ′(ū))>Φ(ū) = 0. Moreover, according to [18, Theorem 5.4], Step 2 of Algorithm 2.1 accepts the unit stepsize when uk is close to a nonsingular solution, thus allowing the globalized algorithm to inherit the superlinear convergence rate of the basic NM. Near singular solutions, the crucial issue of asymptotic acceptance of the unit stepsize is much more involved. The following result was obtained in [6, Theorem 1]. Theorem 2.2: Under the assumptions of Theorem 2.1, a set U in it can be chosen so that for every starting point u0 ∈ U , Algorithm 2.1 with σ ∈ (0, 3/4) uniquely defines the sequence {uk}, and αk = 1 is accepted at Step 2 of this algorithm for all k large enough. Combining Theorems 2.1 and 2.2, one can conclude that once a sequence generated by Algorithm 2.1 enters the domain U specified in those theorems, this sequence converges linearly to ū with the asymptotic ratio of 1/2. Moreover, for the “extrapolated” sequence {ûk} generated according to (2.6), the asymptotic ratio is 1/4. It is important to note that the main iterative sequences {uk} are not affected at all by computing the auxiliary sequences {ûk}, and hence, doing so does not affect the global convergence properties of the algorithm. It should also be mentioned that generating the auxiliary sequence {ûk} costs essentially nothing. Getting back to optimization, we first consider the case when there are no inequality constraints in (1.1) (i.e., m = 0): minimize f(x) subject to h(x) = 0. (2.8) In this case, the KKT system (1.2) can be written as (2.3) with p = n+ l, Φ : Rn × Rl → Rn × Rl, Φ(u) = ( ∂L ∂x (x, λ), h(x) ) , (2.9) where u = (x, λ), and the missing argument µ of L is dropped. Algorithm 2.1, as well as its version with extrapolation, and Theorems 2.1 and 2.2 are certainly applicable to this special case of a nonlinear equation (2.3). The problem, however, is that in the optimization context, using the residual of a first-order optimality system as a merit function can hardly be justified: a merit function should be reflecting the intention to find a solution of the optimization problem rather than just any stationary point of it. One typical choice of an optimization-related merit function for problem (2.8) is the l1 exact penalty function ϕc : Rn → R, ϕc(x) = f(x) + c‖h(x)‖1, Copyright © 2022 ASSA. Adv Syst Sci Appl (2022) ACCELERATING SEQUENTIAL QUADRATIC PROGRAMMING FOR INEQUALITY-CONSTRAINED 77 where c > 0 is a penalty parameter; see [1, Section 17], [23, Section 18.4], [18, Section 6.2]. Observe, however, that after this change of the merit function, Theorem 2.2 is no more applicable, and the issue of asymptotically accepting the unit stepsize becomes even more involved, in particular, due to inherent nonsmoothness of ϕc, and to the presence of a penalty parameter c. Another source of complications is that the primal part ξk of the NM step vk = (ξk, ηk) computed at the current iterate uk = (xk, λk) can be guaranteed to be a direction of descent for ϕc (with an appropriate choice of c) at uk only provided the Hessian of the Lagrangian Hk = ∂2L ∂x2 (xk, λk) (2.10) is positive definite, which is not at all automatic even close to a solution/multiplier pair satisfying the second-order sufficient optimality condition. To that end, the algorithm for problem (2.8), presented next, has Step 3, where the sufficient descent property of the generated primal direction is verified. Algorithm 2.2: Fix the parameters c̄ > 0, c̃ > 0, ρ > 0 and σ, θ ∈ (0, 1). Choose u0 = (x0, λ0) ∈ Rn × Rl, and set c−1 = 0 and k = 0. 1. Compute vk = (ξk, ηk) solving the linear system ∂L ∂x (xk, λk) +Hkξ + (h′(xk))>η = 0, h(xk) + h′(xk)ξ = 0, (2.11) with Hk given by (2.10). If (2.11) cannot be solved, stop with failure. 2. Set ck = max { ck−1, 4(1− σ)‖λk + ηk‖∞ + ‖λk‖∞ 3− 4σ + c̄ } . (2.12) If maximum in (2.12) is attained at the second argument, replace ck by ck + c̃. 3. If the inequality ϕ′ck(xk; ξk) ≤ −ρ‖ξk‖2 (2.13) is violated, stop with failure. 4. Set α = 1. If the inequality ϕck(xk + αξk) ≤ ϕck(xk) + σαϕ′ck(xk; ξk) (2.14) holds, set αk = α and go to Step 5. Otherwise, keep replacing α by θα until (2.14) is satisfied. 5. Define uk+1 = (xk+1, λk+1) as uk+1 = uk + αkv k, increase k by 1, and go to Step 1. Observe that (2.11) with Hk given by (2.10) is nothing else but the NM equation (2.4) with v = (ξ, η), for Φ defined in (2.9). On the other hand, this method can be seen as the SQP algorithm, since (2.11) is the Lagrange optimality system for the equality-constrained quadratic programming subproblem minimize 〈f ′(xk), ξ〉+ 1 2 〈Hkξ, ξ〉 subject to h(xk) + h′(xk)ξ = 0; (2.15) for details see, e.g., [18, Section 4.2]. Asymptotic acceptance of the true Hessian and unit stepsize by Algorithm 2.2 has been investigated very recently in [10]. In order to present the main result of that work, some more terminology will be needed. According to [11, Proposition 1] and [12, Proposition 2], the Copyright © 2022 ASSA. Adv Syst Sci Appl (2022) 78 A. F. IZMAILOV, I. S. RODIN assumptions of Theorem 2.1 for Φ defined in (2.9) may hold at ū = (x̄, λ̄) only if λ̄ is a critical Lagrange multiplier associated with a stationary point x̄ of problem (2.8). Criticality of λ̄ means that the linear subspace Q(x̄, λ̄) = { ξ ∈ kerh′(x̄) ∣∣∣∣ ∂2L ∂x2 (x̄, λ̄)ξ ∈ im(h′(x̄))> } (2.16) is nontrivial; see [18, Definition 1.41]. In particular, the assumptions of Theorem 2.1 do hold under the following requirements: • The multiplier µ̄ is critical of order 1, which means that dimQ(x̄, λ̄) = 1, or, in other terms, Q(x̄, λ̄) = span{ξ̄} (2.17) with some ξ̄ ∈ Rn \ {0}. • It holds that rankh′(x̄) = l − 1 (2.18) and h′′(x̄)[ξ̄, ξ̄] 6∈ imh′(x̄). (2.19) The following counterpart of Theorem 2.2 for equality-constrained optimization problems appears in [10, Theorem 3.1]. Theorem 2.3: Let f : Rn → R and h : Rn → Rl be three times differentiable near x̄ ∈ Rn, with their third derivatives Lipschitz-continuous with respect to x̄. Let x̄ be a stationary point of problem (2.8), with an associated Lagrange multiplier λ̄ ∈ Rl, and assume that (2.17) with Q(x̄, λ̄) defined in (2.16), (2.18), and (2.19) hold with some ξ̄ ∈ Rn \ {0}. Assume also that h′(x̄) 6= 0, and set ū = (x̄, λ̄). Then there exist ρ > 0 and a set U ⊂ Rn × Rl starlike with respect to ū, asymptotically dense at ū, and such that for every starting point u0 ∈ U , Algorithm 2.2 with σ ∈ (0, 3/4) uniquely defines the sequence {uk}, this sequence converges to ū, and αk = 1 is accepted at Step 4 of this algorithm for all k large enough. Getting back to a general optimization problems involving inequality constraints, we recall that according to [18, Definition 7.8], a Lagrange multiplier (λ̄, µ̄) associated with a stationary point x̄ of problem (1.1) is called critical with respect to an index set A ⊂ A(x̄) if µ̄A(x̄)\A = 0 and (λ̄, µ̄A) is a critical Lagrange multiplier associated with the stationary point x̄ of the equality-constrained optimization problem minimize f(x) subject to h(x) = 0, gA(x) = 0. (2.20) As discussed in [15] and [18, Section 7.1.2], the effect of criticality on performance of primal- dual algorithms shows up precisely this way, through stabilization of specific activity patterns. This is confirmed by the numerical experiments discussed in Section 4 below. 3. OPTIMIZATION WITH INEQUALITY CONSTRAINTS Unlike Algorithms 2.1 and 2.2, the algorithm for problem (1.1) that we present next is not just a prototype algorithm in the following sense: it is supplied by a mechanism intended to avoid unnecessary terminations with failure, especially far from solutions, and it is exactly what we used in numerical experiments in Section 4. More precisely, the issue of possible infeasibility of the iteration subproblem (3.23) (a counterpart of (2.15) involving linearized Copyright © 2022 ASSA. Adv Syst Sci Appl (2022) ACCELERATING SEQUENTIAL QUADRATIC PROGRAMMING FOR INEQUALITY-CONSTRAINED 79 inequality constraints) is not addressed in this algorithm, and it still stops with failure in such cases, as this is a general issue concerned with SQP methods, and there are known tools for tackling it; see, e.g., [18, Section 6.2]. However, when the iteration subproblem is feasible but cannot be solved, or when the generated primal direction ξk does not pass the sufficient descent test (3.25) (a counterpart of (2.13)), the algorithm does not stop with failure, but rather sequentially adjusts the basic choice of Hk = ∂2L ∂x2 (xk, λk, µk) (3.21) in such a way that eventually the needed ξk satisfying (3.25) is always found. For every value of the penalty parameter c > 0, we now need to re-define the penalty function ϕc : Rn → R in order to incorporate the inequality constraints: ϕc(x) = f(x) + c(‖h(x)‖1 + ‖max{0, g(x)}‖1). (3.22) Algorithm 3.1: Fix the parameters c̄ > 0, c̃ > 0, ρ > 0 and σ, θ ∈ (0, 1). Choose u0 = (x0, λ0, µ0) ∈ Rn × Rl, and set c−1 = 0 and k = 0. 1. Define Hk according to (3.21). 2. Compute ξk as a stationary point of the quadratic programming problem minimize 〈f ′(xk), ξ〉+ 1 2 〈Hkξ, ξ〉 subject to h(xk) + h′(xk)ξ = 0, g(xk) + g′(xk)ξ ≤ 0, (3.23) and (λk+1, µk+1) as an associated Lagrange multiplier. If (3.23) is infeasible, stop with failure. If (3.23) is feasible, but a stationary point of (3.23) cannot be found, go to Step 5. Otherwise, set ηk = λk+1 − λk, ζk = µk+1 − µk, vk = (ξk, ηk, ζk). 3. Set ck = max { ck−1, 4(1− σ)‖(λk+1, µk+1)‖∞ + ‖(λk, µk)‖∞ 3− 4σ + c̄ } . (3.24) If max in (3.24) is attained at the second argument, replace ck by ck + c̃. Set ∆k = 〈f ′(xk), ξk〉 − ck(‖h(xk)‖1 + ‖max{0, g(xk)}‖1). 4. If the inequality ∆k ≤ −ρ‖ξk‖2 (3.25) is satisfied, go to Step 6. 5. Choose τk > 0, replace Hk by Hk + τkI , and go to Step 2. 6. Set α = 1. If the inequality ϕck(xk + αξk) ≤ ϕck(xk) + σα∆k (3.26) holds with ϕck defined according to (3.22), set αk = α and go to Step 7. Otherwise, keep replacing α by θα until (3.26) is satisfied. 7. Define uk+1 = (xk+1, λk+1, µk+1) as uk + αkv k, increase k by 1, and go to Step 1. In Step 7 of this algorithm, it is quite typical to use the stepsize parameter in primal updates only; see, e.g., [18, Algorithm 6.7]. However, the theory developed in [10] for the equality-constrained case requires using the stepsize parameter in dual updates as well. This Copyright © 2022 ASSA. Adv Syst Sci Appl (2022) 80 A. F. IZMAILOV, I. S. RODIN variant of linesearch SQP algorithms is also quite common; see, e.g., [23, Algorithm 18.3]. Theoretical global convergence properties of this kind of algorithms, as well as their rate of convergence properties under “standard” assumptions were recently discussed in [20]; see also references therein. Here, however, we are concerned with the cases when those “standard” assumptions for fast local convergence may fail to hold, and in particular, with the cases of convergence to multipliers critical with respect to a relevant index set. Motivated by its success in the equality-constrained case, our proposal being tested in this work is to generate an extrapolated auxiliary sequence {ûk} according to (2.6), with expectation that its convergence will be faster than that of the main sequence {uk}. 4. NUMERICAL EXPERIMENTS Algorithm 3.1 will be abbreviated as SQP-mH, while its variant supplied with extrapolation, and hence, generating an auxiliary sequence {ûk} according to (2.6), will be abbreviated as SQP-mH-EP. Note that (2.6) always employs the direction vk obtained by the basic choice of Hk specified in (3.21). If for some k such vk cannot be defined, we put ûk+1 = uk+1. Instead of modifying the Hessian of the Lagrangian, possible lack of positive definiteness of the true Hessian can be tackled by defining Hk as quasi-Newton approximations of the Hessian. A standard choice adopted here consists of using BFGS approximations complemented by Powell’s correction at Step 1 of Algorithm 3.1 instead of (3.21); see [18, Section 4.1] for details. Furthermore, in the resulting algorithm abbreviated as SQP-BFGS, Steps 4 and 5 of Algorithm 3.1 are dropped, and the stepsize parameter is used for primal updates only, i.e., λk+1 = λk + ηk and µk+1 = µk + ζk. The parameter values used in our computations were as follows: c̄ = c̃ = 1, ρ = 10−9, σ = 0.01, θ = 0.5. Furthermore, we adopted the following rule for controlling τk at Step 5 of Algorithm 3.1: for every given k, when Step 5 is invoked for the first time, we set τk = 1, and then multiply it by 2 every next time Algorithm 3.1 invokes Step 5. Runs of algorithms SQP-BFGS and SQP-mH were terminated with success if a newly generated iterate uk satisfied ‖Φ(uk)‖ ≤ 10−6. (4.27) For SQP-mH-EP, for every k = 1, 2, . . . we first compute ûk, and terminate the run with success if ‖Φ(ûk)‖ ≤ 10−6; (4.28) otherwise, we proceed with computing uk and verifying (4.27). When successful termination did not occur after 200 iterations, or the backtracking procedure at Step 6 of Algorithm 3.1 produced a trial value α such that α‖vk‖ ≤ 10−10, the process was terminated with failure. The experiments were performed in Matlab environment, with PATH solver [3, 4] used for quadratic programming subproblems. Problem instances for the experiments were taken from DEGEN test collection [13] that contains optimization problems with constraints violating LICQ. We employed all problems from DEGEN involving inequality-constraints (numerical results for purely equality-constrained problems can be found in [10]), except for those where the solutions of interest is not a stationary point (problems 20112, 20220, 20224, 30212 and 30302). This leaves 59 test problems. For each of those problems, all algorithms being tested were initialized at the same 100 starting points u0 = (x0, λ0, µ0) generated randomly in the l∞-ball of radius 100, centered at the primal solution of interest (reported in DEGEN) for the primal part x0, and at (0, 0) for the dual part (λ0, µ0), with the additional nonnegativity restriction on µ0. Copyright © 2022 ASSA. Adv Syst Sci Appl (2022) ACCELERATING SEQUENTIAL QUADRATIC PROGRAMMING FOR INEQUALITY-CONSTRAINED 81 As measures of efficiency, we used the average number of quadratic programming subproblems solved and the average number of evaluations of the objective function, constraint mappings, and their derivatives, per one successful run. These results are presented in the form of performance profiles [2]. For DEGEN test set, we employed a modified construction of performance profiles, intended for the case of multiple runs for each test problems; see [21] for details. 2 4 6 8 10 12 14 16 18 20 22 0 0.2 0.4 0.6 0.8 1 SQP-mH SQP-mH-EP BFGS-SQP (a) Performance profile 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.2 0.4 0.6 0.8 1 SQP-mH SQP-mH-EP (b) Subproblems per iteration Fig. 4.1. Subproblems solved The performance profile in Figure 4.1a demonstrates that SQP-mH-EP is as robust as SQP-mH, somehow more robust than SQP-BFGS, and by far more efficient than both by the number of quadratic subproblems solved per a successful run. The effect of using extrapolation is achieved because the true Hessian and the unit stepsize are typically asymptotically accepted, and the set of active constraints of subproblems asymptotically stabilizes. In addition, the plots in Figure 4.1b show which portion of problems required solving no more than a given number of quadratic subproblems per iteration on the average. It can be seen that the behavior of SQP-mH and SQP-mH-EP in this respect is quite similar. In particular, for 40% of test problems, only one quadratic subproblem per iteration had to be solved on the average, while solving two or more subproblems was needed for less than 10% of problems. Copyright © 2022 ASSA. Adv Syst Sci Appl (2022) 82 A. F. IZMAILOV, I. S. RODIN 2 4 6 8 10 12 14 16 18 20 22 0 0.2 0.4 0.6 0.8 1 SQP-mH SQP-mH-EP BFGS-SQP (a) Evaluations of f 2 4 6 8 10 12 14 16 18 20 22 0 0.2 0.4 0.6 0.8 1 SQP-mH SQP-mH-EP BFGS-SQP (b) Evaluations of h and g 2 4 6 8 10 12 14 16 18 20 22 0 0.2 0.4 0.6 0.8 1 SQP-mH SQP-mH-EP BFGS-SQP (c) Evaluations of derivatives Fig. 4.2. Performance profiles by evaluations The performance profiles in terms of evaluations of f , of the constraint mappings, and of their derivatives are presented in Figures 4.2a, 4.2b and 4.2c, respectively. According to Copyright © 2022 ASSA. Adv Syst Sci Appl (2022) ACCELERATING SEQUENTIAL QUADRATIC PROGRAMMING FOR INEQUALITY-CONSTRAINED 83 Figure 4.2a, the relative efficiency of SQP-mH-EP by evaluations of f is lower than that by the count of subproblems solved, but the positive effect of using extrapolation is still evident. The picture in Figure 4.2b by evaluations of h and g, and in Figure 4.2c by evaluations of the derivatives is quite different though from that in Figures 4.1a and 4.2a: SQP-mH-EP demonstrated the best results by these measures of efficiency for about 10% of problems only. This behavior is explained by the fact that even if the true Hessian and the full step are typically asymptotically accepted, and the set of active constraints stabilizes, additional evaluations of h, g, and derivatives are required at extrapolated points for verifying the stopping test (4.28). 5. CONCLUDING REMARKS We have discussed the use of extrapolation techniques for acceleration of convergence of the SQP method equipped with linesearch, when convergence is to a critical Lagrange multiplier of an optimization problem involving inequality constraints. A specific algorithmic implementation of these constructions has been developed and numerically tested. It would be interesting to develop a reasonably complete supporting theory, i.e., conditions ensuring asymptotic acceptance of the true Hessian and the unit stepsize, as it is known to be possible for equality-constrained problems, and hopefully can be achieved by some kind of asymptotic reduction to the equality-constrained case. Observe, however, that the key difficulty for such development apparently consists of establishing the relation between the values and directional derivatives of the penalty functions for (1.1) and for (2.20) in the SQP directions. ACKNOWLEDGEMENTS This research was supported in part by the Russian Foundation for Basic Research Grants 19-51-12003 NNIO a and 20-01-00106 a, and by the Volkswagen Foundation. REFERENCES 1. Bonnans, J.F., Gilbert, J.Ch., Lemaréchal, C. & Sagastizábal, C. (2006) Numerical optimization. Theoretical and practical aspects. 2nd edition. Berlin: Springer-Verlag. 2. Dolan, E.D. & More, J.J. (2002) Benchmarking optimization software with performance profiles. Math. Program. 91, 201–213. 3. Ferris, M.C. & Munson, T.S.. PATH. [Online]. Available http://pages.cs.wisc.edu/˜ferris/path.html. 4. Ferris, M.C. & Munson, T.S. (1999) Interfaces to PATH 3.0: design, implementation and usage. Comput. Optim. Appl. 12, 207–227. 5. Fischer, A., Izmailov, A.F. & Scheck, W. (2020) Adjusting dual iterates in the presence of critical Lagrange multipliers. SIAM J. Optim. 30, 1555–1581. 6. Fischer, A., Izmailov, A.F. & M.V. Solodov. (2021) Unit stepsize for the Newton method close to critical solutions. Math. Program. 187, 697–721. 7. Griewank, A. (1980) Analysis and modification of Mewton’s method at singularities. PhD Thesis. Australian National University, Canberra. 8. Griewank, A. (1980) Starlike domains of convergence for Newton’s method at singularities. Numer. Math. 35, 95–111. 9. A. Griewank. (1985) On solving nonlinear equations with simple singularities or nearly singular solutions. SIAM Rev. 27, 537–563. Copyright © 2022 ASSA. Adv Syst Sci Appl (2022) 84 A. F. IZMAILOV, I. S. RODIN 10. Izmailov, A.F. (2021) Accelerating convergence of a globalized sequential quadratic programming method to critical Lagrange multipliers. Comput. Optim. Appl. 80, 943– 978. 11. Izmailov, A.F., Kurennoy, A.S. & Solodov, M.V. (2018) Critical solutions of nonlinear equations: Local attraction for Newton-type methods. Math. Program. 167, 355–379. 12. Izmailov, A.F., Kurennoy, A.S. & Solodov, M.V. (2018) Critical solutions of nonlinear equations: Stability issues. Math. Program. 168, 475–507. 13. Izmailov, A.F. & Solodov, M.V. (2009) DEGEN. [Online]. Available http://w3.impa.br/˜optim/DEGEN collection.zip. 14. Izmailov, A.F. & Solodov, M.V. (2009) On attraction of Newton-type iterates to multipliers violating second-order sufficiency conditions. Math. Program. 117, 271– 304. 15. Izmailov, A.F. & Solodov, M.V. (2009) Examples of dual behaviour of Newton-type methods on optimization problems with degenerate constraints. Comput. Optim. Appl. 42, 231–264. 16. Izmailov, A.F. & Solodov, M.V. (2011) On attraction of linearly constrained Lagrangian methods and of stabilized and quasi-Newton SQP methods to critical multipliers. Math. Program. 126, 231–257. 17. Izmailov, A.F. & Solodov, M.V. (2012) Stabilized SQP revisited. Math. Program. 133, 93–120. 18. Izmailov, A.F. & Solodov, M.V. (2014) Newton-Type Methods for Optimization and Variational Problems. Cham, Switzerland: Springer Series in Operations Research and Financial Engineering, Springer International Publishing. 19. Izmailov, A.F. & Solodov, M.V. (2015) Critical Lagrange multipliers: what we currently know about them, how they spoil our lives, and what we can do about it. TOP 23, 1–26. 20. Izmailov, A.F. & Solodov, M.V. (2016) Some new facts about sequential quadratic programming methods employing second derivatives. Optim. Meth. Software 31, 1111–1131. 21. Izmailov, A.F., Solodov, M.V. & Uskov, E.I. (2015) Combining stabilized SQP with the augmented Lagrangian algorithm Comput. Optim. Appl. 62, 405–429. 22. Li, D.-H. & Qi, L. (2000) Stabilized SQP method via linear equations. Applied Mathematics Technical Reptort AMR00/5. University of New South Wales, Sydney. 23. Nocedal, J. & Wright, S.J. (2006) Numerical Optimization. 2nd edition. New York, Berlin, Heidelberg: Springer-Verlag. 24. Wright, S.J. (1998) Superlinear convergence of a stabilized SQP method to a degenerate solution. Comput. Optim. Appl. 11, 253–275. Copyright © 2022 ASSA. Adv Syst Sci Appl (2022) Introduction Nonlinear equations and equality-constrained optimization Optimization with inequality constraints Numerical experiments Concluding remarks