M Arnst

C Soize

Identification and sampling of Bayesian posteriors of high-dimensional symmetric positive-definite matrices for data-driven updating of computational models

Keywords: data-driven modeling, Bayesian inversion, symmetric positive-definite matrix, identification criterion, nonparametric probabilistic model, Markov Chain Monte Carlo

Identification and sampling of Bayesian posteriors of high-dimensional symmetric positive-definite matrices for data-driven updating of computational models

Maarten Arnst, Christian Soize

Introduction

We propose a probabilistic methodology for data-driven updating of non-Gaussian high-dimensional symmetric positive-definite matrices involved in computational models. For instance, the computational model could be a reduced-order model parameterized by symmetric positive-definite reduced matrices, or it could involve constitutive models with material tensors represented by symmetric positive-definite matrices, among other possibilities. In this methodology, we cast the updating as a Bayesian identification of the symmetric positive-definite matrices. Our aim is not only to improve the accuracy of the predictions but also to learn from the data about the impact of parametric uncertainties and modeling errors on the predictions.

Whereas an inverse problem formulated within a deterministic framework typically takes the form of an optimization problem, the Bayesian approach to inverse problems takes the form of an application of the Bayes formula to update a prior probability distribution by means of a likelihood function to obtain a posterior probability distribution [START_REF] Tarantola | Inverse Problem Theory[END_REF][START_REF] Flath | Fast algorithms for Bayesian uncertainty quantification in large-scale linear inverse problems based on low-rank partial Hessian approximations[END_REF][START_REF] Dashti | The Bayesian approach to inverse problems[END_REF]. The numerical implementation of an inverse problem within this Bayesian approach typically involves a Markov chain Monte Carlo (MCMC) method to sample from the posterior, and descriptors of the posterior are then estimated from the samples thus obtained by means of methods from mathematical statistics [START_REF] Robert | Monte Carlo Statistical Methods[END_REF].

Casting the data-driven updating as a Bayesian identification of the symmetric positive-definite matrices raises challenges, both from the point of view of the formulation of the Bayesian inversion as from the point of view of the numerical implementation of the Bayesian inversion. First, the prior must be a matrix-variate probability distribution whose support is limited to a set of symmetric positive-definite matrices. The sets of symmetric positivedefinite matrices are not linear spaces (vector spaces) and this probability distribution cannot be Gaussian. Second, the prior and the likelihood function must be defined-and their hyperparameters identified-in a manner that is consistent with our aim of obtaining not only improved predictions but also a probabilistic representation of the model uncertainties. Third, an MCMC method is required that is able to explore efficiently the set of symmetric positive-definite matrices in accordance with the posterior. The MCMC method must ensure the symmetry and the positive definiteness of the matrix samples. And the dimensionality of the problem grows with the square of the matrix dimension and hence can be high, so that conventional MCMC methods can be inefficient. For example, Metropolis-Hastings variants of MCMC methods require a kernel to accept and reject proposals from, and choosing an adequate kernel can be difficult in high dimension.

We address these challenges as follows. First, we select the prior from the family of probability distributions for symmetric positive-definite random matrices introduced with a constructive approach based on the maximumentropy principle in [START_REF] Soize | Nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF][START_REF] Soize | Non-Gaussian positive-definite matrix-valued random fields for elliptic stochastic partial differential operators[END_REF]. This family of maximum-entropy probability distributions lends itself well to representing uncertainty in symmetric positivedefinite matrices of computational models [START_REF] Soize | Nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF][START_REF] Soize | Non-Gaussian positive-definite matrix-valued random fields for elliptic stochastic partial differential operators[END_REF][START_REF] Soize | A comprehensive overview of a non-parametric probabilistic approach of model uncertainties for predictive models in structural dynamics[END_REF][START_REF] Chebli | Experimental validation of a nonparametric probabilistic model of nonhomogeneous uncertainties for dynamical systems[END_REF][START_REF] Arnst | A nonparametric probabilistic model for ground-borne vibrations in buildings[END_REF][START_REF] Cottereau | Construction of a probabilistic model for impedance matrices[END_REF][START_REF] Mignolet | Nonparametric stochastic modeling of linear systems with prescribed variance of several natural frequencies[END_REF][START_REF] Soize | Time-domain formulation in computational dynamics for linear viscoelastic media with model uncertainties and stochastic excitation[END_REF][START_REF] Batou | Rigid multibody system dynamics with uncertain rigid bodies[END_REF][START_REF] Guilleminot | Stochastic model and generator for random fields with symmetry properties: application to the mesoscopic modeling of elastic random media[END_REF]. Second, we define the likelihood function as an exponential of minus a weighted distance that represents the model-data misfit. The posterior thus obtained exhibits several hyperparameters that control the dispersion of the prior and the weight of the weighted distance that represents the model-data misfit involved in the likelihood function. Using an identification criterion that quantifies the agreement between the predictions and the data, we identify these hyperparameters so as to obtain not only improved predictions but also a probabilistic representation of the model uncertainties. Third, we introduce a transformation of measure from the set of symmetric positive-definite matrices to a linear space. This transformation of measure allows us to set up an MCMC method that ensures the symmetry and the positive definiteness of the matrix samples while being defined in terms of state variables with values in a linear space. We use the Itô-SDE-based MCMC method introduced in [START_REF] Soize | Construction of probability distributions in high dimension using the maximum entropy principle: Applications to stochastic processes, random fields and random matrices[END_REF], which does not require a kernel to be chosen and inherits computational efficiency in high dimension from leveraging the gradient of the target to know "where to go" in the high-dimensional space [START_REF] Soize | Construction of probability distributions in high dimension using the maximum entropy principle: Applications to stochastic processes, random fields and random matrices[END_REF][START_REF] Arnst | Itô-SDE MCMC method for Bayesian characterization of errors associated with data limitations in stochastic expansion methods for uncertainty quantification[END_REF].

We apply the proposed methodology to a problem of a data-driven updating of a reduced-order model in computational structural dynamics. In this domain of application, the data-driven updating with a maximum-entropy probability distribution used as the prior can also be seen as a data-driven updating of a nonparametric probabilistic model constructed as introduced in [START_REF] Soize | Nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF]. However, the end result is not an identification of hyperparameters of a nonparametric probabilistic model as in [START_REF] Soize | A comprehensive overview of a non-parametric probabilistic approach of model uncertainties for predictive models in structural dynamics[END_REF][START_REF] Chebli | Experimental validation of a nonparametric probabilistic model of nonhomogeneous uncertainties for dynamical systems[END_REF][START_REF] Arnst | A nonparametric probabilistic model for ground-borne vibrations in buildings[END_REF][START_REF] Soize | Probabilistic model identification of uncertainties in computational models for dynamical systems and experimental validation[END_REF][START_REF] Arnst | Inversion of probabilistic structural models using measured transfer functions[END_REF]; indeed, the posterior is in general no longer a member of the family of maximum-entropy probability distributions. Rather, we obtain an updated probabilistic reduced-order model with a posterior for the reduced matrices that is specifically tailored to provide improved predictions of the dynamical behavior of the structure as well as a probabilistic representation of the model uncertainties.

The proposed methodology relates to the broader context of MCMC methods [START_REF] Robert | Monte Carlo Statistical Methods[END_REF][START_REF] Liu | Monte Carlo strategies in scientific computing[END_REF] as follows. MCMC methods for sampling from probability distributions whose support is not a linear space may be obtained by using a rejection mechanism to discard samples that lie outside the support, by using barrier functions to repel the Markov chain from the boundary of the support, by using reparameterizations, as well as in other ways. In the proposed methodology, the use of the transformation of measure may be viewed as the use of a reparameterization. MCMC methods for sampling from complex probability distributions may be accelerated by exploiting local or global structure of the target probability distribution. Such local or global structure may be deduced from analyses of previous samples, from gradients and Hessians of the posterior, as well as from other sources of information; and it may be exploited by adaptivity mechanisms; by multiscale methods; by surrogate modeling, for instance, in low dimension, the method in [START_REF] Marzouk | Stochastic spectral methods for efficient Bayesian solution of inverse problems[END_REF]; by transformations such as transport maps, for instance, in low dimension, the method in [START_REF] Parno | Transport map accelerated Markov Chain Monte Carlo[END_REF]; by deducing transition kernels from stochastic dynamics in a manner that exploits derivative information. In the proposed methodology, the use of the Itô-SDE-based MCMC method may be viewed as a method that may be classified in the latter class of methods.

