C Soize

C Farhat

Probabilistic learning for modeling and quantifying model-form uncertainties in nonlinear computational mechanics

Keywords: probabilistic learning, model-form uncertainties, nonparametric probabilistic method, model reduction, uncertainty quantification, machine learning

Recently, a novel, nonparametric, probabilistic method for modeling and quantifying model-form uncertainties in nonlinear computational mechanics was proposed. Its potential was demonstrated through several uncertainty quantification (UQ) applications in vibration analysis and nonlinear computational structural dynamics. This method, which relies on projection-based model order reduction in order to achieve computational feasibility, exhibits a vector-valued hyperparameter in the probability model of the random reduced-order basis and associated stochastic, projection-based reduced-order model. It identifies this hyperparameter by formulating a statistical inverse problem grounded in target quantities of interest and solving the corresponding nonconvex optimization problem. For many practical applications however, this identification approach is computationally intensive. For this reason, this paper presents a faster, predictor-corrector approach for determining the appropriate value of the vector-valued hyperparameter that is based on a probabilistic learning on manifolds. It also demonstrates the computational advantages of this alternative identification approach through the UQ of two three-dimensional, nonlinear, structural dynamics problems associated with two different configurations of a MEMS device.

Notation, nomenclature and acronyms

Throughout this paper:

A real, deterministic variable is denoted by a lower case letter such as y. A real, deterministic vector is denoted by a boldface, lower case letter such as in y = (y 1 , . . . , y N ). A real, random variable is denoted by an upper case letter such as Y . A real, random vector is denoted by a boldface, upper case letter such as in Y = (Y 1 , . . . , Y N ). A real, deterministic matrix is denoted by an upper (or lower) case letter between brackets such as [A] (or [a]). A real, random matrix is denoted by a boldface, upper case letter between brackets such as [A]. y denotes the Euclidean norm of vector y. E denotes the mathematical expectation. M N,n denotes the set of N ×n real matrices. M n denotes the set of square n×n real matrices.

INTRODUCTION

In general, once a computational model must be considered for any purpose, model-form uncertainties become unavoidable. They can result either from the lack of knowledge of the true physics underlying the problem of interest, or the omission or truncation of modeling details. Their raisons d'être range from the lack of available information, to the inability to discern in some circumstances between important and unimportant modeling details. Such uncertainties affect the ability of even an HDM to deliver predictive results -that is, results that match reasonably well their experimental counterparts (which themselves can also be tainted by other types of uncertainties). Examples in computational structural dynamics include constitutive, multiscale, friction, homogenization, and free-play modeling errors. A parametric PROM constructed for accelerating parametric studies or stochastic computations meant to be carried out using an HDM inherits all model-form uncertainties associated with this HDM, because it is obtained by projecting this same HDM onto a carefully computed ROB. Such a reduced-order model is also tainted by additional uncertainties due to, among other factors: the finite sampling used during the offline training of this PROM where the ROB is computed (for example, see [START_REF] Paul-Dubois-Taine | An adaptive and efficient greedy procedure for the optimal training of parametric reduced-order models[END_REF]); the truncation of the subspace of approximation represented by this ROB using POD or SVD (for example, see [START_REF] Xu | Truncated SVD methods for discrete linear ill-posed problems[END_REF]); and the approximation errors induced by hyperreduction when this process is used to accelerate the projections incurred by the construction of this PROM (for example, see [START_REF] Farhat | Structure-preserving, stability, and accuracy properties of the Energy-Conserving Sampling and Weighting (ECSW) method for the hyper reduction of nonlinear finite element dynamic models[END_REF]).

