Convergence of successive linear programming algorithms for noisy functions
Abstract
EconStor is a publication server for scholarly economic literature, provided as a non-commercial public service by the ZBW.
Full text
Hansknecht, Christoph; Kirches, Christian; Manns, Paul Article — Published Version Convergence of successive linear programming algorithms for noisy functions Computational Optimization and Applications Provided in Cooperation with: Springer Nature Suggested Citation: Hansknecht, Christoph; Kirches, Christian; Manns, Paul (2024) : Convergence of successive linear programming algorithms for noisy functions, Computational Optimization and Applications, ISSN 1573-2894, Springer US, New York, NY, Vol. 88, Iss. 2, pp. 567-601, https://doi.org/10.1007/s10589-024-00564-w This Version is available at: https://hdl.handle.net/10419/315241 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. http://creativecommons.org/licenses/by/4.0/
Computational Optimization and Applications (2024) 88:567–601 https://doi.org/10.1007/s10589-024-00564-w Convergence of successive linear programming algorithms for noisy functions Christoph Hansknecht1·Christian Kirches1·Paul Manns2 Received: 16 February 2023 / Accepted: 30 January 2024 / Published online: 26 February 2024 © The Author(s) 2024 Abstract Gradient-based methods have been highly successful for solving a variety of both unconstrained and constrained nonlinear optimization problems. In real-world applications, such as optimal control or machine learning, the necessary function and derivative information may be corrupted by noise, however. Sun and Nocedal have recently proposed a remedy for smooth unconstrained problems by means of a stabilization of the acceptance criterion for computed iterates, which leads to convergence of the iterates of a trust-region method to a region of criticality (Sun and Nocedal in Math Program 66:1–28, 2023. https://doi.org/10.1007/s10107-023-01941-9). We extend their analysis to the successive linear programming algorithm (Byrd et al. in Math Program 100(1):27–48, 2003. https://doi.org/10.1007/s10107-003-0485-4, SIAM J Optim 16(2):471–489, 2005. https://doi.org/10.1137/S1052623403426532) for unconstrained optimization problems with objectives that can be characterized as the composition of a polyhedral function with a smooth function, where the latter and its gradient may be corrupted by noise. This gives the flexibility to cover, for example, (sub)problems arising in image reconstruction or constrained optimization algorithms. We provide computational examples that illustrate the findings and point to possible strategies for practical determination of the stabilization parameter that balances the size of the critical region with a relaxation of the acceptance criterion (or descent property) of the algorithm. Keywords Nonsmooth optimization ·Global convergence ·Noisy optimization BChristoph Hansknecht [email protected] Christian Kirches [email protected] Paul Manns [email protected] 1Institute for Mathematical Optimization, TU Braunschweig, Braunschweig, Germany 2Chair of Numerical Analysis and Optimization, TU Dortmund University, Dortmund, Germany 123
568 C. Hansknecht et al. 1 Introduction Handling non-smoothness is an ubiquitous research question in nonlinear optimization because it arises naturally in different areas, for example, penalty functions for constrained optimization [4], statistical data analysis and signal processing [5,6], and neural network architectures [7]. In this work we study the convergence properties of successive linear programming algorithms to solve the optimization problem min x∈Rnφ(x):=ω(F(x)), (P) where ω:Rp→Ris convex and Lipschitz continuous with polyhedral epigraph, and F:Rn→Rpis twice continuously differentiable. Moreover, we assume that Fand its Jacobian can only be accessed inexactly so that their evaluations are corrupted by noise. This and similar problems have been studied in the literature, see, for example, [8–12] and the references therein. Many optimization problems can be formulated in terms of problem (P) such as the Lagrangian form min x∈Rny−Ax2 2+βx1 of the famous LASSO problem [5,13] with A∈Rm×n,β>0, and y∈Rmthat is particularly popular among data scientists for sparse parameter identification in overparameterized models. More broadly speaking, a wide class of nonlinear optimization problems fall under (P) as well. As an example, the unconstrained minimization of a smooth objective f:Rn→Rcan be modeled by setting p=1 and ω(x)=xfor x∈R. We can also examine nonlinear programs constrained by smooth functions g:Rn→Rmand h:Rn→Rkyielding problems of the form min x∈Rnf(x)s.t. g(x)≤0,h(x)=0.(NLP) These problems may be solved by minimizing a non-smooth exact penalty function of the form φ(x,ν):= f(x)+ν g(x)T +,h(x)TT 1 ,(1) where y+:= max(y,0)and ν>0. In fact, strict local solutions of (NLP) are local minimizers of φ(x,ν) for a sufficiently large value of νif gand hare smooth and satisfy the Mangasarian–Fromovitz constraint qualification [14, Theorem 4.4], [4, Theorem 17.3] at the respective points. The penalty function φcan be expressed as ω(F(x)), where F(x):=(f(x), g(x)T,h(x)T)Tis smooth and ω(x,y,z):=x+ ν(y+T,zT)T1is convex and polyhedral. Besides problems of type (NLP), a variety of other problems, such as linear or nonlinear fitting problems can be formulated in terms of (P) as well. Noisy functions The combination of unconstrained optimization with noisy observations has recently been examined [1,15–21]. The authors consider the minimization of a smooth function φ:Rn→Rwhile only having access to 123
Convergence of successive linear programming… 569 f(x)=φ(x)+ε(x)and g(x)=∇φ(x)+e(x), (2) where the only assumptions are that both |ε|and eare uniformly bounded. Consequently, it is not generally possible to generate a sequence {xk}of iterates converging to a local optimum or stationary point of φ. Intuitively, while the gradient noise eis small compared to ∇φ, the direction gis a suitable search direction with respect to φ. This allows for the use of an Armijo-like globalization strategy [16,20] or, in case of [1,17], a trust-region method, where the noise is handled by stabilizing the reduction ratio, which is of course closely related to the Armijo condition. As soon as a region is reached where the noise produced by εand ebecomes too large relative to φand ∇φ respectively, no further progress can be expected and the algorithm may stall. However, this critical region is visited infinitely often and once reaching it, the algorithm does not produce objective values much larger than the objective values attained in the critical region. The authors also study the problem of adapting quasi-Newton methods to the noisy setting. Contribution We build on the ideas in [1,15,21] and consider the non-smooth problem (P)ina setting, where function and derivative evaluations are only available as noisy observations. As the authors in [1,15,21], we assume the following noise model: Rather than being able to evaluate Fand its derivative Fdirectly, we only have access to ˜ F(x):=F(x)+δF(x)and G(x):= F(x)+δF(x). These proxies consist of the original functions Fand Fas well as error functions δF: Rn→Rpand δF:Rn→Rp×n. In terms of the problem (NLP), this is tantamount to noise in the objective f, the constraints g,h, and their respective derivatives. Contrary to this, we assume that the function ωdoes not suffer from any noise. What is more, we presume that the structure of ωis well understood in the sense that, for example, its Lipschitz constant is known, which is certainly the case for the penalty function in (1). In order to solve optimization problems of the form (P), we propose a trust-region algorithm leaning on the successive linear programming template proposed in [2] and a convergence analysis that builds on the ideas in [3,15,21]. Specifically, we use a stabilization of the iterate acceptance test in order to assert that a neighborhood of a stationary point is visited infinitely often by the iterates produced by the algorithm. The polyhedral structure of ωis handled by first solving a linear program in order to determine a direction for a subsequent Cauchy point determination. This can also be interpreted as an active set determination for the corresponding kinks of the polyhedral epigraph of ω. We also provide computational examples that illustrate the theoretical results and the practical behavior of the algorithm. Moreover, the results point to open questions and possible approaches regarding the choice of the correct stabilization parameter in the acceptance test. Structure of the Remainder We introduce the successive linear programming algorithm and the modified acceptance test in Sect.2. The asymptotics of the algorithm are analyzed in Sect. 3.We 123
570 C. Hansknecht et al. provide computational examples and the corresponding results in Sect.4.Wedrawa conclusion in Sect.5. 2 A noise-tolerant successive linear programming algorithm In the noisy setting, we cannot expect to find the true optimum or stationary points of φ, since we do not have access to Fand F. Specifically, in a small region around the true optimum x∗,˜ Fand Gmay oscillate, thereby making their evaluations unreliable. This impairs globalization strategies in nonlinear programming because their acceptance tests require reliable evaluations of Fand a model function involving F. In the non-noisy regime, a trust-region method produces a sequence {xk}of iterates by assembling and subsequently optimizing model functions qk:Rn→R, yielding astepdk. The quality of dkis determined according to the reduction ratio ρk:= φ(xk)−φ(xk+dk) φ(xk)−qk(dk), which is used to determine whether or not the step will be accepted. However, in the noisy setting, we only have access to ˜ Fleading to a noisy composite function ˜ φ(x):=ω( ˜ F(x)). While we can build a model ˜qk:Rn→Rpwhich coincides with ˜ Fat xk, we cannot control the numerator ˜ φ(xk)−˜ φ(xk+dk). Indeed, if we reduce the trust region, sending dkto zero, the denominator of ρkwill tend to zero while the numerator will oscillate, making the ratio unreliable. To alleviate this problem, we turn towards a recent adaptation [1] of trust-region methods in order to solve the noisy counterpart of (P). The authors of [1] add a correction term, that is a positive constant ϑ>0, to both the numerator and denominator of the reduction ratio ρkto mitigate the effect of noisy evaluations, yielding a modified ratio ˜ρk:= ˜ φ(xk)−˜ φ(xk+dk)+ϑ ˜ φ(xk)−˜qk(dk)+ϑ. The parameter ϑcan then be chosen appropriately in order to stabilize the ratio. As we will see, this means that for ϑlarge enough, the iterates of the successive linear programming algorithm converge to a critical region around a stationary point. The downside is that this region grows with ϑand the algorithm also accepts steps that do not improve the objective. Apart from this adjustment, we follow the algorithmic approach in [3]. Specifically, we use the following partially linearized and quadratic models at xk∈Rn ˜ (x;d):=ω˜ F(x)+G(x)d,(3) ˜ k(d):= ˜ (xk;d), and (4) ˜qk(d):= ˜ k(d)+1 2d,Bkd,(5) 123
Convergence of successive linear programming… 571 where the Bk∈Rn×nare symmetric (not necessarily positive definite) approximations of the curvature of ω◦F. Here and throughout we let ·,· be the standard scalar product in Rn, inducing the 2-norm, which we simply denote by ·. For matrices, we let ·be the corresponding operator norm, given by its largest singular value. Algorithm Based on the models above, our noise-tolerant approach to solving (P) is laid out in Algorithm 1. In each iteration, an initial step dLP is computed in Line 3by solving the problem min dLP≤LP k˜ k(d), where ·LP is a norm on Rndefining the trust region associated with ˜ k. While we make no further assumptions regarding this norm, it is in practice advantageous to cast this subproblem as a linear program to be solved using state-of-the-art LP solvers [22, 23]. To this end, two conditions should be met. First, the epigraph of ωshould be polyhedral. Second, the feasible region should be polyhedral as well, or, equivalently, ·LP should be a polyhedral norm, such as ·1or ·∞. In any case, due to the equivalence of norms in Rnthere exists a constant γ>0 such that for each d∈Rn it holds that d≤γdLP.(6) The algorithm proceeds to compute a Cauchy step dC kin Lines 4–7. To this end, it employs a line search initialized with a step size sufficiently small to ensure that the Cauchy step falls into the trust region bounded by LP k. During the line search the step size is shortened by a factor of 0 <τ<1 until the quadratic reduction achieved by the Cauchy point is within a factor of 0 <η<1 of its linear reduction. The actual step dkto be taken in Line 8can be different from the Cauchy step dC k, provided that it improves upon the quadratic reduction of dC k. This gives some algorithmic flexibility, allowing for the computation of Newton-type steps in order to achieve local quadratic convergence. Based on the stabilized reduction ratio ˜ρk computed in Line 9, the step is either accepted (Lines 11–12) or rejected (Lines 14– 15) according to an acceptance threshold of ρu>0. Additionally, the trust-region radii k+1and LP k+1are adjusted based on ˜ρk: 1. The value of k+1is increased or decreased based on whether ˜ρkachieves a value of at least ρs. The decrease is such that the new trust-region radius is at most κu<1 times as large as the previous one, thereby ensuring a true reduction, while being at least κldk(with 0 <κ l≤κu) in order to prevent an immediate collapse of the trust region. 2. If ˜ρkachieves at least ρu(note that ρu≤ρs), the LP trust-region radius LP k+1 is increased beyond dC kLP, as long as it does not exceed the upper bound of LP max ≥1. The new LP trust-region radius is also only increased beyond LP k+1 if the full LP step dLP was accepted (i.e., αk=1), indicating that the partially linearized model ˜ kis a good approximation of ˜ φkacross the entire LP trust region. If ˜ρkfalls short of ρu,LP k+1is decreased while being kept within a factor of θ>0 of dkLP. 123
572 C. Hansknecht et al. Remark 2.1 When applied to problem (NLP), Algorithm 1uses the strategies introduced in [2], which form the basis of the active set method in the highly successful Knitro code [24], which combines sequential linear programming with equality constrained quadratic programming approaches in order to achieve robust performance over a range of large-scale nonlinear programming problems. Algorithm 1: A noise-tolerant algorithm to minimize φ(x)=ω(F(x)) Input : Functions ω,˜ F,G, Initial point x0∈Rn, Initial trust region radii 0 < LP 0≤LP max,0< 0 Parameters: Acceptance thresholds 0 <ρ u≤ρs<1, Step adjustments 0 <κ l≤κu<1, θ>0, Cauchy line search parameters 0 <η<1, 0 <τ <1, Ratio stabilizer ϑ>0, Maximum LP trust region radius LP max ≥1 Output : Primal point x∗∈Rn 1k←0 2until Some termination criterion is satisfied 3Compute LP step dLP k←arg mindLP≤LP k˜ k(d) 4αk←min 1, k/dLP k 5while ˜ φ(xk)−˜qk(αkdLP k)<η˜ φ(xk)−˜ k(αkdLP k)do 6αk←ταk 7dC k←αkdLP k 8dk←Step dsuch that d≤kand ˜qk(d)≤˜qk(dC k) 9Compute stabilized reduction ratio ˜ρk←˜ φ(xk)−˜ φ(xk+dk)+ϑ ˜ φ(xk)−˜qk(dk)+ϑ 10 if ˜ρk≥ρuthen Accept step 11 Set xk+1←xk+dk 12 Pick LP k+1∈dC kLP, LP maxsuch that LP k+1≤LP kif αk<1 13 else Reject step 14 Set xk+1←xk 15 Pick LP k+1∈min(θdkLP, LP k), LP k 16 if ˜ρk≥ρsthen 17 Set k+1≥k 18 else 19 Choose k+1∈[κldk,κ uk] 20 k←k+1 21 return x∗=xk 123
Convergence of successive linear programming… 573 3 Convergence analysis of Algorithm 1 We begin our convergence analysis with the introduction of the standing assumptions and a recap of the relevant stationarity concept for (P) in Sect.3.1. We analyze the criticality measure for this notion of stationarity in the noisy setting in Sect.3.2.Weuse these results to prove lower bounds on the trust-region radii that occur in Algorithm 1 in Sect.3.3, which are then used to obtain sufficient decrease and, as a consequence, convergence of the produced iterates to critical regions in Sect.3.4. 3.1 Standing assumptions and stationarity In order to study the convergence properties of Algorithm 1, we make several assumptions regarding the amount of noise, the functions ω,F, and the matrices Bkused in the quadratic models ˜qk. Assumption 1 We assume that the noise is uniformly bounded via δF(x)≤εFand δF(x)≤εFfor all x∈Rn. We refer to εFand εFas the noise levels of the functions ˜ Fand Grespectively. Assumption 2 ωis Lipschitz-continuous with constant Lω, i.e., it holds for all x,y∈ Rnthat |ω(x)−ω(y)|≤Lωx−y. Assumption 3 Fand Fare Lipschitz-continuous with constants LFand LF, i.e., it holds for all x,y∈Rnthat F(x)−F(y)≤LFx−yand F(x)−F(y)≤LFx−y. Assumption 4 The Hessian approximations Bkare bounded, i.e., there exists β>0 such that for all d∈Rn,k>0 it holds that |d,Bkd| ≤ βd2. Our aim in the following is to find a local optimum of φ. A first-order necessary condition (see [25, p. 184]) of optimality for (P) states that x∗∈Rncan only be a local optimum if max λ∈∂ω(F(x∗))λ, F(x∗)d≥0 for all d∈Rn,(7) where ∂ω(z)denotes the subdifferential of ωat z∈Rp. To measure how close a point x∈Rnis to satisfying these conditions, we use the criticality measure introduced in [26]. Specifically, we let (x;) :=φ(x)−min dLP≤(x;d) 123
574 C. Hansknecht et al. and set k() :=(xk;) as before. Clearly, since d=0 is a feasible solution of the inner optimization problem, the value (x;) is always non-negative. On the other hand, the following result establishes that a vanishing reduction over a nontrivial trust region is tantamount to reaching a point satisfying first-order conditions, which we call a critical point: Lemma 3.1 ([26], Lemma 2.1) A point x∗∈Rnsatisfies conditions (7)iff there exists >0such that (x∗;) =0. As a consequence of this result, an algorithm solving (P) should aim at generating a sequence of iterates such that lim infk→∞ k() =0 for some fixed >0 (assumed to be 1 in the following). This ensures the existence of an accumulation point x∗of the iterates satisfying first-order conditions. 3.2 Analysis of model function and criticality measure in the presence of noise Since we do not have access to the values of Fand Frequired to compute k,we define a noisy measure of criticality via ˜ (x;) := ˜ φ(x)−min dLP≤˜ (x;d) and set ˜ k() := ˜ (xk;) as before. This function is also non-negative if the same realization of the function noise δF(xk)is used when computing ˜ φ(xk)and constructing the linear approximation ˜ k(d). We analyze its properties and relationship to kbelow. Since δFcannot be assumed to be continuous, neither can ˜ φ. This differs from the analysis in [3], where the Lipschitz-continuity of φis used to argue that the reduction ratio approaches one if the trust-region radius is driven to zero. We can, however, state that the criticality measures kand ˜ kare related by the following approximation result: when considering a fixed xk, we claim that ˜ k(1)→k(1)for εF→0 and εF→0 and that we also have convergence of the minimizers of the convex programs in the definitions of ˜ k(1)and k(1). This follows from the epi-convergence of the functionals LεF,εF(d):= ˜ k(d)+idLP≤(d)and L0,0(d):=k(d)+idLP≤(d), where iA:Rd→{0,∞} is the indicator function of A⊂Rd, that is iA(x)=∞if x/∈Aand iA(x)=0 else. We recall that the functionals LεF,εFepi-converge to L0,0 if and only if for all d∈Rnthe inequalities L0,0(d)≤lim inf εF,εF→0 LεF,εF(dε)for all sequences dε→dand (8) L0,0(d)≥lim sup εF,εF→0 LεF,εF(dε)for some sequence dε→d(9) hold, see, for example, [27, § 7], which is shown below. 123
Convergence of successive linear programming… 581 Proof The proof is in Appendix A. The following convergence theorem states that when the objective of (P) is bounded below, an application of Algorithm 1will produce one of two possible mutually exclusive outcomes: the algorithm may stop at a critical point after a finite number of iterations as described in Lemma 3.13 or, alternatively, Algorithm 1visits a critical region infinitely often. In terms of the functions ω,˜ F, and G,thecritical region is defined as C(δ) := x˜ (x;1)≤δ. By definition, an iterate xkproduced during the execution of Algorithm 1is contained in C(δ) if ˜ k(1)≤δ. What is more, Proposition 3.2 establishes that ˜ k(1)tends to k(1)as the errors εFand εFapproach zero. These results therefore suggest that the iterate is close to being optimal in the sense of Lemma 3.1. Theorem 3.14 Consider an application of Algorithm 1to the noisy variant of problem (P). Suppose that Assumptions 1to 4and (11)hold. Then either ˜ k(1)=0for some k ≥0, or lim k→∞ ˜ φ(xk)=−∞, or there are infinitely many k ∈Nsuch that xk∈C(δmax), where δmax := max ⎛ ⎝ϑ(1−ρu)γ LP max ρuηB,ϑ(1−ρu)γ LP max ρuηA⎞ ⎠ is given in terms of the constants A,B from Lemma 3.12. Proof If there are only finitely many accepted steps, the result follows from Lemma 3.13, yielding the first possibility. Otherwise, we can assume that during the algorithm, an infinite number of accepted steps occurs. If ˜ φ(xk)tends to −∞, the second possibility occurs, so we can assume in the following that ˜ φ(xk)(and hence φ(xk))is bounded below. Let Kbe the sequence of accepted steps, i.e., consisting of those kwhere xk+1= xk. Clearly, if lim infk→∞ ˜ k(1)=0, then the result follows. So we can assume that there exists a δ>0 such that ˜ k(1)≥δfor all k≥k0. The claim stating that the region C(δmax)is visited infinitely often is tantamount to ensuring that δ≤δmax, which will be the aim of the remainder of this proof. For each k∈K,k≥k0we have that ˜ φ(xk)−˜ φ(xk+1)+ϑ ˜ φ(xk)−˜qk(dk)+ϑ≥ρu>0. 123
582 C. Hansknecht et al. We deduce using Lemma 3.10 that ˜ φ(xk)−˜ φ(xk+1)≥ρu˜ φ(xk)−˜qk(dk)+(ρu−1)ϑ ≥ρuηαkmin(LP k,1)˜ k(1)+(ρu−1)ϑ. It follows that ˜ φ(xk)−˜ φ(xk+1)≥ρuηαkLP kmin 1 LP k ,1˜ k(1)+(ρu−1)ϑ ≥ρuηαkLP k 1 LP max ˜ k(1)+(ρu−1)ϑ. We can now apply Lemma 3.12 to bound αkLP kbelow based on min and the constants Aand B: ˜ φ(xk)−˜ φ(xk+1)≥ρuηmin γ LP max ˜ k(1)+(ρu−1)ϑ =ρuηmin(A,Bδ) γ LP max ˜ k(1)+(ρu−1)ϑ ≥ρuηmin(A,Bδ) γ LP max δ+(ρu−1)ϑ. Let us assume towards a contradiction that δ>δ max. We distinguish two cases with respect to the minimum min(A,Bδ): 1. The minimum is attained at A, implying that ˜ φ(xk)−˜ φ(xk+1)≥ρuηA γ LP max δ+(ρu−1)ϑ ≥ρuηA γ LP max δmax+(ρu−1)ϑ +C1. for a constant C1>0. Using the fact that δmax ≥ϑ(1−ρu)γ LP max ρuηA by definition of δmax, this implies that ˜ φ(xk)−˜ φ(xk+1)≥C1>0. 2. The minimum is attained at Bδ, implying that ˜ φ(xk)−˜ φ(xk+1)≥ρuηBδ2 γ LP max +(ρu−1)ϑ ≥ρuηBδ2 max γ LP max +(ρu−1)ϑ +C2 123
Convergence of successive linear programming… 583 for a constant C2>0. We now use the fact that δmax ≥ϑ(1−ρu)γ LP max ρuηB to deduce that ˜ φ(xk)−˜ φ(xk+1)≥C2>0. In either case ˜ φdecreases by min(C1,C2)>0 from xkto xk+1. Since this decrease is strictly positive and there are infinitely many accepted steps in the sequence K,it follows that ˜ φ(xk)tends to −∞, which is a contradiction. It must therefore hold that δ≤δmax as desired. Interpretation of Theorem 3.14 In theoretical terms, the result in Theorem 3.14 is as expected: the size of the critical region Cdepends on the stabilization parameter ϑ, If we increase ϑ, Algorithm 1can tolerate a larger amount of noise at the cost of a decreased accuracy with respect to the criticality measure ˜ . Of course, problem (P) can generally also be unbounded. The remaining case, where ˜ k(1)=0forsomek∈N, can involve different scenarios. If k(1)=0 holds as well, the iterate xkis a critical point of (P), which is the ideal situation. Otherwise, the noises δF(xk)and δF(xk)attain values such that xkappears to be critical in the noisy model. In case of an unconstrained version of (NLP), this is tantamount to a non-zero gradient that is canceled out by noise. For a function Fthat is not afflicted by noise, i.e., satisfying εF=εF=0, it holds that Mε 0=Mε 1=0, where Mε 0and Mε 1are the constants from Lemma 3.4. This allows us to set ϑ=0, whereby we recover the original algorithm discussed in [3]. If we apply Theorem 3.14 in this situation, it follows from ϑ=0 that δmax =0, implying that the region C(δmax)contains precisely the points critical with respect to ˜ φ, which itself coincides with φin this particular case. Thus, the original convergence result [3, Theorem 3.8] follows for Algorithm 1in the noiseless case. As mentioned in the introduction, we can also apply Algorithm 1to solve smooth unconstrained nonlinear problems affected by noise. To this end, we can set ω(x)=x, achieving a Lipschitz constant of Lω=1. Lemma 3.12, specifically (11), then suggests a stabilization of ϑ≥2εF+εF 1−ρs >2εF+εF,(22) which may be weaker than the stabilization of rεFanalyzed in [1]ifεFbecomes sufficiently large. This is particularly true for all εF>0 for the choice r=2/(1−ρs), which corresponds to the choice r=2/(1−c2)in (8) in [1]. ThereasonistheestimateinLemma3.7 based on the convexity of ω. Conversely, the authors of [1] have a Lipschitz continuous derivative of the objective at hand. In that case, the criticality measure satisfies x∈C(δ) ⇐⇒ −min dLP≤1G(x)d≤δ, 123
584 C. Hansknecht et al. Table 1 Parameters used for numerical experiments Symbol Explanation Value LP 0Initial LP trust-region radius 1 LP max Maximum LP trust-region radius 10 0Initial trust-region radius 1 ρuStep acceptance threshold 0.1 ρsThreshold for increase of 0.5 κlLower bound for adjustment of after failed step 0.1 κuUpper bound for adjustment of after failed step 0.8 θLower bound for adjustment of LP after failed step 0.5 ηFactor of relative decrease for Cauchy step 0.1 τShortening factor for Cauchy line search 0.5 which in turn is equivalent to G(x)being bounded above by a constant. 4 Numerical experiments In order to illustrate the performance and examine the behavior of Algorithm 1, we implement the algorithm1in Python (3.11.5), using numpy [28] (1.26.0), scipy [29] (1.11.3) (including HiGHS [23] (1.5.0) as LP solver), and Ipopt [30] (3.14.13) to solve the respective subproblems. We generally compare the performance of the classical algorithm (i.e., Algorithm 1with a stabilization of ϑ=0), with its stabilized counterpart (where ϑ>0). In terms of termination criteria, we first of all impose an iteration limit, after which the algorithm terminates. Secondly, we monitor the LP trust-region radius, LP.Ifthe radius contracts to a value close to zero (1 ×10−10), we see this as a failure of the algorithm and let it terminate. Lastly, when the noise criticality ˜ k(1)falls below the threshold of 1 ×10−6, we terminate the algorithm, knowing that the current iterate is very close to being optimal. We then examine the final iterate xf, i.e., the iterate xk in Algorithm 1of the iteration at which the termination criterion becomes satisfied. If not indicated otherwise, we choose the parameters of Algorithm 1according to the values in Table 1and the stabilization parameter ϑ∗ ε. In order to obtain inexact evaluations of a given function F:Rn→Rpand its derivative, we inject noise by setting δF(x)=XF∈Rp,XF∼Bn(εF), and δF(x)=XF∈Rp×n,XF∼Bnp(εF), (23) where Bn(s)denotes the uniform distribution on the n-dimensional Euclidean ball centered at the origin with radius s≥0. We can compute a sample from this distribution 1Available at https://github.com/chrhansk/noisy-nonlinear. 123
Convergence of successive linear programming… 585 by first sampling ntimes from a standard normal distribution and scaling the resulting vector [31] by its Euclidean norm to obtain a sample from the uniform distribution on Sn(1). We can then use inverse transform sampling to obtain a suitable radius in [0,s] and rescale the vector accordingly. The injected noise is bounded according to the noise levels εFand εFwhile being sufficiently unpredictable to significantly affect the solution process. In Sect.4.1, we augment an example of a quadratic test problem from [1] with a non-smooth term. We obtain qualitatively similar results in this case. Then we compare and visualize the different behaviors of the unstabilized algorithm and the stabilized algorithm for the Rosenbrock test function in Sect.4.2. In Sect.4.3, we apply the algorithm to an image reconstruction problem with total variation regularization and assess the impact of different choices of the stabilization parameter. Finally, in Sect. 4.4, we apply the algorithm in a penalty method for a small constrained optimization problem from CUTest [32] as motivated in the introduction, which points to future research directions. Sections4.1,4.2 and 4.4 use essentially the same type of a nonsmooth objective function that includes an 1-penalty term and we provide the required estimates of its Lipschitz constant in Appendix B. 4.1 Failure of the classical algorithm To illustrate the difference in performance between the classical algorithm and Algorithm 1, we consider the case of 1-penalized optimization problems of the form x→ f(x)+λx1with a smooth function f:Rn→R. It is clear that these problems are non-smooth due to the presence of the ·1term, while being expressible as problems of type (P) based on suitable choices of ωand F. This problem class also enables us to minimize ˜ kand ˜qkover the trust regions defined in terms of LP and by solving linear or quadratic programs respectively. What is more, the only curvature information in this problem class is due to f, enabling us to either use the Hessian of f(or any quasi-Newton approximation) to obtain the matrices Bk. We specifically examine the case where fis a quadratic of the form f(x)=1 2x,Dx, where Dis the matrix in Rn×nfor n=8 given as D=diag(10−5,10−4.75,10−4.5,...,10−3.25), taken from [1], where this optimization problem has been studied without an 1penalty. It is apparent that the optimal solution of this instance of (P)isx∗=0. We set the parameter λto 1 ×10−2while injecting noise according to (23) with noise levels εF=1×10−1and εF=1×10−5, and initialize Algorithm 1with the initial point x0=(1000,0,...,0)T, limiting the number of iterations to 50, and performing quadratic steps based on the true Hessian D. We show an example of the difference in performance in Fig. 1, where the values of the reduction ratio are clipped to ±5 in order to properly display the results. We 123
586 C. Hansknecht et al. see that the classical algorithm performs dramatically worse than its stabilized counterpart. Indeed, the classical algorithm stalls almost immediately, due to the reduction ratio ρkbecoming unreliable. Consequently, the LP trust region collapses, and the classical algorithm makes no progress towards optimality. Conversely, the addition of a stabilization yields an algorithm rapidly approaching the optimum, both in terms of primal distance and objective value while maintaining a reasonably large LP trustregion radius. Similarly, noisy and noiseless criticality decrease rapidly throughout the iterations of the stabilized algorithm. Unfortunately, the criticality bound established in Theorem 3.14 attains a value of δmax ≈54 000, limiting its use in terms of the criticality actually achieved throughout the iterations. We would like to point out that the failure of the classical algorithm is not guaranteed in this scenario: By running the experiment with 100 different random seeds we found that the classical algorithm stalls in about half (57) of the cases, while performing well in the other half. Conversely, the stabilized algorithm consistently performs well in all cases. Its characteristics are qualitatively similar to the case that is depicted in Fig.1. A key problem of the classical algorithm is therefore its unreliability when applied to noisy functions. 4.2 A variant of the Rosenbrock problem Following the previous experiments conducted based on the quadratic function, we go on to examine the performance on a variant of the famous Rosenbrock function, given by R(x,y):=(a−x)2+b(y−x2)2 with parameters of a=1, b=100. The Rosenbrock function has a unique optimum at (x∗,y∗)=(a,a2), i.e., at (1,1)for our choice of parameters. We modify the problem by adding a penalty of λ(x,y)−(x∗,y∗)1with a value of λ=1×10−1, yielding a problem of type (P) having the same global optimum as R. We show an example of the difference in performance between the classical and stabilized algorithms in Fig.2. The figure shows the trajectories generated by Algorithm 1with and without stabilization starting at (x0,y0)=(−1.5,0), injecting noise according to (23)for different values of εFand a fixed value of εF=1×10−5, performing quadratic steps according to the true Hessian of Rwith an iteration limit of 50. Examining the trajectories of the classical algorithm, shown in Fig. 2a, we find that for different values of εF, the trajectories are initially almost identical, until the algorithm stalls at points with a distances to the optimum increasing with εF. Conversely, the trajectories of the stabilized algorithm, shown in Fig. 2b, vary significantly for different noise levels. However, the stabilization yields trajectories leading significantly closer to the optimum than those of the classical algorithm even for larger noise levels. This is confirmed by the statistics shown in Fig. 3, displaying the distribution of the distance to the optimum for various noise levels for 100 different random seeds, demonstrating that the stabilized algorithm consistently outperforms the classical one, in particular for larger noise levels. 123
Convergence of successive linear programming… 587 Fig. 1 Performance on an 1-penalized quadratic problem over 50 iterations of Algorithm 1with noise levels of εF=1×10−1and εF=1×10−5.Thex-axes always show the iteration count of Algorithm 1, the y-axes show the quantities indicated in the captions below the respective subplots 4.3 Image reconstruction Although this is not the focus of this article, we also provide a computational example that has a meaningful problem size. Specifically, we consider an artificial task of reconstructing an image under noisy observations. That is we seek to recover a matrix Y∈RM×Nwith values normalized to be in [0,1]. In our setting, Yis only available in the form of noisy observations. Specifically, for an input X∈RM×N, the fidelity 123
588 C. Hansknecht et al. Fig. 2 Trajectories for the modified Rosenbrock problem plotted over the shifted criticality 1 +φ(x)− mindLP≤1ω(F(x)+F(x)d), which allows to show it on a logarithmic scale. The markers show the position of the final iterate of X, given by the squared Frobenius norm of X−Y,1 2X−Y2 F,aswellasits derivative with respect to Xcannot be evaluated. Instead, we have access to the map X→ 1 2X−˜ Y2 F, where ˜ Yis a noisy version of Y, redrawn for each guess X.We obtain the term ˜ Yby sampling from a componentwise uniform distribution δY(X)=YF∈RM×N,YF∼UM×N(−εimg,ε img), and setting ˜ Yto Y+δY(X)clipped back to have coefficients in [0,1]. The amount of noise injected to the image is in turn governed by the parameter εimg ≥0. This noise model translates into noise injected into the evaluations of Fand F, which can be 123
Convergence of successive linear programming… 589 Fig. 3 Distribution of the distance between the final iterate xfand noiseless optimum x∗for different values of the noise level εF estimated in terms of εimg,M, and N(see Appendix B) while not conforming to the noise model (23). We also impose an anisotropic total variation (TV) regularization penalty, defined as TV(X):= M−1 i=1 N j=1|Xi+1,j−Xi,j|+ M i=1 N−1 j=1|Xi,j+1−Xi,j|, to our objective, turning the problem non-smooth, balancing off fidelity and regularity by a parameter λ>0. The regularization term can be expressed as AM,NX1with a suitable matrix AM,N. Consequently, we can formulate the reconstruction problem as problem of type (P), consisting of a smooth term (the fidelity), and an 1-penalized linear function. Naturally, the regularization does not suffer from any noise. Based on a regularization parameter of λ=5×10−3we reconstruct the image showninFig.4a. To avoid having to solve large quadratic problems, we do not compute quadratic steps and opt to instead increase the number of iterations to 100 starting from X0=0. As a baseline, Fig. 4b shows the image when we apply Algorithm 1to the original image (i.e., setting εimg =0). The restored image closely resembles the original one. We proceed to study the effect of the value of ϑon the quality of the reconstructed image. In principle, it must hold that ϑ≥ϑ∗ εin order for the criticality to provably converge. It is however unclear whether setting ϑto ϑ∗ εyields the best results in practice. The large value of δmax seen in Sect. 4.1 seems to suggest that (11) is rather pessimistic. We therefore examine the performance of Algorithm 1for values of ϑnot necessarily satisfying the inequality. 123
590 C. Hansknecht et al. Fig. 4 Sample image for the image reconstruction To gauge performance, we record both the original and noisy objective after the iterations. The results, shown in Fig.5, demonstrate the effect of ϑ: For small stabilization values, Algorithm 1stalls early on, as was the case in our previous experiments. As we increase ϑ, there appears to be an optimal choice or small region, where both the noisy and the noiseless evaluation of the final objective are minimized. This effect is more pronounced for higher values of εimg, where the noiseless objective for ϑ=0 is about 4 times as large as that of the optimal choice of ϑ. It is interesting to see that this sweet spot also shows in the noisy objective, suggesting that noisy observations may be sufficient to find it. Lastly, as we increase ϑbeyond the sweet spot, the final objective increases sharply. This is likely due to the case that the algorithm simply accepts too many steps, even when they are in fact disadvantageous in terms of progressing towards an optimum. Ultimately, for a sufficiently large value of ϑ,all steps are accepted, which, as the final objective suggests, leads to poor solutions. The values of ϑ∗ εare given by 4 ×104,2×105, and 4 ×105for the respective noise levels, significantly exceeding the optimal values and beyond the point, where all steps are accepted. We also find that the objectives are consistent with the visual appearance of the reconstructed images, shown in Fig. 6: While setting ϑto zero yields satisfactory results, even though a grainy appearance remains for larger noise levels, a disproportionately large value of ϑproduces a distorted result with visible artifacts. For our best guess of ϑ, the restored images do not suffer from artifacts and closely resemble the original one even for larger noise levels. 4.4 Constrained optimization As a final example and in order to demonstrate the possible use of Algorithm 1as a subproblem solver in constrained optimization algorithms, we study a constrained optimization problem of the type (NLP). Specifically, we examine the behavior of Algorithm 1when applied to the HS71 benchmark problem of the CUTest [32] 123
Convergence of successive linear programming… 597 Since the only difference between the linearized and quadratic model is the quadratic term, we have that 1 2(αk/τ)2dLP k,BkdLP k≥(1−η) ˜ φ(xk)−˜ k(αdLP k/τ). The left hand side can be bounded above by using Assumption 4and relation (6) to yield 1 2(αk/τ)2dLP k,BkdLP k≤1 2(αk/τ)2βγ2dLP k2 LP ≤1 2(αk/τ)2βγ2dLP kLPLP k. Similarly, for the right hand side we can use Lemmas 3.7 and 3.8 to obtain ˜ φ(xk)−˜ k(α/τdLP k)≥αk/τ ˜ φ(xk)−˜ k(dLP k) ≥αk/τ min(1, LP k)˜ k(1), Putting these inequalities together yields the bound dC kLP =αkdLP kLP ≥2(1−η)τ βγ2min 1,1 LP k˜ k(1) required to complete the proof. Proof of Lemma 3.13 Let k0be the index of the last accepted step. Then, xk+1= xk=:x∗for all k>k0. Consequently, after finishing the k0-the iteration, ˜ k(1)stays at a constant value of δ≥0. What is more, due to the rejection of the steps following iteration k0it holds for all k>k0that k+1≤κuk< k(since κu<1). Therefore, ktends to zero. Recall from Lemma 3.12 that if δ>0, then kis bounded away from zero. Therefore, since ktends to zero, it must hold that δ=0. B Estimations Lipschitz constant of the 1-penalty function In the following, we give an estimation for the Lipschitz constant Lωof the penalty function ω:R×Rm→R,ω(x,y)=x+νy1, based on the constant ν>0 and the dimension m∈N. Since we use this function in all of the examples in Sect.4, and since the value of ϑ∗ εdepends on the value of Lω, we make its derivation explicit. To obtain an optimal value of Lω, we solve the optimization problem max x,y,x,y|ω(x,y)−ω(x,y)| s.t. (x−x)2+y−y2≤1, 123
598 C. Hansknecht et al. i.e., we maximize the difference in values of ωwhile controlling the distance between the points (x,y)and (x,y). Observe that |ω(x,y)−ω(x,y)|=|(x−x)+νy−y1|, from which it follows that both the objective and the constraint value only depend on x−xand y−y. We can therefore simplify the problem by setting y=0 and y=0: max x,y|x+νy1| s.t. x2+y2≤1. We can simplify the problem further by realizing that we can assume both xand y to be non-negative, eliminating the absolute value in the objective. The largest ratio of νy1over y2is achieved by setting all entries of yto the same value y0∈R, yielding the problem max x,y0 x+νmy0 s.t. x2+my2 0≤1 x,y0≥0. By setting z:=√my0, we obtain the problem max x,zx+ν√mz s.t. x2+z2≤1 x,z≥0. The optimal solution of this problem is attained at x∗ z∗=1 √1+ν2m1 ν√m, yielding the objective √1+ν2m=:Lω. Image Reconstruction In the following, we provide estimations regarding the noise levels associated with the image fidelity map introduced in Sect.4.3. Recall that the squared Frobenius norm of amatrixA∈RM×Nis given by A2 F:= M i=1N j=1a2 ij. Thus, if |aij|≤εfor a given ε>0, it follows that A2 F≤ M i=1 N j=1 ε2=ε2MN, 123
Convergence of successive linear programming… 599 and therefore, that AF≤ε√MN. The noisy fidelity function ˜ F(X)satisfies the identities ˜ F(X)=1 2X−˜ Y2 F=1 2X−(Y+δY(X))2 F =1 2X−Y2 F+X−Y,δ Y(X)F+1 2δY(X)2 F, where ·,·Fdenotes the inner product that induces the Frobenius norm. Since the entries of Xand Yare in [0,1], and therefore have absolute values bounded by 1, it follows that X−YF≤√MN. The choice of distribution implies that the values in δY(X)are bounded by ±εimg, and therefore that δY(X)F≤εimg√MN,from which it follows that |˜ F(X)−F(X)|≤X−YFδY(X)F+1 2δY(X)2 F≤εimg +1 2ε2 imgMN, by means of the Cauchy–Schwarz inequality for ·,·F. Furthermore, it holds for any matrix A∈Rm×nthat A≤AF. Any estimation with respect to the Frobenius norm therefore also produces an upper bound for our standard norm ·. Specifically, the noise level εimg yields a corresponding value for εF. Similarly, it holds that G(X)=F(X)−δY(X), and therefore G(X)−F(X)F≤δY(X)F≤ εimg√MN, corresponding to a value for εF. References 1. Sun, S., Nocedal, J.: A trust region method for noisy unconstrained optimization. Math. Program. 66, 1–28 (2023). https://doi.org/10.1007/s10107-023-01941-9 2. Byrd, R.H., Gould, N.I., Nocedal, J., Waltz, R.A.: An algorithm for nonlinear optimization using linear programming and equality constrained subproblems. Math. Program. 100(1), 27–48 (2003). https:// doi.org/10.1007/s10107-003-0485-4 3. Byrd, R.H., Gould, N.I., Nocedal, J., Waltz, R.A.: On the convergence of successive linearquadratic programming algorithms. SIAM J. Optim. 16(2), 471–489 (2005). https://doi.org/10.1137/ S1052623403426532 4. Wright, S., Nocedal, J., et al.: Numerical Optimization, 2nd edn. Springer, Berlin (2006). https://doi. org/10.1007/978-0-387-40065-5 5. Tibshirani, R.: Regression shrinkage and selection via the lasso. J. R. Stat. Soc. Ser. B Methodol. 58(1), 267–288 (1996). https://doi.org/10.1111/j.2517-6161.1996.tb02080.x 6. Candes, E.J., Romberg, J.K., Tao, T.: Stable signal recovery from incomplete and inaccurate measurements. Commun. Pure Appl. Math. A J. Issued Courant Inst. Math. Sci. 59(8), 1207–1223 (2006). https://doi.org/10.1002/cpa.20124 7. Fukushima, K.: Cognitron: a self-organizing multilayered neural network. Biol. Cybernet. 20(3), 121– 136 (1975). https://doi.org/10.1007/BF00342633 8. Aspelmeier, T., Charitha, C., Luke, D.R.: Local linear convergence of the admm/douglas-rachford algorithms without strong convexity and application to statistical imaging. SIAM J. Imaging Sci. 9(2), 842–868 (2016). https://doi.org/10.1137/15M103580X 9. Cartis, C., Gould, N.I., Toint, P.L.: On the evaluation complexity of composite function minimization with applications to nonconvex nonlinear programming. SIAM J. Optim. 21(4), 1721–1739 (2011). https://doi.org/10.1137/11082381X 123
600 C. Hansknecht et al. 10. Apkarian, P., Noll, D., Ravanbod, L.: Nonsmooth bundle trust-region algorithm with applications to robust stability. Set-Valued Var. Anal. 24(1), 115–148 (2016). https://doi.org/10.1007/s11228-0150352-5 11. Boyd, S., Parikh, N., Chu, E., Peleato, B., Eckstein, J., et al: Distributed optimization and statistical learning via the alternating direction method of multipliers. Found. Trends®Mach. Learn. 3(1), 1–122 (2011). https://doi.org/10.1561/2200000016 12. Nesterov, Y.: Gradient methods for minimizing composite functions. Math. Program. 140(1), 125–161 (2013). https://doi.org/10.1007/s10107-012-0629-5 13. Santosa, F., Symes, W.W.: Linear inversion of band-limited reflection seismograms. SIAM J. Sci. Stat. Comput. 7(4), 1307–1330 (1986). https://doi.org/10.1137/0907087 14. Han, S.-P., Mangasarian, O.L.: Exact penalty functions in nonlinear programming. Math. Program. 17(1), 251–269 (1979). https://doi.org/10.1007/BF01588250 15. Shi, H.-J.M., Xie, Y., Byrd, R., Nocedal, J.: A noise-tolerant quasi-newton algorithm for unconstrained optimization. SIAM J. Optim. 32(1), 29–55 (2022). https://doi.org/10.1137/20M1373190 16. Irwin, B., Haber, E.: Secant penalized bfgs: a noise robust quasi-newton method via penalizing the secant condition. Comput. Optim. Appl. 84(3), 651–702 (2023). https://doi.org/10.1007/s10589-02200448-x 17. Cao, L., Berahas, A.S., Scheinberg, K.: First-and second-order high probability complexity bounds for trust-region methods with noisy oracles. Math. Program. 66, 1–52 (2023). https://doi.org/10.1007/ s10107-023-01999-5 18. Gal, R., Haber, E., Irwin, B., Saleh, B., Ziv, A.: How to catch a lion in the desert: on the solution of the coverage directed generation (cdg) problem. Optim. Eng. 22(1), 217–245 (2021). https://doi.org/ 10.1007/s11081-020-09507-w 19. Curtis, F.E., Scheinberg, K., Shi, R.: A stochastic trust region algorithm based on careful step normalization. Informs J. Optim. 1(3), 200–220 (2019). https://doi.org/10.1287/ijoo.2018.0010 20. Berahas, A.S., Byrd, R.H., Nocedal, J.: Derivative-free optimization of noisy functions via quasinewton methods. SIAM J. Optim. 29(2), 965–993 (2019). https://doi.org/10.1137/18M1177718 21. Xie, Y., Byrd, R.H., Nocedal, J.: Analysis of the BFGS method with errors. SIAM J. Optim. 30(1), 182–209 (2020). https://doi.org/10.1137/19M1240794 22. Gurobi Optimization, LLC: Gurobi Optimizer Reference Manual (2022). https://www.gurobi.com 23. Huangfu, Q., Hall, J.J.: Parallelizing the dual revised simplex method. Math. Program. Comput. 10(1), 119–142 (2018). https://doi.org/10.1007/s12532-017-0130-5 24. Byrd, R.H., Nocedal, J., Waltz, R.A.: Knitro: an integrated package for nonlinear optimization. In: Large-Sscale Nonlinear Optimization, pp. 35–59. Springer, Berlin (2006). https://doi.org/10.1007/0387-30065-1_4 25. Fletcher, R.: Practical Methods of Optimization: Vol. 2: Constrained Optimization (1981). https://doi. org/10.1002/9781118723203 26. Yuan, Y.-x: Conditions for convergence of trust region algorithms for nonsmooth optimization. Math. Program. 31(2), 220–228 (1985). https://doi.org/10.1007/BF02591750 27. Rockafellar, R.T., Wets, R.J.-B.: Variational Analysis, vol. 317. Springer, Berlin (2009). https://doi. org/10.1007/978-3-642-02431-3 28. Harris, C.R., Millman, K.J., van der Walt, S.J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N.J., Kern, R., Picus, M., Hoyer, S., van Kerkwijk, M.H., Brett, M., Haldane, A., del Río, J.F., Wiebe, M., Peterson, P., Gérard-Marchant, P., Sheppard, K., Reddy, T., Weckesser, W., Abbasi, H., Gohlke, C., Oliphant, T.E.: Array programming with NumPy. Nature 585(7825), 357–362 (2020). https://doi.org/10.1038/s41586-020-2649-2 29. Virtanen, P., Gommers, R., Oliphant, T.E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S.J., Brett, M., Wilson, J., Millman, K.J., Mayorov, N., Nelson, A.R.J., Jones, E., Kern, R., Larson, E., Carey, C.J., Polat, ˙ I, Feng, Y., Moore, E.W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E.A., Harris, C.R., Archibald, A.M., Ribeiro, A.H., Pedregosa, F., van Mulbregt, P.: SciPy 1.0 Contributors: SciPy 1.0: fundamental algorithms for scientific computing in python. Nat. Methods 17, 261–272 (2020). https://doi.org/10.1038/s41592-019-0686-2 30. Wächter, A., Biegler, L.T.: On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Math. Program. 106(1), 25–57 (2006). https://doi.org/10.1007/ s10107-004-0559-y 123
Convergence of successive linear programming… 601 31. Marsaglia, G.: Choosing a point from the surface of a sphere. Ann. Math. Stat. 43(2), 645–646 (1972). https://doi.org/10.1214/aoms/1177692644 32. Gould, N.I., Orban, D., Toint, P.L.: Cutest: a constrained and unconstrained testing environment with safe threads for mathematical optimization. Comput. Optim. Appl. 60(3), 545–557 (2015). https://doi. org/10.1007/s10589-014-9687-3 Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. 123