This paper is organized as follows. In Sec. 2, we describe in a general setting a formulation and a numerical implementation of a Bayesian identification of symmetric positive-definite matrices. In Sec. 3, we provide details relevant to the use of a maximum-entropy probability distribution as the prior. In Sec. 4, we cast a problem of a data-driven updating of a reducedorder model in structural dynamics as a Bayesian identification of symmetric positive-definite matrices. In Sec. 5, we provide numerical results.

Bayesian identification of symmetric positive-definite matrices

We now describe in a general setting a formulation and a numerical implementation of a Bayesian identification of a symmetric positive-definite matrix. Specifically, we consider a Bayesian inversion of a model

y = g([X]), g : M + n (R) → R m , (1)
that solves a forward problem, that is, that associates to a value [X] of the symmetric positive-definite matrix, which we assume to be of dimension n, a corresponding value y of observables, which we assume to be collected in a vector of dimension m; here and in the remainder of this paper, M + n (R) denotes the set of n-dimensional symmetric positive-definite matrices with real entries. We denote by y obs the actually observed value of the observables.

Formulation within the Bayesian approach to inverse problems

The Bayesian approach to inverse problems leads to the Bayes formula that updates a prior with a likelihood to obtain a posterior as follows [START_REF] Tarantola | Inverse Problem Theory[END_REF][START_REF] Flath | Fast algorithms for Bayesian uncertainty quantification in large-scale linear inverse problems based on low-rank partial Hessian approximations[END_REF][START_REF] Dashti | The Bayesian approach to inverse problems[END_REF]:

ρ [X] [X] y obs posterior ∝ f y obs [X] likelihood × ρ 0 [X] prior . (2)
In this Bayesian approach to inverse problems, the prior probability density function (PDF) ρ 0 is construed to represent any information available about the symmetric positive-definite matrix other than that contained in the data; the likelihood function f is construed to represent the information extracted about the symmetric positive-definite matrix from the data; and the posterior PDF ρ [X] is construed to represent the plausibility of values of the symmetric positive-definite matrix as this plausibility is implied by all the available information. In [START_REF] Flath | Fast algorithms for Bayesian uncertainty quantification in large-scale linear inverse problems based on low-rank partial Hessian approximations[END_REF] and in the remainder of this paper, the notation ∝ expresses that the PDF in the left-hand side is proportional to the right-hand side, whereby the proportionality constant is implied by the requirement that the PDF be normalized: it must integrate to 1 over its domain of definition. The determination of a suitable prior is where the symmetry and the positive definiteness of the matrix pose a first difficulty. Indeed, the prior must assign positive weight only to symmetric positive-definite matrices, that is, it must satisfy π 0 (M + n (R)) = 1, where π 0 is the probability distribution that admits ρ 0 as a PDF, written as π 0 (dX) = ρ 0 ([X])dX; for details about the definition of the volume element, dX, on M + n (R), we refer to [START_REF] Soize | Nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF]. Through the Bayes formula, the posterior then inherits the symmetry and positivedefiniteness properties encoded in the prior, so that samples from the posterior are symmetric positive-definite matrices. We address this difficulty by using a prior from the family of PDFs for symmetric positive-definite random matrices constructed by means of the maximum-entropy principle in [START_REF] Soize | Nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF][START_REF] Soize | Non-Gaussian positive-definite matrix-valued random fields for elliptic stochastic partial differential operators[END_REF]. However, as the proposed methodology is applicable even for other types of prior, we proceed with presenting it in a general setting with the prior denoted generically by ρ 0 ; in Sec. 3, we describe the use of the proposed methodology with a maximum-entropy PDF used as the prior.

We use for the likelihood function an expression as an exponential of minus a weighted least-squares model-data misfit represented by a distance:

f y obs [X] ∝ exp - 1 2 y obs -g([X]) 2 [Γ] ; (3)
here,

2 [Γ] = T [Γ] with = y obs -g([X]
). If the data are perturbed by noise, this expression can be interpreted as a representation of the impact of the noise as a centered multivariate Gaussian random variable added to the observables predicted by the forward problem; thus, [Γ] -1 is the covariance matrix of the Gaussian random variable. If the data are free of noise, this expression can be interpreted as a Gaussian regularization of a Dirac impulse concentrated at the data.

Transformation of measure

The posterior does not in general possess a special structure that can be exploited to set up a direct simulation method to sample from it. Thus, to sample from the posterior, an MCMC method must be used.

The construction of such an MCMC method is where the symmetry and the positive definiteness of the matrix pose a second difficulty: the samples produced by the MCMC method must satisfy the symmetry and positivedefiniteness constraint. We address this difficulty by using a transformation of measure. We assume that even if the posterior does not possess a special structure, the prior (is chosen such that it) does possess a special structure that can be exploited to represent the prior on the set of symmetric positivedefinite matrices in a relatively simple way as a deterministic transformation of another PDF on a linear space. Specifically, we assume that there exists a relatively simple deterministic transformation X from the linear space R n(n+1)/2 to the set of of symmetric positive-definite matrices M + n (R),

X : R n(n+1)/2 → M + n (R) : ξ → [X] = [X (ξ)], (4)
where the dimension n(n + 1)/2 reflects the matrix symmetry, such that the image of the standard Gaussian on R n(n+1)/2 under X is the prior on M + n (R):

π 0 = π N • X -1 , π 0 (dX) = ρ 0 ([X])dX, π N (dξ) ∝ exp - 1 2 ξ 2 dξ. ( 5
)
For the maximum-entropy PDF, such a transformation can be deduced immediately from the direct simulation method available to obtain samples from it, as we show in Sec. 3. Given such a transformation, we define the PDF ρ Ξ on R n(n+1)/2 whose image under X is the posterior on M + n (R):

ρ Ξ ξ y obs ∝ exp - 1 2 y obs -g [X (ξ)] 2 [Γ] × exp - 1 2 ξ 2 . (6)
Thus, owing to this transformation of measure, an MCMC method can be set up to sample from ρ Ξ on R n(n+1)/2 where it does not face the symmetry and positive-definiteness constraint.

Change of variables

The first factor in the expression of ρ Ξ represents the information extracted from the data and the second factor represents the prior information. The transformation of measure is such that the second factor is standard Gaussian with a covariance matrix equal to the identity matrix. However, as a consequence of the typical ill-posedness of inverse problems, the data can be expected to provide significant information only on certain components of the symmetric positive-definite matrix and leave other components rather uninformed. Correspondingly, the first factor in the expression of ρ Ξ can be expected to concentrate its weight on narrower ranges of values only in certain directions in R n(n+1)/2 and distribute its weight over wider ranges of values in other directions.

To ensure that an MCMC method works well in such a situation, we introduce an additional change of variables defined as a linear transformation involving a Cholesky factor of an approximation to the covariance matrix. This change of variables is a standard approach related to principal component analysis. Specifically, we describe in Sec. 2.5 how an approximation [ C] to the covariance matrix of ρ Ξ can be obtained via a linearization of the inverse problem. Given such an approximation [ C] to the covariance matrix of ρ Ξ , we use its Cholesky factorization into the product of a lower-triangular matrix and its transpose,

[ C] = [L C ][L C ] T , (7)
to define a transformation from R n(n+1)/2 to R n(n+1)/2 given by

q = [L C ] -1 ξ. (8)
Given this transformation, we define the PDF ρ Q on R n(n+1)/2 whose image under this transformation is ρ Ξ on R n(n+1)/2 :

ρ Q q y obs ∝ exp - 1 2 y obs -g X ([L C ]q) 2 [Γ] × exp - 1 2 [L C ]q 2 . ( 9
)
Provided that [ C] is a good approximation to the covariance matrix of ρ Ξ , the covariance matrix of ρ Q is approximately equal to the identity matrix.

Itô-SDE-based MCMC method

The dimension of R n(n+1)/2 grows with the square of the matrix dimension and hence can be high, so that conventional MCMC methods can be inefficient. We address this difficulty by using the Itô-SDE-based MCMC method introduced in [START_REF] Soize | Construction of probability distributions in high dimension using the maximum entropy principle: Applications to stochastic processes, random fields and random matrices[END_REF]. This MCMC method is based on the Itô SDE

   dQ = P dt, dP = -∇ q ψ(Q)dt - 1 2 ζ 0 P dt + ζ 0 dW , (10)
