{'fraction_non_alphanumeric': 0.04188188846849811, 'fraction_numerical': 0.011449251521631848, 'mean_word_length': 4.5940922057605595, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An efficient, globally convergent method for optimization under uncertainty using adaptive model reduction and sparse grids

This work introduces a new method to efficiently solve optimization problems constrained by partial differential equations (PDEs) with uncertain coefficients. The method leverages two sources of inexactness that trade accuracy for speed: (1) stochastic collocation based on dimension-adaptive sparse grids (SGs), which approximates the stochastic objective function with a limited number of quadrature nodes, and (2) projection-based reduced-order models (ROMs), which generate efficient approximations to PDE solutions. These two sources of inexactness lead to inexact objective function and gradient evaluations, which are managed by a trust-region method that guarantees global convergence by adaptively refining the sparse grid and reduced-order model until a proposed error indicator drops below a tolerance specified by trust-region convergence theory. A key feature of the proposed method is that the error indicator---which accounts for errors incurred by both the sparse grid and reduced-order model---must be only an asymptotic error bound, i.e., a bound that holds up to an arbitrary constant that need not be computed. This enables the method to be applicable to a wide range of problems, including those where sharp, computable error bounds are not available; this distinguishes the proposed method from previous works. Numerical experiments performed on a model problem from optimal flow control under uncertainty verify global convergence of the method and demonstrate the method's ability to outperform previously proposed alternatives.

rules with natural dimension-adaptive refinement and error estimation to approximate the integrals that arise in the computation of risk measures [3]; reduced-order models cheaply approximate the solution to the parametrized nonlinear system at the quadrature nodes. The method solves the resulting optimization problem in a globally convergent manner using the trust-region method proposed in Ref. [17], which allows for inexact objective and gradient evaluation, even at trust-region centers. Global first-order convergence of this trust-region method relies on (1) inexpensive, computable error indicators for the objective and gradient approximations and (2) models that can be refined until these error indicators drop below required tolerances. The only requirement on the error indicators is that they bound the corresponding error up to an arbitrary constant that need not be computed or estimated; this provides substantial flexibility to develop a posteriori residual-based error indicators that are applicable to a wide range of problems. To this end, we derive error indicators for the objective function and its gradient. The gradient error indicator comprises a combination of both primal and adjoint residual-based error indicators for the ROM approximation [29,28] and an approximation of the quadrature truncation error on the forward neighbors of the sparse grid [13].
With the error indicators defined, we introduce two algorithms to construct an anisotropic sparse grid and reduced basis, which together form the trust region approximation models. The first algorithm constructs a sparse grid and reduced basis such that an accuracy condition on the gradient error indicator [16] is satisfied at the trust-region center. The second algorithm constructs a (possibly different) sparse grid and reduced basis such that an accuracy condition on the objective error indicator [17] is satisfied. Both algorithms combine the dimension-adaptive approach proposed in Ref. [13] to construct an anisotropic sparse grid with a variant of the greedy method proposed in Refs. [22,24] to construct a reduced basis into a single nested algorithm that constructs both simultaneously. Here, the outer loop refines the sparse grid by adding the index set from the set of forward neighbors that contributes most to the quadrature truncation error for the ROM gradient. A direct search for the index set that maximizes this error is inexpensive because it requires only ROM evaluations. Once the sparse grid is refined, the inner loop greedily constructs the reduced basis by sampling the HDM at sparse grid collocation points where the primal and adjoint residual-based error indicators are maximized. To ensure the reduced basis can accurately approximate the primal and adjoint solutions, it is constructed from both primal and dual snapshots. A crucial ingredient in this setting is the use of minimum-residual primal and adjoint reduced-order models, as this approach (1) ensures that the reduced-order-model solution minimizes the objective and gradient error indicators over the associated trial subspaces, (2) guarantees the greedy method based on the primal and adjoint residual-based indicators terminates with a reduced basis that satisfies the trust-region global convergence conditions, and (3) ensures the presence of adjoint solutions in the primal basis does not degrade the primal approximation and vice versa, which is not guaranteed in general [2,12].
Chen et al. have related work [6,7,8,9] that also integrates reduced-order models and sparse grids for uncertainty quantification and stochastic optimal control. Unlike the method proposed in this work, however, their approach adheres to an offline-online decomposition that constructs a sparse grid and reduced-order model in a computationally expensive offline stage and subsequently deploys them in a computationally inexpensive online stage. Their methods are equipped with rigorous a posteriori error bounds and convergence guarantees, provided the governing equations are linear elliptic PDEs. By breaking the offline-online decomposition and globalizing the sparse-grid-reduced-order-model (SG-ROM) approximation model with a practical trust-region method [17], we can establish global convergence for a significantly wider range of governing equations, including nonlinear problems.
There are a number of other works that have combined reduced-order, or more generally, surrogate models, to accelerate optimization under uncertainty [19,27,23,31] and others that have used reducedorder models to estimate risk measures, including the conditional value-at-risk [14,30]. Ref. [19] used reduced-order models to efficiently approximate the solution of stochastic linear elasticity problems and parametrized the stochastic space using polynomial chaos. While their method does adaptively re-sample the parameter space throughout the optimization to promote convergence, the adaptation is heuristic and does not necessarily guarantee global convergence. Ref [31] presents a method that uses local reduced-order models and Monte Carlo sampling inside a globally convergent trust region method to minimize conditional value-at-risk. The reduced-order models are local to a Voronoi cell of the stochastic space, which limits the method to rather low-dimensional stochastic spaces, and built using samples of the PDE solution and its gradient with respect to the stochastic variables.
The remainder of the paper is organized as follows. Section 2 formulates the optimization under uncer-tainty problem governed by a large-scale system of nonlinear equations. Section 3 introduces the primary approximation techniques that form the foundation for our method: stochastic collocation via anisotropic sparse grid and minimum-residual, projection-based primal and adjoint reduced-order models. Section 4 presents the new adaptive algorithm for optimization under uncertainty that manages the inexactness introduced from the sparse grid quadrature and reduced-order-model evaluations using the trust-region method introduced in Ref. [17]. Appendix A introduces and derives residual-based error bounds used to define the error indicators required by the trust-region convergence theory. Numerical results presented in Section 5 show the proposed method dramatically reduces the number of PDE queries required to solve two stochastic optimal flow control problems, as compared with the benchmark method proposed in Ref. [17].
2. Problem formulation. Let (Ω, F, P ) be a probability space: Ω denotes the set of outcomes, F ⊆ 2 Ω denotes a σ-algebra of events and P : F → [0, 1] denotes a probability measure. We consider the high-dimensional model (HDM) to be a large-scale system of nonlinear equations parametrized by both a random vector y ∈ Ξ-which results from the mapping Ω ω → y ∈ Ξ ⊆ R ny with Lebesgue density ρ : Ξ → [0, +∞)-and deterministic parameters µ ∈ R nµ . The corresponding problem statement is the following: Given a realization of the random vector y ∈ Ξ and an instance of the parameters µ ∈ R nµ , compute the (primal) solution u satisfying where r : (u, y, µ) → r(u, y, µ) with r : R nu × R ny × R nµ → R nu denotes the residual. We assume that for every (y, µ) pair, there exists a unique (primal) solution u = u (y, µ) satisfying Eq. (1); we further assume that the implicit map (y, µ) → u is continuously differentiable with respect to its arguments. We primarily consider the case where the residual r arises from the high-fidelity discretization of a parametrized, stochastic partial differential equation (SPDE). In most practical applications, the dimension n u is large, which causes Eq. (1) to be computationally expensive to solve. Let f : R nu ×R ny ×R nµ → R denote a scalar-valued quantity of interest (QoI) associated with the system (e.g., the integral of a quantity over a region of the domain in the SPDE setting) and let F : R ny × R nµ → R denote its restriction to the manifold of solutions to Eq. (1) such that (2) F : (y, µ) → f (u (y, µ), y, µ).
This work considers only the risk-neutral measure, i.e., the expectation risk measure R ≡ E; however, the approach can be extended to other smooth risk measures. Because the parameter-space dimension n µ is potentially large, we employ the adjoint method to compute the gradient of the quantity of interest. The problem statement associated with the HDM adjoint problem is: Given a realization of the random vector y ∈ Ξ, an instance of the parameters µ ∈ R nµ , and the primal solution u ∈ R nu , satisfying Eq. (1), compute the adjoint solution λ ∈ R nu satisfying where the adjoint residual is defined as From the primal-adjoint pair (u , λ ) satisfying Eqs. (1) and (5), the gradient of the quantity of interest can be computed as where operator that reconstructs the gradient from the adjoint solution is The gradient of the objective, i.e., ∇J, can be can be computed directly from the gradient ∇ µ F and the analytical form of the risk measure. In the risk-neutral case, this corresponds to the expectation of the gradient of F , i.e., 3. Two sources of inexactness: sparse grids and model reduction. In this section, we introduce the two approximation tools used in this work to accelerate the solution of the stochastic optimization problem (3): dimension-adaptive anisotropic sparse grids and projection-based reduced-order models. Sparse grids provide a relatively small number of quadrature nodes (and associated weights) to enable the efficient evaluation of the risk measure in a moderate-dimensional stochastic space Ξ, while reduced-order models rapidly approximate the quantity of interest at each quadrature node. Section 4 combines these approximation techniques with a trust-region method [17] to yield an efficient, globally convergent method for solving stochastic optimization problems.
3.1. Stochastic collocation based on sparse grids. Assuming the stochastic space Ξ corresponds to a tensor product of one-dimensional intervals, i.e., Ξ = Ξ 1 ⊗ · · · ⊗ Ξ ny , Ξ k ⊂ R, k = 1, . . . , n y , the sparsegrid construction begins with the definition of a one-dimensional quadrature rule of level i ∈ N that will be used in the kth dimension, which we denote by E i k . The level i indicates refinement of the one-dimensional quadrature rule such that where the probability density is a product of one-dimensional probability densities, ρ = ρ 1 · · · ρ ny . Let Ξ i k ⊂ [−1, 1] be the quadrature nodes associated with the quadrature rule E i k . From the one-dimensional quadrature rules, the corresponding difference operators are defined as The requirement in (10) on the quadrature rules implies ∆ i k [h] → 0 as i → ∞. A multi-dimensional difference operator is constructed from a tensor product of one-dimensional difference operators, each possibly at a different level of refinement, i.e., (12) ∆ i := ∆ i1 1 ⊗ · · · ⊗ ∆ in y ny .
A multi-index i ≡ (i 1 , . . . , i ny ) ∈ N ny is used to track the refinement level of each one-dimensional difference operator, i.e., i k is the refinement level of the difference operator in dimension k. From the multi-dimensional difference operator, a quadrature rule E I is defined by summing over all multi-indices in a prescribed multiindex set I ⊂ N ny , i.e., ny .
An index set I must satisfy the standard admissibility requirement [13] for the above expression to be a convergent quadrature approximation. This requirement is: k − e j ∈ I for all k ∈ I, j ∈ { | 1 ≤ ≤ n y , k > 1}, where e j denotes the jth canonical unit vector. The collection of multi-indices that neighbor a multi-index set I is defined as where I c denotes the complement of the multi-index set I in N ny , i.e., I c := {i ∈ N ny | i ∈ I}. Following Refs. [13,16,17], the truncation error, which can be written as an infinite sum, can be approximated as Finally, let Ξ I ⊂ Ξ denote the quadrature nodes associated with the multi-dimensional quadrature rule E I . If nested one-dimensional quadrature rules are used, then Ξ I ⊂ Ξ J for I, J multi-index sets such that I ⊂ J . This can lead to substantial savings since evaluations of h can be recycled as the sparse grid is refined.

Projection-based model reduction.
In the construction of projection-based reduced-order models, the distinction between the stochastic variables, y, and parameters, µ, is unimportant so, for notational brevity, we combine them into a single vector, (16) z where n z := n y + n µ .
3.2.1. Minimum-residual primal reduced-order model. The projection-based reduced-order model construction begins with the ansatz that the state vector lies in a low-dimensional subspace (17) u ≈ Φq, where Φ ∈ R nu×ku with k u n u denotes reduced basis matrix that spans the low-dimensional subspace, q ∈ R ku denotes the generalized coordinates of the state u in the reduced subspace, and R m×n denotes the the set of full-column-rank m × n real-valued matrices (i.e., the noncompact Stiefel manifold).
The governing equations for the reduced-order model are obtained by substituting the ansatz (17) into the HDM governing equations (1) and enforcing orthogonality of the residual to a test basis matrix Ψ : R nu × R nz → R nu×ku , yielding the following problem statement: Given z ∈ R nz , compute the reduced primal solution q satisfying We assume that for each z ∈ R nz there exists a unique reduced solution q = q (z) satisfying Eq. (18) and the implicit map z → q is continuously differentiable with respect to its arguments.
The reduced-order model in (18) is said to possess the minimum-residual property [18,5,28] if it is equivalent to minimizing the HDM residual r in some metric over the low-dimensional trial space Ran(Φ), with Ran(A) denoting the range of matrix A. given z ∈ R nz for some symmetric positive-definite matrix Θ ∈ R nu×nu .
The first-order optimality conditions of the unconstrained optimization problem (19) provide a necessary condition for a reduced-order model to possess the minimum-residual property, namely as this choice yields equivalence between the ROM governing equations (18) and the first-order optimality conditions of problem (19). Minimum-residual reduced-order models have the benefit of minimizing the error bound on the QoI in (25) and possess three key properties defined in Proposition 1 [28]: (1) optimality, (2) monotonicity, and (3) interpolation.
• (Monotonicity) Let q ∈ R ku be the solution for a given z corresponding to any projection-based reduced-order model with reduced basis Φ ∈ R nu×ku such that Ran(Φ ) ⊆ Ran(Φ), then • (Interpolation) Let u satisfy r(u , z) = 0. If u ∈ Ran(Φ), then where a superscript + denotes the Moore-Penrose pseudo-inverse.
Proof. Optimality follows directly from the assumption that the q is the global minimum of (19): since q is the global minimum of the residual over the trial space, any other point in the trial space must have a larger residual. Monotonicity and interpolation follow directly from optimality.
The quantity of interest associated with any reduced-order model characterized by a unique mapping z → q can be computed from the mapping which can be equipped with a residual-based error indicator if the quantity of interest and residual are sufficiently regular, as defined in Appendix A, i.e., where κ > 0 is independent of z, Θ ∈ R nu×nu is a symmetric positive-definite (SPD) matrix, and x 2 Θ := x T Θx for any x ∈ R nu ; see Appendix A.
The problem statement associate with the adjoint of the reduced-order model in (18) for a general test basis Ψ = Ψ(u, z) is the following: Given z ∈ R nz and the reduced primal solution q satisfying Eq. (18), compute the reduced adjoint solution η ∈ R ku satisfying where all quantities in Eq. (26) are evaluated at (Φq , z). From the reduced primal-adjoint pair (q , η ) satisfying Eqs. (18) and (26), the gradient of the reduced quantity of interest can be computed as where again all quantities are evaluated at (Φq , z). If the reduced primal solution q is exact such that r(Φq (z), z) = 0 or the test basis Ψ is independent of the primal state vector u, the reduced adjoint equations (26) simplify to and the gradient of the reduced quantity of interest simplifies to The reduced adjoint equations (26) and the corresponding expression to construct the gradient of the reduced quantity of interest (27) can be challenging to implement, as computing the derivative of the test basis entails computing the second derivative of the residual. Furthermore, the reduced adjoint variable η satisfying Eq. (26) cannot necessarily be reconstructed in the full space such that it provides an accurate approximation of the full-system adjoint variable λ satisfying Eq. (5), nor does the reduced gradient ∇ µ F r accurately approximate its full-system counterpart ∇ µ F . Therefore, instead of adhering to this reducethen-optimize [20] approach, we instead construct a minimum-residual reduced-order model for the fullsystem adjoint equations (5). While this approach forfeits consistency between the reduced quantity of interest F r and its derivative, it gains the notion of residual minimization and its properties (i.e., optimality, monotonicity, interpolation).

3.2.2.
Minimum-residual adjoint reduced-order model. The minimum-residual adjoint construction begins with the ansatz that the adjoint solution can be well approximated in a low-dimensional subspace. For simplicity, we employ the same trial subspace to approximation the primal and adjoint solutions (30) λ ≈ Φη, but this is not required. Reference [28] takes the adjoint trial subspace to be equal to the primal test subspace and shows this to be sufficient for the minimum-residual adjoint and true reduced-order model adjoint to agree whenever the reduced primal solution is exact. The· notation is used to distinguish the reduced-order adjoint model from the true adjoint of the reduced system. We propose to compute the reduced adjoint variable by minimizing the Θ λ -norm of the adjoint residual over the trial subspace such that the associated problem statement becomes: Given z ∈ R nz and the reduced primal solution q satisfying Eq. (18), compute the minimum-residual adjoint solutionη ∈ R ku satisfying The first-order optimality conditions of (31) can be written as the linear system where all quantities in Eq. (32) are evaluated at (Φq , z). The minimum-residual adjoint reduced-order model also possess similar concepts of optimality, monotonicity, and interpolation as the primal [28].
define a minimum-residual adjoint reduced-order model whose solution η ∈ R ku for a given z ∈ R nz , satisfies (31). Then the following properties hold: • (Optimality) For any λ ∈ Ran(Φ), • (Monotonicity) Let (Φ , Θ λ ) define a minimum-residual adjoint reduced-order model such that whereη satisfies • (Interpolation) If the adjoint solution λ satisfying Eq. (5) also satisfies λ ∈ Ran(Φ), then Proof. Optimality follows directly from problem (31): becauseη minimizes the adjoint residual over the trial space and the optimization problem is convex (i.e., it is a linear-least-squares problems characterized by a full-column-rank matrix), any other element of the trial space will yield a larger objective-function value. Monotonicity and interpolation follow directly from optimality.
Minimum-residual sensitivity reduced-order models and the corresponding approximation of the gradient of the QoI can be constructed similarly. We refer to Refs. [28,29] for details.
We propose to approximate the gradient of the QoI using the adjoint-based gradient operator g λ defined in (8) as where the minimum-residual primal q and adjointη solutions satisfy (19) and (31), respectively. We emphasize that ∇ µ F r = ∇ µ F r in general.
From Appendix A, any approximation of this form is equipped with the following residual-based error bound for some constants κ, τ > 0, any q,η ∈ R ku , and any symmetric-positive-definite matrices Θ, Θ λ ∈ R nu×nu . We justify the choice of employing the minimum-residual primal q and adjointη solutions in the gradient estimate (37) by noting that this choice minimizes the gradient-error bound (38); critically, this property is not generally shared with the exact gradient of the reduced quantity of interest, i.e., ∇ µ F r .

Approximation model and managed inexactness via trust-regions.
We have now introduced two approximation techniques equipped with associated error indicators: dimension-adaptive sparse grids, which enable efficient quadrature in moderate-dimensional stochastic spaces, and projection-based reducedorder models, which enable inexpensive PDE solutions. In this section, we combine these into a single approximation model for the risk measure (4) associated with a stochastic, parametrized nonlinear system (1). The level of the sparse grid and dimension of the reduced-order model serve as two refinement mechanisms that trade cost for improved approximation accuracy. With this refinement-equipped approximation model, we develop a trust-region-based approximation model management method [1] to efficiently solve the original stochastic optimization problem (3). Critically, we show that the resulting approach enables global convergence to the solution of the original problem. Section 4.1 reviews details of the trust-region method proposed in Ref. [17], which provides the foundation for the proposed method. This approach allows for approximation models with inexact objective and gradient evaluations, as well as flexible error indicators. The remainder of the section formulates our approximation model and error indicator in the context of the global convergence theory of Ref. [17] and introduces a refinement algorithm for the approximation model to guarantee global convergence of our method.

4.1.
Trust-region method with inexact objective and gradient evaluations. Let J : R nµ → R be a smooth functional satisfying assumptions stated below. At trust-region iteration k characterized by trust-region center µ k , a trust-region method constructs a smooth approximation model m k : R nµ → R such that m k ≈ J within the trust region {µ ∈ R nµ | µ − µ k ≤ ∆ k }, where ∆ k > 0 denotes the trustregion radius. A trust-region method computes a trial iterateμ k by approximately solving the trust-region subproblem The computed trial iterate must satisfy the fraction of Cauchy decrease condition [11].
4.1.1. The gradient condition. Following Refs. [15,17], the gradient of the model m k must sufficiently approximate that of the true objective function J at the trust-region center µ k , i.e., where ξ > 0 is independent of k. Suppose the model m k is equipped with the error bound where ϕ k : R nµ → R denotes an error indicator for the gradient. Then requirement (40) can be restated solely in terms of the error bound as with κ ϕ > 0 a chosen constant, which we refer to as the gradient condition.

The objective condition.
A trust-region method accepts the trial iterateμ k if it produces sufficient decrease in the objective function relative to the decrease predicted by the model [11]. Traditionally, the ratio of actual-to-predicted reduction is used to assess whether or not the trial iterate is successful; this is costly to compute in the present context because it requires evaluating the true objective function at the trial iterate J(μ k ), which in turn requires computing the exact expectation over stochastic space. However, Ref. [17] introduced a method for scrutinizing a trust-region step that does not require evaluating the true objective function yet nonetheless ensures global convergence. This technique introduces another approximation ψ k : R nµ → R for the true objective function J that must be equipped with the error bound where σ > 0 is independent of k and θ k : R nµ → R denotes an error indicator for the decrease in the objective function from µ k . To ensure global convergence [17], the following objective condition must hold for a fixed ω ∈ (0, 1): such that lim k→∞ r k = 0. Also, η 1 and η 2 with 0 < η 1 < η 2 < 1 are the thresholds used to determine if the trial step should be accepted, which occurs if and if the trust-region radius should be increased, which occurs if k ≥ η 2 . The complete algorithm, summarized in Algorithm 1, is globally convergent [17] provided the following assumptions hold: 1. J is twice continuously differentiable and bounded below.
4. There exist κ 1 , κ 2 > 0 such that, for all µ ∈ R nµ and for all k, The remainder of this section introduces approximation models m k and ψ k and their associated error bounds ϕ k and θ k , respectively, for the risk measure associated with the quantity of interest (4) based on sparse grids and reduced-order models. The section also develops an algorithm that ensures global convergence of the method by adaptively refining the sparse grid and reduced-order model such that the approximation models are sufficiently accurate at trust-region centers according to the gradient condition (42) and the objective condition (44).

4.2.
Approximation model based on sparse grids and reduced-order models. In this section, we combine the approximation techniques from Section 3 with the trust-region method described in the Section 4.1 to develop a globally convergent optimization method that leverages and manages inexactness to solve the risk-neutral optimization problem As indicated by problem (48), we return to the notation introduced in Section 2 that distinguishes between stochastic variables y and optimization parameters µ. At the kth trust-region iteration, define the approximation models m k and ψ k discussed in Section 4.1 as respectively, where (I k , Φ k ) and (I k , Φ k ) are sparse-grid/reduced-basis pairs and, for any y ∈ Ξ and µ ∈ R nµ , q k (y, µ) denotes the unique solution satisfying Ψ T k r(Φ k q k , y, µ) = 0 and q k (y, µ) denotes the unique solution satisfying Ψ T k r(Φ k q k , y, µ) = 0. We assume that the projection uniquely determines the test basis matrix Ψ k (resp. Ψ k ) from the trial basis matrix Φ k (resp. Φ k ), e.g., Galerkin projection leads to Ψ k = Φ k , minimum-residual leads to (20). Approximation models (49) comprise the sparse grid approximation of the expectation with high-dimensional-model evaluations replaced by reduced-order-model evaluations. The choice of the sparse-grid/reduced-basis pairs will be driven by accuracy requirements on Algorithm 1 Trust region method with inexact objective evaluations [17] 1: Initialization: Given Model update: Choose a model, m k (µ), and gradient error bound, ϕ k (µ), such that where ξ > 0 is an arbitrary constants 3: Step computation: Approximately solve the trust-region subproblem for a candidate stepμ k that satisfies : Actual-to-predicted reduction: Compute actual-to-predicted reduction ratio approximation where ψ k (µ) and θ k (µ) satisfy where η < min{η 1 , 1 − η 2 } and σ > 0 is an arbitrary constant 5: Step acceptance: the approximation models provided by the gradient condition (42) and the objective condition (44); we defer a complete discussion to Section 4.3. From Eq. (27), the gradient of the approximation model is For this expression to hold, the adjoint η k appearing in ∇ µ F r ( · , µ) (see Eq. (27)) must satisfy the true reduced adjoint equations (26); however, for the case of minimum-residual primal reduced-order models, the reduced adjoint solution is difficult or impractical to compute for the reasons discussed in Section 3.2.1. Therefore, as described in Section 3.2.2, we replace the reduced adjoint solution with the minimum-residual adjoint solution satisfying Eq. (31) and approximate the model gradient according to whereη k (y, µ) uniquely satisfies (52)η k = arg min η∈R ku 4.2.1. Error indicator for the gradient condition. To properly embed the proposed approximation model in the globally convergent, inexact trust-region framework of Section 4.1, the approximation models must be equipped with error bounds of the form (41) and (43). The error indicator for the model gradient takes the form where β 1 , β 3 , β 4 > 0 are constants chosen to balance the contribution of the three terms and the individual error terms account for the primal reduced-order model error (E 1,k ), the adjoint reduced-order model error (E 3,k ), and the quadrature truncation error (E 4,k ), i.e., Each term can be efficiently computed since they only require primal and adjoint reduced-order model solves over the collocation nodes of the current sparse grid, I k , and its forward neighbors, N (I k ).
To demonstrate that the error indicator ϕ k (µ) as prescribed in (53) provides an asymptotic error bound for the model gradient, consider the error between the objective and model gradient: Employing E = E I k + E I c k and the triangle inequality yields From Proposition A.2, under the regularity and boundedness assumptions in Assumptions 1-2, there exist constant κ , τ > 0 such that Finally, because the integrals above cannot be computed exactly for general integrands, we approximate them as E I c k ≈ E N (I k ) and E = E I k ∪I c k ≈ E I k ∪N (I k ) , which yields From this bound and the definition of ϕ k in (53), the required gradient error bound in (41) follows.

Error indicator for the objective condition.
Similarly, the error indicator for the decrease in the objective function is defined as where α 1 , α 2 > 0 are constants chosen to balance the contributions of the two terms and the individual error terms account for the error in the primal solution (E 1,k ) and quadrature truncation error (E 2,k ), i.e., Both terms can be efficiently computed since they only require primal reduced-order model solves over the collocation nodes of the current sparse grid, I k , and its forward neighbors, N (I k ).
To demonstrate that the error indicator θ k (µ) as prescribed in (59) provides an asymptotic error bound for the objective decrease, we use the triangle inequality to bound the error in the objective decrease from µ k to µ by the error in the objective at both points, i.e., The objective error for any µ ∈ R nµ takes the form Again, the simple relation E = E I k +E I c k and triangle inequality bound this by two terms: one that accounts for the error in the reduced-order model, and one that accounts for the quadrature truncation error, i.e., From Proposition A.1, under the regularity and boundedness assumptions in Assumptions 1-2, there exists a constant κ > 0 such that Finally, we approximate the integrals in the above expression as The required objective decrease error bound in (43) follows from the above bound and the definition of θ k in (59)-(60).

Adaptive construction of sparse grid and reduced basis.
With the approximation models based on sparse grid quadrature and reduced-order model evaluations defined in the previous section, along with the associated error indicators, this section defines algorithms to construct the sparse-grid/reducedbasis pairs (I k , Φ k ) and (I k , Φ k ) such that the gradient condition (42) and the objective condition (44), which are required for convergence, are satisfied.
The first algorithm constructs (I k , Φ k ), which define model m k (µ) according to definitions (49), in Step 2 of Algorithm 1 to ensure satisfaction of the gradient condition (42); the second algorithm constructs (I k , Φ k ), which define model ψ k (µ k ) according to definitions (49), in Step 4 of Algorithm 1 (i.e., after a trial stepμ k is available) to ensure satisfaction of the objective condition (44). Both algorithms combine the concepts of dimension-adaptive sparse grid construction [13] and greedy construction of a reduced basis [22,24], wherein (1) the outer loop refines the sparse grid based on the local truncation error in the neighboring index set, and (2) the inner loop greedily samples over the current sparse grid to ensure the reduced-order-model contributions to the error are sufficiently small.

4.3.1.
Refinement algorithm for satisfying the gradient condition. As described in Section 4.1.1, global convergence of the trust-region method in Section 4.1 is predicated on gradient condition (42). A sufficient condition for the gradient condition to hold is that each term of ϕ k (µ k ) as defined in Eq. (53) satisfies an appropriate fraction of the condition, i.e., The purpose of the positive weights β 1 , β 3 , β 4 , introduced in the previous section, is to balance the individual contributions of the error terms such that the uniform split above is justified. This decomposition reduces the task of satisfying gradient condition (42) to ensuring each inequality in (66) is satisfied. While the interplay between the error terms in (66) and refinement of the sparse grid and reduced-order model is highly coupled and fairly complex, the following observations suggest an effective training strategy: (1) for a fixed reducedorder model, E 4,k decreases (possibly non-monotonically) as the sparse grid is refined, and (2) for a fixed sparse grid, E 1,k and E 3,k decrease (possibly non-monotonically) as the reduced-order model is hierarchically refined. Therefore, for a fixed reduced-order model, the proposed method adapts the sparse grid using the anisotropic dimension-adaptive approach proposed in Ref. [13] to reduce the truncation error E 4,k and, given the updated sparse grid, the proposed method employs a variant of the classical greedy method [22,24] to adapt the reduced-order model to reduce the error terms E 1,k and E 3,k . The algorithm performs these steps iteratively until inequalities (66) are satisfied. Sparse-grid construction The proposed construction of I k mimics the dimension-adaptive algorithm introduced by Gerstner and Griebel [13] for constructing a goal-oriented, anisotropic sparse grid. This approach approximates the truncation error associated with sparse grid I k using only the neighbors N (I k ). If this truncation-error approximation is larger than a specified tolerance, the multi-index in the set of neighbors that contributes most to the error is added to the index set, i.e., I k ← I k ∪ {i * } with and the integrand is g : R ny → R. In the context of the proposed SG-ROM approximation, the dimensionadaptive algorithm is applied to the integrand that appears in the gradient truncation error E 4,k (µ), i.e., for a fixed reduced basis matrix Φ k and given µ ∈ R nµ , we set the integrand to g : y → g λ (Φ kηk ( y , µ), Φ k q k ( y , µ), y , µ) .
While convergence is not necessarily monotonic, this term approaches zero as I k → N ny .
Reduced-basis construction The construction of the reduced basis follows the well-studied greedy algorithm [22,24]. The original greedy algorithm aims to improve the parametric robustness of a reducedorder model by iteratively enriching the reduced basis with snapshots of the high-dimensional model at the point in parameter space where the reduced-order-model error is largest. Pragmatically, this is executed by evaluating the reduced-order model and an inexpensive error indicator at a (possibly large) set of candidate points in parameter space, and subsequently performing a direct search for the maximum error-indicator value over this candidate set. Ref. [6] developed a weighted variant of the greedy algorithm for stochastic problems with non-uniform probability distributions. This approach uses the probability density ρ(y) to weight the error indicator while training the reduced-order model over the stochastic space Ξ. This is sensible, as stochastic-space regions with larger probability density yield larger contributions to the expected reduced-order-model error. Ref. [6] also coupled the weighted greedy algorithm with sparse grids by using the quadrature nodes as the candidate set; this is motivated by the observation that the reduced-order model is queried only at these points in stochastic space.
This work adopts a similar weighted greedy algorithm to train the reduced-order model over the quadrature nodes and neighbors Ξ I∪N (I) assuming fixed parameters µ and sparse grid I. Because the gradient condition is required to hold only at the trust-region center, the proposed method performs training solely in stochastic space, with Ξ I∪N (I) comprising the candidate set, for µ = µ k fixed. Unlike traditional greedy methods, the proposed method constructs a reduced basis that accurately represents both the primal and adjoint states over the candidate set. This is required because the same reduced basis Φ is employed for the primal (18) and adjoint (31) reduced-order models, and thus the greedy algorithm is responsible for reducing both the primal E 1,k and adjoint E 3,k error terms that arise in the gradient error indicator ϕ k (µ) defined in (53). The approach achieves this by adding both primal and adjoint snapshots to the reduced basis at each greedy iteration.
From the form of the gradient error indicator ϕ k in (53), the primal error indicator is taken as the primal residual norm weighted by the density, i.e., the integrand in E 1,k evaluated at µ = µ k (68) ρ(y) r(Φ k q k (y, µ k ), y, µ k ) Θ , and the adjoint error indicator is taken as the adjoint residual norm weighted by the density, i.e., the integrand in E 3,k evaluated at µ = µ k Refinement algorithm With these error indicators, the dimension-adaptive algorithm to construct the sparse grid I k and reduced basis Φ k at the kth trust-region iteration takes form. The algorithm initializes these entities to their values from the previous trust-region iteration as I k ← I k−1 and Φ k ← Φ k−1 . Then, the algorithm expands the sparse grid using the quadrature nodes in the neighboring index set that contribute most significantly to the truncation error, i.e., (70) With the updated sparse grid, the greedy algorithm uses the primal error indicator (71) y * = arg max to reduce E 1,k error below its required threshold in (66). If a minimum-residual reduced-order model is employed, the algorithm is guaranteed to terminate due to the monotonicity property in Proposition 1. In the limiting case where snapshots have been added for each y ∈ Ξ I∪N (I) , the primal reduced-order model will be exact for each y ∈ Ξ I∪N (I) and thus E 1,k = 0. Furthermore, the minimum-residual property guarantees E 1,k will decrease monotonically from the monotonicity property in Proposition 1.
After E 1,k has been reduced below its required threshold, the algorithm employs another greedy method with the adjoint error indicator (72) to reduce E 3,k below its required threshold in (66). Even if a minimum-residual adjoint reduced-order model is used, E 3,k is not guaranteed to decrease monotonically because modification of the trial basis matrix Φ k alters the linearization point defining the adjoint residual. This can be seen from Proposition 2 where the three properties hold provided the primal state is fixed. However, in the limit where all quadrature nodes in Ξ I∪N (I) have been sampled, the primal reduced-order model and adjoint reduced-order models will be exact and the algorithm will terminate with E 3,k = 0.
Once the reduced basis is constructed such that E 1,k and E 3,k satisfy their required bounds in (66), the dimension-adaptive greedy algorithm terminates if the truncation error indicator E 4,k also satisfies its bound in (66). Otherwise, another iteration of the algorithm is performed, i.e., the algorithm expands the sparse grid according to (70) and enriches the reduced basis according to (71)-(72). Because the truncation error tends to zero as the sparse grid is refined, this iteration is guaranteed to terminate provided the regularity conditions on the state and adjoint described in Assumptions 1-2 of Appendix A are met. Assuming the algorithm terminates, the resulting sparse grid I k and reduced basis Φ k are guaranteed to satisfy the gradient conditions in (66), thereby ensuring the gradient condition (42) required for global convergence is satisfied.
Algorithm 2 summarizes the refinement algorithm that applies the nested iteration described above to satisfy the individual gradient conditions (66). At trust-region iteration k = 0, the algorithm initializes the sparse grid as the uniform level-one sparse grid I = {(1, . . . , 1)}, which corresponds to a single quadrature node Ξ I = {0} 1 . The algorithm constructs the reduced basis from the primal and adjoint snapshot at this single quadrature node computed at the first trust-region center, i.e., the snapshots are u (0, µ 0 ) and λ (0, µ 0 ). Then, the approach may refine the sparse grid and reduced basis to yield I 0 and Φ 0 , respectively, by applying Algorithm 2 with inputs I −1 = {(1, . . . , 1) For all subsequent trust-region iterations, the sparse grid and reduced basis are initialized from their values at the previous trust-region iteration. In addition to providing a natural way to initialize the dimensionadaptive greedy algorithm, this choice has the following benefit: it expands the sparse grid and enriches the reduced basis only if the choices I k = I k−1 and Φ k = Φ k−1 are not sufficient to guarantee convergence.

Algorithm 2 Refinement algorithm for satisfying the gradient condition
Refine index set: Add index set with largest contribution to truncation error: Evaluate primal error indicator: Greedily select y ∈ Ξ i * with largest error indicator:

6:
Reduced-order model construction: Update reduced basis with new snapshots: Evaluate dual error indicator: Greedily select y ∈ Ξ i * with largest error:

10:
Reduced-order model construction: Update reduced basis with new snapshots: end while 12: end while This completes the discussion of the refinement algorithm used to construct I k and Φ k such that the gradient condition (42) is satisfied. We now turn attention to constructing I k and Φ k such that the objective condition (44) holds; this enables assessment of the trust-region step without requiring queries to J(µ).

4.3.2.
Refinement algorithm for satisfying the objective condition. As described in Section 4.1.2, satisfying the objective condition (44) is also required to preserve global convergence of the trust-region method when the approximate objective function ψ k (µ) is used in place of true objective function J(µ) in the computation of the actual-to-prediction reduction k . Analogously to inequalities (66), a sufficient condition for the objective condition to hold is that each term of the error indicator θ k (µ k ) as defined in Eq. (59) satisfies an appropriate fraction of the condition, i.e., This decomposition reduces the monolithic task of satisfying the objective condition (44) to ensuring each inequality in (73) is satisfied. Analogously to the approach described in Section 4.3.1 to construct I k and Φ k , we propose to employ a weighted greedy algorithm to enforce inequalities (73). While proposed dimension-adaptive greedy algorithm to construct (I k , Φ k ) for the objective decrease condition will be very similar to that used to construct (I k , Φ k ), there will be two critical differences. First, the error terms in (73) involve two points in parameters space: the trust-region center µ k and the candidate stepμ k . In contrast, the error terms in the gradient-condition inequalities (66) impose requirements only at the trust-region center µ k . This has implications for both the dimension-adaptive sparse grid construction and greedy method. Second, comparing definitions (54) and (60) reveals that inequalities (73) impose requirements only on the primal reduced-order model accuracy (through E 1,k ) and truncation error (through E 2,k ), whereas the gradient-condition inequalities (66) also placed requirements on the adjoint accuracy (through E 3,k ). This implies only primal snapshots are required during the greedy construction of the reduced-order model. However, we choose to also include adjoint snapshots because the candidateμ k will become the trust-region center at the next iteration µ k+1 if the trust-region step is successful, and adjoint accuracy is required at trust-region centers to satisfy the gradient condition.
Sparse-grid construction For a given (I k , Φ k ), if the truncation error conditions in inequalities (73), i.e., the requirements on E 2,k , are not satisfied, the algorithm updates the sparse grid according to The integrand in each term is precisely the integrand of E 2,k at the two parameters instances of interest: the trust-region center µ k and the candidate stepμ k . Therefore this refinement process can be repeated iteratively until the conditions on E 2,k in (73) are satisfied.
Reduced-basis construction Following the combined dimension-adaptive greedy method introduced in Section 4.3.1 for the gradient condition, the algorithm interleaves sparse-grid refinement steps with greedy reduced-basis construction. For a fixed (I k , Φ k ), define y * ∈ Ξ I k ∪N (I k ) and µ * ∈ {µ k ,μ k } as the quantities that maximize the weighted residual-based error indicator, i.e., (y * , µ * ) = arg max The reduced basis Φ k is updated according to Φ k ← Φ k u (y * , µ * ) λ (y * , µ * ) ; an optional orthogonalization step is usually used to ensure the reduced basis is full rank and the resulting reduced-order model is well-conditioned. The argument of the maximization problem in (75) is precisely the integrand of E 1,k . Assuming a minimum-residual reduced-order model is used, the terms E 1,k (µ k ) and E 1,k (μ k ) will monotonically decrease with each greedy iteration; the greedy algorithm proceeds until the requirements on E 1,k in inequalities (73) are satisfied. Refinement algorithm The proposed algorithm alternates between sparse-grid refinement and reducedbasis enrichment exactly as in Algorithm 2: for a fixed sparse grid, the greedy method is applied to enrich the reduced basis, then the reduced-order model is fixed and the sparse grid is refined. The combined algorithm terminates when all conditions in (73) are satisfied.
Algorithm 3 summarizes the proposed algorithm. Similarly to Algorithm 2, this algorithm refines a given sparse-grid/reduced-basis pair and thus implicitly requires initialization of each quantity. At any iteration k, the pair (I k , Φ k ) generated by Algorithm 2 to satisfy the gradient-condition inequalities (73) at µ k is used to initialize Algorithm 3. If the gradient-condition inequalities turn out to be more restrictive than that in the objective-condition inequalities (44), the algorithm will not modify the sparse grid or reduced basis, i.e., I k = I k and Φ k = Φ k . In this case, the actual-to-predicted ratio is unity and acceptance of the step is guaranteed.
Algorithm 3 Refine reduced basis and sparse grid for objective condition Input: Refine index set: Add index set with largest contribution to truncation error Evaluate error indicator: Greedily select y ∈ Ξ I k ∪N (I k ) , µ ∈ {µ k ,μ k } with the largest error y * , µ * = arg max Reduced-order model construction: Update reduced basis with new snapshot end while 8: end while 5. Numerical results: optimal boundary control of the incompressible Navier-Stokes equations. In this section, we apply the proposed method based on sparse grids and reduced-order models to solve an optimal control problem governed by the steady-state incompressible Navier-Stokes equations; we also compare the method's computational efficiency to that of existing methods.
Let Ω x ⊂ R 2 be the channel with a backward facing step shown in Figure 1. The goal of the control problem is to minimize the expected vorticity in the region immediately downstream of the step, denoted by Ω * , by controlling the inflow velocity along the boundary Γ c subject to uncertainty both in the fluid viscosity ν and in the inlet velocity along boundary Γ i . The optimization problem takes the form where u(x, y, µ) is the solution of the incompressible Navier-Stokes equations The fluid viscosity and inflow are uncertain and parametrized with a two-dimensional stochastic space with a uniform probability distribution ρ(y)dy = 2 −2 dy. In all numerical experiments, the partial differential equation in (77) is discretized with 232 Q 2 − Q 1 Taylor-Hood elements, yielding a state space of dimension n u = 2034, after application of the essential boundary conditions. The velocity at all nodes along the parametrized inlet Γ c , are taken as optimization variables, yielding a parameter space of dimension n µ = 38, and the finite element basis functions restricted to Γ c define the boundary control g(x, µ). The control regularization parameter is set to α = 0.1.
The flow velocity at the parametrized inlet is initially set to zero (µ = 0), i.e., a viscous wall, which leads to a recirculation in the mean flow and standard-deviation offset defined as respectively (see Figure 2). The local solution of the stochastic optimal control problem shown in Figure 2 effectively eliminates the recirculation region in not only the mean flow, but also in the standard-deviation offsets.
We apply the proposed method to solve problem (76). We construct sparse grids using one-dimensional, nested Clenshaw-Curtis quadrature rules [10], which ensure HDM and ROM evaluations can be recycled between sparse-grid refinement levels because the quadrature nodes at level j are a subset of those at level j + 1. We employ minimum-residual primal and adjoint reduced-order models because the PDE operators are not symmetric positive definite. For simplicity and implementation efficiency, the primal and adjoint reduced-order models are constructed to minimize their respective residuals in the Euclidean norm, i.e., Θ = Θ λ = I. Because we use minimum-residual primal and adjoint reduced-order models, Algorithms 2-3 are guaranteed to return sparse grids and reduced bases that satisfy the gradient condition (42) and objective condition (44), and global convergence is thus ensured, provided the regularity conditions in Appendix A are satisfied. The reduced basis is initialized using the primal and adjoint solution and primal solution sensitivity at y = y 0 = 0 and µ = µ 0 , i.e., to enhance the parametric robustness of the reduced-order model at early trust-region iterations with few HDM samples. Subsequently, the basis is updated using only primal and adjoint solutions as prescribed by Algorithms 2 and 3, because including n µ = 38 sensitivity solutions at every HDM sample would result in a computationally inefficient method. The trust region algorithm parameters, e.g., those used in Algorithm 1, are taken as: η 1 = 0.1, η 2 = 0.75, η = 0.1, γ = 0.5, r k = 1/(k + 1), and ω = 0.9. An approximate solution of the trust region subproblem is computed using the Steihaug-Toint Conjugate Gradient method [26,25,11] to minimize a quadratic approximation of m k (µ) in (49) within the trust region. The constants used to weight the contributions of the individual terms in the gradient error indicator are taken as unity, β 1 = β 3 = β 4 = 1.0, and κ ϕ = 1.0. To equally weight the terms in the objective function error indicator and weaken the error condition on the objective function, i.e., to make step acceptance more likely, the constants are taken as α 1 = α 2 = 0.01. Table 1 reports the convergence history of the proposed method. The method converges to a firstorder critical point ( ∇J(µ k ) → 0); after only 8 trust-region iterations the first-order optimality condition has reduced nearly 4 orders of magnitude from the initial (sub-optimal) control. At early iterations, the approximation model m k (µ) and true objective J(µ) agree only to one digit at the candidate stepμ k ; this is expected because a coarse sparse grid and only a few additional primal and adjoint HDM samples are required ( Figure 3) to meet the gradient condition. However, despite the inexpensive approximation with limited accuracy, the resulting approximation model is sufficiently accurate to significantly reduce the true objective function during the trust-region iteration. As the algorithm approaches a local minimum, the the gradient condition (42) places more stringent requirements on the model error and, as a result, the approximation model consists of a finer sparse grid with more HDM samples used to construct the reduced basis ( Figure 3) and thus provides a better approximation of the objective function.
We assess the computational performance of the proposed method using two baseline methods for comparison. The first method uses a fixed, 5-level isotropic sparse grid to perform integration in stochastic space using only HDM evaluations and the optimization problem is solved with the BFGS algorithm [21]. The second method uses the trust-region dimension-adaptive sparse grid approach proposed in Ref. [17]. For brevity, we abbreviate these methods SG-ISO and SG-TR, respectively. We refer to our proposed trust-region approach with dimension-adaptive sparse grids and reduced-order models as SG-ROM-TR.
The SG-TR method requires a large number of primal and adjoint HDM evaluations as the trust-    Fig. 3. The sparse-grid index set ( top) and corresponding quadrature nodes ( bottom) used by the proposed method at trust-region iterations k = 0, 3, 5, 6 ( left-to-right). Legend: sparse grid I k ( ), sparse grid neighbors N (I k )( ). The points in stochastic space where the primal and adjoint HDM are sampled to construct the reduced basis Φ k are indicated with a ( ). region iterations progress, because the evaluation of the models and their error indicators rely solely on HDM evaluations, albeit on a highly adapted sparse grid. In contrast, our SG-ROM-TR method requires HDM evaluations only to construct the reduced basis; evaluations of the trust-region models and their error indicators rely only on reduced-order model evaluations. As a result, the proposed SG-ROM-TR method requires one to two orders of magnitude fewer queries to the HDM (Figure 4) at the cost of a large number of ROM queries for various reduced-basis dimensions, k u ∈ [39, 165] ( Figure 5).
Because the proposed and baseline methods have different costs (i.e., the proposed method requires only HDM and ROM evaluations, while the baseline methods employ only HDM evaluations), care must be taken when assessing their performance. The ultimate cost metric of interest is wall time; however, it is wellknown that hyperreduction [4] must be incorporated to reduce the complexity associated with evaluating the nonlinear terms. Integration of hyperreduction in the proposed framework is the subject of ongoing work so we turn to another error metric. To assess the speedups that can be realized by this method, the following simple cost model is introduced where C is the total cost associated with a particular method in the units of equivalent number of primal HDM queries, n hp is the number of primal HDM queries, n ha is the number of adjoint HDM queries,n h is the average number of nonlinear iterations required to solve the primal HDM, n rp is the number of primal ROM queries, n ra is the number of adjoint ROM queries,n r is the average number of nonlinear iterations required to solve the primal ROM, and τ is the ratio of the cost of a primal HDM query to a primal ROM query. This cost model assumes a primal HDM (ROM) solve isn h (n r ) times as expensive as an adjoint solve and a primal HDM solve is τ times as expensive as a primal ROM solve. Numerical experiments at the initial and optimal controls suggestn h =n r = 5 and a range of τ values are considered.
As shown in Figure 6, for costly-to-evaluate reduced-order models (τ = 1), the proposed SG-ROM-TR method demonstrates faster convergence, in terms of the cost C, than the brute-force approach SG-ISO and similar convergence to the state-of-the-art SG-TR method. This is expected because SG-TR and SG-ROM-TR are adaptive algorithms with refinement at each iteration tailored to the convergence requirements, while the SG-ISO method relies solely on HDM evaluations on a fine sparse grid for all iterations. For a modest ROM speedup (τ = 10), the computational cost required for the SG-ROM-TR method to reach a given value of the objective function is nearly an order of magnitude lower than that of the state-of-the-art SG-TR method. The hypothetical situation of a free reduced-order model (τ = ∞) shows the speedup attainable by our method in comparison to SG-TR is at most 500. The τ -scenarios are not realistic in the sense that they assume the speedup of the ROM is constant across all iterations, which cannot be the case since the size of the ROM varies by more than a factor of 4 between the smallest and largest ROM used ( Figure 5). However, the scenarios do provide useful information because they bound the true cost, i.e., if the smallest and largest ROMs considered have speedup factors of τ 0 and τ 1 , respectively, the true cost-objective curve will lie between the τ 0 and τ 1 curves. For example, in the modest scenario of τ 0 = 100 and τ 1 = 10, the true cost-objective curve will lie between the ( ) and ( ) curves in Figure 6. This will lead to a significant speedup of the proposed method over both the SG-ISO and SG-TR methods.
6. Conclusion. This work proposes and demonstrates the merits of an efficient method for solving optimization problems constrained by large-scale nonlinear systems of equations with uncertain parameters. The method approximates the objective function, which comprises an integral over the stochastic space where the integrand depends on the solution of the governing system of equations, by combining two approximation techniques: dimension-adaptive sparse grids and projection-based reduced-order models. The reduced-order models generate inexpensive approximate solutions of the governing system, and sparse grids produce efficient quadrature schemes to integrate the objective function and its gradient. The proposed approximation model is used to define the trust-region model problem in the context of a trust-region method that allows for inexact objective and gradient evaluations, provided the model objective and gradient are equipped with weak error estimators. We derive error estimators for our proposed sparse-grid/reduced-basis model and introduce dimension-adaptive, greedy algorithms to refine the sparse grid and reduced basis until the trust-region convergence requirements are satisfied. Minimum-residual primal and adjoint reduced-order models ensure these algorithms will terminate and return a sparse grid and reduced basis that ensure convergence. Unlike existing methods that combine reduced-order models and sparse grids for stochastic optimization in the context of linear elliptic PDEs [7], our approach applies to general nonlinear systems and general quantities of interest with sufficient regularity; further, the method is equipped guaranteed to be globally convergent to the solution of the true optimization problem under mild regularity assumptions. Numerical experiments on a model problem from optimal stochastic flow control demonstrated the ability of the proposed method to significantly outperform alternative approaches in terms convergence as a function of computational cost. Promising future research directions include the extension of the proposed method to non-smooth risk