{'fraction_non_alphanumeric': 0.03863909527972028, 'fraction_numerical': 0.013234812062937062, 'mean_word_length': 4.193431692438644, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A globally convergent method to accelerate topology optimization using on-the-fly model reduction

We present a globally convergent method to accelerate density-based topology optimization using projection-based reduced-order models (ROMs) and trust-region methods. To accelerate topology optimization, we replace the large-scale finite element simulation, which dominates the computational cost, with ROMs that reduce the cost of objective function and sensitivity evaluations by orders of magnitude. To guarantee convergence, we first introduce a trust-region method that employs generalized trust-region constraints and prove it is globally convergent. We then devise a class of globally convergent ROM-accelerated topology optimization methods informed by two theories: the aforementioned trust-region theory, which identifies the ROM accuracy conditions required to guarantee the method converges to a critical point of the original topology optimization problem; a posteriori error estimation theory for project-based ROMs, which informs ROM construction procedure to meet the accuracy conditions. This leads to trust-region methods that construct and update the ROM on-the-fly during optimization; the methods are guaranteed to converge to a critical point of the original, unreduced topology optimization problem, regardless of starting point. Numerical experiments on three different structural topology optimization problems demonstrate the proposed reduced topology optimization methods accelerate convergence to the optimal design by a factor of at least two.


Introduction
Topology optimization enables the discovery of new, intricate, and non-intuitive structural designs, which can now be manufactured due to recent advances in additive manufacturing techniques [44]. However, topology optimization remains computationally expensive; for instance, density-based methods require billions of computational voxels to accurately represent complex geometries, and optimization requires hundreds of thousands of computing hours [1]. The computational cost increases further for robust optimization that accounts for uncertainties in loading conditions or material properties, often rendering it intractable in practical engineering settings.
A typical topology optimization approach uses first-order optimization methods and requires hundreds of design iterations for convergence. The dominant cost in each design iteration is the solution of the finite element problem associated with the (linear) elasticity equations. Hence, one class of approaches to improve the efficiency of topology optimization targets the (iterative) solver for the linear(ized) finite element system. Amir et al. [5] and Choi et al. [12] use inexact linear solves to efficiently compute approximate finite element solutions, while others have developed improved preconditioners, e.g., based on substructuring [14], multigrid [4], and Krylov methods with subspace recycling [37]. Alternatively, Rojas-Labanda et al. [28] investigate the use of second-order optimization methods to reduce the number of design iterations. Nguyen et al. [24] introduce a multi-resolution approach to topology optimization that uses two different meshes for finite element analysis and geometry representation to reduce the analysis cost.
Another class of approaches [40,15,23,12,39] to reduce the computational cost of topology optimization is based on the concept of reanalysis [19,45,35,31], i.e., they recognize that the optimization problem is a "many-query" problem that requires the analysis of many closely related designs along the optimization path and use previous finite element solutions to reduce the cost of subsequent solutions. Most of these methods are based on projection-based reduced-order models (ROMs), which approximate the solution to the finite element problem using a low-dimensional reduced basis. Yoon [40] uses ROMs in the frequency domain, with eigenmodes as a reduced basis, to reduce the cost of optimizing the frequency response to structures. Gogu [15] constructs ROMs on-the-fly, i.e., during the topology optimization procedure, using previously computed finite element solutions; the reduced basis is adapted when the finite element residual exceeds a predefined tolerance. Choi et al. [12] use a similar approach that adapts the residual tolerance based on the first-order optimality criteria to ensure more stringent accuracy requirements on the ROM are used near convergence; they also investigate its efficient implementation for large-scale problems. More recently, Xiao et al. [39] introduced a similar method to [15] that constructs a reduced basis on-the-fly using proper orthogonal decomposition (POD) [34] rather than Gram-Schmidt orthogonalization. A commonality of these approaches is they have been shown to effectively reduce the overall cost of topology optimization; however, to our knowledge, they make no attempt to rigorously assess the impact of ROM approximation on the convergence or to design algorithms that guarantee convergence to an optimal solution of the original (unreduced) problem. This work addresses the need for rigorous convergence theories and algorithms to accelerate topology optimization using ROMs constructed on-the-fly.
The first contribution of this work is the development of a model management framework [3] and global convergence theory for optimization problems with convex constraints based on generalized, error-aware trust regions [42]. Similar to (traditional) trust-region methods, error-aware trust-region methods (approximately) solve an inexpensive optimization problem associated with an approximate model within a trust region to advance the trust-region center toward the solution of the original optimization problem. The key difference is that error-aware trust-region methods consider a trust region based on the sublevel sets of an error indicator for the approximation model, whereas traditional trust regions are based on a notion of distance in the design space [13]. In this work, we generalize the error-aware trust-region method in [42], developed for unconstrained optimization, to problems with inexpensive convex constraints, which commonly arise in topology optimization, e.g., box constraints on the design variables and volume fraction bounds. We prove global convergence of the method assuming (i) the value and gradient of the approximation model at any iteration match the optimization objective function at the corresponding trust-region center and (ii) the trust-region constraint is the sublevel set of an asymptotic error bound for the approximation model. While assumptions (i) on the accuracy of the approximation model at the trust-region center can be weakened considerably by appealing to more general trust-region theory [11,17,20,16,21], this turns out to be unnecessary; the assumptions are easy to satisfy if projection-based ROMs are used as the approximation model. This error-aware trust-region method provides a general framework to develop globally convergent solvers for optimization problems constrained by partial differential equations using ROMs constructed onthe-fly, which leads to the second contribution of this work.
The second key contribution of this work is the development of an efficient, globally convergent topology optimization method that uses ROMs constructed on-the-fly as the approximation model in the proposed error-aware trust-region framework. At a given iteration, the ROM is updated on-the-fly based on the solution and adjoint snapshots at all trust-region centers encountered and a POD compression procedure that guarantees the ROM satisfies the aforementioned accuracy requirements (i) and (ii) for the trust-region approximation model. We will consider two trust-region constraints: the traditional trust region based on distance from the center and the sublevel sets of the residual-based error indicator of the ROM. The latter choice follows the work in [43,42] for general PDE-constrained optimization problems as well as the work in [15,12,39] that use the residual as an indicator to update the ROM. We appeal to a posteriori error estimation theory for projection-based ROMs (see, e.g., review paper [29] and textbook [18]) to prove that the procedure satisfies the criteria (i) and (ii) for global convergence in the error-aware trust-region framework; therefore, the method is guaranteed to converge to a local minimum of the original optimization problem from an arbitrary starting point. We use three benchmark topology optimization problems to demonstrate global convergence of the method and demonstrate speedups of at least a factor of two relative to a standard topology optimization approach that applies the method of moving asymptotes [36] to the original (unreduced) topology optimization problem. While this is not the first work to embed on-the-fly ROMs in a trust-region setting, e.g., [8,9,41,2,30,43,42,27], to our knowledge it is the first to do so in the context of topology optimization and, other than [42], use the notion of an error-aware trust region.
The remainder of the paper is organized as follows. Section 2 introduces a density-based topology optimization formulation with Helmholtz filtering and the adjoint method to compute gradients with respect to the design variables. Section 3 reviews projection-based reduced-order models for the linear elasticity finite element system and provides the associated error analysis that is required to establish global convergence of the proposed ROM-based topology optimization method. Section 4 introduces the error-aware trust-region method that utilizes generalized trust-region constraints and establishes its global convergence theory. Then we use projection-based reduced-order models and the error analysis from Section 3, and introduce a procedure to construct a reduced basis from from state and adjoint snapshots, to define a class of globally convergent ROM-based topology optimization methods. These new methods are tested on a suite of benchmark topology optimization problems in Section 5. Finally, Section 6 offers conclusions.

Density-based topology optimization with Helmholtz filtering
In this section, we present a computational procedure for topology optimization where we seek a structure that minimizes a given objective function (e.g., compliance) under a volume constraint. Our presentation of topology optimization follows [32,7]. We review a linear elasticity formulation (Section 2.1), review a Helmholtz filtering technique (Section 2.2), state the topology optimization problem (Section 2.3), derive the adjoint method to efficiently evaluate the objective function gradient (Section 2.4), and analyze the cost of the optimization procedure (Section 2.5).

Linear elasticity
We first review the finite element approximation of the linear elasticity equation for density-based topology optimization. We represent the geometry of the structure by a density distribution over a d-dimensional design domain, where the (near) zero and unity density indicate the absence and presence, respectively, of the material in a given region. To this end, we partition our design domain Ω into N e non-overlapping polygonal elements, tΩ e u Ne e"1 . This partition, or mesh, is delineated by N v vertices. The element-wise discontinuous density is represented as a vector ρ P P Ă R Ne with entries ρ " pρ 1 , . . . , ρ Ne q T , where ρ e is the density of element e and P :" rρ l , 1s Ne Ă R Ne is the space of admissible densities for a lower bound ρ l " 10´3. The volume constraint is expressed as where |Ω e | is the volume of element e. We model the response of the structure under the specified load condition using linear elasticity. To this end, we introduce a density-dependent stiffness matrix Kpρq " R dNvˆdNv and a (fixed) load vector f P R dNv associated with a linear (i.e., polynomial degree one) finite element approximation of the linear elasticity equations. Formally, K : P Ñ R dNvˆdNv maps the density vector to the stiffness matrix. To show the explicit dependence of the stiffness matrix on the density, we introduce the assembly operator P e P R dNvˆdN e v for element e " 1, . . . , N e . The assembly operator P e is a subset of the columns of the dN vˆd N v identity matrix such that v e " P T e v, where v e P R dN e v is the elemental displacement vector and v P R dNv is the global displacement vector. The stiffness matrix and load vector are represented in terms of their element contributions as where K e P R dN e vˆd N e v is the element stiffness matrix for unity density and f e P R dN e v is the element load vector, both of which are independent of ρ. The exponent p is the penalization parameter; we choose p " 3 which ensures good convergence to 0-1 solutions [33]. If all elements are the same size and shape, all element stiffness matrices are identical; this is the case for typical density-based topology optimization problems [33].
The displacement field of the linear finite element approximation is represented by a displacement vector u P R dNv associated with the N v mesh vertices. To impose essential (or Dirichlet) boundary conditions, we partition the global displacement vector into N u u unconstrained degrees of freedom u u P R N u u and N c u constrained degrees of freedom u c P R N c u such that Notice this partition assumes a particular ordering of the global degrees of freedom. In this work, the constrained degrees of freedom correspond to clamped displacements so we have u c " 0. Given ρ P P , the unconstrained degrees of freedom of the solution uůpρq P R N u u must satisfy r u puůpρq; ρq " 0; here r u : R N u uˆP Ñ R N u u is the residual operator associated with the force-equilibrium condition given by where K uu pρq P R N u uˆN u u and f u P R N u u are the portions of the stiffness matrix and load vector associated with unconstrained degrees of freedom.

Helmholtz density filter
As noted in numerous works and summarized in a review paper [33], the topology optimization problem is ill-posed if we choose the density field as the design variable ρ P P , which defines the elemental density of the stiffness matrix Kpρq. To stabilize the formulation, we choose a separate "unfiltered" density field ψ as the design variable, which implicitly defines the "filtered" density ρ P P . The unfiltered density ψ P Ψ Ă R Ne is expressed as ψ " pψ 1 , . . . , ψ Ne q T , where ψ e is the unfiltered density associated with element e and Ψ :" rρ l , 1s Ne is the space of admissible unfiltered densities. There exists many filtering procedures [33]; in this work we consider the so-called Helmholtz filter [22].
We provide a brief description of the Helmholtz filter; we refer to [22] for details. The main idea behind the Helmholtz filter is to relate the (function representation of) the unfiltered density field ψ and the filtered density field ρ by the Helmholtz equation´r 2 ∆φ`φ " ψ with homogeneous Neumann boundary condition, where r " 2 ? 3R and R is the characteristic radius of the filter. To this end, given an elemental unfiltered density vector ψ P Ψ Ă R Ne , we first compute nodal filtered density φ ‹ pψq P R Nv that satisfies a finite element approximation of the Helmholtz equation: here R φ : R NvˆΨ Ñ R Nv is the Helmholtz residual operator given by where H P R NvˆNv is the assembled Helmholtz "stiffness matrix", and b : Ψ Ñ R Nv is an operator that maps the unfiltered density vector to the Helmholtz "load vector". To show the explicit dependence of the load vector on the density, we introduce the Helmholtz assembly operator Q e P R NvˆN e v for element e " 1, . . . , N e . The assembly operator is a subset of the columns of the N vˆNv identity matrix such that w e " Q T e w, where w e is the elemental vector and w P R Nv is the global vector. The stiffness matrix and load vector are represented in terms of their element contributions as where H e P R N e vˆN e v is the element Helmholtz stiffness matrix and b e P R N e v is the element Helmholtz load vector, both of which are independent of ψ. If all elements are the same size and shape, all element stiffness matrices are identical. Given the nodal filtered density φ P R Nv , the elemental filtered density ρ ‹ pφq P P is given by an element-wise nodal averaging operator ρ˚: R Nv Ñ P given by We also introduce the filtering operator ρ ‹ : Ψ Ñ P from elemental unfiltered density ψ P Ψ to elemental filtered density ρ P P given by ρ ‹ e pψq :" ρe pφ ‹ pψqq " Before we conclude this section, we recall that we wish to impose volume constraint (1) on the filtered density ρ, which describes the physically relevant density distribution. Owing to the volume preservation property of the Green operator of the Helmholtz equation [22], this volume constraint on the filtered density field ρ is satisfied as long as the unfiltered density field ψ satisfies the identical volume constraint Hence, the volume constraint can be imposed as a linear constraint on the unfiltered density field.

Optimization problem
We now introduce an objective function j : R N u uˆP Ñ R. Our formulation will treat general objective functions, but we assume that the function is twice continuously differentiable to develop convergence proofs. One objective function that is of particular interest is the compliance output j c : R N u uˆP Ñ R given by where f is the load vector given by (2). For notational convenience, we introduce the mapping u ‹ u : Ψ Ñ R N u u from the unfiltered density ψ to the unconstrained elasticity displacements u ‹ u pψq given by u ‹ u pψq :" uůpφ ‹ pψqq. Note that u ‹ u pψq satisfies r u pu ‹ u pψq, ρ ‹ pψqq " 0 for any ψ P Ψ. In the remainder, we will refer to the mapping u ‹ u as the (primal) high-dimensional model (HDM). We also introduce the topology optimization objective function J : Ψ Ñ R that maps the unfiltered density ψ P Ψ to the output of the objective function j : R N u uˆP Ñ R that properly accounts for the solution of the elasticity and Helmholtz equations: i.e., The final form of the topology optimization problem is as follows:

Derivative computations via the adjoint method
Given the large number of design variables in topology optimization problems, its efficient solution requires gradient-based optimization methods. That is, in addition to the objective function Jpψq, we also require its gradient ∇J : Ψ Ñ R Ne . As Jpψq " jpu ‹ u pψq, ρ ‹ pψqq, the gradient can be expanded as (Throughout this work, we adhere to the convention where the partial derivative of scalar-valued and vectorvalued functions yield a row vector and a matrix, respectively.) The gradient ∇J requires matrix-vector products of two sensitivities, Bu ‹ u Bψ and Bρ ‹ Bψ , with the partial derivatives of j : R N u uˆR Ne Ñ R, which can be derived analytically from the expression for jpu u , ρq. For compliance output, we have In the remainder of this section, we derive adjoint expressions to efficiently evaluate We first consider the evaluation of Bρ ‹ Bψ T v for an arbitrary v P R Ne . From the definition of φ ‹ , the Helmholtz residual is zero irrespective of variations in ψ P R Ne , which in turn implies the total derivative of the Helmholtz residual with respect to ψ is zero, i.e., Expanding the total derivative in (14) via the chain rule and dropping arguments, we have which can be rearranged to obtain the following expression for the product of From the definition of the Helmholtz residual (6) and (7), we have BR φ Bφ pφ; ψq " H, BR φ Bψ e pφ; ψq "´Q e b e , e " 1, . . . , N e where summation is not implied over the repeated index, which yields where µpvq P R Nv is the adjoint variable that satisfies Finally, from the definition of ρ ‹ e : R Nv Ñ R in (9), we obtain where qpvq " ř Ne We next consider the evaluation of Bu ‹ u Bψ T v for an arbitrary v P R N u u . From the definition of u ‹ u , the elasticity residual is zero irrespective of variations in ψ P R Ne , which in turn implies the total derivative of the elasticity residual with respect to ψ is zero, i.e., Expanding the total derivative in (18) via the chain rule and dropping arguments, we have which can be rearranged to obtain the following expression for the product of here λ ‹ : R Ne Ñ R N u u is the adjoint map such that the adjoint variable λ ‹ pψq satisfies where we have appealed to the definition of the elasticity residual (4) to obtain Bru Buu pu u ; ρq " K uu pρq. Combining (12), (17), (20), the gradient ∇J can be written compactly as where λ ‹ is defined in (21). This shows the primary operations required to compute the gradient are an elasticity adjoint solve (to compute λ ‹ ) and a Helmholtz adjoint solve (to apply Bρ ‹ Bψ to a vector).

Remark 2.
In the special case of compliance minimization, the gradient computation simplifies significantly.
Due to the special form of the compliance function in (13) where Bj Bu u T " f u and the symmetry of the matrix K uu , the adjoint system in (21) is identical to the primal system in (4) and therefore we have Due to this equivalence between the primal and adjoint elasticity state in this special case, the adjoint solve would be redundant. Furthermore, the compliance does not directly depend on the density field (13) and the gradient reduces to Thus the gradient computation only involves a Helmholtz adjoint solve to apply Bρ ‹ Bψ to a vector.

Computational cost
To close this section, we analyze the cost of topology optimization. For simplicity, we focus on compliance minimization problems for which λ ‹ pψq " u ‹ u pψq and hence no explicit computation of the adjoint is required; see Remark 2. We decompose the mapping u ‹ u : Ψ Ñ R N u u from the (unfiltered) density ψ P Ψ to the linear elasticity displacement u ‹ u pψq into five steps, and analyze the cost from both the theoretical and practical perspectives; the practical assessment is based on the absolute and fractional run-time result shown in Figure 1 for a typical two-dimensional topology optimization problem (MBB beam; Section 5.1).
1. The computation of bpψq P R Nv by (7). This vector assembly requires OpN e q operations. Figure 1 shows that in practice this cost is negligible. 2. The solution of the Helmholtz equation (5) for φ P R Nv . Because the Helmholtz stiffness matrix H P R NvˆNv does not depend on the design variable ψ P Ψ, the Cholesky factorization of the matrix can be computed for once and for all. Given the Cholesky factors associated with a minimal-fill ordering, the marginal cost of each Helmholtz solve is OpN e logpN e qq and OpN 4{3 e q in two-and three-dimensions, respectively. Figure 1 shows that in practice the cost grows superlinearly but is negligible.
3. The computation of ρ˚pφq P P by (8). The vector assembly again requires OpN e q operations. Figure 1 shows that in practice this cost is negligible.
4. The assembly of the linear elasticity stiffness matrix Kpρq by (2). The assembly is done in two steps: the evaluation of the density-weighted elemental stiffness matrices ρ p e K e ; the assembly of the matrices into a sparse matrix. While both operations are OpN e q in theory, the latter operation dominates the run time on a modern computational platform as the operation requires unstructured memory access. Figure 1 shows that this is one of the dominant costs of the map u ‹ u : Ψ Ñ R N u u in practice.
5. The solution of the linear elasticity system (3) for uůpρq P R N u u . The operation is carried out using a sparse Cholesky solver. Using a minimal-fill ordering, the cost of the linear solve is OpN 3{2 e q and OpN 2 e q in two-and three-dimensions, respectively. Figure 1 shows that this is the dominant cost of the map u ‹ u : Ψ Ñ R N u u in practice, especially as the problem size increases.
The cost analysis shows that the assembly and solution of the linear elasticity system accounts for roughly 98% of the cost in computing u ‹ u : Ψ Ñ R N u u and the cost of Helmholtz filtering is negligible. In Section 3 we introduced a surrogate model to accelerate the solution of the linear elasticity problem.

Projection-based model reduction for inexpensive approximation of linear elasticity equations
In this section we introduce the projection-based reduced-order model (ROM) which will be used as the approximation model in our trust-region method. We present the Galerkin ROMs to approximate the solution u ‹ pψq and adjoint λ ‹ pψq (Sections 3.1 and 3.2), residual-based error estimates for ROM solutions (Section 3.3), and the computational cost associated with the construction and evaluation of the ROM (Section 3.4).

Galerkin reduced-order model for the primal solution
We first discuss the construction of the Galerkin ROM for the primal solution u ‹ pψq. To this end, we introduce a reduced basis (matrix) Φ k P R dNvˆk ; here, k ď dN v formally, but k ! dN v in practice. (We will detail the construction of Φ k in Section 4.2; for now we assume it is given.) We then represent the ROM solution u ‹ k pψq P R dNv as a linear combination of the reduced basis: u ‹ k pψq " Φ kû ‹ k pψq. The coefficients u ‹ k pψq are given by the following residual statement: given ψ P R Ne , find the ROM coefficientsû ‹ k pψq P R k such thatr where the ROM residualr k : R kˆN e Ñ R k is given bŷ and the ROM stiffness matrixK k : R Ne Ñ R kˆk and load vectorf k P R k are given bŷ The solution to (25) satisfies a linear system: given ψ, findû ‹ k pψq such that Finally, the ROM objective function J k : Ψ Ñ R is given by We make a few remarks: Remark 3. Assuming the columns of Φ k are linearly independent, the ROM stiffness matrixK k pρq is symmetric positive definite and the ROM problem is well-posed.

Remark 4.
The Galerkin approximation is optimal in the energy norm in the sense that, for any given where the energy norm, which depends on ρ P R Ne , is defined by }w} Kpρq :" a w T Kpρqw @w P R dNv , @ρ P R Ne .
Remark 5. Suppose we are given a basis matrix Φ k P R dNvˆk , and we construct a new basis matrix Φ k`1 " rΦ k , ys P R dNvˆpk`1q by augmenting Φ k by a vector y P R dNv . If u ‹ k pψq and u ‹ k`1 pψq are the ROM approximations associated with Φ k and Φ k`1 respectively, then the energy norm of the error is a non-increasing function under basis enrichment.
Remark 7. The expressionK k pρq " ř Ne e"1 ρ p e pP T e Φ k q T K e pP T e Φ k q emphasizes that the ROM stiffness matrix K k pρq can be assembled using only dense operations. The elemental reduced basis P T e Φ k P R dN e vˆk comprises the rows of Φ k associated with the element e; we then perform the dense operation pP T e Φ k q T K e pP T e Φ k q for each element. The use of dense operations enables an efficient assembly of the ROM residual using BLAS operations. The same remark applies to the construction off k .

Galerkin reduced-order model for the adjoint solution
We now discuss the construction of the ROM for the approximation of the adjoint. We again assume that the reduced basis (matrix) Φ k P R dNvˆk for k ! dN v is given. As before, given ψ P R Ne , we construct an approximation of the form λ ‹ k pψq " Φ kλ ‹ k pψq P R dNv to the solution λ ‹ pψq P R dNv using Galerkin projection. The coefficientsλ ‹ k pψq P R k are given by the following adjoint problem: given ψ P R Ne , find the ROM coefficientsλ ‹ k pψq P R k such that where the ROM stiffness matrixK k : R Ne Ñ R kˆk is given by (26) and the adjoint load vectorĝ k : R kˆRNe Ñ R k is given bŷ For the compliance output, the right hand side simplifies toĝ k pu ‹ k pψq, ρ ‹ pψqq "f k and hence λ ‹ k pψq " u ‹ k pψq. Remark 8. Remarks 4-7 also hold for the adjoint system. Namely, the Galerkin approximation is optimal in the energy norm in the sense that, for any given ψ P R Ne , which is the counterpart of Remark 4. The optimality implies that the energy norm of the error in the adjoint is a non-increasing function under basis enrichment (Remark 5). In addition, if Φ k is chosen such that u ‹ pψq P ImgpΦ k q and λ ‹ pψq P ImgpΦ k q, then the ROM reproduces the adjoint λ ‹ pψq; i.e., λ ‹ k pψq " λ ‹ pψq (Remark 6). Note that we require both u ‹ pψq and λ ‹ pψq to be in ImgpΦ k q, because the adjoint equation itself depends on u ‹ k pψq and the adjoint would not be exact if the primal solution u ‹ k pψq is not exact. Finally, the adjoint system can be efficiently assembled using dense linear algebra operations (Remark 7).

Error estimation
We now introduce error estimates for our ROM approximation of the solution, adjoint, objective function, and objective function sensitivity. To begin, we introduce the key ingredients of our (a posteriori ) error estimates. The (primal) residual, which has been defined earlier, is given by The adjoint residual is given by the adjoint residual depends on both the adjoint approximation v and the linearization point of the ROM adjoint problem u ‹ k pψq. In addition, given any arbitrary matrix A P R mˆn , we denote its maximum and minimum singular values by σ max pAq and σ min pAq, respectively. Finally, given a symmetric positive matrix A P R mˆm , the symmetric positive definite matrix A 1{2 is the matrix that satisfies A " A 1{2 A 1{2 .
We now state a series of theorems that relates various errors to the residuals; the associated proofs are provided in Appendix C.
Theorem 1. For any given ψ P R Ne , the energy norm of the error in the solution is bounded by Theorem 2. For any given ψ P R Ne , the energy norm of the error in the adjoint is bounded by with all variables and forms evaluated about ψ and ρ ‹ pψq.
Theorem 3. For any given ψ P R Ne , the error in the objective function is bounded by where B is given by with all variables and forms evaluated about ψ and ρ ‹ pψq.
Theorem 4. For any given ψ P R Ne , the error in the shape sensitivity is bounded by where the entries of matrices C P R NeˆdNv and D P R NeˆdNv are given by where B is given by (29), all variables and forms evaluated about ψ and ρ ‹ pψq, and the summation on repeated indices is implied.
Remark 9. If the objective function j : R dNv Ñ R is linear, then the matrix A given by (28) and consequently B vanishes. As a result, the adjoint and output error estimates in Theorems 2 and 3 simplify to In addition, in Theorem 4, the last term of C in (30) vanishes.
Remark 10. For the compliance output, the output and sensitivity error estimates further simplify to because λ ‹ k pψq " u ‹ k pψq and rp¨;¨q " r adj p¨;¨q. We may also obtain a "looser-version" of Theorem 3, which does not involve the adjoint residual: Theorem 5. For any given ψ P R Ne , the error in the shape sensitivity is bounded by where B is given by (29).
Theorems 1-5 show that all quantities that are relevant in topology optimization-and in particular the objective function value and the associated shape sensitivity-are well approximated by the ROM as long as the primal and adjoint residuals are small. In other words, the primal and adjoint residuals serve as an indicator of the errors in the ROM approximations (up to a constant). We will leverage this observation to devise residual-aware ROM-accelerated topology optimization strategies in Section 4.

Computational cost
We now assess the (online) computational cost of the ROM solution and the associated residual.
Solution evaluation. We first analyze the cost of ROM solution evaluation: ρ ‹ pψq Þ Ñû k pψq. We consider the evaluation of the primal solution; the cost for the adjoint solution can be analyzed in a similar manner. Given the density distribution ρ ‹ pψq P R Ne and a reduced basis matrix Φ k P R dNvˆk , we decompose the computation of the ROM solution u k pρq P R dNv into three steps and assess the associated costs: 1. Assembly of the ROM stiffness matrix and load vector: given ρ ‹ pψq P R Ne , assemblê The computation ofK k pρ ‹ pψqq dominates the cost of this step. The operation count is N e p2kpdN e v q 22 k 2 pdN e v qq; assuming k Á dN e v , the cost is hence « 2pdN e v qk 2 N e . For Q 1 finite elements in two and three dimensions, the cost is « 16k 2 N e and 48k 2 N e , respectively. The cost of this step scales linearly with the number of finite elements N e and quadratically with the dimension of the reduced basis space k. We also note thatf k is independent of ρ ‹ pψq, and hencef k needs to be recomputed only when the reduced basis Φ k is updated; this is unlikeK k pρ ‹ pψqq which must be recomputed for every new density ρ ‹ pψq encountered during the optimization procedure.
2. Solution of the ROM linear system: findû k pψq P R k such that The solution of the linear system requires « k 3 {3 operations using Cholesky factorization.

Representation of the solution in the original vector space
The multiplication requires 2dN v k operations. For Q 1 finite element in two and three dimensions, N v « N e and hence the cost of this step is « 4kN e and 6kN e , respectively. The cost of this step scales linearly with both N e and k.
The first step dominates the overall computational cost, and hence the cost of finding the ROM solution is approximately 2pdN e v qk 2 N e (assuming k Á dN e v q; for Q 1 finite element in two and dimensions, the cost is « 16k 2 N e and 48k 2 N e , respectively.
Residual evaluation. We next analyze the cost of evaluating the 2-norm of the residual:û ‹ k pψqˆρ ‹ pψq Þ Ñ }rpΦ kû ‹ k pψq; ρ ‹ pψqq} 2 . In particular, our interest is in the marginal cost of the residual evaluation for many different densities ρ ‹ pψq and ROM coefficientsû k pψq for a fixed reduced basis Φ k . To minimize the marginal cost, we first precompute the density-independent quantity for each e " 1, . . . , N e . This computation requires pdN e v q 2 kN e operations, but can be performed once and for all for a fixed Φ k . Given the precomputed quantity A e , the computation of the residual for each ρ ‹ pψq and the associated solutionû k pψq is performed in two steps: 1. Construction of the residual vector rpu ‹ k pψq; ρ ‹ pψqq: compute  Figure 2: A decomposition of the runtime to perform the mappings uk : ρ Þ Ñ R dNv and pu, ρq Þ Ñ }rpu, ρq} 2 for a typical two-dimensional problem as a function of reduced basis size for a fixed mesh with 600ˆ200 Q 1 elements (left) and as a function of mesh size for reduced basis sizes of k P r5, 20s (right): evaluation and assembly ofK k pρq ( ); solution of reduced elasticity equations ( ); reconstruction of u k " Φ kûk ( ); and evaluation of }rpu, ρq} 2 ( ). In the right plot, for a fixed color, the lower line corresponds to a reduced basis of size k " 5, the upper line corresponds to k " 20, and the shaded region corresponds to k P p5, 20q.
The two major costs of the evaluation are associated with (i) the computation of rpρ ‹ e pψqq p A e u k pψqs, which requires « dN e v kN e operations, and (ii) the sparse assembly ř Ne e"1 P e r¨s, which is OpdN e v N e q (though the actual cost is highly dependent on the implementation and computer architecture). For a sufficiently large k, the overall cost is « dN e v kN e ; for Q 1 finite element in two and three dimensions, the cost is « 8kN e and 24kN e .
2. Evaluation of the 2-norm }rpu ‹ k pψq; ρ ‹ pψqq} 2 . This cost of the norm evaluation is « 2dN v . For Q 1 finite elements in two and three dimensions, the cost is « 4N e and 6N e , respectively.
The first step dominates the overall computational cost, and hence the marginal cost of evaluating the residual is approximately « dN e v kN e , which evaluates to « 8kN e and 24kN e for Q 1 finite element in two and three dimensions. For a sufficiently large k, the marginal cost of residual evaluation is insignificant compared to the cost of the solution evaluation.
We compare the cost of the ROM analysis to the finite element analysis. For the finite element analysis, we consider the MBB beam configuration in Section 5.1, which has three times as many Q 1 elements in the first dimension as in the second, and we employ a direct solver with an asymptotically optimal (i.e., nested dissection) ordering. We recall from Section 2.5 that the computational cost for the finite element analysis in two and three dimensions are OpN  To provide a concrete cost assessment of the ROM approximation in topology optimization, we show in Figure 2 a breakdown of the ROM runtime as a fraction of the HDM runtime. Consistent with the theory, for a sufficiently large k the dominant cost of the ROM solve and residual evaluation is associated with the construction of the ROM stiffness matrix; nevertheless, the runtime for the assembly is a small fraction of the HDM solve, varying from « 0.3% for k " 5 to « 1% for k " 20. In other words, we can perform « 100 ROM solves and residual evaluations in the time it takes to complete a single HDM solve.