where ζ 0 > 0 is a free parameter and {W (t), t ≥ 0} is a normalized Wiener stochastic process. Under certain regularity conditions, this Itô SDE is such that when completed with an initial condition, its solution {(Q(t), P (t)), t ≥ 0} is ergodic for the so-called Gibbs PDF

ρ (Q,P ) (q, p) ∝ exp -ψ(q) - p 2 2 , (11)
where the ergodicity property means that time averages over trajectories of the solution can be used to approximate probabilistic averages with respect to the Gibbs PDF. With reference to this ergodicity property, the Itô-SDEbased MCMC method is obtained by setting the so-called potential ψ in the Itô SDE in [START_REF] Cottereau | Construction of a probabilistic model for impedance matrices[END_REF] equal to the negative logarithm of the target PDF, that is, the PDF that must be sampled from, and then discretizing this Itô SDE in time. Thus, in the Itô-SDE-based MCMC method, it is the time discretization of an Itô SDE that provides the transition kernel of the Markov chain. Through the gradient of the potential involved in the Itô SDE, the transition kernel is able to leverage the gradient of the target PDF to guide the Markov chain to regions to which the target PDF attaches larger weight. We thus set ψ equal to the negative logarithm of PDF ρ Q defined in (9):

ψ(q) = 1 2 y obs -g X ([L C ]q) 2 [Γ] + 1 2 [L C ]q 2 , (12)
where we can use an equality sign because the multiplicative normalization constant in [START_REF] Arnst | A nonparametric probabilistic model for ground-borne vibrations in buildings[END_REF] need not be known, because the logarithm transforms it into an additive constant that the gradient in [START_REF] Cottereau | Construction of a probabilistic model for impedance matrices[END_REF] then absorbs. Before discretizing the Itô SDE in time, we introduce in it the transformation in [START_REF] Chebli | Experimental validation of a nonparametric probabilistic model of nonhomogeneous uncertainties for dynamical systems[END_REF] as well as the transformation

p = [L C ] T η (13)
to obtain

   dΞ = [ C]Hdt, dH = -∇ ξ φ(Ξ)dt - 1 2 ζ 0 Hdt + ζ 0 [L C ] -T dW , (14)
with φ equal to the negative logarithm of PDF ρ Ξ defined in (6):

φ(ξ) = 1 2 y obs -g [X (ξ)] 2 [Γ] + 1 2 ξ 2 . ( 15
)
The purpose of ( 13) is to arrive in (14) at a system of equations that is a Hamiltonian equation for the Hamiltonian h(ξ, η) = φ(ξ) + 1 2 η T [ C]η, in which there is still the additional friction term and the additional stochastic forcing term. Owing to the transformations in ( 8) and ( 13), the Itô-SDEbased MCMC method will sample directly from ρ Ξ , while still benefitting from the change of variables, which still manifests itself in ( 16) by the presence of the matrices [ C] and [L C ] -T . The effect of the presence of these matrices is that of a normalization of the dynamics of the Itô SDE. One way of seeing this is as follows: injecting the second equation into the first equation in [START_REF] Guilleminot | Stochastic model and generator for random fields with symmetry properties: application to the mesoscopic modeling of elastic random media[END_REF] reveals that the square roots of the eigenvalues of [ C]D ξ ∇ ξ φ are indicative of the time scales involved in the dynamics of the Itô SDE; considering that the inverse of the Hessian of the negative logarithm of a PDF provides an approximation to the covariance matrix of this PDF-an approximation that is exact for a Gaussian PDF-the effect of the presence of the matrices [ C] and [L C ] -T is that of a clustering of the eigenvalues of [ C]D ξ ∇ ξ φ close to 1 and thus that of a normalization of the time scales involved in the dynamics of the Itô SDE.

To discretize in time, we use a Störmer-Verlet method [START_REF] Hairer | Geometric numerical integration illustrated by the Störmer-Verlet method[END_REF] that had already been used for the Itô-SDE-based MCMC method in [START_REF] Soize | Data-driven probability concentration and sampling on manifold[END_REF], leading to

               Ξ ( + 1 2 ) = Ξ ( ) + t 2 [ C]H ( ) , H ( +1) = 1 -b 1 + b H ( ) - t 1 + b ∇ ξ φ(Ξ ( + 1 2 ) ) + √ ζ 0 1 + b [L C ] -T W ( +1) , Ξ ( +1) = Ξ ( + 1 2 ) + t 2 [ C]H ( +1) , (16)
where

b = t ζ 0 /4, (17)
in which t denotes the time step and W (1) , W (2) , . . . are mutually independent centered multivariate Gaussian random variables with covariance matrix t[I], in which [I] is here the n(n + 1)/2-dimensional identity matrix. The change of variables and its effect of a normalization of the dynamics of the Itô SDE facilitate the choice of the time step and the control of the error entailed by the time discretization. The end result of the change of variables is that it allows the Markov chain to make faster progress in directions that are less informed by the data, without being constrained by a stabilityrelated or error-related requirement of slower progress in directions that are more informed by the data. When completed with an initial condition, the time-discrete equations in ( 16) define a Markov chain {(Ξ ( ) , H ( ) ), ∈ N}. The transformation of measure associates the chain of random matrices {[X ( ) ], ∈ N} with the component chain of random variables {Ξ ( ) , ∈ N}:

[X ( ) ] = [X (Ξ ( ) )], (18)
Owing to the transformation of measure, owing to the ergodicity property of the solution to the Itô SDE, and owing to the choice of the potential, but up to an error entailed by the time discretization, trajectories of this chain of random matrices can be construed as samples drawn from the posterior. The gradient of the potential in the transition kernel is the sum of two terms, with the first term stemming from the weighted distance that represents the model-data misfit and the second term stemming from the prior:

∇ ξ φ(ξ) = [J(ξ)] T [Γ] g [X (ξ)] -y obs + ξ; (19)
here, [J(ξ)] denotes the gradient of the composition of X and g, which, with the chain rule, can be written, in terms of entries, as follows:

J ij (ξ) = ∂g i ∂ξ j [X (ξ)] = n =1 n =1 ∂g i ∂X [X (ξ)] ∂X ∂ξ j (ξ). ( 20
)
Thus, for the proposed methodology to work well, the forward problem must afford an efficient computation of not only the observables but also their gradient with respect to the matrix. For the maximum-entropy PDF, the expression of the gradient of the transformation can be readily obtained, as we show in Sec. 3. In the problem of a data-driven updating of a reducedorder model in computational structural dynamics, we use an adjoint-based computation of the gradient, as we discuss in Sec. 4.

Obtaining [ C] from a linearization of the inverse problem

We will conclude the section with discussing how a linearization of the inverse problem allows an approximation [ C] to the covariance matrix of ρ Ξ to be obtained. Let the composition of the transformation X and the forward problem g be linearized about a reference value ξ 0 :

g [X (ξ)] ≈ g [X (ξ 0 )] + [J(ξ 0 )](ξ -ξ 0 ); (21)
here, among other choices, ξ 0 can be chosen equal to a mode of the prior or a mode of the posterior, such as a value that maximizes ρ Ξ , which can be computed with an optimization method such as a quasi-Newton method:

[J(ξ ( ) )] T [Γ][J(ξ ( ) )] + [I] ξ ( +1) = -[J(ξ ( ) )] T [Γ] g [X (ξ ( ) )] -[J(ξ ( ) )]ξ ( ) -y obs . (22)
It can be seen by injecting ( 21) into ( 6) that the approximation of the composition of the transformation X and the forward problem g by its linearized expression leads to the approximation of the weighted distance that represents the model-data misfit by a quadratic expression and ultimately the approximation of ρ Ξ by a Gaussian PDF as follows:

ρ Ξ (ξ|y obs ) ≈ ρΞ (ξ|y obs ) ∝ exp - 1 2 ξ -ξ 2 [ C] -1 , (23)
with mean value given by

ξ = [J(ξ 0 )] T [Γ][J(ξ 0 )] + [I] -1 [J(ξ 0 )] T [Γ] y obs -g [X (ξ 0 )] + [J(ξ 0 )]ξ 0 (24)
and covariance matrix given by

[ C] = [J(ξ 0 )] T [Γ][J(ξ 0 )] + [I] -1 . ( 25
)
It is this matrix that can be used in the change of variables as the approximation to the covariance matrix of ρ Ξ .

