Funding: This material is based upon work supported by the National Science Foundation [Awards CCF-1740796 and CCF-2139735] and the Office of Naval Research [Award N00014-21-1-2532].
]]>We propose, analyze, and present experimental results with a stochastic inexact sequential quadratic optimization (SISQO) algorithm for minimizing an objective function subject to (s.t.) nonlinear equality constraints. Specifically, our algorithm is designed to solve problems of the form
where f : R n → R and c : R n → R m are continuously differentiable, ω is a random variable with probability space (Ω, F , P), F : R n × Ω → R, and E ω [•] represents expectation taken with respect to the distribution of ω. Problems of this type arise in numerous important application areas. A partial list is the following: (i) learning a deep convolutional neural network for image recognition that imposes properties (e.g., smoothness) of the systems of partial differential equations (PDEs) that the convolutional layers are meant to interpret (Ruthotto and Haber 2020); (ii) multiple deep learning problems (see Márquez-Neila et al. 2017), including physics-constrained deep learning for high-dimensional surrogate modeling and uncertainty quantification without labeled data (Zhu et al. 2019), natural language processing with constraints on output labels (Nandwani et al. 2019), image classification, detection, and localization (Ravi et al. 2019), deep reinforcement learning (Achiam et al. 2017), deep network compression (Chen et al. 2018), and manifold-regularized deep learning (Tomar andRose 2014, Kumar Roy et al. 2018); (iii) accelerating the solution of PDE-constrained inverse problems by using a reduced-order model in place of a full-order model coupled with techniques to learn the discrepancy between the reduced-and full-order models (Sheriffdeen et al. 2019); (iv) multistage modeling (Shapiro et al. 2014); (v) portfolio selection (Shapiro et al. 2014); (vi) optimal power flow (Summers et al. 2015); and (vii) statistical problems, such as maximum likelihood estimation with constraints (Geyer 1991, Chatterjee et al. 2016).
Let R denote the set of real numbers, R ≥p (respectively, R >p ) denote the set of real numbers greater than or equal to (respectively, strictly greater than) p ∈ R, and N :� {0, 1, 2, : : : } denote the set of natural numbers. Let R n denote the set of n-dimensional real vectors, R m×n denote the set of m-by n-dimensional real matrices, and S n denote the set of n-by n-dimensional symmetric real matrices. For any p ∈ N \ {0}, let [p] :� {1, : : : , p}. The ℓ 2 -norm is written simply as ‖ • ‖. Any run of our algorithm generates a sequence of iterates {x k }, where x k ∈ R n for all k ∈ N. For all k ∈ N, we append the subscript k to other quantities defined with respect to the kth iteration, and for brevity, we define ∇f k :� ∇f (x k ), c k :� c(x k ), and ∇c k � ∇c(x k ). We refer to the range space of ∇c k as Range(∇c k ) and refer to the null space of ∇c T k as Null(∇c T k ). The fundamental theorem of linear algebra provides that these spaces are orthogonal and Range(∇c k ) + Null(∇c T k ) � R n . Finally, recall (see, e.g., Nocedal and Wright 2006) that a primal point x ∈ R n and a dual point y ∈ R m constitute a first-order stationary point for Problem (1) if and only if c(x) � 0 and ∇f (x) + ∇c(x)y � 0. These conditions are necessary for x to be a local minimizer when the constraint functions satisfy a constraint qualification as is assumed in the paper; see Assumption 1.
Our algorithm is presented in Section 2, with a convergence analysis in Section 3. The results of numerical experiments are presented in Section 4, and concluding remarks are presented in Section 5.
A run of our proposed SISQO algorithm generates a sequence
where for all k ∈ N, (x k , y k ) ∈ R n × R m is a primal-dual iterate pair; g k ∈ R n is a stochastic gradient estimate; v k ∈ R n is a normal primal direction that aims to reduced infeasibility by reducing a local derivative-based model of the ℓ 2 -norm constraint violation measure; u k ∈ R n is a tangential primal direction that aims to maintain the reduction in linearized infeasibility achieved by the normal primal direction while also aiming to reduce the objective by reducing a stochastic gradient-based quadratic approximation of the objective; d k :� v k + u k ∈ R n is a full primal direction; δ k ∈ R m is a dual direction; (ρ k , r k ) ∈ R n × R m is a primal-dual linear system residual pair; (τ trial k , τ k ) ∈ R >0 ∪ {∞} × R >0 is a pair of trial and employed merit parameter values; (ξ trial k , ξ k ) ∈ R >0 × R >0 is a pair of trial and employed ratio parameter values; and (α k, min , α k, max , α k ) ∈ R >0 × R >0 × R >0 is a tuple of minimum, maximum, and employed step size values, the last of which aims to produce a subsequent iterate x k+1 ← x k + α k d k yielding sufficient reduction in the ℓ 2 -norm merit function. We present our algorithm in the context of the generation of a realization of the sequence (2), although our analysis ultimately considers the stochastic process defined by the algorithm, namely
of which ( 2) is a realization. In the rest of this section, we present our algorithm in the context of the generation of a realization (2) toward our complete statement of Algorithm 1.
For the remainder of the paper, we make the following assumption.
Assumption 1. Let X ⊆ R n be an open convex set containing {x k } generated by every run of Algorithm 1. The objective f : R n → R is continuously differentiable and bounded below over X , and its gradient ∇f : R n → R n is Lipschitz continuous with constant L ∈ R >0 (with respect to the ℓ 2 -norm) and bounded over X . The constraint c : R n → R m (with m ≤ n) is continuously differentiable and bounded over X, and its Jacobian ∇c(•) T : R n → R m×n is Lipschitz continuous with constant Γ ∈ R >0 (with respect to the vector-induced ℓ 2 -norm) and bounded over X . In addition, for all x ∈ X , the singular values of ∇c(x) T are bounded uniformly below by a positive real number.
The elements of this assumption are standard in the continuous constrained optimization literature. Note that it does not include an assumption that X is bounded. One could relax the assumption to say that X contains {X k } almost surely, but because such a relaxation would only make it necessary to remark constantly on the probabilityzero event that {X k } ⊄ X without adding much value to our ultimate results, we employ Assumption 1 as it is stated.
Motivated by the success of numerous line-search-SQO methods for solving deterministic equality constrained optimization problems, our algorithm employs an exact penalty function as a merit function; in particular, it
Curtis, Robinson, and Zhou: Stochastic Inexact Sequential Quadratic Optimization INFORMS Journal on Optimization, Articles in Advance, pp. 1-23, © 2024 INFORMS
employs the ℓ 2 -norm merit function φ : R n × R >0 → R defined by φ(x, τ) � τf (x) + ‖c(x)‖, where τ is a positive merit parameter that is updated adaptively. (The choice of the ℓ 2 -norm in φ is not essential. Another norm could be used. The choice of the ℓ 2 -norm merely makes certain calculations simpler for our presentation and analysis.) A model l : T d‖, with which we define the model reduction function ∆l :
The merit and model reduction functions play critical roles in our inexactness conditions.
During the kth iteration, the algorithm computes a normal direction
Instead of solving (5) exactly, the algorithm allows for an inexact solution to be employed by only requiring the computation of v k ∈ Range(∇c k ) satisfying the Cauchy decrease condition
where ɛ c ∈ (0, 1] is user defined. In ( 6), v c k :� �∇c k c k is the negative-gradient direction for the objective of ( 5) at v � 0, and
otherwise, if c k � 0, then ∇c k c k � 0, and v k � 0 is the unique solution to (5). Popular choices for computing a normal direction satisfying the aforementioned conditions include Krylov subspace methods, such as the linear conjugate gradient (CG) method; see, for example, Nocedal and Wright (2006). For describing the tangential direction computation, let us first describe what would be the computation in a deterministic variant of our approach. Given (x k , y k ), ∇f k , v k ∈ Range(∇c k ), and H k ∈ S n satisfying Assumption 2, consider the quadratic optimization subproblem
which has the unique solution
We make the following assumption pertaining to {H k } throughout the paper.
Assumption 2. There exists ζ ∈ R >0 and κ H ∈ R >ζ such that for all k ∈ N in any run of the algorithm, ‖H k ‖ ≤ κ H and u T H k u ≥ ζ‖u‖ 2 for all u ∈ Null(∇c T k ). The introduction of ( 8) and ( 9) allows us to define, for the purposes of our analysis only (i.e., not for actual computation), the true and exact primal-dual search direction conditioned on the behavior of the algorithm up to the kth iteration as
Because our algorithm only presumes access to a stochastic gradient estimate g k of ∇f k , the corresponding exact, but not true, primal-dual search direction is
(For the description of our algorithm and our initial analysis, it suffices that g k ∈ R n . Our ultimate required assumption about the stochastic gradient estimators is Assumption 6.) Our algorithm, to avoid having to form or factorize the matrix in (10), computes a tangential direction by computing (u k , δ k ) through iterative linear algebra techniques applied to the symmetric indefinite system (10). In particular, our algorithm computes (u k , δ k ) such that the full primal search direction d k :� v k + u k , dual direction δ k , and residual
satisfy at least one of two sets of conditions. Next, we describe these conditions that the algorithm employs to determine what constitutes an acceptable search direction and corresponding residuals.
In the deterministic setting, line-search-SQO methods commonly combine the search direction with an updating strategy for the merit parameter in a manner that ensures that the computed direction is one of sufficient descent for the merit function. The required descent condition is guaranteed to be satisfied by choosing the merit parameter to be sufficiently small so that the reduction in a model of the merit function (recall (4)) is sufficiently large; see, for example, Byrd et al. (2008, lemma 3.1). Following such an approach, our algorithm requires that (u k , δ k ) (yielding d k :� v k + u k ) be computed and τ be set such that the model reduction condition
holds for some user-defined σ u ∈ (0, 1), ɛ u ∈ (0, ζ) (see Assumption 2), and σ c ∈ (0, 1). The value for τ for which ( 12) is required to hold depends on one of two different situations as described next. Condition (12) plays a central role in the conditions that we require (u k , δ k ) to satisfy. We define these in the context of termination tests (TTs) because they dictate conditions that once satisfied, can cause termination of an iterative linear system solver applied to (10). (The tests are inspired by the sufficient merit approximation reduction termination tests developed in Byrd et al. (2008Byrd et al. ( , 2010) ) and Curtis et al. (2009) for a deterministic SQO method.) Our first termination test states that an inexact solution is acceptable if ( 12) is satisfied with the current merit parameter value (i.e., τ ≡ τ k ← τ k�1 ), the norms of the residuals satisfy certain upper bounds, and either the tangential direction is sufficiently small in norm compared with the normal direction or the tangential direction is one of sufficiently positive curvature for H k and yields a sufficiently small objective value for (8) (with g k in place of ∇f k ). The test makes use of a sequence {β k } that will also play a critical role in our step size selection scheme that is described in the next subsection.
, and v k ∈ Range(∇c k ) computed to satisfy (6), the pair (u k , δ k ) satisfies TT1 if with the pair (ρ k , r k ) defined in (11), it holds that
and ( 12) is satisfied with τ ≡ τ k�1 . (In this case, τ k ← τ k�1 , so (12) holds with τ ≡ τ k .) TT1 cannot be enforced in every iteration, even in the deterministic setting, because there may exist points in the search space at which all of the conditions required cannot be satisfied simultaneously, even if the linear system (10) is solved accurately. In short, the algorithm needs to allow for the computation of a search direction for which (12) can only be satisfied with a decrease of the merit parameter. That said, the algorithm needs to be careful in terms of the situations in which such a decrease is allowed, or else, the merit parameter sequence may behave in a manner that ruins a convergence guarantee for solving the original constrained optimization problem. For our algorithm, we employ the following termination test for this situation.
), and v k ∈ Range(∇c k ) computed to satisfy (6), the pair (u k , δ k ) satisfies Termination Test 2 if with the pair (ρ k , r k ) defined in (11), ( 13)-( 15) and
Curtis, Robinson, and Zhou: Stochastic Inexact Sequential Quadratic Optimization INFORMS Journal on Optimization, Articles in Advance, pp. 1-23, © 2024 INFORMS
hold. (In this case, for user-defined ɛ τ ∈ (0, 1), the algorithm will set
where
(18) so ( 12) is satisfied with τ ≡ τ k . See Lemma 3 for a proof.)
In Lemma 1, we show under a loose assumption about the iterative linear system solver that for all k ∈ N, the algorithm can compute a pair (u k , δ k ) satisfying at least one of TT1 or TT2. Therefore, the index of each iteration of each realization of our method is contained in one of two index sets:
It is worthwhile to emphasize that in terms of the stochastic process defined by the algorithm, the index sets K 1 and K 2 are random (i.e., they may contain different indices in different runs). This randomness is handled as part of our convergence analysis.
Upon computation of d k ← v k + u k , our algorithm computes a positive step size α k to set x k+1 . Given positive Lipschitz constants L and Γ (recall Assumption 1), it follows for all α ∈ R >0 that
Combining these with (4) yields
This derivation provides a convex piecewise quadratic upper-bounding function for the change in the merit function corresponding to a step from x k to x k + αd k . Given user-defined η ∈ (0, 1) and the aforementioned sequence {β k } ⊂ (0, 1], our algorithm's step size selection scheme makes use of
The definition of α suff k can be motivated as follows. Its value, when β k � 1, is the largest in [0, 1] such that for all α ∈ [0, α suff k ], the right-hand side of (20) (with ∇f k replaced by g k ) is less than or equal to �ηα∆l(x k , τ k , g k , d k ). Such an inequality is representative of one enforced in deterministic line-search-SQO methods. Otherwise, with β k ∈ (0, 1] introduced and not necessarily equal to one, the value of α suff k can be diminished during the optimization, which allows for step size control as is required for convergence guarantees for certain stochastic gradient-based methods; see, for example, Bottou et al. (2018). The first term inside the min appearing in ( 21) is important for the convergence guarantees that we prove for our method, but it can behave erratically because of the algorithm's use of stochastic gradients. To account for this, given user-defined ɛ ξ ∈ (0, 1), our algorithm defines
Combining this inequality with (21), the monotonically nonincreasing behaviors of {ξ k } and {τ k }, and assuming that
where ξ �1 and τ �1 initialize the sequences {ξ k } and {τ k }, respectively, one finds that
The value α min k serves as a minimum value (i.e., a lower bound) for our choice of step size. In our analysis, we show that the sequence {ξ k } is bounded below and away from zero by a positive real number that is common to all runs of the algorithm (see Lemma 9).
Next, let us derive a maximum value (i.e., an upper bound) for our algorithm's choice of step size. If d k � 0, then without loss of generality, the algorithm can set α k ← α k, min . Otherwise, if d k ≠ 0, consider the strongly convex function φ : R → R defined by
Notice that when β k � 1, it holds that φ(α) ≤ 0 for all α ∈ R ≥0 if and only if the right-hand side of (20) with ∇f k replaced by g k is less than or equal to �ηα∆l(x k , τ k , g k , d k ). Thus, following a similar argument, one can be motivated as to the fact that our algorithm never allows a step size larger than α
Finally, again to mitigate adverse effects caused by the use of stochastic gradients, our algorithm employs the maximum step size
where θ ∈ R ≥0 is user defined. Overall, our algorithm allows any step size with
Lemma 4 in our analysis shows that this interval is nonempty.
In the primal space, our algorithm employs the iterate update x k+1 ← x k + α k d k . However, in the dual space, it allows additional flexibility. For consistency with the deterministic setting (see, e.g., Curtis et al. 2009, equation 2.19), our algorithm is stated to require y k+1 to satisfy
Clearly, choosing y k+1 ← y k + δ k is one particular option satisfying (27), although other choices, such as leastsquares multipliers, could also be used.
Algorithm 1 (Stochastic Inexact Sequential Quadratic Optimization (SISQO))
Require:
using the definitions in ( 24) and ( 26) 15: set x k+1 ← x k + α k d k , and choose y k+1 satisfying ( 27) 16: end for
Our analysis is presented in three parts. In Section 3.1, we show that Algorithm 1 is well posed, which is followed by Section 3.2, in which a set of lemmas is proved that hold for every run of the algorithm. The analysis in Sections 3.1 and 3.2 is written generically in terms of a realization of a run of the algorithm. Then, in Section 3.3, we prove convergence properties for the algorithm that are written in terms of the stochastic process defined by the algorithm. This analysis focuses on an event that presumes certain behavior of the merit and ratio parameter sequences.
Our aim in this subsection is to prove that each procedure in each iteration of every run of Algorithm 1 is performed in a manner that terminates finitely under mild assumptions. Along the way, we establish other useful properties. In this subsection and the following one, our analysis merely requires for all k ∈ N that g k ∈ R n and that H k ∈ S n satisfies Assumption 2.
For the sake of generality, we make the following assumption about the iterative linear system solver employed in line 6 of Algorithm 1, which merely requires that the residual of the linear system solve vanishes asymptotically over any run of the solver. We emphasize, however, that there exist approaches, such as MINRES (Paige and Saunders 1975), that guarantee that an exact solution-which would satisfy our termination tests-can be computed in a number of linear system solver iterations that is bounded uniformly for all linear systems arising throughout any realization of the algorithm. However, we merely make the following assumption because it is all that is required by our analysis, and it offers more flexibility in the choice of linear system solver. Assumption 3. For all k ∈ N in any run, the iterative linear system solver employed in line 7 of Algorithm 1 to compute
We also make the following assumption concerning the algorithm iterates and corresponding stochastic gradient estimates computed in each iteration. Assumption 4. For all k ∈ N in any run, it holds that c k ≠ 0 or g k ∉ Range(∇c k ).
We justify Assumption 4 in the following manner. In the deterministic setting, the algorithm encounters a point x k such that c k � 0 and ∇f k ∈ Range(∇c k ) if and only if there exists y k such that (x k , y k ) is first-order stationary for Problem (1). In such a scenario, it is reasonable to require that an exact solution of ( 10) is computed or at least that a sufficiently accurate solution is computed such that a practical termination condition for (10) is triggered and the algorithm terminates. In the stochastic setting, the algorithm encounters c k � 0 and g k ∈ Range(∇c k ) if and only if x k is exactly feasible and the stochastic gradient lies exactly in the range space of ∇c k . Because g k is a stochastic gradient, we contend that it is unlikely that it will lie exactly in Range(∇c k ), except in special circumstances. Thus, for simplicity in our analysis, we impose Assumption 4. Note that if Assumption 4 did not hold, then one of the following could be employed in a practical implementation. (i) If a sufficiently accurate solution of (10) satisfies neither TT1 nor TT2, then a new stochastic gradient could be sampled, perhaps following a procedure to ensure that if multiple new stochastic gradients are computed, then each is computed with lower variance, or (ii) random (e.g., Gaussian) noise could be added to g k for all k ∈ N so that Assumption 4 holds almost surely in all iterations, in which case the convergence result that we prove holds almost surely.
We can now show that the search direction computation is well posed. We remark in passing that if one was to employ a linear system solver, such as MINRES, that would produce an exact solution of the linear system within a uniformly bounded number of iterations, then the arguments in the proof of the following lemma would show that the linear system solver computes (u k , δ k ), satisfying at least one of TT1 or TT2 in a uniformly bounded number of iterations.
Lemma 1. For all k ∈ N in any run, the iterative linear system solver computes (u k , δ k ) satisfying at least one of TT1 or TT2 in a finite number of iterations.
Proof. We prove the result by considering two cases.
Case 1. c k ≠ 0. For this case, we show that (u k , δ k ) ≡ (u k, t , δ k, t ) satisfies TT2 for sufficiently large t ∈ N. Let us first observe that it follows from Assumption 3, Assumption 4, and β k ∈ (0, 1] that ( 13) and ( 14) hold with (ρ k , r k ) ≡ (ρ k, t , r k, t ) for all sufficiently large t ∈ N.
Let us now show that ( 15) holds for all sufficiently large t ∈ N. Because c k ≠ 0, it follows under Assumption 1 that v k ≠ 0. If u k, * � 0, then Assumption 3 implies {u k, t } → u k, * � 0, in which case it follows from κ u ∈ R >0 that the former condition in (15) holds with u k ≡ u k, t for all sufficiently large t ∈ N. On the other hand, if u k, * ≠ 0, then (10) and Assumption 2 imply
Combining this inequality with ɛ u ∈ (0, ζ), κ v ∈ R >0 , v k ≠ 0, and Assumptions 2 and 3, it follows that the latter set of conditions in (15) holds with u k ≡ u k, t for all sufficiently large t ∈ N.
Finally, let us show that ( 16) holds for all sufficiently large t ∈ N, which combined with the previous conclusions, shows that TT2 is satisfied by (u k , δ k ) ≡ (u k, t , δ k, t ) for all sufficiently large t ∈ N. By Assumption 3, ( 6), and
, which shows that ( 16) holds with r k ≡ r k, t for all sufficiently large t ∈ N.
Case 2. c k � 0. For this case, we show that (u k , δ k ) ≡ (u k, t , δ k, t ) satisfies TT1 for all sufficiently large t ∈ N. First, recall that c k � 0 implies that v k � 0. We also claim that u k, * ≠ 0. To prove this by contradiction, suppose that u k, * � 0. Combining this fact with v k � 0 and (10), it follows that g k + ∇c k (y k + δ k, * ) � 0, which with c k � 0, violates Assumption 4. Thus, u k, * ≠ 0.
Next, notice that the argument used in the beginning of case (1) still applies in this case, which allows us to conclude that both ( 13) and ( 14) hold with (ρ k , r k ) � (ρ k, t , r k, t ) for all sufficiently large t ∈ N. Combining (29) with Assumptions 2 and 3 and ɛ u ∈ (0, ζ) allows us to deduce that the second set of conditions in (15) holds with u k ≡ u k, t for all sufficiently large t ∈ N. Next, from the fact that v k � 0 and (10), it follows that ∇c 10), and Assumption 2 shows that 12) holds with τ ≡ τ k�1 for all sufficiently large t ∈ N. In summary, we have shown that for all sufficiently large t ∈ N, the pair (u k , δ k ) ≡ (u k, t , δ k, t ) satisfies TT1. w Next, we prove that every full primal direction is nonzero.
Lemma 2. For all k ∈ N in any run, it holds that d k ≠ 0.
then this shows that the inequality in (13) cannot hold, meaning that (u k , δ k ) satisfies neither TT1 nor TT2, which contradicts Lemma 1. Hence, the only possibility is that c k ≠ 0, which we assume for the rest of the proof.
Notice that from 16) is not satisfied; thus, (u k , δ k ) does not satisfy TT2. Also, observe from v k ≠ 0 (which follows from c k ≠ 0 and Assumption 1), d k � 0, and ( 6 12) is not satisfied with τ � τ k�1 ; thus, (u k , δ k ) does not satisfy TT1. Overall, we have reached a contradiction to Lemma 1, and because we have reached a contradiction in all cases, the original supposition that d k � 0 cannot be true. w
We now show that our update strategy for the merit parameter ensures that the model reduction condition (12) always holds for τ ≡ τ k . We also show another important property of {τ k }. Lemma 3. For all k ∈ N in any run, the inequality in (12) holds with τ ≡ τ k . In addition, for all k ∈ N such that
Proof. The desired conclusion follows for k ∈ K 1 because of the manner in which TT1 is defined and the fact that the algorithm sets τ k ← τ k�1 for all k ∈ K 1 . Hence, let us proceed under the assumption that k ∈ K 2 . The inequality in (12) holds for τ ≡ τ k with
We now proceed to show that this inequality holds by considering two cases.
Curtis, Robinson, and Zhou: Stochastic Inexact Sequential Quadratic Optimization INFORMS Journal on Optimization, Articles in Advance, pp. 1-23, © 2024 INFORMS Case 1. g
In this case, the algorithm sets τ k ← τ k�1 . Combining this with ( 16), ∇c T k u k � r k , and ɛ r ∈ (σ c , 1) yields
The update (17) yields τ k ≤ τ trial k , which combined with ( 16), ( 18), ∇c T k u k � r k , and ɛ r ∈ (σ c , 1) yields
as desired. Moreover, from (17), we have
We conclude this subsection by showing that the interval defining our step size selection scheme (i.e., [α min k , α max k ]) is positive and nonempty for all k ∈ N. We also show a useful property of the computed step size that is needed in our analysis.
Proof. It follows from ( 24) and the fact that {β k }, {ξ k }, and {τ k } are positive sequences that α min k > 0 for all k ∈ N. Hence, considering ( 24) and ( 26), to prove that 0
, which with (25), shows that
For this case, it follows from (25), α suff k ∈ (0, 1], and the triangle inequality that 4) and ( 25), one finds (as previously mentioned) that φ is strongly convex. In addition, one finds that φ(0
In this subsection, we prove general results about the behavior of Algorithm 1. As in the previous subsection, our analysis here merely requires for all k ∈ N that g k ∈ R n and H k ∈ S n satisfy Assumption 2. The next lemma gives a lower bound on ‖c k ‖ � ‖c k + ∇c T k v k ‖ relative to ‖c k ‖. Lemma 5. There exists κ 1 ∈ R >0 (a constant common to all runs of the algorithm) such that for all k ∈ N in any run, one has that
Proof. This result follows as in Curtis et al. (2009, lemma 3.5) but with small straightforward modifications to account for the fact that in our analysis here, the singular values of {∇c T k } are bounded away from zero uniformly over all runs as part of Assumption 1. w The next lemma shows that ‖v k ‖ is bounded below and above proportionally to ‖c k ‖. Lemma 6. There exists (κ 2 , κ 3 ) ∈ R >0 × R >0 (constants common to all runs of the algorithm) such that for all k ∈ N in any run, one has that κ 2 ‖c k ‖ ≤ ‖v k ‖ ≤ κ 3 ‖c k ‖.
Proof. Assumption 1 ensures the existence of λ min ∈ R >0 such that ∇c T k ∇c k ≽ λ min I for all k ∈ N in any run. We now prove each desired inequality. First, consider the former inequality. Because this holds trivially when c k � 0, let us proceed under the assumption that c k ≠ 0. One finds
It follows from this inequality, the triangle inequality, and ( 6) that
Substituting in for the value of α c k (recall ( 7)) and
Again, substituting the value of α c k and using the definition of λ min , one finds
It follows from these inequalities and Assumption 1 that there exists κ 2 ∈ >0 as claimed.
Let us now turn to the latter inequality. It follows from the normal direction computation that
Putting these facts together shows that
which combined with Assumption 1-namely, that the Jacobian function ∇c(•) T is bounded over the set X containing the iterates-establishes the existence of κ 3 ∈ R >0 as claimed. w
The next result gives a useful bound on the size of the search direction.
Lemma 7. There exists κ 4 ∈ R ≥2 (a constant common to all runs of the algorithm) such that for all k ∈ N in any run, one finds that
, the triangle inequality, and Lemma 6, it follows that
The existence of the required κ 4 ∈ R ≥2 now follows from Assumption 1 because max{2, 2κ 2 3 ‖c k ‖} is uniformly bounded for all k ∈ N in any run, which completes the proof. w
The next lemma shows that the model reduction ∆l(x k , τ k , g k , v k + u k ) is bounded below by a similar quantity as the upper bound for ‖d k ‖ 2 in the previous lemma.
Lemma 8. There exists κ 5 ∈ R >0 (a constant common to all runs of the algorithm) such that for all k ∈ N in any run, one has that ∆l
Lemma 3 shows that (12) holds with τ ≡ τ k . Combining this fact with Lemma 5 and the monotonically nonincreasing behavior of {τ k } shows that
Curtis, Robinson, and Zhou: Stochastic Inexact Sequential Quadratic Optimization INFORMS Journal on Optimization, Articles in Advance, pp. 1-23, © 2024 INFORMS which proves the existence of the claimed κ 5 ∈ R >0 because σ u , ɛ u , σ c , κ 1 , and τ �1 are positive real numbers. The remaining inequalities follow from Lemmas 2 and 7. w
We next prove a lower bound for {ξ k } that is uniform over every run of the algorithm.
Lemma 9. There exists ξ min ∈ R >0 (a constant common to all runs of the algorithm) such that in any run, there exists k ξ ∈ N and
Proof. For all k ∈ N, it follows from ( 22) and Lemmas 7 and 8 that
Now, consider any iteration such that ξ k < ξ k�1 . For such iterations, it follows from ( 22) and ( 30) that
Combining this fact with the initial choice of ξ �1 shows that ξ k ≥ ξ min :� min{(1 � ɛ ξ )κ 5 =κ 4 , ξ �1 } for all k ∈ N. Combining this result with the that ξ k < ξ k�1 implies that ξ k ≤ (1 � ɛ ξ )ξ k�1 (it decreases by at least a factor of 1 � ɛ ξ ) gives the desired result. w
The next lemma gives a bound on the change in the merit function in each iteration.
Lemma 10. For all k ∈ N in any run, one has that
Proof. By Lemma 4, one has φ(α k ) ≤ 0. Hence, starting as in (20); adding and subtracting the terms
; using the definition of φ(•); and using the fact that ∇c
Our goal in this subsection is to prove a convergence result for our algorithm. In general, in a run of the algorithm, one of three possible events can occur with respect to the merit parameter sequence. One possible event is that the merit parameter sequence eventually remains constant at a value that is sufficiently small. This is the event that we consider in our analysis here, where the meaning of sufficiently small is defined formally in the event E that is introduced shortly. The other two possible events are that the merit parameter sequence vanishes or eventually remains constant at a value that is too large. As discussed in Berahas et al. (2021, section 3.2.2), the former of these two events does not occur for the algorithm in that paper if the differences between the stochastic gradient estimates and the true gradients of the objective are uniformly bounded in norm; in particular, see Berahas et al. (2021, proposition 3.18). It is straightforward to show that such a conclusion also holds for Algorithm 1 in this paper because the merit parameter update strategy follows the same kind of approach as for the algorithm in Berahas et al. (2021); in particular, see the consistency between Berahas et al. (2021, equations (3.3) and (3.4)) and in this paper, ( 17) and (18) as well as Byrd et al. (2008, lemma 4.7), which considers the setting of inexact linear system solutions using inexactness tolerance conditions of the same type as in this paper. Moreover, in Berahas et al. (2021, section 3.2.2), it is shown that the latter type of event (namely, that the merit parameter remains constant at a value that is too large) occurs with probability of zero if one makes a reasonable assumption about the influence of the stochastic gradient estimates on the computed search directions; see also Berahas et al. (2023, section 4.3) for additional discussion of this case in the context of an algorithm that employs a step decomposition approach, like in Algorithm 1. Again, it is straightforward to see that such a conclusion also holds for Algorithm 1 because the merit parameter update strategy is of the same form. Consequently, for our purposes here, we do not consider these latter events because we contend that for practical purposes, one can focus on the first event for the same reasons as in Berahas et al. (2021Berahas et al. ( , 2023)).
Our main convergence result for Algorithm 1 considers an assumption that combines all of the assumptions required for our analysis until this point and assumes certain behavior of the merit parameter sequence through an event denoted as E. For this event and the subsequent analysis, recall the stochastic process (3) defined by the algorithm. Consider for each k ∈ N the condition
similar to the one appearing in ( 18). (For the sake of brevity in our notation, we overload the meaning of H k ; here, it may be a random variable satisfying Assumption 6 introduced shortly, which is consistent with the previously introduced and employed Assumption 2.) With this condition, let us define the following trial value of the merit parameter that would be computed in iteration k ∈ N if the algorithm was to employ ∇f (X k ) in place of G k and solve (8) exactly:
(To be clear, the quantity T trial, true k never needs to be computed by our algorithm; it is only used in our analysis.) Using this quantity, we define our event of interest, namely E, as the following.
For some (k min , τ min , f sup ) ∈ N × R >0 × R, the event E :� E(k min , τ min , f sup ) occurs if and only if f (X k min ) ≤ f sup , and there exists (K, T ′ , Ξ ′ ) ∈ N × R >0 × R >0 with K ′ ≤ k min , T ′ ≥ τ min , and Ξ ′ ≥ ξ min (see Lemma 9) such that
In other words, event E is the event in which by iteration k min , the merit and ratio parameter sequences become constant at values at least τ min and ξ min , respectively, and the objective value at iteration k min is bounded. With respect to this event, we make the following assumption for the rest of our analysis, our ultimate focus of which will be on the behavior of the algorithm starting in iteration k min , at which point the adaptive merit and ratio parameters are constant.
Assumption 5. For some (k min , τ min , f sup ) ∈ N × R >0 × R, the event E :� E(k min , τ min , f sup ) occurs, and conditioned on the occurrence of E, Assumptions 1, 2, 3, and 4 hold. In addition, along with the restrictions that {β k } ⊂ (0, 1] and (23) holds for all k ∈ N, the sequence {β k } k≥k min is chosen in a manner that is F k min -measurable.
A few remarks are in order with respect to Assumption 5. First, that the event E includes that the merit parameter sequence is bounded below can, as previously mentioned, be justified for the same reasons as in Berahas et al. (2021); the additional requirement that it eventually remains constant can be justified by Lemma 3, which shows that if the merit parameter is decreased, then it is decreased by at least a constant factor. Note that it is the inequality T ′ ≤ T trial, true k that represents the aforementioned notion of the merit parameter ultimately being sufficiently small for all large k. Second, that E includes that the ratio parameter sequence is bounded below and eventually remains constant is not a strong assumption; it follows under our prior assumptions (that are carried forward in Assumption 5) because of Lemma 9. That said, the critical aspect here is that E requires that this sequence has become constant by iteration k min . Third, that E includes that f (X k min ) is bounded above is a relatively weak assumption but
Curtis, Robinson, and Zhou: Stochastic Inexact Sequential Quadratic Optimization INFORMS Journal on Optimization, Articles in Advance, pp. 1-23, © 2024 INFORMS
necessary for our purposes of ultimately showing that a sequence of stationarity measures vanishes in expectation. Finally, with respect to {β k } k≥k min , we state in Theorem 1 particular choices satisfying Assumption 5 for which our convergence guarantees hold. Precise strategies for setting these values that are consistent with our convergence guarantees are stated after Theorem 1.
Let G 0 be the σ-algebra defined by the initial conditions of Algorithm 1, and for all k ∈ N, let G k be the σ-algebra generated by the initial conditions and {G 0 , : : : , G k�1 }. In addition, for all k ∈ N, let the trace σ-algebra of E on G k be F k :� G k ∩ E. Hence, {F k } is a filtration. For brevity, let
where P ω denotes probability with respect to the distribution of ω (and as for ( 1), E ω denotes expectation with respect to the distribution of ω). Observe that conditioned on E, one has that
and one has that the random variables T ′ and Ξ ′ are F k -measurable for k � k min ≥ K ′ . We make Assumption 6 about {G k } and {H k }. That the stochastic gradient estimators are unbiased is standard for algorithms based on stochastic approximation. One may be able to relax the so-called bounded-variance assumption introduced here, but we contend that this assumption is sufficient for showing the general type of convergence guarantee that our algorithm offers. Hence, we make a bounded-variance assumption here so as not to obfuscate the other details. Assumption 6. There exists M g ∈ R >0 such that for all k ∈ N, the gradient estimator
In addition, for all k ∈ N, the matrix H k (satisfying Assumption 5; i.e., the bounds in Assumption 2) is F k -measurable.
Combining Assumption 6 with Jensen's inequality, it holds for all k ∈ N that
ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi
We now return to our analysis. First, let us derive bounds on the expected difference between U k and U true k . To that end, let us define Z k ∈ R n× (n�m) as an (F k -measurable) matrix whose columns form an orthonormal basis for Null(∇c(X k ) T ), which implies that Z T k Z k � I and ∇c(X k ) T Z k � 0. Under Assumption 5 (namely, Assumption 1), let U k, 1 ∈ R m and U k, 2 ∈ R n�m be vectors forming the orthogonal decomposition of U k into Range(∇c(X k )) and Null(∇c
The corresponding values for D k and D true k are found to be
In the proof of our next lemma, we use the fact that
which can be seen as follows. The nonzero eigenvalues of AB are equal to those of BA when the products are valid, meaning that the nonzero eigenvalues of
, which are all one; hence, the bound in (37) holds.
Lemma 11. There exists κ 6 ∈ R >0 such that for all k ∈ N with k ≥ k min , one finds
ffi ffi ffi ffi ffi ffi ffi M g q + κ 6 β k :
Proof. It follows from ( 35) that
which combined with Assumption 6, shows that
Combining this equation with the triangle inequality, Assumption 5 (specifically, Assumptions 1 and 2), ( 14), and (37) ensures the existence of κ 6 ∈ R >0 such that for all k ∈ N,
which is the first desired result. Next, to derive the desired bound on
one can combine the expression for U k � U true k with the triangle inequality to obtain
Taking conditional expectation and using Assumption 6, (34), (37), and ( 14), one finds
where κ 6 is the same value as used, which completes the proof. w
We now bound the difference in expectation between ∇f (X k
For the first term on the right-hand side, the Cauchy-Schwarz inequality,
Lemma 11, and Assumption 5 (i.e., Assumption 1) imply that there exists
Now, for the second term, first observe from Assumption 6 that
Combining this fact with (35), the Cauchy-Schwarz inequality, Assumption 5 (namely, Assumptions 1 and 2),
Curtis, Robinson, and Zhou: Stochastic Inexact Sequential Quadratic Optimization INFORMS Journal on Optimization, Articles in Advance, pp. 1-23, © 2024 INFORMS (37), and (34) shows that there exists (κ 8 , κ 8
Combining the results gives the desired result. w
We now proceed to bound in expectation the last few terms appearing in the right-hand side of the inequality proved in Lemma 10. The next lemma considers the last pair of terms.
Lemma 13. There exists κ 9 ∈ R >0 such that for all k ∈ N with k ≥ k min , one finds
Proof. From ( 11), ( 14), the fact that
, ( 24), ( 23), and the monotonically nonincreasing behavior of {T k } and {Ξ k }, it follows that there exists κ 9 ∈ R >0 such that
which gives the desired conclusion. w
Our next result provides an upper bound in expectation for the second term appearing on the right-hand side of the inequality in Lemma 10. Lemma 14. There exists κ 10 ∈ R >0 such that for all k ∈ N with k ≥ k min , one finds
Proof. Let I k be the event that ∇f (X k ) T (D k � D true k ) ≥ 0, and let I c k be its complementary event. It follows from (32), the definition of I k , the fact that
and the law of total expectation that for all k ≥ k min , one finds Combining this with the fact that (26
and the law of total expectation shows for all k ≥ k min that
Combining this with Lemma 11, (23
and Assumption 1 shows that there exists κ 10 ∈ R >0 such that for all k ≥ k min , one finds
which is the desired conclusion. w
We now use the model reduction based on the true step D true k to show an upper bound on the expected reduction in the model based on the step D k .
Proof. It follows from Lemma 12; (4); the fact that
, and D true k are all F k -measurable for k ≥ k min ; (9); and ( 14) that for all k ≥ k min ,
which is the desired result. w
We now prove our main result. In the result, the quantity ∆l(X k , T k , ∇f (X k ), D true k ) serves as a measure of stationarity with respect to (1); after all, the proof for Lemma 8 shows, with (∇f
Thus, in a run, if there exists infinite K ⊆ N with lim k∈K, k→∞ ∆l(x k , τ k , ∇f k , d true k ) � 0, then (38) and Lemma 6 imply that lim k∈K, k→∞ ‖c k ‖ � lim k∈K, k→∞ ‖u true k ‖ � lim k∈K, k→∞ ‖v k ‖ � 0, which combined with (9), shows that any limit of {(x k , y k + δ true k )} is a first-order stationary point for (1). In our stochastic setting, we prove for two different choices of {β k } k≥k min that an expected average measure of stationarity exhibits desirable properties. These properties match those ensured by a stochastic gradient method in the unconstrained setting (where ‖∇f (X k )‖ 2 is the measure of stationarity).
Then, the following results hold.
Curtis, Robinson, and Zhou: Stochastic Inexact Sequential Quadratic Optimization INFORMS Journal on Optimization, Articles in Advance, pp. 1-23, © 2024 INFORMS i. If β k � β � ψA ′
(1�η)(A ′ +θ) for some ψ ∈ (0, 1) for all k ≥ k min , then
where φ min ∈ R is a lower bound for φ(•, T ′ ) over X by Assumption 5 (namely, Assumption 1).
ii.
for some ψ ∈ (0, 1) for all k ≥ k min , then
In either case, if in a run, there exists K ⊆ N with |K | � ∞ and lim k∈K, k→∞ ∆l(x k , τ k , ∇f k , d true k ) � 0, then any limit point of {(x k , y k + δ true k )} is a first-order stationary point for (1). Proof. By the definition of A ′ , {β k } ⊂ (0, 1], and line 15 of Algorithm 1, it follows that 38)); Lemmas 10, 14, 8, 13, and 15; and the fact that {β k } ⊂ (0, 1] that for all k ≥ k min , one finds
where
Then, taking total expectation conditioned only on the event E, it follows for k ∈ N that
After rearrangement, one finds that (39) holds, where the limit as k → ∞ holds because of the aforementioned fact that E[φ(X k min , T ′ ) | E] is bounded under Assumption 5. Case ii. By the conditions on {β k } k≥k min , it follows in a similar manner as in case (i) that
which after rearrangement and taking limits as k → ∞, proves that (40) holds.
The final conclusion of the theorem follows by the arguments provided before the theorem. w
We close this section by remarking that as described in Berahas et al. (2021), the elements of {β k } can be chosen to satisfy Assumption 5 and the conditions of Theorem 1. Specifically, to obtain the convergence guarantees in case (i) of Theorem 1, the algorithm can set
for all k ∈ N, which clearly ensures that β k ∈ (0, 1] and that (23) holds. In addition, assuming that event E occurs, the value α ′ k for sufficiently large k becomes the realization of A ′ stated in the theorem, in which case β k � β for all subsequent k satisfies the condition stated in the theorem. To obtain the convergence guarantees in case (ii) of Theorem 1, the algorithm can "reset" a diminishing sequence after each iteration, in which the merit parameter and/or the ratio parameter are decreased. Specifically, in any iteration, say k ∈ N, in which the merit parameter and/or the ratio parameter were reduced from the prior iteration, one can set
. If the value for k is reset after every time the merit parameter and/or the ratio parameter are decreased, under event E, the value for k eventually will not be reset, meaning that in subsequent iterates, {β k } will satisfy the conditions of the theorem.
In this section, we demonstrate the performance of a MATLAB implementation of Algorithm 1 for solving (i) a subset of the constrained and unconstrained testing environment with safe threads (CUTEst) set (Gould et al. 2015) and (ii) two optimal control problems from Hinterm üller et al. (2003). The goal of our testing is to demonstrate the computational benefits of using inexact subproblem solutions obtained based on our termination tests from Section 2.2.
To obtain the normal direction v k as an inexact solution of (5), we applied CG to ∇c k ∇c T k v � �∇c k c k . Denoting the tth CG iterate as v k, t , where v k, 0 � 0, the method sets v k ← v k, t , where t is the first CG iteration such that
. The properties of CG as a Krylov subspace method ensure that v k, t ∈ Range(∇c k ) for all t ∈ N; hence, v k ∈ Range(∇c k ).
To obtain the tangential direction u k and associated dual search direction δ k , we applied the MINRES method (Paige andSaunders 1975, Choi et al. 2011) to (10). (We discuss our choice of H k along with each set of experiments.) Letting (u k, t , δ k, t ) denote the tth MINRES iterate, where (u k, 0 , δ k, 0 ) � (0, 0), the method sets (u k , δ k ) ← (u k, t , δ k, t ), where t is the first MINRES iteration such that for some κ ∈ (0, 1) (recalling the definition of (ρ k, t , r k, t ) in ( 28)),
and TT1 and/or TT2 hold. The choice of κ ∈ (0, 1) is discussed with each experiment.
Curtis, Robinson, and Zhou: Stochastic Inexact Sequential Quadratic Optimization INFORMS Journal on Optimization, Articles in Advance, pp. 1-23, © 2024 INFORMS 4.2. Choosing the Step Size Algorithm 1 (see line 15) stipulates that the step size α k chosen for the kth iteration satisfies α k ∈ [α min k , α max k ]. Keeping in mind that α min k ≤ α suff k ≤ min{α φ k , 1} (see Lemma 4), we take advantage of this flexibility in choosing the step size by defining
(1:1)
t k α min k otherwise, 8 > < > : where t k is the largest value of t ∈ N such that (1:1) t α min k ≤ min{α φ k , α min k + θβ 2 k , 1}. We do not explicitly compute α φ k in our code. Instead, we can verify whether (1:1) t α min k ≤ α φ k (as needed) because it is equivalent to verifying whether φ((1:1) t α min k ) ≤ 0, which is computable.
To test the utility of using inexact subproblem solutions in Algorithm 1, we consider two algorithm variants: SISQO and SISQO_exact. SISQO is Algorithm 1 with inexact solutions computed as described in Section 4.1 with a relatively large value for κ in (42). On the other hand, SISQO_exact is identical to SISQO with the exception that it uses a relatively small value for κ in ( 42). (Because of the similarities of the algorithms, SISQO_exact acts as a proxy for the stochastic sequential quadratic programming (SQP) algorithm from Berahas et al. (2021), although because it employs iterative linear algebra techniques, we are able to compare SISQO_exact with SISQO more readily.) We specify the values of κ ∈ (0, 1) used along with each of our tests in Sections 4.5 and 4.6. Our reason for comparing these two variants is to focus attention on the numerical gains obtained as a result of using inexact subproblem solutions. Both variants use the same computation for the normal step, so the performance difference can be attributed directly to the inexact tangential step computation. Additionally, we compare SISQO with a stochastic subgradient method employed to minimize the merit function φ directly (for various fixed values of τ). We refer to our implementation of this algorithm as Subgrad. Because H k is a diagonal matrix for all k ∈ N in all of our experiments, one CG or MINRES iteration is comparable computationally with two iterations of Subgrad.
Our metrics of interest are infeasibility and stationarity. Given any iterate x k in a run of SISQO, we consider the termination conditions ‖c(x k )‖ ∞ ≤ 10 �6 and ‖∇f k + ∇c k y k, ls ‖ ∞ ≤ 10 �2 , where y k, ls is the least-squares multiplier at x k . If an iterate satisfying these conditions is found in the first 1,000 iterations, then SISQO terminates and returns x SISQO ← x k . Otherwise, SISQO terminates after the 1,000th iteration and sets k ′ ← arg min i∈{0}∪ [1,000]
This allows us to associate with each run of SISQO the two measures err feasibility � ‖c(x SISQO )‖ ∞ and err stationarity � ‖∇f (x SISQO ) + J(x SISQO ) T y SISQO ‖ ∞ , where y SISQO ∈ R m is the least-squares multiplier at x SISQO . We use the total number of CG and MINRES iterations performed by SISQO as a budget for the total number of CG and MINRES iterations performed by SISQO_exact; no other termination condition is used for SISQO_exact. Upon termination of SISQO_exact, we define x exact -the iterate with which we define the feasibility and stationarity errors-using the same strategy as for setting x SISQO . Finally, we ran Subgrad for multiple instances of τ. (Further details on the choices of τ and the iteration budget for Subgrad are given with each experiment.) Upon termination of Subgrad, we define x subgrad -the iterate with which we define the feasibility and stationarity errors-using the same strategy as for setting x SISQO . In all cases, we define the KKT error as the maximum of the feasibility and stationarity errors.
In the CUTEst set (Gould et al. 2015), there are 138 equality constrained problems with m ≤ n. From these, we selected those such that (i) (n + m) ∈ [500, 10, 000], (ii) the objective function is not constant, (iii) the objective function remained above �10 50 over the sequences of iterates generated by runs of our algorithm, and (iv) the linear independence constraint qualification was satisfied at all iterates encountered in each run of our algorithm. This process of elimination resulted in the following 11 test problems: ELEC, LCH, LUKVLE1, LUKVLE3, LUKVLE4, LUKVLE6, LUKVLE7, LUKVLE9, LUKVLE10, LUKVLE13, and ORTHREGC.
The function and derivative evaluations from CUTEst are deterministic, and for the purpose of these experiments, we exploited this fact to compute values as needed by our algorithm, including using function evaluations to estimate Lipschitz constants. However, we introduced noise into the computation of the objective gradients for each application of our stochastic algorithm. In particular, we generated stochastic gradients as
� , where for testing purposes, we considered the three noise levels ɛ N ∈ {10 �4 , 10 �2 , 10 �1 }. This particular choice for defining the stochastic gradients ensured that an appropriate value for M g as indicated in Assumption 6 would be given by M g � {10 �8 , 10 �4 , 10 �2 }, corresponding to the values for ɛ N .
We set κ � 0:1 for SISQO and κ � 10 �7 for SISQO_exact. All of the remaining parameters were set identically: τ �1 � σ � κ v � 0:1, η � 0:5, ɛ uv � 1, 000, ξ �1 � ɛ c � 1, ɛ τ � ɛ ξ � 0:01, κ ρ � κ r � 100, ɛ 2 � 0:9, κ u � 10 �8 , χ � 1 � 10 �8 , θ � 10 4 , and β k � 1 for all k ∈ N. For all k ∈ N, we randomly generated a sample point near x k , and then, we estimated L k and Γ k using finite differences of the objective gradients and constraint Jacobians between x k and the sampled point. These values were used in place of L and Γ, respectively, in our step size selection. Here, H k � I for all k ∈ N in all runs.
For each test problem, we ran SISQO, SISQO_exact, and Subgrad with five different random seeds. As previously justified, the iteration budget for Subgrad was set to be twice the total numbers of CG and MINRES iterations used by SISQO. Also, for Subgrad, we ran the algorithm with the 11 merit parameter values in τ ∈ {10 0 , 10 �1 , : : : , 10 �10 } with step sizes set as α k � τ τL k +Γ k for all k ∈ N, and then, we selected the best iterate over all of these runs. We computed the feasibility and KKT errors for all algorithms as described in Section 4.4; see Figure 1.
From Figure 1, one finds that SISQO performs better than SISQO_exact and Subgrad in terms of both feasibility and KKT errors. SISQO achieves smaller errors for smaller noise levels, which may be expected because of the fact that these experiments are run with constant {β k }.
In our second set of experiments, we considered two optimal control problems motivated by those in Hinterm üller et al. (2003). In particular, we modified the problems to have equality constraints only and finite sum objective functions. Specifically, given a domain Ξ ∈ R 2 , a constant N ∈ N >0 , reference functions w ij ∈ L 2 (Ξ) and z ij ∈ L 2 (Ξ) for (i, j) ∈ {1, : : : , N} × {1, : : : , N}, and a regularization parameter λ ∈ R >0 , we first considered the problem min w, z
� � s:t: � ∆w � z in Ξ and w � 0 on ∂Ξ: (43)
Second, with the same notation but z ij ∈ L 2 (∂Ξ), we also considered min w, z
∂Ξ) � � s:t: � ∆w + w � 0 in Ξ and ∂w ∂p � z on ∂Ξ, (44) where p represents the unit outer normal to Ξ along ∂Ξ. As reference functions for both problems, we chose z ij � 0 and
x 2 � for all (i, j) ∈ {1, : : : , N} × {1, : : : , N} for some Note. Subgrad, stochastic subgradient method. Curtis, Robinson, and Zhou: Stochastic Inexact Sequential Quadratic Optimization
(ɛ S , ɛ N ) ∈ R >0 × R >0 . We selected the following values for the constants: N � 3, λ � 10 �5 , ɛ S � 50, and ɛ N ∈ {10 �4 , 10 �2 , 10 �1 }. Because the objective functions of ( 43) and ( 44) are finite sums, to generate stochastic gradients as unbiased estimates of the true gradient, we first uniformly generated random (i, j) ∈ {1, : : : , N} × {1, : : : , N}, and then, we computed the gradient corresponding to the (i, j)th term in the objective function. We note that with the choice of parameters, it follows that an appropriate value for M g in Assumption 6 is given by M g ≈ {10 �8 , 10 �4 , 10 �2 } to correspond, respectively, to the values for ɛ N .
Because the optimal control problems have a quadratic objective function and linear constraints, we used the exact second derivative matrix H k � diag(I, λI) for all k ∈ N. For this choice, the curvature condition on H k in Assumption 2 is trivially satisfied.
In terms of algorithm parameters, we set κ � 10 �4 for SISQO and κ � 10 �7 for SISQO_exact. All of the remaining parameters were set identically for the two variants in the same manner as in the previous section with the following exceptions: τ �1 � 10 �4 , L � 1, and Γ � 0, where the latter choice is valid because the objectives are quadratic and the constraints are linear.
For each of the two optimal control problems in ( 43) and (44), we ran SISQO, SISQO_exact, and Subgrad with five different random seeds, and then, we computed their average feasibility and KKT errors as described in Section 4.4. We observed that {τ k } was constant in all runs of SISQO and SISQO_exact. Therefore, we ran Subgrad with only three merit parameter values, namely τ ∈ {10 �2 , 10 �4 , 10 �6 }, and we choose step sizes as α k ← τ τL+Γ for all k ∈ N. (In these experiments, the budget for Subgrad iterations was set to the total numbers of CG and MINRES iterations used by SISQO because the constraint function evaluations, required in each iteration of Subgrad, are as expensive computationally as each CG and MINRES iteration.) In Table 1, we report average feasibility and KKT errors as well as the average number of iterations performed by Algorithm 1 before termination ("iter.") and the number of CG and MINRES iterations ("C + M iter."), with the latter discussed in Section 4.1. The results are given in Table 1. One can observe that SISQO performs better than the others in terms of average feasibility and KKT errors.
We have proposed, analyzed, and tested an inexact stochastic SQP algorithm for solving stochastic optimization problems involving deterministic, smooth, nonlinear equality constraints. We proved a convergence guarantee (in expectation) for our algorithm that is comparable with that proved for the exact stochastic SQP method recently presented by Berahas et al. (2021), which in turn, is comparable with that known for the stochastic gradient in unconstrained settings (Bottou et al. 2018). Our MATLAB implementation, SISQO, illustrated the benefits of allowing inexact step computation for solving problems from the CUTEst set (Gould et al. 2015) and two optimal control problems.
Downloaded from informs.org by [98.7.209.7] on 24 May 2024, at 09:40 . For personal use only, all rights reserved.
Curtis, Robinson, and Zhou: Stochastic Inexact Sequential Quadratic Optimization
INFORMS Journal on Optimization, Articles in Advance, pp. 1-23, © 2024 INFORMS Downloaded from informs.org by [98.7.209.7] on 24 May 2024, at 09:40 . For personal use only, all rights reserved.
Curtis, Robinson, and Zhou: Stochastic Inexact Sequential Quadratic Optimization INFORMS Journal on Optimization, Articles in Advance, pp. 1-23, © 2024 INFORMS