Globally convergent, error-aware trust-region method for efficient topology optimization
In this section, we introduce an extension of the error-aware trust-region method developed in [42] for unconstrained problems to problems with convex constraints. We specialize the method to the present topology optimization setting and use the projection-based ROM introduced in Section 3 as the trust-region model and the residual as the error-aware trust-region constraint. Global convergence theory for the proposed method is provided in Appendix A.

Error-aware trust-region method
We begin by introducing the proposed error-aware trust-region method as a general method to solve an optimization problem with convex constraints. To this end, consider the following abstract optimization problem on a convex set C Ă R N : minimize where F : C Ñ R and C satisfy the following assumptions: Assumption 1 (Objective function). For a domain C Ă R N that satisfies Assumption 2, the objective function F : C Ñ R satisfies the following conditions: 1) F is twice continuously differentiable on C; 2) F is bounded from below on C; and Assumption 2 (Constraints). The constrained set C Ă R N satisfies the following conditions: 2) C is nonempty, closed, convex; and 3) a first-order constraint qualification holds at any critical point of (32).
To solve (32), we introduce Algorithm 1 that produces a sequence of points tx k u 8 k"1 , which we call trust-region centers, that converges to a first-order critical point. We use the following criticality measure χpxq :"ˇˇˇˇˇˇmin where χpx˚q " 0 implies x˚is a first-order critical point of (32) [13]. At each trust-region center x k , we introduce an inexpensive model m k : C Ñ R intended to approximate the objective function F in the error-aware trust region where ϑ k : C Ñ R ě0 is the trust-region constraint and ∆ k P R ą0 is the trust-region radius. Standard trustregion algorithms take ϑ k pxq " }x´x k } 2 ; however, we generalize the notion of a trust-region constraint to allow for trust regions that take into account the error in the approximation model. We require the approximation models and trust-region constraints to satisfy the following assumptions: Assumption 3 (Approximation model). For a domain C Ă R N that satisfies Assumption 2, approximation models m k : C Ñ R associated with trust-region centers x k P C, k P N, satisfy the following conditions: 1) m k is twice continuously differentiable on C; 2) m k px k q " F px k q; 3) ∇m k px k q " ∇F px k q; 4) there exists β ą 0 (independent of k) such that β k ď β for all k P N, where Assumption 4 (Trust-region constraint). For a domain C Ă R N that satisfies Assumption 2, the trustregion constraints ϑ k : C Ñ R ě0 associated with the trust-region centers x k P C, k P N, satisfy the following conditions: 1) ϑ k is twice continuously differentiable on C; 2) there exists κ ∇ϑ ą 0 (independent of k) such that max xPC }∇ϑ k pxq} ď κ ∇ϑ for all k P N; 3) ϑ k px k q " 0; 4) there exist ζ ą 0 and ν ą 1 (independent of k) such that |F pxq´m k pxq| ď ζϑ k pxq ν for all x P B k .
Algorithm 1 Error-aware trust-region method with convex constraints 1: Initialization: Given x 0 , ∆ 0 , ∆ max , 0 ă γ 1 ă γ 2 ă 1, 0 ă η 1 ă η 2 ă 1 2: Model and constraint update: Choose a model m k : C Ñ R and trust-region constraint ϑ k : C Ñ R that satisfy Assumptions 3 and 4. 3: Step computation: Approximately solve the trust-region subproblem min xPB k m k pxq for a candidate stepx k that satisfies the fraction of Cauchy decrease condition (37). 4: Actual-to-predicted reduction: Compute actual-to-predicted reduction ratio approximation k :" F px k q´F px k q m k px k q´m k px k q .