With reference to the change of variables, it is interesting to note that as in [START_REF] Flath | Fast algorithms for Bayesian uncertainty quantification in large-scale linear inverse problems based on low-rank partial Hessian approximations[END_REF], insight into the information provided by the data can be gained from an eigendecomposition. Indeed, with [Λ r ] the diagonal matrix with the r dominant eigenvalues and [V r ] collecting in its columns the corresponding eigenvectors of

[J(ξ 0 )] T [Γ][J(ξ 0 )], [J(ξ 0 )] T [Γ][J(ξ 0 )] ≈ [V r ][Λ r ][V r ] T , [V r ][V r ] T = [I], (26)
the Sherman-Morrison-Woodbury formula provides the approximation

[ C] = [J(ξ 0 )] T [Γ][J(ξ 0 )] + [I] -1 ≈ [I] -[V r ][D r ][V r ] T , (27) with [
D r ] = [Λ r ] + [I] -1 [Λ r ], (28)
thus suggesting that the eigenvectors corresponding to the dominant eigenvalues of [J(ξ 0 )] T [Γ][J(ξ 0 )] are indicative of the directions in R n(n+1)/2 that the data provide significant information on. As a consequence of the typical ill-posedness of inverse problems, the spectrum of [J(ξ 0 )] T [Γ][J(ξ 0 )] can be expected to decay rapidly. Later in this paper, we will also find this result helpful for the identification of hyperparameters of the posterior.

Using a maximum-entropy PDF as the prior

We now describe the use of the proposed methodology with a prior taken from the family of PDFs for symmetric positive-definite random matrices constructed by means of the maximum-entropy principle in [START_REF] Soize | Nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF][START_REF] Soize | Non-Gaussian positive-definite matrix-valued random fields for elliptic stochastic partial differential operators[END_REF].

Concise review of the family of maximum-entropy PDFs

The family of PDFs for symmetric positive-definite random matrices constructed by means of the maximum-entropy principle in [START_REF] Soize | Nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF][START_REF] Soize | Non-Gaussian positive-definite matrix-valued random fields for elliptic stochastic partial differential operators[END_REF] has two hyperparameters: the mean value and a dispersion parameter. A symmetric positive-definite random matrix, say, [X], has such a maximum-entropy PDF with mean value [X], which must be a symmetric positive definite matrix, and dispersion parameter δ if given the Cholesky factorization of [X] into the product of a lower-triangular matrix and its transpose,

[X] = [L X ][L X ] T , (29)
the normalized symmetric positive-definite random matrix [G] in

[X] = [L X ][G][L X ] T (30) has a PDF ρ [G]
that is obtained by maximizing entropy, max

ρ [G] - M + n (R) log ρ [G] ([G])ρ [G] ([G])dG, (31)
under the following constraints:

M + n (R) ρ [G] ([G])dG = 1, (32)
M + n (R) [G]ρ [G] ([G])dG = [I], (33)
M + n (R) log(det[G])ρ [G] ([G])dG < +∞. ( 34
)
The solution of the constrained optimization problem (31) with (32)-(34) by means of the method of Lagrange multipliers leads to

ρ [G] ([G]) ∝ 1 M + n (R) ([G])(det[G]) n+1 2δ 2 + 1-n 2 exp - n + 1 2δ 2 tr[G] , (35)
where

1 M + n (R) ([G]) = 1 if [G]
is an n-dimensional symmetric positive-definite matrix with real entries and 1

M + n (R) ([G]) = 0 otherwise, det[G] is the determi- nant of [G],

and tr[G] is the trace of [G]

. The Lagrange multiplier associated with the constraint in (34) results in the presence of the hyperparameter in (35), namely, the dispersion parameter, δ, defined such that integrate to 1 over M + n (R) so that samples of [G] and therefore also samples of [X] are symmetric and positive definite. The constraint in (33) imposes that the mean value of [G] be [I] and therefore the mean value of [X] be [X]. The constraint in (34) serves to ensure the existence of the secondorder moment of the inverse of [G] and therefore of [X], as required when such symmetric positive-definite random matrices are used to represent uncertain symmetric positive-definite matrices of computational models, see, for instance, [START_REF] Soize | Nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF][START_REF] Soize | Non-Gaussian positive-definite matrix-valued random fields for elliptic stochastic partial differential operators[END_REF][START_REF] Soize | A comprehensive overview of a non-parametric probabilistic approach of model uncertainties for predictive models in structural dynamics[END_REF][START_REF] Chebli | Experimental validation of a nonparametric probabilistic model of nonhomogeneous uncertainties for dynamical systems[END_REF][START_REF] Arnst | A nonparametric probabilistic model for ground-borne vibrations in buildings[END_REF][START_REF] Cottereau | Construction of a probabilistic model for impedance matrices[END_REF][START_REF] Mignolet | Nonparametric stochastic modeling of linear systems with prescribed variance of several natural frequencies[END_REF][START_REF] Soize | Time-domain formulation in computational dynamics for linear viscoelastic media with model uncertainties and stochastic excitation[END_REF][START_REF] Batou | Rigid multibody system dynamics with uncertain rigid bodies[END_REF][START_REF] Guilleminot | Stochastic model and generator for random fields with symmetry properties: application to the mesoscopic modeling of elastic random media[END_REF].

δ 2 = M + n (R) [G] -[I] 2 F ρ [G] ([G])dG [I] 2 F , (36)
The PDF ρ [G] possesses a special structure that can be exploited to set up a direct simulation method based on the representation of [G] as the product of a lower-triangular random matrix and its transpose,

[G] = [L][L] T , (37)
with the entries {L , n ≥ ≥ ≥ 1} of [L] mutually statistically independent random variables distributed as follows:

for = :

L = √ 2 δ √ n + 1 Y , with Y ∼ Γ(α , 1), α = n + 1 2δ 2 + 1 - 2 ; (38)
for > :

L ∼ N (0, σ 2 ), σ = δ √ n + 1 ; (39)
here, we denote by Y ∼ Γ(α , 1) that Y is a gamma random variable with shape parameter α and scale parameter 1, and we denote by L ∼ N (0, σ 2 ) that L is a Gaussian random variable with mean 0 and variance σ 2 .

Using a maximum-entropy PDF as the prior in the proposed methodology

The use of a maximum-entropy PDF with mean value [X] and dispersion parameter δ as the prior in the proposed methodology amounts to setting

ρ 0 ([X]) ∝ ρ [G] ([L X ] -1 [X][L X ] -T ). ( 40
)
The transformation of measure requires a representation of the prior on the set of symmetric positive-definite matrices M + n (R) as the image under a transformation X of a standard Gaussian PDF on the linear space R n(n+1)/2 . Such a transformation X immediately follows from the representation used in the direct simulation method in (37) with ( 38)-(39) by setting

[X (ξ)] = [L X ][L(ξ)][L(ξ)] T [L X ] T , ( 41
)
where L is the transformation from the linear space R n(n+1)/2 to the set of lower-triangular matrices with positive diagonal entries L L ,

L : R n(n+1)/2 → L L : ξ → [L] = [L(ξ)], (42)
such that, with the so-called half-vectorization ξ = (ξ 11 , . . . , ξ n1 , ξ 22 , . . . , ξ n2 , . . . , ξ (n-1)(n-1) , ξ n(n-1) , ξ nn ),

the entries {L (ξ), n ≥ ≥ ≥ 1} of [L(ξ)] are as follows:

for = :

L (ξ) = √ 2 δ √ n + 1 √ y , y = c -1 Γ c N (ξ ; 0, 1); α , 1 ; (44)
for > :

L (ξ) = δ √ n + 1 ξ ; (45)
here, c Γ (•; α , 1) is the cumulative distribution function of a gamma random variable with shape parameter α and scale parameter 1, and c N (•; 0, 1) is the cumulative distribution function of a Gaussian random variable with mean 0 and variance 1. It is owing to the isoprobability transform that the gamma PDF with with shape parameter α and scale parameter 1 is the image under the transformation c -1 Γ c N (•; 0, 1); α , 1 of a standard Gaussian PDF. The approximation of the covariance matrix in the change of variables and the computation of the gradient in the Itô-SDE-based MCMC method require the gradient of the transformation X . From (41) with ( 44)-(45), it immediately follows that the components of the gradient of the transformation X are given by