In [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and highdimensional nonlinear models[END_REF], a nonparametric, probabilistic method was developed to model and quantify modelform uncertainties in deterministic computational models. This method recognizes that data, () Prepared using nmeauth.cls whether experimental or numerical and high-dimensional, may contain fundamental information or knowledge that is not captured by deterministic computational models. Therefore, it exploits available data to adapt the subspace in which the solution of the problem formulated using the computational model is searched. The method is innovative as it is a significant departure from the dominant parametric approach to UQ where model parameters are adapted instead, in order to discover some form of the information or knowledge encapsulated in data. To this end, the nonparametric, probabilistic method developed in [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and highdimensional nonlinear models[END_REF] parameterizes the solution subspace by randomizing its representative basis using a vector-valued hyperparameter. This transforms the deterministic computational model into a hyperparametric stochastic one. In order to determine the hyperparameter, the method constructs and solves an inverse statistical problem designed to minimize the discrepancies between the mean values and statistical fluctuations of the QoIs predicted using the hyperparametric stochastic model, and target values obtained from the available data. If the computational model is a deterministic HDM, the method implements all of the above ideas into an SPROM in order to make them computationally feasible. It follows that the method developed in [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and highdimensional nonlinear models[END_REF] for modeling and quantifying model-form uncertainties essentially extracts the information missing in a given deterministic computational model from data, and infuses this information into a stochastic, lower-dimensional counterpart. For this reason, this method has also been described in [START_REF] Farhat | Modeling and quantification of model-form uncertainties in eigenvalue computations using a stochastic reduced model[END_REF] as a data-driven, probabilistic, model adaptation method.

As mentioned above, a PROM constructed using a Galerkin [START_REF] Carlberg | Efficient nonlinear model reduction via a least-squares Petrov-Galerkin projection and compressive tensor approximations[END_REF] or Petrov-Galerkin [START_REF] Carlberg | The GNAT method for nonlinear model reduction: effective implementation and application to computational fluid dynamics and turbulent flows[END_REF] projection of an HDM inherits all uncertainties present in this HDM. For this reason, and because it is low-dimensional, such a reduced-order model can be used to effectively model and quantify the uncertainties of the HDM of interest. More importantly however -and unlike most if not all other UQ methods proposed in the literature that mention or rely on model reduction for accelerating the stochastic computations they incur -the nonparametric, probabilistic method developed in [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and highdimensional nonlinear models[END_REF] accounts for the modeling errors introduced by model reduction to construct its underlying SPROM. If experimental data are not available, it leads to the modeling and quantification of the modeling errors due to model reduction only. Otherwise, it leads to the modeling and quantification of the modeling errors of the given HDM, while accounting for those introduced during the process by model reduction.

The potential of this nonparametric, probabilistic method for the UQ of model-form uncertainties in nonlinear structural dynamics problems was demonstrated in [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and highdimensional nonlinear models[END_REF]. The method was extended to model-form uncertainties arising in generalized eigenvalue computations in [START_REF] Farhat | Modeling and quantification of model-form uncertainties in eigenvalue computations using a stochastic reduced model[END_REF], where its performance for such applications was demonstrated through the vibration analysis of two different real-world aircaft for which experimental data is available [START_REF] Gupta | Ground vibration tests on a flexible flying wing aircraft[END_REF]. In [START_REF] Farhat | A stochastic projection-based hyperreduced order model for model-form uncertainties in vibration analysis[END_REF], the computational aspects of this method were refined by hyperreducing the SPROM to accelerate its execution time. The potential of the resulting enhanced approach for quantifying model-form uncertainties in eigenvalue computations, while accounting for modeling errors due to hyperreduction, was demonstrated through what-if vibration analysis scenarios. These were associated with shape changes of a supersonic jet engine nozzle with inner walls, outer walls, stringers, and baffles [START_REF] Farhat | A stochastic projection-based hyperreduced order model for model-form uncertainties in vibration analysis[END_REF].

Given an HDM of interest of dimension N , the family of nonparametric, probabilistic methods for UQ described in [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and highdimensional nonlinear models[END_REF], [START_REF] Farhat | Modeling and quantification of model-form uncertainties in eigenvalue computations using a stochastic reduced model[END_REF], and [START_REF] Farhat | A stochastic projection-based hyperreduced order model for model-form uncertainties in vibration analysis[END_REF] consists in substituting the ROB used for constructing a nonlinear PROM of dimension n ≪ N by an adapted SROB, which leads to the construction of a nonlinear SPROM of dimension n. The probability model of the SROB depends on the aforementioned vector-valued hyperparameter, α, which belongs to an admissible set C α ⊂ R mα , where m α = 2 + n(n + 1)/2. This family of UQ methods performs the identification of α by solving a nonconvex optimization problem formulated using an objective (cost) function related to a certain distance between some given targets for the QoIs, and the corresponding predictions carried out using the nonlinear SPROM. Again, the target values for the QoIs can be extracted from experimental data when these are available, or corresponding HDM predictions when experimental data are not available. If n is greater than a few units, m α = n(n + 1)/2 is sufficiently large to challenge the fast solution of this optimization problem and therefore to compromise the execution in near-real-time of this family of nonparametric, probabilistic UQ methods. To this end, this paper presents an alternative, more economical approach for identifying the vector-valued hyperparameter () Prepared using nmeauth.cls α for practical values of the low-dimension n of the SPROM. First, this approach constructs a low-dimensional predictor α 0 for α of dimension m 0 α ≪ 2 + n(n + 1)/2 whose computational cost is relatively small. If it judges the predictor to be insufficiently accurate, it constructs next a corrector α opt of the same dimension using a probabilistic learning on manifolds (PLM) based on ideas that are similar to those recently described in [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF][START_REF] Soize | Polynomial chaos representation of databases on manifolds[END_REF] and demonstrated in [START_REF] Ghanem | Probabilistic nonconvex constrained optimization with fixed number of function evaluations[END_REF]. By avoiding to call a large number of times the nonlinear SPROM, this alternative approach achieves a lower computational complexity than the random search optimizer originally proposed in [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and highdimensional nonlinear models[END_REF] for identifying the hyperparameter α.

The remainder of this paper is organized as follows. Section 2 sets the context of this paper to that of nonlinear computational structural dynamics problems in the time domain. It also overviews the nonparametric probabilistic method proposed in [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and highdimensional nonlinear models[END_REF] for modeling and quantifying model-form uncertainties for two purposes: to keep this paper as self-contained as possible; and more importantly, to highlight the presence, role, and significance of the hyperparameter α of this method. Section 3 briefly recalls the nonconvex optimization problem formulated in [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and highdimensional nonlinear models[END_REF] for determining this vector-valued hyperparameter and outlines the potential difficulties associated with the solution of this problem using a conventional approach. Section 4 introduces the main contribution of this paper. It starts by describing the construction of a low-dimensional predictor α 0 for α that is based on the Fourier transform of the stochastic solution of the nonlinear problem of interest. This stochastic solution is computed in near-real-time using the SPROM. Next, Section 4 describes the construction of the proposed corrector based on the aforementioned probabilistic learning on manifolds. Section 5 discusses the performance of the nonparametric probabilistic method for UQ equipped with the proposed predictor-corrector approach for two 3D, nonlinear, structural dynamics problems associated with two different configurations of a MEMS device. Finally, Section 6 concludes this paper. REMARK 1. PROMs are usually designed to accelerate the repeated processing of parametric HDMs known as µ-HDMs. In this context, µ denotes a vector-valued parameter representing, for example, a parameteric shape, material property, boundary condition, or source term. Consequently, the PROMs and SPROMs mentioned above have been referred to in [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and highdimensional nonlinear models[END_REF][START_REF] Farhat | Modeling and quantification of model-form uncertainties in eigenvalue computations using a stochastic reduced model[END_REF][START_REF] Farhat | A stochastic projection-based hyperreduced order model for model-form uncertainties in vibration analysis[END_REF] as µ-PROMs and µ-SPROMs, respectively. Here, it is noted however that the vector-valued parameter µ is unrelated to the vector-valued hyperparameter α discussed above, and that the contribution of this paper is independent from and agnostic about how the parametric aspect of a µ-PROM is treated. The problem discussed in this paper and the proposed solution to this problem are also independent of the concept of hyperreduction and agnostic about how such an approximation is performed. For these reasons, and in order to avoid any unnecessary distractions, the contribution of this paper is presented without any reference to, or discussion of, the parametric and hyperreduction aspects of a PROM or SPROM. For a demonstration of how all these concepts perform well together, the reader is referred to [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and highdimensional nonlinear models[END_REF][START_REF] Farhat | Modeling and quantification of model-form uncertainties in eigenvalue computations using a stochastic reduced model[END_REF][START_REF] Farhat | A stochastic projection-based hyperreduced order model for model-form uncertainties in vibration analysis[END_REF].

NONPARAMETRIC PROBABILISTIC METHOD FOR MODEL-FORM UNCERTAINTIES IN NONLINEAR COMPUTATIONAL MECHANICS

In principle, the family of nonparametric, probabilistic methods for UQ described in [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and highdimensional nonlinear models[END_REF], [START_REF] Farhat | Modeling and quantification of model-form uncertainties in eigenvalue computations using a stochastic reduced model[END_REF], and [START_REF] Farhat | A stochastic projection-based hyperreduced order model for model-form uncertainties in vibration analysis[END_REF] can be extended to any linear or nonlinear computational problem for which a PROM can be constructed. However, the details associated with any extension can be problem dependent. For this reason, and for the sake of clarity and simplicity, the focus of this entire paper is set here on nonlinear transient dynamics problems.

Nonlinear high-dimensional computational model

After semi-discretization by the FE method, a large class of nonlinear transient dynamics problems can be written as

[M ] ÿ(t) + g(y(t), ẏ(t)) = f(t) , t ∈ ]t 0 , T ] , (2.1)
() Prepared using nmeauth.cls where y(t) = (y 1 (t), . . . , y N (t)) is the R N vector of the N displacement dofs before the application of the boundary conditions, ẏ(t) = dy(t)/dt and ÿ(t) = d 2 y(t)/dt 2 are the corresponding velocity and acceleration vectors, [M ] is the FE mass matrix belonging to M + N and is assumed here to be independent of t, g(y(t), ẏ(t)) is the R N vector of internal forces at time t and is assumed here to depend on y(t) and ẏ(t), f(t) is the R N vector of external forces at time t, and finally, t 0 and T define the time interval of interest and satisfy -∞ < t 0 < T < +∞.

The above equation defines an HDM whose initial conditions can be written as

y(t 0 ) = y 0 , ẏ(t 0 ) = y 1 , (2.2)
and whose boundary conditions and any other governing linear constraints can be written altogether as

[B] T y(t) = 0 NCD , t ∈ [t 0 , T ] . (2.3) Here, [B] is a time independent matrix in M N,NCD satisfying [B] T [B] = [I NCD ], N CD < N
, and y 0 and y 1 are two given vectors in R N that satisfy the constraint (2.3). The QoIs (also called observations) at time t associated with the above HDM are represented by the vector o(t) = (o 1 (t), . . . , o mo (t)) with values in R mo . This vector is written here as

o(t) = h time (y(t)) , t ∈ [t 0 , T ] , (2.4)
where h time is a given mapping with values in R mo (for example, 2 , where

h time (y(t)) = [O] d 2 y(t)/dt
[O] is a given constant matrix in M mo,N ).

Construction of a nonlinear projection-based reduced-order model

Let [V ] ∈ M N,n be a ROB of dimension n ≪ N constructed for approximating the solution {y(t), t ∈ [t 0 , T ]}, and satisfying the constraint equation

[B] T [V ] = [0 NCD,n ] .
(2.5)

In [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and highdimensional nonlinear models[END_REF], [V ] was furthermore constructed to satisfy the orthonormality condition

[V ] T [M ] [V ] = [I n ].
In this work, [V ] is constructed to satisfy instead the orthonormality condition

[V ] T [V ] = [I n ] . (2.6)
The Galerkin projection of the HDM defined in Equations (2.1) to (2.4) onto the ROB [V ] yields the following nonlinear PROM

y (n) (t) = [V ] q(t) , t ∈ [t 0 , T ] , (2.7)
[V ] T [M ] [V ] q(t) + [V ] T g([V ] q(t), [V ] q(t)) = [V ] T f(t) , t ∈ ]t 0 , T ] , (2.8)
whose initial conditions are

q(t 0 ) = [V ] T y 0 , q(t 0 ) = [V ] T y 1 , (2.9)
where {y (n) (t), t ∈ [t 0 , T ]} is the n-order approximation of {y(t), t ∈ [t 0 , T ]}. Concerning the QoIs, the corresponding approximation o (n) of o is given by

o (n) (t) = h time (y (n) (t)) , t ∈ [t 0 , T ] . (2.10)
For a given n ≪ N , the approximation error associated with the prediction of the QoIs using the nonlinear PROM instead of the nonlinear HDM can be estimated a priori as

ε(n) = T t0 o(t) -o (n) (t) 2 dt .
() Prepared using nmeauth.cls

Construction of a nonlinear stochastic projection-based reduced-order model

The details of the construction of the SPROM associated with the HDM defined in Equations (2.1) to (2.4) can be found in [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and highdimensional nonlinear models[END_REF]. Here, a short summary of this construction is adapted to the othonormalization condition (2.6) of the underlying ROB.

The nonparametric probabilistic method for the modeling and quantification of model-form uncertainties described here consists in substituting the deterministic ROB [V ] by a stochastic counterpart [W], referred to as the SROB. In view of the constraint Equation (2.5) and the orthonormality condition (2.6), this SROB must verify the following properties:

• [W] : θ → [W(θ)
] is a random variable (random matrix), defined on a probability space (Θ, T , P), with values in M N,n . • The support of its probability distribution, P [W] , is in the subset of M N,n defined by the following constraints which are verified almost surely

[B] T [W] = [0 NCD,n ] ,
(2.11) Consequently, y (n) , q, and o (n) are transformed into the stochastic processes Y (n) , Q, and O (n) , and the nonlinear SPROM is given by

[W] T [W] = [I n ] . ( 2
Y (n) (t) = [W] Q(t) , t ∈ [t 0 , T ] , (2.13)
[W] T [M ] [W] Q(t) + [W] T g([W] Q(t), [W] Q(t)) = [W] T f(t) , t ∈ ]t 0 , T ] , (2.14)
and the initial conditions

Q(t 0 ) = [W] T y 0 , Q(t 0 ) = [W] T y 1 .
(2.15)

The R n -valued stochastic solution (2.12), is a second-order random variable defined on a probability space (Θ, T , P). Specifically, it is built as follows:

[W] = ([V ] + s [Z]) [H s (Z)] ,
(2.17)

[H s (Z)] = ([I n ] + s 2 [Z] T [Z]) -1/2 , (2.18) [Z] = [A] -[V ] [D] ,
(2.19)

[D] = ([V ] T [A] + [A] T [V ])/2 , (2.20)
[A] = [U] -[B] {[B] T [U]} , (2.21) [U] = [G(β)] [σ] , (2.22)
() Prepared using nmeauth.cls where [G(β)] is a second-order, centered, M N,n -valued random matrix defined on (Θ, T , P) and [σ] is a given upper triangular matrix in M u n (positive or zeros diagonal entries). The probability distribution of [W] and a generator of this random matrix are described in Appendix D of [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and highdimensional nonlinear models[END_REF]. The m α = 2 + n(n + 1)/2 hyperparameters of the stochastic model of [W] are as follows:

-A deterministic, real parameter s satisfying 0 ≤ s L ≤ s ≤ s U ≤ 1, where the lower bound s L and upper bound s U are fixed; note that if s L = s = 0, then [W] = [V ] is deterministic and there are no statistical fluctuations. -A deterministic, real parameter β satisfying 0 < β L ≤ β ≤ β U < +∞, where the lower bound β L and upper bound β U are fixed. -An upper triangular matrix [σ] in M u n (positive or zero diagonal entries) parameterized by m σ = n(n + 1)/2 and characterized by the following property: for all

1 ≤ k ≤ k ′ ≤ n, [σ L ] kk ′ ≤ [σ] kk ′ ≤ [σ U ] kk ′ .
The lower bounds are the entries of the upper triangular matrix [σ L ]: they are fixed in M u n . The upper bounds are the entries of the upper triangular matrix [σ U ]: they are fixed in M u n . Hence, the hyperparameter of the nonparametric probabilistic method reviewed here is the

vector α = (s, β, {[σ] kk ′ , 1 ≤ k ≤ k ′ ≤ n}) of dimension m α = 2 + n(n + 1)/2.
It can be simply re-written as α = (s, β, [σ]). It belongs to the admissible set C α defined by the lower and upper bounds of every one of its components. For any fixed value of α in C α , the generator of the random matrix [G(β)] described in [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and highdimensional nonlinear models[END_REF] enables the computation of any realization [G(β; θ)] of [G(β)] for θ in Θ. It also enables the obtention of the corresponding realization [W(θ)] of the random matrix [W], using Equations (2.17) to (2.22).

Note that even though the SROB [W] depends on the hyperparameter α (see Equations (2.17) to (2.22)), this dependency is not highlighted in the remainder of this paper by writing [W(α)], for the sake of simplicity. REMARK 2. Regarding the interpretation of the individual hyperparameters constituting the vectorvalued hyperparameter α, the following is noteworthy:

-From Eqs. (2.17) to (2.22), it follows that the random matrix [W] depends on the hyperparameters s and [σ], and s is a scaling parameter of the matrix-valued hyperparameter

[σ].
-From the construction of the random matrix [G(β)] described in details in [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and highdimensional nonlinear models[END_REF], it follows that the hyperparameter β enables the control of the statistical correlations between the random components of a same column of the random matrix [G(β)].

-In [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and highdimensional nonlinear models[END_REF], [σ] was introduced as the upper triangular matrix resulting from the Cholesky factorization of an unknown but positive-definite square matrix [c n ] that participates to the construction of the correlation tensor of the random matrix The identification of the hyperparameter α in C α ⊂ R mα can be performed using various approaches for formulating a statistical inverse problem related to o target , including the maximum likelihood and nonlinear Least-Squares methods [START_REF] Spall | Introduction to Stochastic Search and Optimization[END_REF][START_REF] Serfling | Approximation Theorems of Mathematical Statistics[END_REF][START_REF] Soize | An Accelerated Course with Advanced Applications in Computational Engineering[END_REF]. For example, in the context of nonlinear transient dynamics problems, the vector-valued hyperparameter α was determined in [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and highdimensional nonlinear models[END_REF] by:

• Defining a total cost function of the form

J(α) = w J J mean (α) + (1 -w J ) J std (α) , (3.1)
where:

-J mean (α) and J std (α) are partial cost functions formulated using the deterministic, time-dependent functions {o target (t), t ∈ [t 0 , T ]} and {o (n) (t; α), t ∈ [t 0 , T ]}, and the stochastic, time-dependent functions {O (n) (t; α), t ∈ [t 0 , T ]}. These functions are designed to control the discrepancies between the mean values and standard deviations of O (n) (t) and o target (t), respectively. -w J is a real scalar statisfying 0 ≤ w J ≤ 1, so that the total cost function J is a convex combination of the partial cost functions J mean and J std . • Solving the associated minimization problem

α opt = min α∈Cα J(α) .
(3.2)

In [START_REF] Farhat | Modeling and quantification of model-form uncertainties in eigenvalue computations using a stochastic reduced model[END_REF] and [START_REF] Farhat | A stochastic projection-based hyperreduced order model for model-form uncertainties in vibration analysis[END_REF], this approach was extended to the context of eigenvalue computations where the observed targets are algebraic rather than time-dependent quantities, by tailoring the definitions of J mean (α) and J std (α) to this context. In any case, the reader is referred to [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and highdimensional nonlinear models[END_REF] for the detailed development of the expressions of these two partial cost functions in the case of applications in nonlinear structural dynamics, and to [START_REF] Farhat | Modeling and quantification of model-form uncertainties in eigenvalue computations using a stochastic reduced model[END_REF] and [START_REF] Farhat | A stochastic projection-based hyperreduced order model for model-form uncertainties in vibration analysis[END_REF] for that of applications in vibration analysis. These expressions are not repeated here because an alternative approach to (3.1) and (3.2) is proposed in this paper for identifying the hyperparameter α, for the main reason outlined below. Indeed, the optimization problem defined by (3.1) and (3.2) can be solved using either a deterministic method such as the interior point algorithm (with potential problems due to the non convexity of the cost function J(α)), or a probabilistic technique such as a genetic algorithm. For each candidate/iterate solution α = (s, β, [σ]) proposed by the optimizer, the evaluation of J(α) is performed in either case using a Monte Carlo simulation method equipped with ν sim independent realizations {[G(β; θ ℓ )], ℓ = 1, . . . , ν sim } of the random matrix [G(β)] as the stochastic solver for Equations (2.13) to (2.22), and the mean-square convergence with respect to ν sim is controlled. If the chosen optimizer incurs ν J evaluations of the cost function J(α), it incurs ν J ×ν sim repeated solutions of the governing forward problem using the nonlinear SPROM (2.13, 2.14). During the application of the nonparametric probabilistic method outlined above to the UQ of nonlinear transient dynamics problems [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and highdimensional nonlinear models[END_REF] and generalized eigenvalue problems [START_REF] Farhat | Modeling and quantification of model-form uncertainties in eigenvalue computations using a stochastic reduced model[END_REF][START_REF] Farhat | A stochastic projection-based hyperreduced order model for model-form uncertainties in vibration analysis[END_REF], it was found that as soon as the dimension of the SPROM n dictated by accuracy requirements becomes greater than a few units, m α becomes sufficiently large to induce a very large number of optimization iterations ν J . In this case, the numerical cost of the identification approach of the hyperparameter α defined by (3.1) and (3.2) becomes significant, if not prohibitive, particularly in the context of model reduction. For this reason, a faster alternative approach for identifying α is presented next.

PREDICTOR-CORRECTOR METHOD FOR THE ESTIMATION OF THE HYPERPARAMETER

The alternative approach presented here for identifying the vector-valued hyperparamter α of the probability model of the SROB operates on Fourier transforms of target values of time-dependent QoIs, and Fourier transforms of counterpart values predicted using the constructed nonlinear SPROM. The choice of working in the frequency domain instead of the time domain is made here () Prepared using nmeauth.cls because in general, the phase between a time-dependent target value of a QoI and its counterpart value predicted using a PROM is generally unknown, and that between a time-dependent target value of a QoI and its counterpart value predicted using an SPROM is random and therefore never known. Hence, it is easier to formulate the discrepancy (error) in the frequency domain by focusing on the modulus of the complex-valued Fourier transforms of measured and computed time domain responses.

The method presented here for identifying α is a predictor-corrector approach based on a probabilistic learning on manifolds. Its predictor component can be written as

α 0 = (s 0 , β 0 , [σ 0 ]) in C α , where [σ 0 ] is a predictor for the matrix-valued hyperparameter [σ] ∈ M u
n that is algebraically constructed through adapted approximations, β 0 = (β L + β U )/2, where the subscripts L and U designate here and throughout the remainder of this paper lower and upper bounds of the quantity they designate, and s 0 is identified by constructing an optimization problem and solving it using a fast grid search in algorithm.

If the predictor must be improved, the proposed corrector is based on solving another optimization problem where the cost function is constructed using the extreme values statistics of a positivevalued random variable Q that quantifies the error between the frequency-dependent target values of some QoIs and their frequency-dependent predicted values. The independent realizations of the extreme values statistics are extracted from a large number of additional realizations that are computed by a probabilistic learning on manifolds for Q without exercising the nonlinear SPROM. In this case, the optimizer is the interior point algorithm. The PLM requires the knowledge of an initial dataset of a relatively small number of points {α ℓ ′ d } ℓ ′ localized in the neighborhood of the predictor α 0 , and the corresponding realizations

{q ℓ ′ d } ℓ ′ of {Q(α ℓ ′ d )} ℓ ′ computed using the nonlinear SPROM.

Observation in the frequency domain and stochastic error

Let B o = [0 , ω o ]
denote a frequency band of observation that is sampled by m ω points using a constant frequency increment ∆ω, which implies that ω 0 = m ω ∆ω. For any vector-valued time function {z(t), t ∈ [t 0 , T ]}, the associated complex vector-valued frequency function

{ z(ω), ω ∈ B o } is defined as z(ω) = T t0
e -iω t z(t) dt. This function can be viewed as the restriction to B o of the Fourier transform of {z(t) t ∈ [t 0 , T ]}, if z(t) = 0 for all t in R\[t 0 , T ] (which will be the case, up to a small relative error, for the quantities to which this Fourier transform will be applied to in the remainder of this paper).

Let also o(ω) = ( o 1 (ω), . . . , o mo (ω)) ∈ C mo , with o(ω) = h freq (-ω 2 y(ω)), denote the complex vector of frequency-dependent QoIs, { o(ω), ω ∈ B o }, associated with the time-dependent QoIs, {o(t), t ∈ [t 0 , T ]}. The mapping h freq from C N into C mo is linear and represented by the matrix [O] in M mo,N . It extracts the observed dofs from the vector -ω 2 y(ω) ∈ C N . Such a definition of frequency-dependent QoIs is applied here to {o [START_REF] Soize | An Accelerated Course with Advanced Applications in Computational Engineering[END_REF]), and to the target data

(n) (t), t ∈ [t 0 , T ]} defined in Equation (2.10), the stochastic observation {O (n) (t, α), t ∈ [t 0 , T ]} defined in Equation (2.
{o target (t), t ∈ [t 0 , T ]} defined in Section 3. It follows that for all ω in B o , o (n) (ω) = -ω 2 [O] y (n) (ω)) , O (n) (ω; α) = -ω 2 [O] Y (n) (ω; α)) , o target (ω) = -ω 2 [O] y target (ω)) . Now, let {dB (n) (ω; α) = (dB (n) 1 (ω; α), . . . , dB (n) mo (ω; α)) , ω ∈ B o }
denote the stochastic process defined on (Θ, T , P) and indexed by B o , with values in R mo such that, for all j = 1, . . . , m o , dB

(n) j (ω; α) = log 10 (| O (n) j (ω; α)|) .
The following functions ω → db target (ω) = (db target 1 (ω), . . . , db target mo (ω)) and ω → db (n) (ω) = (db

(n) 1 (ω), . . . , db (n) mo (ω))
with values in R mo are similarly defined on B o -that is, such that for all j = 1, . . . , m o ,

db target j (ω) = log 10 (| o target j (ω)|) , db (n) j (ω) = log 10 (| o (n) j (ω)|) . () Prepared using nmeauth.cls
For a fixed value of the dimension n of the SPROM, the stochastic discrepancy (error) between the target values of the QoIs and their counterpart values predicted using the nonlinear SPROM is denoted by the following positive-valued random variable Q(α) = {Q(α; θ), θ ∈ Θ} defined on (Θ, T , P) as follows

Q(α; θ) = Bo db target (ω) -dB (n) (ω; α; θ) 2 µ dω , (4.1)
where, for any function ω → r(ω) = (r 1 (ω), . . . , r mo (ω)) defined on B o with values in R mo and for all ω ∈ B o , r(ω) µ is defined as

r(ω) µ = { mo j=1 (r j (ω)) 2 µ j (ω)} 1/2 , (4.2)
with

µ j (ω) = e j (ω) max ω∈Bo {e j (ω)} and e j (ω) = 10 db target j (ω) |db target j (ω) -db (n) j (ω)| . (4.3)
For a fixed hyperparameter α, the realization dB (n) (ω; α; θ) of positive-valued random variable dB (n) (ω; α) (for which ω given is fixed in B o ) is computed by solving the realization θ ∈ Θ of the nonlinear SPROM. Next, two real numbers q target and q error associated with the random variable Q(α) are defined as follows

q target = Bo db target (ω) 2 µ dω , q error = Bo db target (ω) -db (n) (ω) 2 µ dω . (4.4)
These two positive numbers are used in the remainder of this paper to normalize the cost function of an optimization problem devoted to the identification of the hyperparameter α.

REMARK 3. For the purpose of the work presented in this paper, the function µ(ω) = (µ 1 (ω), . . . , µ mo (ω)) defined above was designed by numerical experiments to privilege those frequencies in the band B o for which the error is the greatest and make it possible to increase the sensitivity of the random variable Q(α) to variations in the hyperparameter α.

Construction of the predictor

The construction of the predictor is performed in two stages as described below.

Stage 1. The first stage consists in constructing a predictor

[σ 0 ] for [σ].
Examining the algebraic structure of Equation (2.13) and Equations (2.17) to (2.22), the idea here is to estimate [σ 0 ] by using an approximation of the error between the generalized response obtained using the nonlinear SPROM and the projection of the target on the SROB. Let a (n) (ω) = -ω 2 q(ω) ∈ C n be the frequency-dependent generalized acceleration associated with the generalized coordinate q computed using the nonlinear SPROM (see Section 2.2) and let a target (ω) be its counterpart target. From Equation (2.6) and the least-squares inversion of Equation (2.7), it follows that

a (n) (ω) = [V ] T (-ω 2 y (n) (ω)) , a target (ω) = [V ] T (-ω 2 y target (ω)) . (4.5)
In order to focus the construction of the error on the QoIs, the computation of a (n) (ω) = (a

(n) 1 (ω), . . . , a (n)
n (ω)) and a target (ω) = (a target 1 (ω), . . . , a target n (ω) can be performed only for the QoIs rather than all dofs of the HDM. This is well adapted to the case where the target consists of experimental data for the QoIs. In any case, this error can be written as a error (ω) = (a error 1 (ω), . . . , a error n (ω)) where, for all ω in B o ,

a error k (ω) = | (|a target k (ω)| -|a (n) k (ω)|) | , k = 1, . . . , n . () Prepared using nmeauth.cls Define [g] = Bo a error (ω) a error (ω) T dω .
This positive matrix is assumed to be definite ([g] ∈ M + n ) and is scaled as follows to obtain the matrix

[g s ] ∈ M + n [g s ] = 1 g max [g] , g max = max k,k ′ { [g] kk ′ } . Let [ σ 0 f ]/ √
m ω be the upper triangular (n×n) real matrix resulting from the Cholesky factorization of the matrix [g s ], where the normalization constant, m ω , is the number of sampling points of the frequency band B o (see Section 4.1). For constructing [σ 0 ], the matrix [ σ 0 f ] is first transformed into a sparse upper matrix belonging to M u n in order to keep only its entries that have a major influence on controlling the error between the target and the generalized response computed using the nonlinear SPROM. This construction allows the reduction of the number of hyperparameters, m α , whose maximum value is n(n + 1)/2. It is carried out in two steps as follows. Let τ rel be a relative tolerance (for instance τ rel = 0.15), and τ = τ rel max k [ σ 0 f ] kk be the corresponding absolute tolerance. Define the upper triangular matrix

[σ 0 f ] in M u n as follows [σ 0 f ] kk ′ = 0 if [ σ 0 f ] kk ′ ≤ τ and [σ 0 f ] kk ′ = [ σ 0 f ] kk ′ if [ σ 0 f ] kk ′ > τ .
Then, a predictor [σ 0 ] for [σ] is obtained by keeping in each column k ′ of the matrix [σ 0 f ] only the entry corresponding to the maximum over

k of |[σ 0 f ] kk ′ |. For k ′ fixed in {1, . . . , n}, let k 0 (k ′ ) = arg max k {|[σ 0 f ] kk ′ |}. Thus, for 1 ≤ k ≤ k ′ ≤ n, the upper triangular matrix [σ 0 ] is defined as [σ 0 ] kk ′ = 0 if k = k 0 (k ′ ) and [σ 0 ] kk ′ = [σ 0 f ] kk ′ if k = k 0 (k ′ ) . (4.6)
Now, let m 0 σ denote the number of non-zero entries in the matrix

[σ 0 ]. Given that 1 ≤ m 0 σ ≤ n, it follows that m 0 α = 2 + m 0 σ ≤ 2 + n,
which highlights the substantial reduction of the number of hyperparameters in the matrix [σ] achieved by the above construction of a predictor for this matrix.

Stage 2. The second stage consists in constructing a predictor s 0 for s. Let β 0 = (β L + β U )/2 denote a fixed value of the hyperparameter β, and recall the predictor [σ 0 ] for [σ] defined in Equation (4.6). The predictor s 0 for the hyperparameter s ∈ [s L , s U ] ⊂ R + is proposed here as the solution of the following optimization problem

s 0 = arg min s∈[sL,sU ] J freq (s) , (4.7)
where

J freq (s) = w J q target mean Q (s) + 1 -w J γ J q target |std Q (s) -γ J q error | , (4.8) mean Q (s) = E{Q(s, β 0 , [σ 0 ])} , (4.9) std Q (s) = {E{Q(s, β 0 , [σ 0 ]) 2 } -mean Q (s) 2 } 1/2 , (4.10) α = (s, β 0 , [σ 0 ]) , (4.11)
w J is fixed in [0, 1] (for example, w J = 0.5), q target and q error are defined in Equation (4.4), and γ J is a constant fixed in [0, 1]. The value of γ J is application-dependent and must be estimated as explained in Section 5. The nonconvex optimization problem defined in (4.7) is solved using a grid search algorithm associated with a grid s L = s 1 < s 2 < . . . < s ms = s U constituted of m s sampling points of the interval [s L , s U ]. For a given s j , the mean value and standard deviation of the random variable Q(s j , β 0 , [σ 0 ]) are estimated using the classical statistical estimators and ν s independent realizations of the random variable Q(s j , β 0 , [σ 0 ]) computed by solving ν s realizations of the nonlinear SPROM. Therefore, the computational cost of the proposed construction of the predictor s 0 is dominated by that of m s ×ν s realizations of the nonlinear SPROM.

() Prepared using nmeauth.cls This completes the description of the construction of the predictor α 0 for the hyperparameter α, which can be written as

α 0 = (s 0 , β 0 , [σ 0 ]) ∈ R m 0 α with m 0 α = 2 + m 0 σ . (4.12) REMARK 4.
In Stage 2 of the construction of the predictor α 0 , the value of β is fixed to β 0 = (β L + β U )/2 and the hyperparameter s is determined by solving the optimization problem (4.7), where β 0 and [σ 0 ] are given. Alternatively, [σ] could have been fixed to some value [σ 0 ] and a different optimization problem could have been formulated to determine s and β. However, numerical experiments carried out for assessing the performance of the predictor α 0 have revealed that the random matrix [W] is not sensitive with respect to different values of β. For this reason, the first aforementioned approach consisting in fixing β to the midpoint value of the admissible interval

[β L , β U ],
where appropriate values of β L and β U are determined experimentally, was chosen in formulating Stage 2 of the construction of α 0 (note however that in the construction of the corrector α opt described in Section 4.4, β is not fixed to some value a priori).

Quality assessment of the predictor

The quality of the predictor α 0 = (s 0 , β 0 , [σ 0 ]) of dimension m 0 α can be assessed a posteriori, using the following empirical indicator that is adapted to the logarithmic scale,

I = I + + 1 I -, I ± = mo j=1 I ± j , (4.13)
where

I ± j = Bo |db target j (ω) -db ± j (ω)| dω Bo |db target j (ω) + db ± j (ω)| dω , (4.14)
and {db + j (ω), ω ∈ B o } and {db - j (ω), ω ∈ B o } are the upper and lower envelopes of the confidence region of the family of random variables {dB (n) j (ω, α 0 ), ω ∈ B o } for a probability level p c (for instance, p c = 0.98). This means that, for all ω in B o ,

Proba{dB (n) j (ω, α 0 ) ≤ db + j (ω)} ≥ p c , Proba{dB (n) j (ω, α 0 ) ≥ db - j (ω)} ≥ 1 -p c . (4.15)
The dimensionless indicator I defined in Equation (4.13) is a measure of a membership of the target

{db target j (ω, α 0 ), ω ∈ B o } to the confidence region {[db - j (ω), db + j (ω)], ω ∈ B o }.
The smaller the I, the better the prediction of the probabilistic model of model-form uncertainties to represent the target.

Construction of the corrector

The computation of the corrector α opt associated with the predictor α 0 for the hyperparameter α is based on:

• The following hypotheses:

-The hyperparameter α = (s, β, [σ]) estimated by the corrector is an update of the predictor value α 0 = (s 0 , β 0 , [σ 0 ]) defined in Equation (4.12); therefore, the dimension of α remains m 0 α = 2 + m 0 σ . -The m 0 σ non-zero entries of the upper triangular matrix [σ 0 ] are represented by the vector

σ 0 = (σ 0 1 , . . . , σ 0 m 0 σ
). -The corrector [σ] for [σ 0 ] has the same sparsity pattern as [σ 0 ] and therefore m 0 σ non-zero entries; these are similarly represented by the vector σ = (σ 1 , . . . , σ m 0 σ ). -The bounds of the hyperparameters s, β, and σ, which have to be specified for the corrector algorithm, are defined in terms of the known values s 0 , β 0 , and σ 0 of the predictor. These bounds enable the definition of the admissible set C 0 α of values for the hyperparameter α, where the superscript 0 emphasizes the dependence of C 0 α on α 0 .

() Prepared using nmeauth.cls

• The PLM for solving nonconvex optimization problems described in [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF], and whose advantages relative to alternatives are discussed in [START_REF] Soize | Polynomial chaos representation of databases on manifolds[END_REF] and [START_REF] Ghanem | Probabilistic nonconvex constrained optimization with fixed number of function evaluations[END_REF].

REMARK 5. The probabilistic model of [W] being an informative prior probabilistic model, some of its parameters can be chosen arbitrarily. Here, the sparsity pattern of [σ] is chosen a a priori to be the same as that of [σ 0 ] -which was designed in Section 4.2 for the sole purpose of decreasing the number of individual hyperparameters to identify. This design is relevant to the capacity of this probabilistic model to generate, for given observations, the most narrow confidence regions containing the targets.

4.4.1.

Lower bounds for the corrector algorithm. As already stated, the lower and upper bounds of the hyperparameters s, β, and σ, which depend on α 0 , define the admissible set

C 0 α ⊂ R m 0 α for candidate values of α.
The values of the constants appearing in the bounds presented below result from recent and previously performed numerical experiments [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and highdimensional nonlinear models[END_REF][START_REF] Farhat | Modeling and quantification of model-form uncertainties in eigenvalue computations using a stochastic reduced model[END_REF][START_REF] Farhat | A stochastic projection-based hyperreduced order model for model-form uncertainties in vibration analysis[END_REF]. These bounds are constructed as follows.

-The hyperparameter s belongs to [s L , s U ], where s L = s 0 (1 -0.5) and s U = s 0 (1 + 0.5). This interval is centered at the predictor s 0 and has been identified by numerical experiments.

Since the predictor is assumed to be sufficiently good, s should not deviate significantly from s 0 . -The hyperparameter β belongs to [β L , β U ], where β L = 0.01 and β U = 0.03. These proposed values for the bounds are consistent with the optimal values obtained in all previous applications. However, they are not critical and can be slightly modified.

-The lower bound σ L = (σ L 1 , . . . , σ L m 0 σ
) and the upper bound (for further details, see [START_REF] Soize | A nonparametric probabilistic approach for quantifying uncertainties in low-and highdimensional nonlinear models[END_REF]).

σ U = (σ U 1 , . . . , σ U m 0 σ ) of the hyperparameter σ = (σ 1 , . . . , σ m 0 σ ) are such that, for all k = 1, . . . , m 0 σ , σ k ∈ [σ L k , σ U k ]: -If σ k corresponds to a diagonal entry of [σ], the lower bound is σ L k = max{c dm ×σ 0 k , τ } and the upper bound is σ U k = c dM ×σ 0 k ,
Hence, the admissible set for α is defined as

C 0 α = {α = (s, β, σ) , s ∈ [s L , s U ] , η ∈ [β L , β U ] , σ ∈ m 0 σ k=1 [σ L k , σ U k ] } . ( 4
α opt = arg min α∈C 0 α J(α) , (4.17)
where C 0 α is defined in Equation (4.16). A candidate choice for the cost function J could be the same as that in Equation (4.8) -that is,

J(α) = w J q target mean Q (α) + 1 -w J γ J q target |std Q (α) -γ J q error | , (4.18)
where mean Q (α) and std Q (α) are the mean value and the standard deviation of the random variable Q(α) defined in Equation (4.1). However, this cost function is not the most appropriate here for the reasons given below.

() Prepared using nmeauth.cls

The random variable Q(α) is positive almost surely. If its standard deviation were to tend to zero, its mean value would also have to tend to zero by virtue of the Tchebychev inequality. Since in principle the target does not coincide with the response of the nonlinear SPROM, the confidence region of the prediction computed using the nonlinear SPROM does not have a zero width. This means that the standard deviation of Q(α) must not tend to zero and therefore its mean value cannot a priori be zero for the optimal value α opt . In other words, in order to reduce the mean value of Q(α), it is necessary to generate statistical fluctuations in the nonlinear SPROM. This is the reason why, thanks to the presence of its second term, the cost function (4.8) (or (4.18)) allows the minimization of the mean value while maintaining a significant standard deviation. Nevertheless, since Q(α) is a positive-valued random variable, solving the optimization (4.17) problem using the cost function defined in (4.18) consists in searching for the minimum of a positive-valued cost function. From a statistical viewpoint, the efficiency of the minimization process can be increased by using the minimum statistics of Q(α) instead of all its possible values. For this reason, the cost function J(α) defined in (4.18) is replaced here by

J(α) = w J q target mean Q min (α) + 1 -w J γ J q target |std Q min (α) -γ J q error | , (4.19)
where mean Q min (α) and std Q min (α) are the mean value and standard deviation of the random variable Q min = min{Q}, which depends on α and is defined as the minimum statistics for the random variable Q(α). The mean value and the standard deviation of Q min (α) can be written as

mean Q min (α) = E{Q min (α)} , std Q min (α) = {E{Q min (α) 2 } -mean Q min (α) 2 } 1/2 . (4.20)

Probabilistic learning on manifolds for estimating the cost function.

Recall from the technical arguments presented at the end of Section 3 that whichever algorithm is used to solve the optimization problem (4.17) equipped with a cost function such as that defined in Equations (4. [START_REF] Goldberg | Genetic Algorithms in Search. Optimization & Machine Learning[END_REF]) and (4.20), this algorithm is bound to experience technical difficulties that can be summarized as follows. For each value of α proposed by the optimizer, the mean value and the standard deviation of Q min (α) defined in Equation (4.20) should be estimated using a large number, ν sim , of independent realizations of the random variable Q(α). These are typically obtained using the Monte Carlo technique applied to the nonlinear SPROM. Hence, if the chosen optimization algorithm requires ν α evaluations of the cost function J(α), the nonlinear SPROM must be exercised ν α ×ν sim times.

If n is not sufficiently small, ν α can be expected to be large and the total solution time to be cost prohibitive.

On the other hand, the PLM described in this section will incur only ν d calls to the nonlinear SPROM, where ν d ≪ ν α ×ν sim . Therefore, it will drastically reduce the number of times this computational model is exercised and thereby drastically reduce the total computational cost associated with the identification of α. This method has two key components, namely:

• A PLM that enables the generation of ν ar additional realizations concentrated on a manifold identified using an initial set of ν d data points (with ν d ≪ ν ar ). These additional realizations are determined without performing any additional exploitation of the nonlinear SPROM. • A smoothing technique for estimating J(α g ) at any trial point α g in C 0 α generated for the optimization algorithm, using only the ν ar additional realizations computed from the initial dataset.

The method consists of 4 steps:

Step 1. The first step consists in constructing an initial dataset by exercising ν d times the nonlinear SPROM. To this end, consider a fixed number ν d of values of the hyperparameter α, α 1 d , . . . , α ν d d in C 0 α , obtained using a uniform random sampling of ν d points in the compact set C 0 α defined in Equation (4.16). For each ℓ ′ fixed in {1, . . . , ν d }, a realization θ ℓ ′ of the nonlinear SPROM is constructed with α = α ℓ ′ and excercised to obtain the realization

q ℓ ′ d = Q(α ℓ ′ ; θ ℓ ′ ) of the stochastic error. Then, ν d data points x 1 d , . . . , x ν d d in R m 0 α +1 are () Prepared using nmeauth.cls introduced such that, for ℓ ′ = 1, . . . , ν d , x ℓ ′ d = (α ℓ ′ d , q ℓ ′ d ) ∈ C 0 α ×R + ⊂ R m 0 α +1 .
Step 2. The second step consists in constructing a diffusion maps basis and in generating additional independent realizations. For this purpose, the following random variable is introduced

X = (A, Q) ,
defined on (Θ, T , P), with values in

C 0 α ×R + ⊂ R m 0 α +1 such that {x ℓ ′ d , ℓ ′ = 1, . . . , ν d } are ν d independent realizations of X.
The diffusion maps basis is constructed using only the initial dataset {x ℓ ′ d , ℓ ′ = 1, . . . , ν d }. Then, ν ar ≫ ν d additional independent realizations, {x ℓ ar , ℓ = 1, . . . , ν ar }, of the random vector X are generated using the PLM without using the nonlinear SPROM (for further details, the reader is referred to [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF] as the details of this part of the method are similar to those described in this reference). Once these additional realizations have been computed, the ν ar additional independent realizations {α ℓ ar , ℓ = 1, . . . , ν ar } of the random variable A and ν ar additional independent realizations {q ℓ ar , ℓ = 1, . . . , ν ar } of the random variable Q are obtained by performing the following extraction (α ℓ ar , q ℓ ar ) = x ℓ ar , ℓ = 1, . . . , ν ar .

If ν d is sufficiently large, the learning can be expected to be successful and the realizations {α ℓ ar } ℓ and {q ℓ ar } ℓ will belong to C 0 α and R + , respectively.

Step 3. The third step pertains to the computation of independent realizations of the minimum statistics, represented by the positive-valued random variable Q min on (Θ, T , P), for the random variable Q. Let ν r and ν e be two integers such that ν ar = ν e ×ν r . Consequently, ν r = ν ar /ν e . It is assumed here that ν e and ν r are sufficiently large -say, ν e = 100, ν r = 10 000, and then ν ar = 1 000 000. Consider the ν ar realizations {α ℓ ar } ℓ ∈ C 0 α of A and the ν ar realizations {q ℓ ar } ℓ ∈ R + of Q, all of which have been computed in Step 2. For a sufficiently large value of ν e , the minimum statistics of the positive-valued random variable Q is typically defined as the positive-valued random variable Q min such that Q min = min{Q (1) , . . . , Q (νe) }, where Q (1) , . . . , Q (νe) are ν e independent copies of the positive-valued random variable Q. The random variable Q min depends on ν e . However, the dependence on ν e is not highlighted here in order to keep the notation simple. The ν r realizations {q 1 min , . . . , q νr min } of Q min , which are assumed to be independent, are then computed as q r min = min ℓ ∈ {νe(r-1) +1, ... , rνe}

q ℓ ar , r = 1, . . . , ν r , (4.21)
where the minimum q r min = q ℓr ar is reached for ℓ = ℓ r ∈ {ν e (r -1) + 1, . . . , rν e }. Now, let α r min designate α ℓr ar -that is, α r min = α ℓr ar -and let A min denote the random vector associated with Q min and for which {α r min , r = 1, . . . , ν r } are ν r independent realizations. Note that, here, A min is simply a notation: in particular, it does not represent the minimum statistics of A. It is a random variable for which the realizations correspond to those of Q min .

Step 4. The last step consists in evaluating the cost function J(α g ) at any point α g in C 0 α using only the additional independent realizations {(α r min , q r min ), r = 1, . . . , ν r } of the random variable (A min , Q min ). Taking into account Equations (4. [START_REF] Goldberg | Genetic Algorithms in Search. Optimization & Machine Learning[END_REF]) and (4.20), the objective is to compute E{Q min (α g )} and E{Q min (α g ) 2 }, where for α g given in C 0 α , the random variable Q min (α g ) is the minimum statistics for the random variable Q(α g ) (the stochastic error). For κ = 1, 2,

E{(Q min (α g )) κ } = E{(Q min ) κ | A = α g } , (4.22)
() Prepared using nmeauth.cls where the conditional mathematical expectation appearing in the right-hand side of Equation (4.22) is carried out using the kernel density estimation method and the additional realizations {(α r min , q r min ), r = 1, . . . , ν r } computed in Step 3 (the details of the computation of the conditional expectation can be found in [START_REF] Ghanem | Probabilistic nonconvex constrained optimization with fixed number of function evaluations[END_REF]). REMARK 6. There is no theoretical limit for the value of the dimension m 0 α of the hyperparameter α in the proposed PLM. However, the larger m 0 α the larger the overall computational cost, which is true for any optimization algorithm adapted to the solution of a nonconvex problem. In any case, for the PLM to be effective, the number of data points ν d in the initial dataset for which each point x ℓ ′ d is computed using the nonlinear SPROM must be adapted to the dimension m 0 α of the hyperparameter α, so that this method can actually achieve learning. In general, if m 0 α is increased, the number of data points ν d should also be increased. [START_REF] Bergstra | Random search for hyper-parameter optimization[END_REF] is a candidate choice for the solution of the nonconvex optimization problem defined in Equations (4.17), (4.19), and (4.20); however, this choice would be strongly limited to the case where m 0 α is of the order of a few units and therefore is not necessarily appropriate here. Other candidate choices are any Bayesian optimization algorithm [START_REF] Snoek | Practical Bayesian optimization of machine learning algorithms[END_REF], or any random search [START_REF] Bergstra | Random search for hyper-parameter optimization[END_REF] or genetic algorithm [START_REF] Goldberg | Genetic Algorithms in Search. Optimization & Machine Learning[END_REF]. In general, such algorithms tend to generate a very large number of evaluations of the cost function and therefore are relatively computationally intensive. However, thanks to the PLM introduced for evaluating the cost function, this shortcoming is avoided here. Therefore, such algorithms are appropriate in this context.

Optimization algorithms. A priori, a grid search algorithm

In principle, gradient-based algorithms [START_REF] Bengio | Gradient-based optimization of hyperparameters[END_REF] are not suitable for the solution of nonconvex optimization problems. However, they may be considered here for two reasons:

• A good estimation of the solution α opt is provided by the predictor α 0 and can be used as an advantageous initial condition. In this case, α opt can be searched in C 0 α defined in Equation (4.16), which is localized in a neighborhood of α 0 .

• Thanks to the PLM, a very large number ν ar of independent realizations of the stochastic error is available for constructing a good estimate of the minimum statistics (see Step 3 of the PLM described in Section 4.4.3).

These two specificities justify the use of a deterministic algorithm such as the interiorpoint algorithm (for example, see [START_REF] Byrd | An interior point algorithm for large-scale nonlinear programming[END_REF]) for solving the optimization problem defined in Equations (4.17), (4.19), and (4.20). Furthermore, it is noted that in several numerical experiments associated with applications such as those discussed in Section 5, the deterministic interior-point algorithm outperformed various genetic algorithms at solving the nonconvex optimization problem underlying the construction of the corrector for α 0 .

Initialization of the interior-point algorithm.

Because it was computed using the cost function defined by Equation (4.19), the predictor α 0 cannot be applied "as is" to initialize the solution of the minimization problem (4.17) formulated using the cost function defined by Equation (4.19). Instead, the initial condition of the latter problem, α init , should be chosen in C 0 α (4.16) which, as explained in Section 4.4, accounts for α 0 . Furthermore, since the cost function (4.19, 4.20) is formulated in terms of the minimum statistics, Q min , of the random variable Q, α init can be selected by considering the corresponding Q min . Recalling that the minimum statistics are constructed from the ν ar independent realizations {q ℓ ar , ℓ = 1, . . . , ν ar } of Q associated with the ν ar realizations {α ℓ ar , ℓ = 1, . . . , ν ar } of A, it follows that a candidate choice for α init could be α init = α ℓ init ar , where ℓ init = arg min ℓ∈{1,...,νar} q ℓ ar . Unfortunately, this minimum value is selected almost surely and is too sensitive to the statistical fluctuations. For this reason, a more appropriate choice for α init is to first compute q init such that

q init = arg max q∈R + p Q min (q) , () Prepared using nmeauth.cls
where q → p Q min (q) on R + is the PDF of Q min , then choose α init = α ℓ init ar , where ℓ init is such that ℓ init = arg min ℓ∈{1,...,νar} |q initq ℓ ar |. The PDF p Q min (q) can be estimated using the Gaussian kernel-density estimation method (for example, see [START_REF] Bowman | Applied Smoothing Techniques for Data Analysis[END_REF][START_REF] Scott | Multivariate Density Estimation: Theory, Practice, and Visualization, Second Edition[END_REF]) and the ν r independent realizations {q 1 min , . . . , q νr min } computed in Step 3 of the PLM described in Section 4.4.3.

4.4.6. Quality assessment of the corrector. Once the corrector α opt has been determined by solving problem (4.17), its quality can be assessed as explained in Section 4.3 but after substituting α 0 by α opt .

APPLICATIONS

The performance of the predictor-corrector approach outlined in Section 4 for identifying the vectorvalued hyperparameter α of the nonparametric probabilistic approach for modeling and quantifying model form uncertainties outlined in Section 2 is illustrated here for two 3D, nonlinear, structural dynamics problems associated with two different configurations of the same MEMS device. This device is made of two parts: a mobile part, and a suspended part attached to the mobile part. It is excited by a transient acceleration prescribed at the base of the mobile part. The two considered configurations of this device differ only by the type of nonlinearity that characterizes them. In the first considered configuration (Config-1), the structural nonlinearity is induced by nonlinear stiffnesses between the suspended and mobile parts of the device. In the second one (Config-2), the structural nonlinearity is induced by the presence of barriers that induce shocks between the mobile and suspended parts. For Config-1, the level of model-form uncertainties is relatively low; for Config-2, it is high. Hence, these two different configurations are also chosen for the additional purpose of demonstrating the performance of the nonparametric probabilistic approach for modeling and quantifying model-form uncertainties for two very different levels of uncertainty.

For both problems considered here, the solution of the underlying nonlinear PROM (2.8) is timeintegrated using the midpoint rule equipped with adaptive time-stepping. This, in order to guarantee the convergence at each time-step of the fixed point method applied to the solution of the nonlinear discrete equations of motion with a relative precision of 10 -8 .

Three-dimensional MEMS device

A schematic of the 3D MEMS device considered here is shown in Figure 1(a). The mobile part of this device is constituted of a square frame with a vertical beam attached to it. Its suspended part is constituted of a parallelepipedic solid with two attached vertical beams. The suspended part is attached to the mobile part by a 3D suspension made of 20 springs. Its geometry is described in a Cartesian coordinate system Ox 1 x 2 x 3 that is attached to its mobile part and adopted as a reference frame. The origin O of this frame is located at the bottom left corner of the device. The axis Ox 1 is horizontal and oriented positively from left to right. The axis Ox 2 is vertical and oriented positively from bottom to top. The axis Ox 3 is perpendicular to the plane Ox 1 x 2 ; it is also oriented positively from bottom to top. The MEMS itself is a cylinder of noncircular cross section whose major axis is Ox 3 ; the dimensions of its plane section are given in Figure 1(a). The external width of the square frame is 30×10 -6 m, its external height is 31×10 -6 m, and its depth is 4×10 -6 m. The vertical surface at the bottom of the MEMS is described by

(x 1 ∈ [0 , 30×10 -6 ], x 2 ∈ [0 , 31×10 -6 ], x 3 = -4×10 -6
). The suspended and mobile parts are made of a homogeneous, orthotropic, linear elastic material whose mechanical properties in the aforementioned reference frame are those of a standard (100) silicon wafer [START_REF] Hopcroft | What is the Young's modulus of silicon[END_REF]. Specifically, the Young moduli of this material are E 11 = E 22 = 169×10 9 N/m 2 and E 33 = 130×10 9 N/m 2 ; its Poisson ratios are ν 23 = 0.36, ν 31 = 0.28, and ν 12 = 0.064; its shear moduli are G 23 = G 31 = 79.6×10 9 N/m 2 and G 12 = 50.9×10 9 N/m 2 . The mass density of the silicon material is 2 330 Kg/m 3 . The stiffnesses of the springs and suspension depend on the axis along which they act: k s1 = 4 N/m (along Ox 1 ), k s2 = 6 N/m (along Ox 2 ), and k s3 = 1.5 N/m (along Ox 3 ). Structural damping is represented () Prepared using nmeauth.cls using the global damping model described in Appendix A, which enables the control of the modal damping rate of the structural model in the frequency band of analysis.

As stated above, two different configurations of this device are considered here. They differ only by the type of nonlinearity introduced between the two vertical beams of the suspended part and the vertical beam of the mobile part:

L (ε R ). Initially, ε L = ε R = 0.5×10 -6 m.
For both configurations, zero x 1 -, x 2 -, and x 3 -displacement boundary conditions with respect to the reference frame Ox 1 x 2 x 3 are prescribed at the base of the mobile part of the devicespecifically, at the points {x 1 ∈ [0 , 30×10 -6 ] , x 2 = 0 , x 3 ∈ [0 , -4×10 -6 ]}. The following timedependent, square integrable, and real-valued x 1 -acceleration is also prescribed at these same points

Γ(t) = Γ 0 π t {sin(t(ω c + ∆ω c /2)) -sin(t(ω c -∆ω c /2)} , t ∈ [t 0 , T ] ,
where Γ 0 = 120 m/s 2 , ω c = 2 π×13×10 6 rad/s is the central angular frequency, and ∆ω c = 2 π×10×10 6 rad/s is the angular frequency bandwidth. This acceleration is graphically depicted in Figure 2 and ω e = ω c + ∆ω c /2 = 2π×18×10 6 rad/s. At time t 0 , the device is at rest (its displacement and velocity fields are zero). In all analyses, t 0 < 0 is written as t 0 = -m 0 π/ω e = -2.7778×10 -5 s, where m 0 = 1 000, and the time-interval of analysis is set to [t 0 , 1.8559×10 -4 s]. The latter setting is chosen because at T = 1.8559×10 -4 s, the device is returned to its zero equilibrium with a small relative error. The frequency band of observation is chosen as

B o = [0, ω o ],
where ω o = 2 π×70×10 6 rad/s.

() Prepared using nmeauth.cls Time t (s) ×10 -6 -5 0 5 Γ(t) ×10 7 -2 -1 0 1 2 (a) Time-dependent acceleration t → Γ(t)
Frequency ν (Hz) ×10 For both configurations described above, the x 1 -acceleration of the response of the system is observed at 744 spatial points distributed across various locations of the MEMS device, in particular on the boundaries of the vertical beams of the suspended and mobile parts. Hence, the dimension of each of the observation vector o(t) and its frequency domain representation o(ω) is m o = 744. All of these observations are used to identify the hyperparameter α using the predictor-corrector approach presented in Section 4. However, in order to keep the number of figures within a reasonable limit, two observations are selected among all 744 possible ones and denoted by Obs 1 1 and Obs 1 6 for Config-1, and Obs 2 1 and Obs 2 4 for Config-2. If the dynamical system associated with either considered configuration was linear, the energy of the response signal would be concentrated in the same frequency band [ω min , ω e ] as that of the energy of the excitation signal (see Figure 2(b)). Due to the nonlinearity however, part of the energy of the excitation signal is transferred outside its frequency band and consequently, the frequency band of the response is not [ω min , ω e ] but [0, ω max ], where ω max > ω e is to be determined. To this end: the sampling time-step is chosen as ∆t = π/ω max = 6.9444×10 -9 s, where ω max = 4 ω e = 2π×72×10 6 rad/s; and the frequency band of analysis is set to [-ω max , ω max ], where ω e < ω max . This yields n time = 30 725 time points in the time-interval [t 0 , T ]. The sampling frequency step is set to ∆ω = 2π×4686.7 rad/s, yielding also n freq = 30 725 frequency points in the frequency band of analysis [-ω max , ω max ].

Finite element based nonlinear high-dimensional models

First, a linear, FE-based HDM is constructed for the neutral configuration of the 3D MEMS device defined here as the configuration where no material or barriers are inserted between the two vertical beams of the suspended part of this device, and the vertical beam of its mobile part. This HDM has 7 328 8-noded solid elements, 10 675 nodes, and N = 32 025 dofs (see Figure 1(b)). There are 205 of these nodes, which belong to the base of the mobile part of the device: at each of these nodes, all displacement dofs are constrained to zero in the moving reference frame Ox 1 x 2 x 3 , due to the boundary conditions described in Section 5.2. Hence, N C D = 615. Next, two nonlinear variants of this HDM are generated to account for the nonlinear specificities of Config-1 and Config-2. Both variants have the same dimension. Furthermore, this dimension is the same as that of the linear HDM as both nonlinear variants can be constructed without generating additional dofs. For this reason, for the sake of simplicity, and unless otherwise specified, HDM is used to refer to the linear HDM associated with the linear, neutral configuration of the device, or to the nonlinear HDM of the same dimension associated with Config-1 or Config-2. In all cases, the displacement vector y(t) associated with the HDM is measured in the aforementioned reference frame: therefore, it represents () Prepared using nmeauth.cls a relative displacement vector. The HDM is used to analyze the dynamic response of the neutral configuration, configuration Config-1, or configuration Config-2 of the device to the prescribed excitation Γ(t). In all cases, the associated governing equations for the relative displacement vector, y(t), are those given in Equations (2.1) to (2.3), where:

-The vector of external forces f(t) is a function of the mass matrix [M ], the x 1 -dofs where zero relative displacement boundary conditions are applied, and the relative acceleration Γ(t). -The vector of internal forces g(y(t), ẏ(t)) can be written as

g(y(t), ẏ(t)) = [D] ẏ(t) + [K] y(t) + f NL (y(t)) , (5.1)
where:

-As stated in Section5.1, the damping matrix [D] is constructed using the global damping model described in Appendix A, which enables the control of the modal damping rate of the modeled device in the frequency band of analysis. -The stiffness matrix [K] belongs to M +0 N and has a null space of dimension 6 before applying the zero boundary conditions.

-The nonlinear vector function f NL (y(t)) is constructed to account for the structural nonlinearity introduced in Config-1 or Config-2, as applies.

-The initial conditions are y(t 0 ) = ẏ(t 0 ) = 0.

-The matrix

Projection-based reduced-order models

First, a deterministic ROB [V ] ∈ M N,n of dimension n ≪ N is constructed using the classical method of snapshot collection and compression. For each of Config-1 and Config-2, the nonlinear HDM is exercised to compute the time-sampled response {y(j ∆t), j = 1, . . . For n time = 30 725 and µ snp = 10 (n snp = 3 072), the aforementioned convergence rate analysis yields n = 8 for Config-1 and n = 20 for Config-2 (note that in both cases, a ROB based on the natural modes of the linear, neutral configuration of the device requires n mode = 40 modes for convergence and therefore is less efficient).

Next, each constructed ROB [V ] ∈ M N,n is randomized following the procedure given by Equations (2.17-2.22) to obtain the counterpart SROB [W]. Then, the corresponding nonlinear PROM and SPROM are constructed by Galerkin projection of the nonlinear HDM as in Equations (2.7-2.9) and (2.13-2.15), respectively.

Model-form uncertainty quantification

For each nonlinear PROM constructed in Section 5.3, a modified version is designed to generate targets for both selected observations in lieu of experimental data. Specifically, each modified nonlinear PROM constructed for Config-1 or Config-2 is designed to generate model-form uncertainties.

Figure 3 (Config-1) and Figure 4 (Config-2) display the graphs of the function ν → db target (2πν) of the observed targets, the graphs of the function ν → db (n) (2πν) of the counterpart values predicted using the nonlinear PROM, and the graphs of the function ν → db lin (2πν) of the counterpart results obtained using the linear HDM associated with the linear, neutral configuration of the device (black thin line).

• The effect of the structural nonlinearity is very important outside the frequency band [ω min , ω e ] of the excitation, as there are large differences between the numerical predictions obtained using the nonlinear PROM and those obtained using the linear HDM associated with the neutral, linear configuration of the device. • The model-form uncertainties -which can be analyzed by estimating the difference between the target and numerically predicted values -are not homogeneous neither as a function of the frequency, nor as a function of the selected observation. • The effects of the model-form uncertainties are small in the frequency band [ω min , ω e ] of the excitation, but can be very large outside this band. • For observations Obs 1 1 and Obs 1 2 , the effects of the model-form uncertainties are smaller than for the observations Obs 1 6 and Obs 2 4 .

Identification of the hyperparameter via the predictor.

For each of Config-1 and Config-2, the predictor α 0 for the hyperparameter α is computed using the method presented in Section 4.2. For Config-1, n = 8 and therefore the number of hyperparameters associated with the matrix m σ = 36. For Config-2, n = 20 and therefore m σ = 210. Specifically, Stage 1 is used to compute the predictor [σ 0 ]: it leads to m 0 σ = 4 for Config-1 and to m 0 σ = 13 for Config-2. Hence, it achieves a significant reduction in the number of hyperparameters associated with the matrix [σ 0 ]. In Stage 2, the optimization problem defined in Equations (4.7) and (4.8) is solved using the grid search algorithm and the grid {1.0, 1.3, 1.5, 1.8, 2.0, 2.2, 2.5} × 10 -5 with m s = 7 points for Config-1, and the grid {1.0, 2.0, 3.0, 4.0, 5.0, 6.0} × 10 -6 with m s = 6 points for Config-2. For each value of α = (s, β 0 , [σ 0 ]), the nonlinear SPROM is solved for ν s = 100 independent realizations.

Figure 5(a) and Figure 5(b) show the graph of the function s → J freq (s) for Config-1 and Config-2, respectively. The predictor s 0 for s is equal to 2.0×10 -5 for Config-1 and to 4.0×10 -6 for Config-2.

The quality of the computed predictor α 0 is assessed by plotting each target and the corresponding confidence region constructed using the nonlinear SPROM for α = α 0 and computed using Equation (4.15) with p c = 0.98. For each selected observation, Figure 6 (Config-1) and Figure 7 (Config-2) display the graph of the function ν → db target (2πν) of the observed target, the graph of the function ν → db (n) (2πν) of the observation predicted using the nonlinear PROM, and the graph of the confidence region of the random function ν → dB (n) (2πν) of the observation computed using the nonlinear SPROM. A quality criterion is that the target belongs to the confidence region with a probability of 0.98. The reader can observer that for this criterion, the predictor delivers a reasonable performance. Nevertheless, Figure 6 shows that the peak of the target at 48.76×10 6 Hz does not belong to the confidence region and Figure 7 shows that the upper envelope of the confidence region is somehow high. For these reasons, the corrector for α 0 is considered next.

Identification of the hyperparameter via the corrector.

For both Config-1 and Config-2, the parameters w J and γ J of the cost function J(α) defined by Equation (4.19) are set here to w J = 0.5 and γ J = 1.0. A convergence analysis of the PLM is performed with respect to the dimension ν d of the initial dataset computed in its Step 1, using the following values of ν d : 100, 300, 500, 1 000, 1 500, and 2 000. For each of these, the convergence of the quality indicator I is analyzed with respect to the number of additional realizations ν ar using the following set of values for this variable: {200 000, 1 000 000, 2 100 000}. For each value of ν d in this set, the value of I using ν ar = 1 000 000 is found to be close to its counterpart computed using ν ar = 2 100 000, which suggests that convergence is reached for ν ar = 1 000 000. Still, the convergence analysis with respect to ν d and the QA of the corrector are performed using the larger value of ν ar = 2 100 000. To this end, the parameters introduced in Equation (4.21) are chosen such that ν e = 100 and ν r = 21 000. The results of the convergence analysis of the PLM with respect to ν d are reported in Table I. They show that for Config-1, the corrector improves the value of the quality indicator obtained for the predictor as soon as ν d ≥ 1 000; the best value is reached for ν d = 2 000. For Config-2, the corrector improves the value of the quality indicator obtained for the predictor as soon as ν d ≥ 1 500, and the best value of I is obtained for ν d = 1 500. The quality of the corrector for α is further assessed by plotting each target and the corresponding confidence region constructed using the nonlinear SPROM for α = α opt and computed using Equation (4.15) with p c = 0.98. Attention is focused on the results obtained when the PLM is converged with respect to ν d -that is, for ν d = 2 000 for Config-1 and ν d = 1 500 for Config-2 () Prepared using nmeauth.cls -and with respect to the number of additional realizations -that is, for ν ar = 2 100 000 for both configurations. The optimal values of the hyperparameters s and β of α opt are found to be s opt = 1.01×10 -5 and β opt = 0.026 for Config-1, and s opt = 3.56×10 -6 and β opt = 0.016 for Config-2. Figure 8 (Config-1) and Figure 9 (Config-2) display for each selected observation the graph of the function ν → db target (2πν) of the observed target, the graph of the function ν → db (n) (2πν) of the observation computed using the nonlinear PROM, and the graph of the confidence region of the random function ν → dB (n) (2πν) of the observation predicted using the nonlinear SPROM. Again, a quality criterion is that the target belongs to the confidence region with a probability of 0.98. Comparing these figures with Figures 6 and7 shows that the corrector for α improves the confidence region obtained using the predictor for this vector-valued hyperparameter. In particular, Figure 8 shows that the peak of the target at 48.76×10 6 Hz is now inside the computed confidence region and Figure 9 shows that the upper envelope of the computed confidence region is now lower. Figure 8 also shows that for some frequencies, the response of the nonlinear PROM is outside the computed confidence region, but the target is always well within the computed confidence region. These results highlight the ability of the overall nonparametric probabilistic method for modeling and quantifying model-form uncertainties to reach a given target.

5.4.3.

Wall-clock time performance analysis. All computations described above were performed on a Linux cluster in double precision arithmetic. The main steps of the model-form UQ analysis were performed in parallel execution mode on 70 cores of this cluster as described below:

-Predictor (Stage 2): this stage was performed with m s = 6 for Config-1 and m s = 7 for Config-2; for both configurations, the nonlinear SPROM was exercised for ν s = 100 independent realizations. -Predictor (QA): for this purpose, the nonlinear SPROM was exercised for ν sim = 1 000 independent realizations. -Corrector (QA): for this purpose, the nonlinear SPROM was exercised for ν sim = 1 000 independent realizations.

The wall-clock time elapsed in the execution of each of these main steps and the total wall-clock time are reported in Table II for each of Config-1 and Config-2. Two observations are noteworthy:

• In each case, the total wall-clock time of the predictor-corrector approach is dominated by the computational cost of the corrector component: specifically, 94% of the total wall-clock time in the case of Config-1 and 80% in the case of Config-2 is due to the corrector component. • The model-form UQ analysis of Config-2 is more computationally intensive than that of Config-1, due to the presence of shocks in the second configuration. However, most of the difference between the computational costs of both UQ analyses is due to the predictor component for which the wall-clock time for Config-2 is roughly 4.6 times higher than for Config-1.

Table II. Wall-clock time in seconds (and in hours between parentheses) on 70 cores of a Linux cluster for the main steps of the model-form UQ analysis. As shown above, the solution of the optimization problem (4.17) -that is, the computation of the corrector α opt -dominates the total cost of the identification of the vector-valued hyperparameter α. The total cost of this identification itself dominates the overall computational cost of the () Prepared using nmeauth.cls nonparameteric probabilistic method for model-form UQ analysis at the heart of this paper. The computer implementation of the standard solution of problem (4.17) can be described as an overall algorithm with three main loops: an outer-loop where each instance is associated with one iteration of the optimizer and incurs a few evaluations of the cost function to be minimized; an innerloop for performing the ν sim independent realizations (here ν sim = 1 000); and an innermost-loop for performing the SPROM-based simulation in the time domain. On the other hand, once training has been performed, the computer implementation of the solution of problem (4.17) using the PLM described in this paper is organized around two main loops only: the same aforementioned outerloop, and the same aforementioned innermost-loop. Hence, the gain in computational speed enabled by the proposed predictor-corrector approach for identifying α can be easily estimated using the information given above and in Table II. To this end, it is first noted that ν sim = 1 000 executions of the nonlinear SPROM in the time domain using 70 cores and the algorithmic parameters described in the introduction of Section 5 and in Section 5.1 require on average 181 seconds wall-clock time for Config-1, and 1 448 seconds for Config-2. It follows that a good estimate of the wall-clock time on 70 cores for the standard approach for determining α opt is 141 × 181 = 25 521 seconds for Config-1, and 950 × 1 448 = 1 375 600 seconds for Config-2. Therefore, the speedup factors enabled by the PLM proposed in this paper are equal to 1.33 in the case of Config-1, and 53.61 in the case of Config-2 where the effect of the nonlinearity on the solution of the optimization problem (4.17) is much stronger. In any case, it is emphasized here that these speedup factors are only estimates as the analysis given here assumes that the number of outer-loop instances associated with the interior point algorithm is the same in the presence and absence of the PLM.

Main step

SUMMARY AND CONCLUSIONS

A novel, predictor-corrector approach for identifying the vector-valued hyperparameter of a nonparametric probabilistic method for modeling and quantifying model-form uncertainties is presented in this paper. Its key component is a probabilistic learning on manifolds for constructing a low-cost, surrogate model for the Monte Carlo simulations of interest. Specifically, this model enables the economical generation, within the solution of an inverse statistical problem, of stochastic realizations that are concentrated on a manifold identified using an initial set of data points. The computational advantages of this proposed approach are demonstrated in this paper for the modelform uncertainty quantification analysis, using the aforementioned nonparametric probabilistic method, of two three-dimensional, nonlinear, structural dynamics problems associated with two different configurations of a MEMS device. For this application, a speedup factor -with respect to a standard approach for solving an inverse statistical problem -as high as 56 is achieved.

  where c dM = 4 and c dm = 1/c dM . -If σ k corresponds to an upper extradiagonal entry of [σ], there are two cases to consider depending on the sign of σ 0 k : * If σ 0 k > 0, the lower bound is σ L k = c em ×σ 0 k and the upper bound is σ U k = c eM ×σ 0 k , where c eM = 20 and c em = 1/c eM . * If σ 0 k < 0, the lower bound is σ L k = c eM ×σ 0 k and the upper bound is σ U k = c em ×σ 0 k . Note that the proposed values for c dM and c eM are sufficiently large, given the normalization of the matrix [σ]

Figure 1 .

 1 Figure 1. Schematic and FE-based HDM of a 3D MEMS device with a mobile part (yellow) and a suspended part (blue) attached to the mobile part via a suspension system.

Config- 1 .

 1 In this configuration, a nonlinear elastic material is inserted between the aforementioned beams. Its constitutive equation corresponds to a cubic, elastic, restoring force whose elastic constant is k b = 2×10 12 N/m. Config-2. In this alternative configuration, two continuous elastic barriers are introduced at the left and right sides of the vertical beam attached to the mobile part. The constitutive equation of each elastic barrier is linear and characterized by the elastic constant k b = 2×10 12 N/m. The horizontal gap between the left (right) beam of the suspended part and the beam of the mobile part is denoted by ε

  (a). The graph of ω → | Γ(ω)| defined on 2 π×[0 , 70×10 6 ] rad/s, where Γ(ω) = T t0 e -iωt Γ(t) dt, is plotted in Figure 2(b). The energy of the excitation signal is mainly concentrated in the frequency band [-ω e , -ω min ] ∪ [ω min , ω e ], where ω min = ω c -∆ω c /2 = 2π×8×10 6 rad/s

Figure 2 .

 2 Figure 2. Prescribed acceleration at the base of the mobile part of the MEMS device. (a): graph of Γ(t) defined in [-2.78×10 -5 , 1.86×10 -4 ] s and represented in [-5×10 -6 , 5×10 -6 ] s. (b): graph of | Γ(2πν)| (in the log 10 -scale) defined in [-72×10 6 , 72×10 6 ] Hz and represented in the frequency band of observation [0 , 70×10 6 ] Hz.

  [B] ∈ M N,NCD represents the N C D zero boundary conditions and is constructed such that [B] T [B] = [I NCD ].

  , n time }. The solution snapshots {y k snp , k = 1, . . . , n snp }, where y k snp = y(k µ snp ∆t), µ snp denotes the snapshot-step, and n snp = n time /µ snp , are extracted from this response. Next, the thin singular value decomposition of the matrix [y 1 snp . . . y nsnp snp ] is computed, and the singular values s 1 ≥ . . . ≥ s nsnp are ordered in descending order. Then, [V ] is defined as the matrix whose columns are the n left orthonormal vectors associated with the first n ordered singular values s 1 ≥ . . . ≥ s n , where n is determined from the convergence rate of the compression measured by

6 Figure 3 .

 63 Figure 3. Config-1: graphs of the function ν → db target (2πν) of the observed targets (red thick lines), graphs of the function ν → db (n) (2πν) of the counterpart results predicted using the nonlinear PROM (blue mid solid lines), and graphs of the function ν → db (n) lin (2πν) of the counterpart results obtained using the linear HDM associated with the linear, neutral configuration of the device (black thin line).

Frequency 4 Figure 4 .

 44 Figure 4. Config-2: graphs of the function ν → db target (2πν) of the observed targets (red thick lines), graphs of the function ν → db (n) (2πν) of the counterpart results predicted using the nonlinear PROM (blue mid solid lines), and graphs of the function ν → db (n)

2 Figure 5 .

 25 Figure 5. Graph of the function s → J freq (s) for the grid used by the grid search algorithm.

Figure 6 .

 6 Figure 6. QA of the predictor for Config-1: graphs of the function ν → db target (2πν) of the observed target (red thick line), graphs of the function ν → db (n) (2πν) of the observation computed using the nonlinear PROM (blue mid solid line), and graphs of the confidence region of the random function ν → dB (n) (2πν) of the observation predicted using the nonlinear SPROM (yellow region with green lines for the upper and lower envelopes).

(a) Observation Obs 2 1 (b) Observation Obs 2 4 Figure 7 .

 147 Figure 7. QA of the predictor for Config-2: graphs of the function ν → db target (2πν) of the observed target (red thick line), graphs of the function ν → db (n) (2πν) of the observation computed using the nonlinear PROM (blue mid solid line), and graphs of the confidence region of the random function ν → dB (n) (2πν) of the observation predicted using the nonlinear SPROM (yellow region with green lines for the upper and lower envelopes).

Figure 8 .

 8 Figure 8. QA of the corrector for Config-1: graphs of the function ν → db target (2πν) of the observed target (red thick line), graphs of the function ν → db (n) (2πν) of the observation computed using the nonlinear PROM (blue mid solid line), and graphs of the confidence region of the random function ν → dB (n) (2πν) of the observation predicted using the nonlinear SPROM (yellow region with green lines for the upper and lower envelopes).

-

  Corrector (Step 1 PLM): Step 1 of the PLM for constructing the initial dataset was performed with ν d = 2 000 for Config-1 and ν d = 1 500 for Config-2. -Corrector (Step 2 PLM): Step 1 of the PLM for computing additional realizations was performed with ν ar = 1 050×ν d for Config-1 and ν ar = 1 400×ν d for Config-2. -Corrector (Step 3 & Step 4 PLM): the interior point algorithm converged in 19 iterations during which it performed 141 evaluations of the cost function (4.19) in the case of Config-1; it converged in 58 iterations during which it performed 950 evaluations of the cost function in the case of Config-2. () Prepared using nmeauth.cls

Figure 9 .

 9 Figure 9. QA of the corrector for Config-2: graphs of the function ν → db target (2πν) of the observed target (red thick line), graphs of the function ν → db (n) (2πν) of the observation computed using the nonlinear PROM (blue mid solid line), and graphs of the confidence region of the random function ν → dB (n) (2πν) of the observation predicted using the nonlinear SPROM (yellow region with green lines for the upper and lower envelopes).

  .12)• The probability distribution P [W] of [W] depends on a vector-valued hyperparameter α = (α 1 , . . . , α mα ) belonging to a subset C α of R mα .

	The construction of such an SROB is summarized in Section 2.4. The nonlinear SPROM associated
	with the nonlinear PROM described by Equations (2.7) to (2.10) is obtained by substituting [V ] with
	the random matrix [W].

  , T ]} with values in R mo be the vector of observed targets associated with {y target (t), t ∈ [t 0 , T ]} with values in R N . If the model-form uncertainties and/or modeling errors tainting the SPROM defined by (2.13) and(2.14) are due only to model reduction, then o target = o. On the other hand, if they are due to both of the model reduction process and the model-form uncertainties and/or modeling errors existing in the HDM, then o target = o exp , where o exp is based on experimental data.

	[U] (specifically, see Eq. (D.25)
	on page 875 of [4]). Hence, identifying [c n ] = [σ] T [σ] can be performed by identifying
	instead the upper triangular matrix [σ]. More generally, one can identify an upper triangular
	matrix [σ] with positive or zero diagonal entries and define [c n ] = [σ] T [σ]. In both cases, the
	matrix-valued hyperparameter [σ] enables the control of the statistical correlations between
	the columns of the random matrix [G(β)] using roughly half the number of individual
	hyperparameters that would be otherwise needed to identify [c n ] directly.
	3. IDENTIFICATION OF THE HYPERPARAMETER
	Let {o target (t), t ∈ [t 0

() Prepared using nmeauth.cls

  Similarly to the construction of the predictor s 0 for s (see Stage 2 in Section 4.2), the construction of the corrector for the hyperparameter α consists in solving an optimization problem of the form

.16) 4.4.2. Optimization problem.

Table I .

 I Convergence analysis of the PLM with respect to ν d .

		Predictor			Corrector
	ν d Config-1: I Config-2: I	-3.18 2.01	100 -4.30	300 -4.60	500 3.96 2.80	1 000 3.02 4.21	1 500 3.07 1.90	2 000 2.89 1.83

33 hrs) (7.13 hrs)

	Config-1	Config-2
	Predictor (Stage 2)	934	3 991
	Predictor (QA)	185	1 134
	Corrector (Step 1 PLM)	1 850	1 833
	Corrector (Step 2 PLM)	16 018	17 243
	Corrector (Step 3 & Step 4 PLM)	11	13
	Corrector (QA)	178	1 448
	Total (Predictor + Corrector)	19 176	25 662
	(5.

ACKNOWLEDGEMENTS

The second author acknowledges partial support by the Army Research Laboratory through the Army High Performance Computing Research Center under Cooperative Agreement W911NF-07-2-0027, partial support by The Boeing Company under Contract Sponsor Ref. 45047, and partial support by DARPA under the Enabling Quantification of Uncertainty in Physical Systems (EQUiPS) program. This document does not necessarily reflect the position of these institutions, and no official endorsement should be inferred.

APPENDIX A: CONSTRUCTION OF THE DAMPING MATRICES [D]

The structural damping model chosen for constructing the matrix [D] of the FE-based HDMs of the 3D MEMS device described in Section 5 is based on the neutral configuration of this device defined in Section 5.2 -that is, the configuration where no material or barriers are inserted between the two vertical beams of the suspended part of this device and the vertical beam of its mobile part.

() Prepared using nmeauth.cls Hence, this model leads to the same matrix [D] for both Config-1 and Config-2. It corresponds to the following adaptation of the structural damping representation presented in [START_REF] Ohayon | Structural Acoustics and Vibration[END_REF] (Chapter VI, Equation ( 45)).

Let [Φ] = [ϕ 1 . . . ϕ n mode ] denote the matrix of the first n mode natural vibration modes associated with the first n mode lowest eigenvalues 0 < λ 1 ≤ . . . ≤ λ n mode of the undamped neutral configuration of the 3D MEMS device. In the time domain, the dynamics of this undamped configuration are governed by the linear equations of motion [M ] ÿ(t) + [K] y(t) = 0. Then, the deterministic matrix

N is constructed here as follows:

The reader can verify that the above definition of [D] satisfies

where the generalized damping matrix [D] ∈ M + n mode is diagonal and its diagonal entries are

It follows that:

• The structural damping model defined by Equation (A.1) yields a constant damping rate ξ d for the first n mode natural modes of the system, and yields the damping rates

• This model is a global damping model that enables the control of the modal damping rate in the frequency band of analysis.

By construction, n mode is chosen such that n mode ≥ n. Then, the deterministic and stochastic reduced-order damping matrices

are nonsingular. They can also be written as

where χ k = [V ] T [M ] ϕ k is a deterministic vector with values in R n , and

where Ψ k = [W] T [M ] ϕ k is a random vector with values in R n .