5:
Step acceptance: The error-aware trust-region algorithm is described in Algorithm 1. The algorithm is initialized from an initial guess for the optimization solution x 0 , an initial trust-region radius ∆ 0 , and a number of other (standard) constants [13]. At iteration k, an approximation model and trust-region constraint are constructed that satisfy Assumptions 3 and 4 to generate the trust-region subproblem min xPB k m k pxq.
We (approximately) solve this subproblem whose solution, denotedx k , is a candidate for the next trustregion center x k`1 . It is important to note that the trust-region subproblem does not need to be solved exactly; it only needs to find a point that satisfies the fraction of Cauchy decrease condition: where κ P p0, 1q and κ 1 ą 0 are constants (independent of k), and β k is given by (35). Theorem 6 establishes the existence of such a point within the trust region.
Proof. See Appendix A.
Finally, we compute the actual-to-predicted reduction, and use it to assess whether to accept or reject the candidate step as the next trust-region center. If ρ k ą η 1 for some constant η 1 P p0, 1q, the candidate is accepted; otherwise it is rejected. In addition, the actualto-predicted reduction is used to modify the trust-region radius for the next iteration. If ρ k ă η 1 , the step is called unsuccessful and the radius is decreased to ∆ k`1 P rγ 1 ∆ k , γ 2 ∆ k s, where 0 ă γ 1 ă γ 2 ă 1. If ρ k P rη 1 , η 2 q, where η 2 P pη 1 , 1q is a constant, the step is considered successful and the radius is updated to lie in the range ∆ k`1 " rγ 2 ∆ k , ∆ k s (usually ∆ k`1 " ∆ k ). Finally, if ρ k ě η 2 , the step is called very successful and the radius is increased or kept the same. This error-aware trust-region algorithm (Algorithm 1) is globally convergent: for any x 0 P C and ∆ 0 P R ą0 , the algorithm generates a sequence tx k u 8 k"1 whose any convergent subsequence converges to a first-order critical point of (32). The result is summarized in the following theorem: where χ : C Ñ R is the criticality measure (33).
Proof. See Appendix A.