∂ ξ [X (ξ)] = [L X ] ∂ ξ [L(ξ)][L(ξ)] T + [L(ξ)]∂ ξ [L(ξ)] T [L X ] T , (46)
where

∂ ξ [L(ξ)] = ∂ ξ L (ξ) e ⊗ e , (47)
with for = :

∂ ξ L (ξ) = √ 2 δ √ n + 1 1 2 1 √ y ρ N (ξ ; 0, 1) ρ Γ (y ; α , 1) ; ( 48
)
for > :

∂ ξ L (ξ) = δ √ n + 1 ; (49)
here, ρ Γ (•; α , 1) is the PDF of a gamma random variable with shape parameter α and scale parameter 1, and ρ N (•; 0, 1) is the standard Gaussian PDF.

Data-driven update of a reduced-order model

We now describe in the domain of computational structural dynamics a problem of a data-driven updating of a reduced-order model.

Reduced-order model

We consider a deterministic finite element model for the linear dynamical behavior of a dissipative structure about a static equilibrium in the frequency domain. Let this model be projected onto a reduced-order basis to obtain a reduced-order model of the following form:

(-ω 2 [M ] + iω[D] + [K])z(ω) = [Φ] T f (ω), u(ω) = [Φ]z(ω); (50)
here, ω denotes the (circular) frequency; [Φ] the matrix whose columns collect the reduced-order basis; [M ], [D], and [K] the reduced mass, damping, and stiffness matrices; f (ω) the vector of external forces; z(ω) the vector of n reduced degrees of freedom; and u(ω) the vector of displacement degrees of freedom. Although extensions of the proposed methodology to other types of dissipation and boundary conditions can be readily obtained, we assume that the dissipation and the boundary conditions are such that [D] and [K] are symmetric and positive definite.

Data-driven update of the reduced-order model

In computational structural dynamics, data-driven updating most often involves adjusting values of mechanical properties, such as stiffness, dissipation, and inertial properties, so that predictions better fit data relevant to the dynamical behavior of the structure. However, it is not such an adjustment of values of mechanical properties that we pursue. Rather, we seek to adjust directly the reduced matrices of the reduced-order model, and we seek to do so with a probabilistic approach that allows us to obtain not only improved predictions but also a probabilistic representation of the model uncertainties.

Thus, we embed the reduced-order model in (50) in a class of reducedorder models indexed by the values assigned to the reduced matrices, now treated as (generalized) parameters of the reduced-order model:

(-ω 2 [M ] + iω[D] + [K])z(ω; [M ], [D], [K]) = [Φ] T f (ω), u(ω; [M ], [D], [K]) = [Φ]z(ω; [M ], [D], [K]), (51)
and we cast the data-driven update as a Bayesian identification of the (generalized) parameters of the reduced-order model, that is, the reduced matrices.

The specification of the posterior on the set of possible values of the reduced matrices then provides the updated probabilistic reduced-order model.

Using a maximum-entropy PDF as the prior

We use the product of maximum-entropy PDFs with mean values [M ], [D], and [K] and dispersion parameters δ M , δ D , and δ K as the prior:

ρ 0 ([M ], [D], [K]) ∝ A=M,D,K ρ [G A ] ([L A ] -1 [A][L A ] -T ); (52

Interpretation of the prior as a nonparametric probabilistic model

The probabilistic reduced-order model obtained by specifying the prior on the set of possible values of the reduced matrices can be viewed as a nonparametric probabilistic model constructed as in [START_REF] Soize | Nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF]. Indeed, within the present problem setting, the nonparametric probabilistic approach of [START_REF] Soize | Nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF] involves representing the reduced matrices in the reduced-order model in (50) with random matrices [M ], [D], and [K] that are statistically independent and have maximum-entropy PDFs with mean values [M ], [D], and [K] and dispersion parameters δ M , δ D , and δ K , respectively, thus leading to a probabilistic reduced-order model of the following form:

(-ω 2 [M ] + iω[D] + [K])Z(ω) = [Φ] T f (ω), U (ω) = [Φ]Z(ω). ( 53
)
Here, and with reference to Sec.

[A] = [L A ][L A ] T , A = M, D, K, (54)
the normalized symmetric positive-definite matrices

[G M ], [G D ], and [G K ] in [A] = [L A ][G A ][L A ] T , A = M, D, K, (55)
δ 2 A = M + n (R) [G A ] -[I] 2 F ρ [G A ] ([G A ])dG A [I] 2 F , A = M, D, K. (56)
In the entropy maximization, the constraint in (34) serves to ensure the existence of the second-order moments of the inverses of [G M ], [G D ], and [G K ] and therefore of [M ], [D], and [K], thus ensuring that the probabilistic reduced-order model has a second-order solution. In conclusion, the joint PDF of the random matrices [M ], [D], and [K] is the PDF in (52), thus confirming that the probabilistic reduced-order model obtained by specifying the prior on the set of possible values of the reduced matrices is a nonparametric probabilistic model. In the context of the nonparametric probabilistic approach, several methods have been proposed for the identification of the dispersion parameters. In [START_REF] Soize | A comprehensive overview of a non-parametric probabilistic approach of model uncertainties for predictive models in structural dynamics[END_REF][START_REF] Chebli | Experimental validation of a nonparametric probabilistic model of nonhomogeneous uncertainties for dynamical systems[END_REF][START_REF] Arnst | A nonparametric probabilistic model for ground-borne vibrations in buildings[END_REF], a method was proposed based on an identification criterion that gauges the agreement between data and the predictions of the nonparametric probabilistic model. In this method, the data are plotted on the same axes as high-probability confidence regions for the predictions, and a combination of ideally small values of the dispersion parameters is sought for which the data are mostly contained within the high-probability confidence regions for the predictions. In [START_REF] Soize | Probabilistic model identification of uncertainties in computational models for dynamical systems and experimental validation[END_REF][START_REF] Arnst | Inversion of probabilistic structural models using measured transfer functions[END_REF], approaches were proposed which allow the identification to be cast as an optimization problem.

Dynamic test data

Although extensions to other types of data can be readily obtained, we present a data-driven updating of the reduced-order model based on using frequency response functions (FRFs) sampled at frequencies in a frequency band of interest. Specifically, we assume that the data consist of a finite number, say, d, of FRFs sampled at a finite number, say, m, of frequencies ω 1 , . . . , ω m in the frequency band of interest:

y obs = h obs (ω 1 ), . . . , h obs (ω m ) , (57)
with h obs (ω k ) = h obs 1 (ω k ), . . . , h obs d (ω k ) . (58)

Forward problem

We thus define the forward problem as the mapping that represents the use of the reduced-order model to associate to a value of the reduced matrices a corresponding value of the observables:

y = g([M ], [D], [K]) = h(ω 1 ; [M ], [D], [K]), . . . , h(ω m ; [M ], [D], [K]) , (59)
where

h(ω k ; [M ], [D], [K]) = h 1 (ω k ; [M ], [D], [K]), . . . , h d (ω k ; [M ], [D], [K]) , (60)
with

h j (ω k ; [M ], [D], [K]) = c T j [Φ](-ω 2 k [M ] + iω k [D] + [K]) -1 [Φ] T b j , (61)
in which for the j-th FRF, the so-called input-shape vector b j and the socalled output-shape vector c j link the finite element degrees of freedom to the position and the orientation of the actuator and the sensor, respectively.

Bayesian identification of the reduced matrices

Casting the data-driven update as a Bayesian identification of the reduced matrices with the maximum-entropy PDF as the prior leads to the posterior

ρ ([M ],[D],[K]) [M ], [D], [K] y obs ∝ f y obs [M ], [D], [K] × ρ 0 ([M ], [D], [K]), ( 62
)
in which f is the likelihood function. With reference to Sec. 2, we use for f an expression of the form of an exponential of minus a weighted distance that represents the model-data misfit. Although extensions to other types of weighted distance can be readily obtained, we use a logarithmic least-squares model-data misfit, which has the benefits of being as sensitive to resonances as anti-resonances and being insensitive to phase errors:

f y obs [M ], [D], [K] ∝ exp -λ 1 2 d j=1 m k=1 µ jk log |h obs j (ω k )| -log |h j (ω k ; [M ], [D], [K])| 2 . ( 63
)
A comparison with the expression in (3) indicates that the weights λ and µ jk , 1 ≤ j ≤ d, 1 ≤ k ≤ m, can be interpreted as the use of a diagonal covariance matrix (the matrix [Γ] -1 in (3)) whose diagonal elements are the inverses of the products λµ jk . Because f attaches greater weight to values of the reduced matrices that correspond to predictions that fit the data better, the updated probabilistic reduced-order model with the posterior ρ ([M ],[D],[K]) specified on the set of possible values of the reduced matrices can be expected to be specifically tailored to provide improved predictions. The posterior is parameterized by the dispersion parameters δ M , δ D , and δ K and the weights λ and µ jk , 1 ≤ j ≤ d, 1 ≤ k ≤ m, and we comment on the identification of these hyperparameters later in Sec. 4.9.

Numerical implementation of the Bayesian inversion

Because the prior is a product PDF, it suffices to construct separate representations of the maximum-entropy PDFs with mean values [M ], [D], and [K] and dispersion parameters δ M , δ D , and δ K as images under transformations of Gaussian PDFs. We construct these transformations as in (41) with ( 42)-(45) and denote the transformations thus obtained by M, D, and K:

A : R n(n+1)/2 → M + n (R) : ξ A → [A] = [A(ξ A )], A = M, D, K. ( 64
)
Given these transformations, we define the PDF ρ Ξ on R n(n+1)/2 ×R n(n+1)/2 × R n(n+1)/2 whose image under the compound transformation is the posterior:

ρ Ξ ξ y obs ∝ exp -λ 1 2 d j=1 m k=1 µ jk log |h obs j (ω k )| -log h j (ω k ; ξ)| 2 × exp - 1 2 ξ 2 , ( 65
)
in which

h j (ω k ; ξ) = h j ω k ; [M(ξ M )], [D(ξ D )], [K(ξ K )] , (66)
with ξ = (ξ M , ξ D , ξ K ). ( 67
)
The linearization of the composition of the transformations M, D, and K and the forward problem g about a reference value ξ 0 leads as in ( 21) and ( 24) with ( 25) to an approximation [ C] to the covariance matrix of ρ Ξ :

[ C] = λ d j=1 m k=1 µ jk j j (ω k ; ξ 0 ) ⊗ j j (ω k ; ξ 0 ) + [I] -1 , (68)
in which

j j (ω k ; ξ) = ∇ ξ log |h j (ω k ; ξ)|. ( 69
)
We use for ξ 0 a mode of the posterior, namely, a value that maximizes ρ Ξ , which we compute with a quasi-Newton method:

λ d j=1 m k=1 µ jk j j (ω k ; ξ ( ) ) ⊗ j j (ω k ; ξ ( ) ) + [I] ξ ( +1) = -λ d j=1 m k=1 µ jk j j (ω k ; ξ ( ) ) log h j (ω k ; ξ ( ) )| -j j (ω k ; ξ ( ) ) • ξ ( ) -log |h obs j (ω k )| . ( 70
)
The Itô-SDE-based MCMC method then follows from setting the potential φ equal to the negative logarithm of the PDF ρ Ξ defined in (65), injecting this potential φ into the Itô SDE in [START_REF] Guilleminot | Stochastic model and generator for random fields with symmetry properties: application to the mesoscopic modeling of elastic random media[END_REF], and discretizing this Itô SDE in time, thus leading to the time-discrete equations in [START_REF] Arnst | Itô-SDE MCMC method for Bayesian characterization of errors associated with data limitations in stochastic expansion methods for uncertainty quantification[END_REF] with

φ(ξ) = λ 1 2 d j=1 m k=1 µ jk log |h obs j (ω k )| -log h j (ω k ; ξ)| 2 + 1 2 ξ 2 (71)
and

∇ ξ φ(ξ) = λ d j=1 m k=1 µ jk j j (ω k ; ξ) log h j (ω k ; ξ)| -log |h obs j (ω k )| 2 + ξ. ( 72
)
When completed with an initial condition, these time-discrete equations define a Markov chain {(Ξ ( ) , H ( ) ), ∈ N}. The transformation of measure associates the chain of random matrices

{[M ( ) ], [D ( ) ], [K ( ) ],
∈ N} with the component chain of random variables {Ξ ( ) , ∈ N}:

[M ( ) ] = [M(Ξ ( ) M )], [D ( ) ] = [D(Ξ ( ) D )], [K ( ) ] = [K(Ξ ( ) K )]. ( 73
)
As for [START_REF] Arnst | Inversion of probabilistic structural models using measured transfer functions[END_REF], trajectories of this chain of random matrices can be construed as samples drawn from the posterior and hence used as samples from the random matrices in the Monte Carlo simulation with the updated probabilistic reduced-order model.

Identification of the hyperparameters

The posterior is parameterized by the dispersion parameters δ M , δ D , and δ K and the weights λ and µ jk , 1 ≤ j ≤ d, 1 ≤ k ≤ m. These hyperparameters control the dispersion of the prior and the weights of the weighted distance that represents the model-data misfit in the likelihood function. The smaller the values of δ M , δ D , and δ K , the more the posterior concentrates its weight predominantly in the vicinity of the values of the reduced matrices in the reduced-order model in (50). The larger the values of the weights λ and µ jk , 1 ≤ j ≤ d, 1 ≤ k ≤ m, the more the posterior concentrates its weight predominantly in the vicinity of values of the reduced matrices that minimize the weighted distance that represents the model-data misfit. Thus, especially if there are no values of the reduced matrices that result in a perfect representation of the dynamical behavior of the structure and some amount of dispersion should persist to serve as a representation of the model uncertainties, the dispersion parameters δ M , δ D , and δ K cannot be too small and the weights λ and µ jk , 1 ≤ j ≤ d, 1 ≤ k ≤ m cannot be too large.

If the data are perturbed by noise, estimates of the significance of the noise can play a role in the identification of the weights λ and µ jk , 1 ≤ j ≤ d, 1 ≤ k ≤ m. For instance, smaller weights µ jk can be given to FRF samples that are perturbed more significantly by noise.

In a manner that builds on the method that had already been proposed for the identification of the dispersion parameters in the context of the nonparametric probabilistic approach, we use for the identification of the hyperparameters of the posterior a method based on an identification criterion. In this method, we plot the data on the same axes as high-probability confidence regions for the predictions of the updated probabilistic reduced-order model, and we seek a combination of ideally small values of the dispersion parameters δ M , δ D , and δ K and ideally large values of the weights λ and µ jk , 1 ≤ j ≤ 1 ≤ k ≤ m for which the data are mostly contained within the high-probability confidence regions for the predictions.

Adjoint-based computation of the requisite gradients

The approximation of the covariance matrix in (68) in the change of variables requires the computation of gradients as in (69), which are also required for the computation of the gradient of the potential in (72) in the Itô-SDE-based MCMC method. We use for these computations of these gradients an adjoint-based approach that requires for each frequency and for each FRF one direct solve of the reduced-order model and three adjoint solves associated with the reduced-order model. While the computations of the partial derivatives with respect to components of ξ M and ξ D can be carried out analogously, we provide below details about the computation of the partial derivatives with respect to components ξ K of ξ K .

One way of obtaining the adjoint-based computation starts from writing the partial derivative of log |h j (ω

k ; [M ], [D], [K])| with respect to ξ K as ∂ ξ K log |h j (ω k ; [M ], [D], [K])| = 1 |h j (ω k ; [M ], [D], [K])| ∂ ξ K |h j (ω k ; [M ], [D], [K])| = 1 |h j (ω k ; [M ], [D], [K])| 2 Re ∂ ξ K h j (ω k ; [M ], [D], [K])h j (ω k ; [M ], [D], [K]) , (74)
in which we denote by Re(a + bi) = a the real part and by a + bi = a -bi the complex conjugate of a complex number a + bi. The partial derivative of

h j (ω k ; [M ], [D], [K]
) with respect to ξ K reads as

∂ ξ K h j (ω k ; [M ], [D], [K]) = -c T j [Φ](-ω 2 k [M ] + iω k [D] + [K]) -1 (∂ ξ K [K])(-ω 2 k [M ] + iω k [D] + [K]) -1 [Φ] T b j . (75)
Hence, with the chain rule and with (46) with (47)-(49), the partial deriva-tives of log |h j (ω k ; ξ)| with respect to the components ξ K of ξ K read as