Error-aware trust-region method for topology optimization
We propose an error-aware trust-region method (Section 4.1) to efficiently solve the topology optimization problem (11) using projection-based ROMs introduced in Section 3 as the approximation model. The topology optimization problem (11) exactly fits the form of the general optimization problem (32), where the optimization variables are the unfiltered element densities x :" ψ P R N -Ne , the objective function is the topology optimization objective function (10) and the constraint set is the intersection of the simple bounds and volume constraint As we will see in Theorem 8, the objective function and constraint satisfy Assumptions 1 and 2. We now outline the procedure to construct an approximation model based using projection-based reducedorder models (Section 3) that satisfies the requirements for global convergence in Assumption 3. Suppose k trust-region iterations have been completed and we collected the density fields defining the trust-region centers Ψ k :" pψ p0q , . . . , ψ pkq q.
We collect the associated primal and adjoint solutions in matrices of the form U k :" ru ‹ pψ p0q q,¨¨¨, u ‹ pψ pkq qs P R dNvˆpk`1q , Λ k :" rλ ‹ pψ p0q q,¨¨¨, λ ‹ pψ pkq qs P R dNvˆpk`1q , and apply POD to the first k snapshots to obtain "compressed" matrices Φ U k´1 ,POD n k :" POD n k pU k´1 q P R dNvˆn k , where POD n : R dNvˆpk´1q Ñ R dNvˆn applies the singular value decomposition to the input matrix U (resp. Λ) and extract the first n left singular vectors to form Φ U ,POD n (resp. Φ Λ,POD n ). We combine the POD bases and the primal and adjoint solution at the trust-region center, u ‹ pψ pkq q and λ ‹ pψ pkq q, and orthonormalize the resulting matrix using the Gram-Schmidt procedure to obtain the reduced basis (matrix): where n k P t0, 1, . . . , k´1u is the number of vectors retained after POD is applied to the primal and adjoint snapshot matrices and j k :" 2pn k`1 q is the size of the reduced basis. Then the approximation model is the ROM objective function (27): where u ‹ j k pψq is the solution of the ROM associated with the basis Φ j k . The choice of basis in (43) guarantees that u ‹ pψ pkq q, λ ‹ pψ pkq q P ImgpΦ j k q regardless of n k . Hence by Remarks 6 and 8, the primal and adjoint solution of the ROM will be exact at the trust-region center ψ pkq : As a result, the approximation model and its gradient agree with the original topology optimization objective at the trust-region center (Theorem 3), which establishes Assumptions 3.2-3.3, regardless of n k .
Remark 11. For compliance output, the primal and adjoint solutions are identical: λ ‹ pψ piq q " u ‹ pψ piq q, i " 0, . . . , k. Hence, we need to construct only the primal snapshot matrix U k´1 and the associated POD matrix Φ U k´1 ,POD n k . The reduced basis (matrix) is then given by where j k " n k`1 . Note that, because the columns of Φ U k´1 ,POD n k are orthonormal, we need to perform just one step of Gram-Schmidt to orthonormalize u ‹ pψ pkq q with respect to Φ U k´1 ,POD n k .
In this work, we consider two instances of the trust-region constraint. The first is the traditional trustregion constraint of the form ϑ k pψq :" This constraint is simple to implement, efficient to evaluate, and will lead to a globally convergent method [13]; however, it does not account for the error in the approximation model. The other trust-region constraint we consider is the norm of the primal residual evaluated at the reconstructed ROM state since it bounds the error in the approximation model (Theorem 3) and provides a notion of an error-aware trust region. Theorem 8 will verify (46) satisfies Assumption 4, thus leading to a globally convergent method. Furthermore, this choice of trust-region constraint is relevant given that the residual norm has proven to be a popular indicator to trigger adaptation to a reduced basis constructed on-the-fly, both in the context of topology optimization [15,12,39] and more generally [43]; however, these methods are heuristic and cannot guarantee global convergence. The proposed trust-region framework provides a rigorous setting that can be used to adapt existing methods such that they are guaranteed to converge to a local minima from any starting point.