∂ ξ K log |h j (ω k ; ξ)| = -1 |h j (ω k ;ξ)| 2 Re b j (ω k ;ξ)∂ ξ K L (ξ K )f j (ω k ;ξ)+c j (ω k ;ξ)∂ ξ K L (ξ K )e j (ω k ;ξ) , (76) in which q j (ω k ; ξ) = -ω 2 k [M(ξ M )] + iω k [D(ξ D )] + [K(ξ K )] -1 [Φ] T b j , (direct), ( 77
) b j (ω k ; ξ) = [L K ] T q j (ω k ; ξ), (78)
c j (ω k ; ξ) = [L(ξ K )] T b j (ω k ; ξ), (79)
d j (ω k ; ξ) = -ω 2 k [M(ξ M )] + iω k [D(ξ D )] + [K(ξ K )] -1 [Φ] T c j h j (ω k ; ξ), (adjoint), (80)
e j (ω k ; ξ) = [L K ] T d j (ω k ; ξ), (81)
f j (ω k ; ξ) = [L(ξ K )] T e j (ω k ; ξ). (82)

Numerical illustration

We now present a numerical illustration of a data-driven updating of a reduced-order model as described in Sec. 4.

Problem setting and data

We considered the dynamical behavior of a clamped-clamped beam with a length of 10 m, a width of 1.2 m, and a height of 1.6 m. We assumed its constitutive behavior to be isotropic, linear, and dissipative. We assumed a Young modulus of 10 GPa, a Poisson ratio of 0.2, a mass density of 1600 kg m -3 , and a modal damping ratio of 0.01. The frequency band of analysis is [0, 1200] Hz.

We generated (synthetic) data with a three-dimensional finite element model with a mesh of 100 × 12 × 16 three-dimensional 8-node elements of equal size. For the purpose of the data-driven updating of the reduced-order model, we used the three-dimensional finite element model to compute the 15 FRFs that link the force driven by an actuator oriented along the transverse direction (along the height of the beam) and located at 4.75 m (from the leftmost cross section) to the displacements gauged by sensors oriented along the transverse direction and located at equally spaced locations from 0.588 m to 4.116 m with a step of 0.588 m and from 5.294 m to 9.412 m with a step of 0.588 m. This sensor grid corresponds to setting up 16 equally spaced sensors and then leaving out the one closest to the actuator. In addition, for the purpose of enabling a cross-validation, that is, to evaluate the predictive capabilities of the updated probabilistic reduced-order model, we used the three-dimensional finite element model to compute the FRF that links the force driven by an actuator oriented along the transverse direction and located at 5.6 m to the displacement at a sensor oriented along the transverse direction and located at 8.2 m. We sampled all these FRFs with a step of 0.1 Hz from 0 to 200 Hz and with a step of 1 Hz from 201 Hz to 1200 Hz.

Reduced-order model

We set up a one-dimensional finite element model based on the Timoshenko beam theory with a mesh of 200 elements of equal size. We used a shear correction factor of 5/6. We purposely set up such a simple deterministic model in order for there to clearly be modeling errors. Figure 2 compares the data with the corresponding FRFs predicted by the deterministic model. We can observe that at lower frequencies, the FRFs predicted by the deterministic model fit the data well; however, at higher frequencies, there are significant differences, which can be explained as consequences of modeling errors entailed by the use of the Timoshenko beam theory, whose representation of the mechanical behavior involves simplifying hypotheses, such as the rigidity of the cross section.

[Figure 1

about here.]

We projected the deterministic model onto a reduced-order basis constituted of the dynamical eigenmodes associated with the n = 50 lowest eigenfrequencies of the deterministic model. The 50-th eigenpair corresponds to an eigenfrequency of about 2355 Hz (Fig. 1).

Using a maximum-entropy PDF as the prior

[Figure 2 about here.]

[Figure 3 about here.]

We built a nonparametric probabilistic model by representing the reduced stiffness matrix of the reduced-order model with a random matrix and leaving the reduced mass and damping matrices deterministic. The nonparametric probabilistic model thus obtained has one hyperparameter, namely, the dispersion parameter, δ K , which must be identified. To do so, we applied the aforementioned method based on an identification criterion that quantifies the agreement between the data and the predictions of the nonparametric probabilistic model. We set the value of the dispersion parameter equal to the smallest value that suffices for the data to be mostly contained within high-probability confidence regions for the predictions of the nonparametric probabilistic model. This approach led us to set the dispersion parameter equal to δK = 0.25. Figures 2 and3 compare the data with the predictions of the nonparametric probabilistic model for δ K = 0.05 and δK = 0.25, respectively. All results obtained with the nonparametric probabilistic model were obtained with 5000 Monte Carlo samples.

Dynamic test data

The predictions of the deterministic model and those of the nonparametric probabilistic model fit the data more closely at low frequencies than at higher frequencies (Fig. 3). On the one hand, the closer fit at lower frequencies suggests the feasibility of adjusting the reduced-order model so as to improve the predictions at lower frequencies. On the other hand, the looser fit at higher frequencies suggests that the identification of the value of δK = 0.25 followed mostly from the need to increase this value sufficiently for the data at higher frequencies to be mostly contained within the high-probability confidence regions for the predictions of the nonparametric probabilistic model, with the data at lower frequencies playing only a lesser role. These considerations motivated us to seek to update the reduced-order model based on using as data the FRFs sampled at lower frequencies so as to improve the predictions at lower frequencies. Specifically, we used as data the 15 FRFs in the frequency range between 0 and 900 Hz, sampled with the aforementioned steps of 0.1 Hz below 200 Hz and 1 Hz above 200 Hz.

Because the reduced stiffness matrix is defined via a projection onto dynamical eigenmodes, the data must be sufficiently rich to allow modal contributions to be distinguished. We used the modal-analysis criterion (MAC) to verify that the 15 sensor locations suffice to distinguish modes that contribute significantly to predictions in the frequency range between 0 and 900 Hz.

Bayesian identification of the reduced stiffness matrix

Because the dimension of the reduced stiffness matrix is n = 50, casting the data-driven updating as a Bayesian identification of the reduced stiffness matrix is a Bayesian inverse problem of dimension n(n + 1)/2 = 1275.

In the formulation of the likelihood function, we set µ jk = 2π × 0.1 for 1 ≤ k ≤ 2000 and µ jk = 2π × 1 for 2001 ≤ k ≤ m = 3001 for 1 ≤ j ≤ µ = 15, thus assigning equal weight to all FRFs but compensating for the finer sampling below 200 Hz than above 200 Hz. Hence, the posterior has one hyperparameter left, namely, λ. To gain some insight into the range of values of λ of relevance in this particular problem, we carried out a preliminary study in which we linearized the inverse problem about the mode ξ = 0 of the prior and then computed as in (68) for several values of λ approximations to the covariance matrix of ρ Ξ . This preliminary study indicated that for the model-data misfit to result in a meaningful reduction in the variance of certain components, the relevant range of values of λ should be between about 1 and 100. In the following, we first provide results for λ = 12.5 and then address the identification of the hyperparameter λ.

Numerical implementation of the Bayesian inversion

[Figure 4 about here.]

For the change of variables, we used the quasi-Newton method as in (70) to determine a mode of the posterior to serve as reference value ξ 0 (Fig. 4(a)). We linearized the inverse problem about ξ 0 to obtain first the sensitivity matrix d j=1 m k=1 µ jk j j (ω k ; ξ 0 ) ⊗ j j (ω k ; ξ 0 ) and then as in (68) an approximation [ C] to the covariance matrix of ρ Ξ . The rapid decay of the spectrum of the sensitivity matrix (Fig. 4(b)), as well as certain diagonal entries of [ C] being significantly reduced as compared with their unit values under the prior and other diagonal entries of [ C] being rather unchanged (Fig. 4(c)), indicate that the data provide information on certain components but little or no information on others, thus highlighting the importance of the change of variables. The components that undergo the strongest reduction in variance are those with indices such as 1, 51, 100, 148, and so forth (Fig. 4(c)); the transformation of measure links these components to diagonal entries of the Cholesky factor of the reduced stiffness matrix that contribute significantly to diagonal entries of the reduced stiffness matrix and hence have a significant effect on the locations of resonances and anti-resonances of the FRFs. Components with indices larger than about 750 undergo no change in variance (Fig. 4(c)); these components are linked to entries of the Cholesky factor of the reduced stiffness matrix that have an effect mostly on the behavior of the FRFs at frequencies higher than 900 Hz and hence have only little or no impact on the model-data misfit.

[Figure 5 about here.] For the Itô-SDE-based MCMC method, we obtained all results to follow with ζ 0 = 4, ∆t = 0.1, and 5000 steps in the Markov chain. Figure 5 shows the trajectories of the 51-st component and the 1000-th component of the component chain of random variables {Ξ ( ) K , 1 ≤ ≤ 5000}. We can observe that the trajectory of the 51-st component, which Fig. 4(c) indicated the data to provide information on, evolves about a central value that is larger than that under the prior and with a level of variability that is reduced as compared with that under the prior; by contrast, the trajectory of the 1000-th component, which Fig. 4(c) indicated the data to provide little or no information on, evolves about a central value and with a level of variability that are unchanged as compared with those under the prior.

[Figure 6 about here.]

We determined as in (73) the trajectory of the chain of random reduced stiffness matrices {K ( ) , 1 ≤ ≤ 5000} corresponding to the component chain of random variables {Ξ ( ) K , 1 ≤ ≤ 5000} and used it as samples from the random reduced stiffness matrix in the Monte Carlo simulation with the updated probabilistic reduced-order model. Figure 6 compares the predictions of the updated probabilistic reduced-order model for δK = 0.25 and λ = 12.5 with the data. We can observe that as compared with the predictions of the nonparametric probabilistic model in Fig. 3, the predictions of the updated probabilistic reduced-order model in Fig. 6 fit the data significantly more closely at the frequencies in the frequency range between 0 and 900 Hz that we addressed in the updating.

Identification of the hyperparameter

[Figure 7 about here.]

To identify the hyperparameter λ, we applied the aforementioned method based on an identification criterion that quantifies the agreement between the data and the predictions of the updated probabilistic reduced-order model. Whereas increasing δ K can be expected to widen the high-probability confidence regions for the predictions, increasing λ can be expected to narrow the high-probability confidence regions for the predictions. We thus sought to set λ equal to the largest value for which the data remain mostly contained within high-probability confidence regions for the predictions. This approach led us to set λ = 37.5. Figure 7 compares the data with the predictions of the updated probabilistic reduced-order model for δK = 0.25 and λ = 37.5.

Cross-validation

[Figure 8 about here.]

As a cross-validation, we used the nonparametric probabilistic model and the updated probabilistic reduced-order model to predict the FRF that links the force driven by the actuator oriented along the transverse direction and located at 5.6 m to the displacement gauged by the sensor oriented along the transverse direction and located at 8.2 m, which we hadn't used yet (Fig. 8). We can observe that the predictions of both probabilistic models agree with the data and that as compared with the predictions of the nonparametric probabilistic model, the predictions of the updated probabilistic reducedorder model fit the data significantly more closely at the frequencies in the frequency range between 0 and 900 Hz that we addressed in the updating.

Random reduced mass and stiffness matrices

[Figure 9 about here.]

We repeated the previous study, representing, this time, however, not only the reduced stiffness matrix but also the reduced mass matrix with a random matrix. Figure 9 presents the cross-validation result of this study. The comparison with Fig. 8 indicates that the representation of the reduced mass matrix with a random matrix does not have a significant benefit in terms of further improving the predictions in this particular problem.

Conclusion

We proposed a probabilistic methodology for data-driven updating of non-Gaussian high-dimensional symmetric positive-definite matrices involved in parameterized computational models. This methodology is based on the application of the Bayes formula to update the prior probability model by using data in order to obtain the posterior for the high-dimensional symmetric positive-definite random matrices. The hyperparameters of the posterior are identified by using an identification criterion that quantifies the agreement between the predictions and the data. Samples from the posterior are generated with an Itô-SDE-based MCMC method that leverages the gradient of the posterior for computational efficiency.

We applied the proposed methodology to a problem of a data-driven updating of a reduced-order model in computational structural dynamics. The proposed methodology allowed a direct identification of the reduced matrices from data with an informative prior. The data-driven update resulted in an updated probabilistic reduced-order model that provided improved predictions as well as a probabilistic representation of model uncertainties. In this application, the proposed methodology could also be seen as a data-driven updating of a nonparametric probabilistic model. Taking into account the difficulty of the problem that has been solved, the proposed method yields very good results.

There are many interesting directions for future work. The proposed methodology could be extended to other types of computational model, to other types of informative prior, to other types of likelihood function by changing the distance that represents the model-data misfit and/or by changing the Gaussian model of the noise on the observations, as well as to other approaches for the identification of the hyperparameters. In particular, it would be interesting to revisit the approaches that have already been proposed for the identification of the dispersion parameters of the prior probability distributions in the nonparametric probabilistic approach. In addition, whereas the proposed methodology has been applied in computational structural dynamics with synthetic data, it would be interesting to explore the application to more complex computational models with experimental data.

Cross-validation: data (green solid line) and 98%-confidence region for the corresponding FRF predicted by the nonparametric probabilistic model for δK = 0.25 (light blue filled region) and the updated probabilistic reduced-order model with δK = 0.25 and λ = 37.5 (blue filled region

  in which • F denotes the Frobenius norm; thus, δ ascertains the dispersion of [G] and therefore also of [X]. The constraint in (32) imposes that ρ[G]

  ) here, [M ], [D], and [K] are the reduced mass, damping, and stiffness matrices of the reduced-order model in (50), and the PDFs ρ [G M ] , ρ [G D ] , and ρ [G K ] are PDFs as in (35) with dispersion parameters δ M , δ D , and δ K , respectively.

  3, by [M ], [D], and [K] having maximumentropy PDFs with mean values [M ], [D], and [K] and dispersion parameters δ M , δ D , and δ K , we mean that given Cholesky factorizations of [M ], [D], and [K] into products of lower-triangular matrices and their transposes,

have

  PDFs ρ [G M ] , ρ [G D ] , and ρ [G K ] that are obtained by maximizing entropy as in (31) with (32)-(34) and thus have a form as in (35) with dispersion parameters δ M , δ D , and δ K , respectively. Thus, δ M , δ D , and δ K ascertain the dispersion of [G M ], [G D ], and [G K ] and therefore also of [M ], [D], and [K]:

  ). Please refer to the online version for color figures. . . . . . . . . . . . 43 9 Random reduced mass and stiffness matrices: data (green solid line) and 98%-confidence region for the corresponding FRF predicted by the nonparametric probabilistic model for δK = 0.25 and δM = 0.10 (light blue filled region) and the updated probabilistic reduced-order model with δK = 0.25 and δM = 0.10 and λ = 50 (blue filled region). Please refer to the online version for color figures. . . . . . . . . . . . . . . . . . . . . . 44

Figure 1 :

 1 Figure 1: Reduced-order model: 50 lowest eigenfrequencies of the deterministic model.

Figure 2 :

 2 Figure2: Using a maximum-entropy PDF as the prior: data (green solid line), corresponding FRFs predicted by the deterministic model (red solid line), and 98%-confidence regions for the corresponding FRFs predicted by the nonparametric probabilistic model for δ K = 0.05 (light blue filled region). Please refer to the online version for color figures.

  ] Receptance [m/N/Hz] (m) Sensor at 8.253 m. ] Receptance [m/N/Hz] (n) Sensor at 8.824 m. ] Receptance [m/N/Hz] (o) Sensor at 9.412 m.

Figure 3 :

 3 Figure 3: Using a maximum-entropy PDF as the prior: data (green solid line), corresponding FRFs predicted by the deterministic model (red solid line), and 98%-confidence regions for the corresponding FRFs predicted by the nonparametric probabilistic model for δK = 0.25 (light blue filled region). Please refer to the online version for color figures.

  Diagonal entries of [ C].

Figure 4 :

 4 Figure 4: Numerical implementation of the Bayesian inversion: (a) reference value ξ 0 , (b) spectrum of the sensitivity matrix given by d j=1 m k=1 µ jk j j (ω k ; ξ 0 ) ⊗ j j (ω k ; ξ 0 ), and (c) diagonal entries of [ C] = λ d j=1 m k=1 µ jk j j (ω k ; ξ 0 ) ⊗ j j (ω k ; ξ 0 ) + [I] -1 for δK = 0.25

Acknowledgements

Part of this work was carried out while the first author was visiting the second author at the Laboratoire Modélisation et Simulation Multi Echelle of the Université Paris-Est Marne-la-Vallée. Financial support of the Fund for Scientific Research (FNRS) is gratefully acknowledged.

List of Figures