5:
Step acceptance: if k ě η 1 then ψ pk`1q "ψ pkq else ψ pk`1q " ψ pkq end if 6: Trust region update: The final error-aware trust-region method based on ROMs for efficient topology optimization is summarized in Algorithm 2.
Remark 12. The ROM-based topology optimization method that uses traditional trust regions, i.e., the trust-region constraint in (45), is given by Algorithm 2 with the definition of ϑ k in Step 2 replaced with (45).
The algorithm produces a sequence of trust-region centers that converges (liminf sense) to a critical point of (11); see Theorem 8. Lemma 1. Suppose the output functional j : R dNvˆP Ñ R is twice continuously differentiable, and we choose ϑ k pψq :" for P p0, 1q. The error-aware trust-region method based on ROMs for topology optimization satisfy Assumptions 1-4. For " 0, all conditions except the fourth condition of Assumption 4 are satisfied.
Remark 13. In practice we set " 0 in Lemma 1. As a result, the fourth condition of Assumption 4 is not satisfied.
Theorem 8. Consider the optimization problem in (32) with the objective function (41) where j : R N u uˆP Ñ R is twice continuously differentiable with respect to both of its arguments and the feasible set is defined in (42) with V ě ř Ne e"1 ρ l |Ω e |. Then the sequence of trust-region centers produced by Algorithm 2 satisfies lim inf kÑ8 χpψ pkq q " 0 (47) regardless of n k .
Proof. The result is a direct consequence of Lemma 1 and Theorem 7.
To close, we discuss two pertinent details for the implementation of Algorithm 2 in practice. First, the global convergence theory is independent of n k ; however, the behavior of the algorithm does depend on n k . For example, if we take n k " 0, each ROM solve will be extremely fast (since j k " 2 for all k), but the ROM will have little predictive capability and a large number of trust-region center updates will be required. On the other end of the spectrum, choosing n k " k´1 will lead to expensive ROM solves when k becomes large; however, fewer HDM solves will be required because the ROM basis will be as rich as possible. The approach we take is to choose n k " mintk´1, n max u, where n max is the maximum truncation size for the primal/adjoint reduced-order bases, which guarantees the ROM size will be bounded j k ď 2pn max`1 q. Thus n max should be chosen as large as possible such that individual ROM solves are sufficiently fast.

Remark 14.
In practice, we only retain the L most recent trust-region centers in U k´1 and Λ k´1 , rather than all previous centers. This ensures that POD is applied to at most L snapshots and hence controls the POD cost. Furthermore, since the matrix U k´1 (resp. Λ k´1 ) is a low-rank update to U k´2 (resp. Λ k´2 ), efficient algorithms for updating the singular value decomposition [10,12], adapted to the truncated singular value decomposition (POD) in [6,38,26], can be used to ensure the POD cost is negligible in comparison to a HDM solve.
The other practical issue that must be addressed is how to efficiently solve the trust-region subproblem to ensure the fraction of Cauchy decrease condition holds. Given Theorem 6 that establishes the existence of a point in the trust region that satisfies (37), the (global) minimizer of (36) is guaranteed to satisfy (37). Recall from Figure 2 that ROM evaluations are « 100ˆcheaper than HDM evaluations; however, they are not free so trust-region subproblem solvers that require hundreds or thousands of ROM evaluations will not be competitive. This immediately eliminates the use of "global" optimization methods due to the large number of objective evaluations required. In addition, our experiments using gradient-based optimization methods to find a local solution to (36) also required several hundred iterations to even converge to relatively weak optimality tolerances. The approach that proved most effective in practice was to use a gradient-based optimization procedure (Method of Moving Asymptotes [36] in this work) to solve min ψPC m k pψq, i.e., the trust-region subproblem without the trust-region constraint, and terminate the iterations once the trust-region constraint was violated ϑ k pψq ą ∆ k . While this is not guaranteed to produce a step that satisfies the fraction of Cauchy decrease, it works well in practice as we show in the next section.

Numerical experiments
In this section, we study the performance of the proposed ROM-based trust-region methods for topology optimization on three standard benchmark problems: MBB beam, cantilever beam, and simply supported beam. We will consider two instances of the method, one where the trust region is taken as the sublevel sets of the HDM residual (46), subsequently called ROM-TR-RES, and one with traditional trust regions ϑ k pψq " › › ψ´ψ pkq › › 2 [13], to be called ROM-TR-DIST. Given its popularity, we use MMA applied to the original (unreduced) topology optimization, referred to as HDM-MMA, as the standard for comparison. In the remainder of this section, we use the terminology major iteration, to refer to an update from ψ pkq to ψ pk`1q . In all methods considered, a major iteration involves a single HDM evaluation; the ROM-based methods additionally require a number of ROM evaluations.
We use the computational cost required to drive the objective function to within a prescribed tolerance of its optimal value to quantify the performance of the ROM-based optimization algorithms. A more rigorous convergence metric is based on the Karush-Kuhn-Tucker (KKT) conditions; however, it is usually not practical to satisfy these conditions to high precision using MMA as it would require a large number of iterations. Furthermore, the constraints are satisfied in essentially all optimization iterations for all methods considered, and hence the objective function (alone) serves as a reliable surrogate of KKT. For all benchmark problems, we take the optimal value of the objective function, denoted J˚, to be the 2000th iteration of the HDM-MMA algorithm; convergence is achieved for the smallest n P N such that where ą 0 is the (relative) convergence tolerance. The computational cost required to converge to a tolerance of , denoted C , will be measured in units of equivalent HDM solves, i.e., where T is the time required for the method under consideration to converge the topology optimization problem to the tolerance , and t HDM is the time required for a single HDM solve. Because the dominant cost of topology optimization comes from the HDM and ROM solves, T can be expanded as where N HDM and N ROM are the number of HDM and ROM solves, respectively, required for the method to converge the objective function to a tolerance , and t ROM is the time required for a single ROM solve. From Figure 2, t ROM depends on the size of the reduced basis. In our numerical experiments, we fix the maximum basis size to 20 (n k " 19) and therefore use the conservative estimate of t ROM " νt HDM , where ν " 0.01 from Figure 2. Therefore, the computational cost takes the form There are a number of user-defined parameters required by the ROM-based trust-region method; however, most of these are standard trust-region parameters, and we fix them at reasonable values: γ 1 " 0.5, γ 2 " 1, η 1 " 0.1, η 2 " 0.75, ∆ max " 100∆ 0 . The remaining parameter, ∆ 0 , is more delicate and potentially problem-specific so we study its impact on the ROM-based trust-region algorithms in the following sections. To ensure ∆ 0 is appropriately scaled, we consider initial trust-region radii of the form for the TR-ROM-RES and TR-ROM-DIST methods, respectively, where ψ p0q is the (feasible) initial design and τ P R ą0 controls the (scaled) initial radius.

MBB beam
We begin with the most canonical topology optimization problem: the MBB beam ( Figure 3). We use a finite element mesh consisting of 180ˆ60 bilinear quadrilateral elements, minimum filtering length scale of R " 0.12, and maximum volume V " 1 2 |Ω|. The optimal design is shown in Figure 3. For the studies in this section, all algorithms are initialized from the same feasible design: ψ p0q " 0.45¨1.
First, we demonstrate the significance of embedding the reduced topology optimization problem in the adaptive trust-region framework. To this end, we compare two ROM-accelerated methods (with HDM-MMA as the benchmark): ROM-TR-RES, the residual-based trust-region method introduced in Section 4.2; ROM-FIX-RES, an alternative that uses a fixed residual-based "trust-region" radius and does not scrutinize the trust-region step, i.e., accepts any candidate point produced by the trust-region subproblem. ROM-TR-RES (τ " 0.1) requires fewer HDM solves to drive the objective function to a fixed tolerance of the optimal value and converges to " 0.01 in 22 major iterations; HDM-MMA requires 39 iterations (Table 1). Furthermore, ROM-TR-RES is rather insensitive to the initial trust-region radius; we consider three orders of magnitude variation in τ P t0.01, 0.1, 1, 10u and performance degradation only emerges for the largest value τ " 10. The fastest convergence is obtained with τ " 0.1. In contrast, without the trust region adaptivity and candidate scrutiny, ROM-FIX-RES convergence is usually slower than the HDM-MMA algorithm, severely degrades as τ increases, and, in some cases, may not be achieved (Figure 4).
To assess the overall cost C of the ROM-based optimization methods, we account for the cost of the ROM solves as well as the HDM solves according to (52). We note only Op10 2 q ROM evaluations are required to drive the objective function to within the relative error of " 0.01 for reasonable values of τ P t0.01, 0.1, 1u. As a result, the ROM evaluations only add an equivalent of 1 to 2 HDM evaluations, and the ROM-based trust-region methods reduces the computational cost by nearly a factor of two over HDM-MMA (Table 1). Finally, Figure 4 shows the evolution of the trust-region radius; for the non-adaptive methods, the radius is fixed at its initial value, whereas the trust-region method adapts it. It is interesting to note that, for the adaptive method, the trust-region radius converges to nearly the same value, regardless of its starting point.
Remark 15. We repeated this experiment for a number of mesh resolutions, both coarser and finer than the 180ˆ60 quadrilateral mesh used here, and did not observe meaningful differences in the relative performance of the methods.  To highlight the progress of the ROM-based trust-region algorithm (ROM-TR-RES) relative to the HDM-MMA algorithm over the first 30 major iterations, we compare the evolution of the filtered density field ρ in Figure 5. After only 20 major iterations, the ROM-based method has visually converged to the optimal solution, whereas the HDM-MMA method still has volume to remove and boundaries to sharpen. Furthermore, by comparing the two algorithms at iteration 12, it is clear the ROM algorithm makes substantial progress toward the optimal design in the early iterations.
We close this section with a comparison of the two ROM-based trust-region methods: ROM-TR-RES and ROM-TR-DIST. Both algorithms perform similarly and are competitive relative to the HDM-MMA approach ( Figure 6, Table 1) with the ROM-TR-DIST method converging faster for weaker tolerances and ROM-TR-RES converging faster for tight tolerances.
Remark 16. We have extensively tested the ROM-based trust-region methods where the trust-region subproblem is solved exactly using nonlinear programming methods for this benchmark. We found a large number of ROM evaluations (Op10 3 q) are required to solve each subproblem, at least for the suite of optimization methods tested. Therefore, these methods are not competitive in terms of the overall cost. See also the discussion at the end of Section 4.2.   Figure 7: Setup and optimal design of cantilever beam. The point load F in " 3 is implemented as a distributed load of magnitude q in " 1 applied to a segment of length 3.

Cantilever beam
Next, we consider compliance minimization of a cantilever beam ( Figure 7) using a finite element mesh consisting of 160ˆ100 bilinear quadrilateral elements, minimum filtering length scale of R " 2, and maximum volume of V " 1 2 |Ω|. The optimal design is shown in Figure 7. All algorithms studied in this section are initialized from the same feasible design: ψ p0q " 0.45¨1.
Having established the superiority of the (adaptive) trust-region methods over the fixed-radius methods for the MBB beam in Section 5.1, we consider only the (adaptive) trust-region methods for the cantilever beam ( Figure 8, Table 2). The HDM-MMA algorithm requires 31 iterations to converge to a relative objective function tolerance of " 0.01, and significantly more iterations to converge to tighter tolerances; e.g., 298 iterations for " 0.001. In contrast, ROM-TR-DIST with an initial trust-region radius multiplier of τ " 0.1 requires only 11 HDM evaluations and 48 ROM evaluations to converge to " 0.01. As the cost for all ROM evaluations is less than a single HDM evaluation, the overall cost is reduced by a factor of nearly three. Even larger reductions are observed for tighter tolerances; e.g., to achieve " 0.001, ROM-TR-DIST requires 59 HDM evaluations and 971 ROM evaluations (cost of « 10 HDM evaluations) and achieves an overall computational reduction of nearly four. Figure 9 shows the filtered density fields for the ROM-based trust-region method (ROM-TR-RES) and   HDM-MMA at select major iterations. The comparison shows that ROM-TR-RES obtains an excellent approximation to the optimal design after only 14 iterations; the HDM-MMA design contains suboptimal artifacts even after 30 iterations. Furthermore, after only 10 iterations, the ROM algorithm has essentially discovered the optimal topology and uses the next several iterations to refine it; the internal structure has not emerged at this point from the HDM-MMA method.

Simply supported beam
Finally, we consider compliance minimization of a simply supported beam ( Figure 10) using a finite element mesh consisting of 180ˆ90 bilinear quadrilateral elements, minimum length scale of R " 0.5, and maximum volume of V " 2 5 |Ω|. The optimal design is shown in Figure 10. All algorithms studied in this section are initialized from the same feasible design: ψ p0q " 0.36¨1.
We consider the performance of the ROM-based trust-region methods relative to the HDM-MMA algorithm ( Figure 11, Table 3). Unlike the previous problems, HDM-MMA converges to relatively tight tolerances rather quickly, e.g., 52 iterations for " 0.01 and 74 iterations for " 0.001. However, the ROM-based trust-region methods are faster; the relative speedup is similar to the previous problems. which translates to a speedup of greater than three. For the tighter tolerance " 0.001, ROM-TR-DIST only requires 19 HDM evaluations and 101 ROM evaluations, which translates to a speedup of nearly four. Figure 12 shows the filtered density fields for the ROM-based trust-region method (ROM-TR-RES) and HDM-MMA at select major iterations. The comparison reveals ROM-TR-RES obtains an excellent approximation to the optimal design after only 10 iterations; the HDM-MMA design contains suboptimal artifacts even after 30 iterations.

Conclusion
The contributions of this work are twofold: (1) the introduction and analysis of an generalized trustregion method for nonlinear optimization problems with convex constraints and (2) a class of efficient, globally convergent topology optimization methods based on projection-based ROMs and generalized trust regions. We consider two different types of trust regions for the reduced topology optimization problem: the traditional ball (two-norm) in the design space and the sublevel sets of the norm of the HDM residual evaluated at the ROM solution. We provide detailed analysis to show the trust-region method that uses the ROM as the trust-region model and sublevel sets of the HDM residual norm as the trust-region constraint is globally convergent. To ensure the ROM satisfies the accuracy constraints of the trust-region theory, we sample the HDM solution and adjoint at each trust-region center, which guarantees the ROM objective function and gradient match those of the HDM at the trust-region center. This leads to a ROM that is constructed on-the-fly and specialized for the solutions obtained on the optimization path. We apply the ROM-accelerated methods to three benchmark topology optimization problems; the methods converge to the optimal design at least two times faster than the standard MMA method that does not incorporate a ROM.
Having demonstrated the potential of the ROM-based trust-region method for topology optimization on several benchmark problems, there are a number of important avenues for future research. While the ROM evaluations in this work are much cheaper than HDM evaluations, the cost still scale with the large dimension (N e ) due to the Helmholtz filter. The ROMs could be further accelerated by constructing a separate reduced-order model for the Helmholtz system, which would also require an adaptation strategy to ensure global convergence. In addition, the cost could be further reduced by relaxing the requirement that the value and gradient of the reduced objective function be exact at trust-region centers. This would require the development of a more general error-aware trust-region theory based on the work in [20,21,42]. Another interesting avenue of research is to extend the reduced topology optimization approach to more complex topology optimization problems including non-compliance optimization, nonlinear and multiphysics problems, and problems with solution-dependent and non-convex constraints.
Appendix A. Proof of global convergence of error-aware trust-region method The proofs in this section directly follow those setforth in [20].
Proof of Theorem 6. Let We will show that (i) D k Ă B k , (ii) there exists x P D k that satisfies the fraction of (generalized) Cauchy decrease condition, and hence (iii) there exists x P B k that satisfies the condition. First we show that D k Ă B k . Take any x P D k , which can be written as x " x k`pk for }p k } 2 ď δ k . Then ϑ k pxq " ϑ k px k`pk q " ϑ k px k q`∇ϑ k pξq T p k ďˇˇ∇ϑ k pξq T p kˇď κ ∇ϑ δ k ď ∆ k (A.2) for some ξ " tx`p1´tqx k , t P r0, 1s; here, the second equality follows from the continuous differentiability of ϑ and the mean value theorem, the first inequality follows from Assumption 4.3 (ϑ k px k q " 0), the second inequality follows from the Cauchy-Schwartz inequality, and the last inequality follows from the definition of δ k . Since x P C and ϑ k pxq ď ∆ k , we have x P B k , which establishes the inclusion D k Ă B k . Next we recall a standard result from trust-region theory (Theorem 12.2.2 of [13]) that establishes the existence of a point x P D k such that for some κ P p0, 1q. By D k Ă B k and the definition of δ k , there exists x P B k such that m k px k q´m k pxq ě κχpx k q min such that χpx k q ě for all k ą K. Then the sequence of trust-region radii t∆ k u produced by Algorithm 1 satisfies Proof. We first consider the case with a finite number of successful iterations. In this case, there exists K 1 ą 0 such that all iterations k ą K 1 are unsuccessful. Then, above by a geometric series, implying the infinite sum is finite. Therefore, the result holds if there are a finite number of successful iterations. We now consider the case with an infinite sequence of successful iterations tk i u. In this case, for all i such that k i ą K, F px ki q´F px ki`1 q ě F px ki q´F px ki`1 q " F px ki q´F px ki q ě η 1 pm ki px ki q´m ki px ki qq for some constant κ P p0, 1q. (Note that the subscript for the second x in the first and second expressions are k i`1 and k i`1 , respectively.) The first inequality holds because the sequence tF px k qu is non-increasing owing to the step acceptance condition in Algorithm 1, and the first equality holds because iteration k i successful (by construction). The remaining inequalities follow from the step acceptance condition in Algorithm 1, the the fraction of Cauchy decrease (37), and the assumption that χpx ki q ě for all k i ą K. Summing over all i ą I for k I " K, where the equality follows from the telescoping series, and the finiteness of the limit follows from F being bounded below. Since {β and 1 are bounded away from zero, the inequality above implies that ř 8 i"1 ∆ ki ă 8.
Let S Ă N be the ordered set of indices of successful iterations. For every k R S, ∆ k ď γ k´jpkq 2 ∆ jpkq , where jpkq P S is the largest index such that jpkq ă k, i.e. jpkq represents the last successful iteration before the unsuccessful iteration k. Summing over all k R S, This proves the desired result. note that the second derivative vanishes because bpψq is linear in ψ. The twice differentiability of (iii) j : R N u uˆP Ñ R follows from the assumption of the lemma. We next consider (C2). It suffices to show that J : C Ñ R is continuous and C is compact. The function J is continuous because it is twice continuously differentiable as proven for (C1). The domain C is compact by Lemma 5. We finally note that (C3) follows from (C1) and the fact C is compact by Lemma 5.
Lemma 5 (Verification of Assumption 2 for topology optimization.). The constraints of the topology optimization problem (11) satisfies Assumption 2.
Proof. We wish to show the following: (C1) C " X m i"1 C i , where C i " tx P R N | c i pxq ě 0u and each c i : R N Ñ R is twice continuously differentiable on R N ; (C2) C is nonempty, closed, convex; (C3) a firstorder constraint qualification holds at any critical point of (32). The condition (C1) is satisfied because the constraint set (42) is the intersection of linear constraints, which comprise simple bounds and a single general linear constraint. The condition (C2) is satisfied because (i) C is convex and closed because it is the intersection of linear and non-strict inequality constraints and (ii) C is nonempty because ρ l 1 P C since V ě ř Ne e"1 ρ l |Ω e |. Finally, the condition (C3) is satisfied since all constraints are linear functions, a first-order constraint qualification holds (Lemma 12.7 of [25]).
Lemma 6 (Verification of Assumption 3 for topology optimization.). Suppose the objective function j : R N u uˆP Ñ R is twice continuously differentiable. The approximation model introduced in Section 4.2 for the topology optimization problem (11) satisfies Assumption 3.
Proof. We wish to show that the ROM approximation of the objective function J k : R Ne Ñ R (C1) is twice continuously differentiable, (C2) satisfies J k pψ pkq q " Jpψ pkq q where ψ pkq is a trust-region center, (C3) satisfies ∇J k pψ pkq q " ∇Jpψ pkq q where ψ pkq is a trust-region center, and (C4) yields β k :" 1`max ψPC }∇ 2 J k pψq} that is uniformly bounded from the above.
We first note that the conditions (C2) and (C3) are satisfied as a consequence of Theorems 3 and 4, respectively, and our choice of reduced basis which includes the state at the trust-region center.
We next consider condition (C1). The verification of this condition is similar to the verification of condition (C1) of Lemma 4. Given J k pψq " jpuk pρ ‹ pψqq, ρ ‹ pψqq, we need to show each of the following is twice continuously differentiable: (i) uk : P Ñ R dNv ; (ii) ρ ‹ : R Ne Ñ P ; (iii) j : R dNvˆP Ñ R. We have shown (ii) and (iii) are satisfied in the verification of Lemma 4. To analyze the differentiability of uk : P Ñ R dNv , we differentiateK k pρqûk pρq "f twice and appeal to uk pρq " Φ kûk pρq to obtain for p, q " 1, . . . , N e , where all terms are evaluated at ρ. BecauseKpρq :" Φ T k KpρqΦ k , Kpρq is symmetric positive definite for all ρ P P , and Φ k P R dNvˆk has orthonormal columns, the maximum singular value of Kpρq´1 satisfies σ max pKpρq´1q ď σ max pKpρq´1q ă 8. In addition,Kpρq is polynomial in ρ and hence is smooth in ρ. Consequently, BρpBρq is continuous. Because (ii) and (iii) are also twice continuously differentiable as shown for condition (C1) of Lemma 4, it follows that J k pψq is twice continuously differentiable for any k.
We finally consider condition (C4). It suffices to show that t}∇ 2 J k pψq}u kPN,ψPC is uniformly bounded from the above. The boundedness follows from the twice continuous differentiability of J k pψq (i.e., (C1)) and the compactness of C. To see that this bound is uniform, we first note that J k pψq " jpuk pρ ‹ pψqq, ρ ‹ pψqq and the only term that depends on k is uk . We next note thatKpρq :" Φ T k KpρqΦ k , Kpρq is symmetric positive definite for all ρ P P , and Φ k P R dNvˆk has orthonormal columns. It follows that Proof. We wish to show that, for ϑ k pψq :" }rpuk pρ ‹ pψqq; ρ ‹ pψqq} 1´ and ą 0 but ! 1, (C1) ϑ k p¨q is twice continuously differentiable, (C2) max ψPC }∇ϑ k pψq} is uniformly bounded from the above, (C3) ϑ k pψ pkq q " 0 for the trust-region center ψ pkq , and (C4) there exist ζ ą 0 and ν ą 1 independent of k such that |Jpψq´J k pψq| ď ζϑ k pψq ν for all k P N. We first consider (C1). Given (C1) of Lemma 6 is satisfied, we only need to show that r : R dNvˆP Ñ R dNv is twice continuously differentiable. We recall that rpu, ρq " Kpρqu´f . We observe that rpu, ρq is twice continuously differentiable in u and ρ because rpu, ρq is linear in u and Kpρq is polynomial in ρ.
We next note that (C2) is satisfied because rpuk pρ ‹ pψqq; ρ ‹ pψqq is twice continuously differentiable in ψ and C is compact. The uniform boundedness follows from the same argument as Lemma 6 (C4).
We then note that (C3) is satisfied because u ‹ k pψ pkq q " upψ pkq q for any trust-region center ψ pkq as a consequence of Theorem 1 and our choice of reduced basis which includes the state at the trust-region center.
It remains to show that ζ is bounded independent of k. To this end, we first note that max ψPB k }λ ‹ pψq} 2 is bounded because the adjoint problem (21) is well-posed for all ψ P Ψ and B k Ă Ψ is compact. We second note that max ψPB k σ max pBpψqq∆ max is bounded independent of k because all terms of Bpψq in (B.2) are bounded: }Kpρ ‹ pψqq´1} 2 is bounded because Kpρq is nonsingular for all ρ P B k Ă P and B k is compact; the term in the square bracket is bounded because (i) j is twice continuously differentiable by assumption, (ii) u ‹ pψq and u ‹ k pψq are continuous in ψ P B k Ă P , and (iii) ψ belongs to a compact set B k . We note that for " 0 the condition (C4) is not satisfied but conditions (C1)-(C3) are satisfied. that, by (22),