arXiv CSopen access

Adaptive Meta-Learning Stochastic Gradient Hamiltonian Monte Carlo Simulation for Bayesian Updating of Structural Dynamic Models

Adaptive Meta-Learning Stochastic Gradient Hamiltonian Monte Carlo Simulation for Bayesian Updating of Structural Dynamic Models Xianghao Meng1,2, James L. Beck3, Yong Huang1,2,*,Hui Li1,2 1

Key Lab of Smart Prevention and Mitigation of Civil Engineering Disasters of the Ministry of Industry and Information Technology, Harbin Institute of Technology, Harbin, China 2

Key Lab of Structures Dynamic Behavior and Control of the Ministry of Education, Harbin Institute of Technology, Harbin, China

3

Division of Engineering and Applied Science, California Institute of Technology, CA, USA

Abstract In the last few decades, Markov chain Monte Carlo (MCMC) methods have been widely applied to Bayesian updating of structural dynamic models in the field of structural health monitoring. Recently, several MCMC algorithms have been developed that incorporate neural networks to enhance their performance for specific Bayesian model updating problems. However, a common challenge with these approaches lies in the fact that the embedded neural networks often necessitate retraining when faced with new tasks, a process that is time-consuming and significantly undermines the competitiveness of these methods. This paper introduces a newly developed adaptive meta-learning stochastic gradient Hamiltonian Monte Carlo (AM-SGHMC) algorithm. The idea behind AM-SGHMC is to optimize the sampling strategy by training adaptive neural networks, and due to the adaptive design of the network inputs and outputs, the trained sampler can be directly applied to various Bayesian updating problems of the

same type of structure without further training, thereby achieving meta-learning. Additionally, practical issues for the feasibility of the AM-SGHMC algorithm for structural dynamic model updating are addressed, and two examples involving Bayesian updating of multi-story building models with different model fidelity are used to demonstrate the effectiveness and generalization ability of the proposed method. Keywords: Bayesian inference, structural dynamics, Markov chain Monte Carlo, metalearning, neural networks.

1. INTRODUCTION Bayesian inference [1,2] is a general, rational, and robust tool which is widely used for structural model updating in structural health monitoring (SHM) [3-7]. The model updating problem is treated using Bayes’ theorem to determine the posterior distribution of the parameter vector based on the available data. However, multidimensional integrals are usually encountered when calculating the posterior distribution, making analytical calculation intractable [1,2]. As a powerful computational tool, Markov chain Monte Carlo (MCMC) methods are widely employed for posterior distribution numerical simulation [8-15]. The MCMC method generates samples as the states of a special Markov chain whose stationary distribution is equal to the posterior distribution. Under the assumption of ergodicity, a sufficient number of posterior model parameter samples can be obtained from the Markov chain simulation. Many classical MCMC algorithms have been developed based on the theory of the Metropolis–Hastings (MH) algorithm. Two examples are the transitional Markov chain Monte Carlo (TMCMC) algorithm [12,14-15], which can efficiently deal with the large difference between the prior and posterior distributions in practical problems, and the Hamiltonian Monte Carlo (HMC) method [8,16-17]

which enables fast exploration of the state space by simply introducing auxiliary

"momentum" variables and replacing the simple random-walk proposal candidates used in MH algorithms with Hamiltonian dynamics proposals. However, these generic MCMC methods are not tailored to any specific Bayesian updating problem. Further extension and optimization of generic MCMC methods [18-20] are useful for exploring their potential in addressing specific problems, but their effectiveness is still correlated with and limited by their flexibility. However, the more algorithmic parameters that need to be tuned, the greater the difficulty for practitioners. In recent years, the application of neural networks has made the extension and optimization of generic MCMC methods easier and tractable. A large number of MCMC algorithms, enhanced by neural networks, have been proposed [21-23]. These algorithms employ neural networks to replace and enhance specific components within MCMC algorithms to get better proposals for targeted problems, with the goal of achieving fast convergence and efficient exploration of the sampling processes. With the improvement of the performance for specific problems after training, the generalization ability of the methods usually decreases. As a consequent drawback, the embedded neural networks often need to be retrained when faced with new tasks, which is time-consuming and significantly weakens the competitiveness of the methods. Meta-learning techniques [24,25] efficiently acquire the common knowledge contained in a set of similar tasks through various ingenious designs of learning algorithms, aiming to guarantee the generalization ability of learners, thereby enabling them to tackle a range of similar tasks by learning from only one or a few of them. Based on the idea of meta-learning, once trained, the MCMC algorithms can be employed directly for a bunch of similar Bayesian updating tasks, thereby significantly reducing training time and greatly enhancing the competitiveness of the neural networkbased algorithms. This paper achieves this goal by combining two strategies [26,27]:

theoretical improvement and network training. A neural network enhanced SGHMC (NN-SGHMC) algorithm [23] selects key potential adjustable parts of the stochastic gradient Hamiltonian Monte Carlo (SGHMC) algorithm and replaces them with neural networks for training; its meta-learning potential is shown for some Bayesian deep learning tasks. However, its generalization ability is limited, and it cannot be employed directly for structural dynamic model updating due to the variety and complexity of dynamic models. This paper presents a newly developed adaptive meta-learning stochastic gradient Hamiltonian Monte Carlo (AM-SGHMC) algorithm. In AM-SGHMC, the replaced key adjustable parts are based on NN-SGHMC, but a newly developed adaptive metalearning technique is introduced to guarantee the scale-invariance of the trained sampler with respect to the posterior PDF, which greatly improves its generalization ability. In addition, the loss function and its back-propagation path are carefully designed, making them more suitable for network training in a Markov chain environment. As for model updating problems, a prior boundary treatment approach is proposed, which enables the free exploration in the early stage of network training. The above efforts successfully enable meta-learning for Bayesian updating of structural dynamic models. The remainder of this paper is organized as follows. In Section 2, the main challenges and a blueprint of the AM-SGHMC are briefly introduced. In Section 3, the main contributions of the AM-SGHMC are further detailed. In Section 4, the AMSGHMC algorithm is summarized in a pseudocode. In Section 5, two examples involving Bayesian updating of a multi-story shear-building model and a multi-story braced-frame building model, are used to demonstrate the effectiveness and generalization ability of the proposed method. Finally, we conclude in Section 6 with a summary of this study.

2. AM-SGHMC: ADAPTIVE META-LEARNING STOCHASTIC GRADIENT HAMILTONIAN MONTE CARLO 2.1 Problem definition and framework for SG-MCMC AM-SGHMC is based on a complete framework for stochastic gradient MCMC (SG-MCMC) [28]. Here, “complete” signifies that the SG-MCMC framework is capable of accommodating any continuous Markov process that is designed to generate samples from a target distribution. In the sample space 𝜽 ∈ ℝ𝐷 , SG-MCMC treats the samples as the positions of particles and performs stochastic dynamic simulation of an Itô diffusion process to generate Markov chains. For the Itô diffusion process, the stationary distribution of particles in a potential energy field 𝑈(𝜽) is given by 𝜋(𝜽) ∝ exp(−𝑈(𝜽)). In the context of the Bayesian approach, the stationary distribution should correspond to the posterior probability density function (PDF), i.e., 𝜋(𝜽) = 𝑝(𝜽|𝓓), where the posterior PDF is derived from Bayes’ theorem: 𝑝(𝜽|𝓓) = 𝑝(𝓓|𝜽)𝑝(𝜽)⁄𝑝(𝓓) ∝ 𝑝(𝓓|𝜽)𝑝(𝜽)

(1)

Here, 𝜽 represents the model parameter vector, 𝑝(𝜽) is the prior PDF, and 𝑝(𝓓|𝜽) is the likelihood function, which quantifies the probability of observing the data 𝓓 given different values of the model parameters 𝜽 . The evidence function 𝑝(𝓓) = ∫ 𝑝(𝓓|𝜽)𝑝(𝜽) 𝑑𝜽 is a normalizing constant in the denominator that ensures the posterior PDF integrates to one. Consequently, the potential energy field can be set as: 𝑈(𝜽) = − log(𝑝(𝜽|𝓓)) + c ∗∗ = − log(𝑝(𝓓|𝜽)) − log(𝑝(𝜽)) + c ∗

(2)

where c ∗ can be set as an arbitrary constant to facilitate calculation, thereby c ∗∗ = − log(𝑝(𝓓)) + c ∗ . It is worth noting that, as with other Monte Carlo methods, the

evidence

function

𝑝(𝓓)

in

Equation

(1),

which

contains an

intractable

multidimensional integral, is not involved here. In SG-MCMC, the augmented state space is defined as 𝒛 = (𝜽, 𝒑) , where 𝒑 denotes a set of auxiliary variables, and the Hamiltonian is defined as 𝐻(𝒛) = 𝑈(𝜽) + ℊ(𝒑|𝜽), based on the designed PDF of auxiliary variables 𝜋(𝒑|𝜽) ∝ exp(−ℊ(𝒑|𝜽)), such that the augmented stationary distribution is given by 𝜋(𝒛) ∝ exp(−𝐻(𝒛)). The dynamic simulation within SG-MCMC is founded on an Itô diffusion process, which can be expressed by the continuous-time stochastic differential equation (SDE): 𝑑𝒛 = 𝒇(𝒛)𝑑𝑡 + √2𝑫(𝒛)𝑑𝑾(𝑡)

(3)

where 𝒇(𝒛) , 𝑾(𝑡) and 𝑫(𝒛) are the deterministic drift, Wiener process and diffusion matrix, respectively. The drift 𝒇(𝒛) can be explicitly parameterized as: 𝒇(𝒛) = −[𝑫(𝒛) + 𝑸(𝒛)]𝛁𝒛 𝐻(𝒛) + 𝚪(𝒛) 𝑑

𝚪𝑖 (𝒛) = ∑ 𝑗=1

𝜕 [𝑫 (𝒛) + 𝑸𝑖𝑗 (𝒛)] 𝜕𝒛𝑗 𝑖𝑗

(4)

(5)

where 𝑸(𝒛) and 𝚪(𝒛) represent the curl matrix and a correction term, respectively. T

⃖ For brevity, we define 𝑨(𝒛) ∙ 𝛁𝒛 ≜ (𝛁𝒛 ∙ 𝑨T (𝒛)) for any vector 𝒛 and matrix 𝑨(𝒛), such that Equation (5) can be represented efficiently by a single matrix equation 𝚪(𝒛) = ⃖[𝑫(𝒛) + 𝑸(𝒛)] ∙ 𝛁𝒛 . The completeness of the above framework has been proved by [28]: (1) 𝜋(𝒛) ∝ exp(−𝐻(𝒛)) is a stationary distribution of the SDEs (3)-(5) for any pair of skew-symmetric matrix 𝑸(𝒛) and positive semi-definite matrix 𝑫(𝒛); (2) for any Itô diffusion process that possesses a unique stationary distribution 𝜋(𝒛), under mild conditions there exist matrices 𝑸(𝒛) and 𝑫(𝒛) such that the process

is governed by SDEs (3)-(5). The completeness of this framework for SG-MCMC and the flexibility in choosing 𝑸(𝒛) and 𝑫(𝒛) are very suitable for the embedding of neural networks.

2.2 Neural network enhancement of SGHMC and its challenges SG-MCMC is a framework that encompasses all sampling methods with continuous Markov processes, including SGHMC, the predecessor of our AM-SGHMC, based on specific choices of 𝐻(𝒛), 𝑸(𝒛) and 𝑫(𝒛). During the stochastic dynamic simulation of SGHMC, 𝜽 and the auxiliary vector 𝒑 are regarded as the position and momentum of the particles respectively. Assuming 1

an identity mass matrix 𝑴 = 𝐈𝐷 , the Hamiltonian is defined as 𝐻(𝒛) = 𝑈(𝜽) + 2 𝒑T 𝒑, 1

where 𝑈(𝜽) and ℊ(𝒑) = 2 𝒑T 𝒑 can be regarded as potential and kinetic energy terms, respectively. Consequently, 𝒑 has the same dimension as 𝜽 , and in the stationary distribution, it follows the standard Gaussian distribution 𝒑~𝒩(𝟎𝐷 , 𝐈𝐷 ). In the SG-MCMC framework, the curl matrix 𝑸(𝒛) and diffusion matrix 𝑫(𝒛) for SGHMC are chosen as: 𝟎 𝑮

−𝑮 ] 𝟎

(6)

𝟎 𝑫(𝒛) = [ 𝟎

𝟎 ] 𝑪

(7)

𝑸(𝒛) = [

where 𝑪 is a positive definite matrix called here the damping matrix and G is called here the gyroscopic-coupling matrix. According to the chosen 𝐻(𝒛), 𝑸(𝒛) and 𝑫(𝒛), the continuous-time SDEs governing the dynamics of SGHMC can be derived as: 𝑑𝜽 = 𝑮𝒑𝑑𝑡

(8)

𝑑𝒑 = −𝑪𝒑𝑑𝑡 − 𝑮𝛁𝜽 𝑈(𝜽)𝑑𝑡 + √2𝑪𝑑𝑾(𝑡)

(9)

where the potential energy gradient 𝛁𝜽 𝑈(𝜽) can be calculated using the auto-

differentiation technique [8]. Then, the corresponding discretized update rule with modified forward Euler discretization of SGHMC with step-size 𝜂 is: 𝒑𝑡+1 = (1 − 𝜂𝑪)𝒑𝑡 − 𝜂𝑮𝛁𝜽𝑡 𝑈(𝜽𝑡 ) + 𝝐𝑡 ,

𝝐𝑡 ~𝒩(𝟎𝐷 , 2𝜂𝑪)

𝜽𝑡+1 = 𝜽𝑡 + 𝜂𝑮𝒑𝑡+1

(10) (11)

It can be seen that, in SGHMC, the gyroscopic-coupling matrix G in 𝑸(𝒛) mainly controls the interaction between 𝜽 and 𝒑, while the role of the damping matrix 𝑪 in 𝑫(𝒛) is to control the intensity of damping to counteract the influence of the corresponding intensity of the random impulse 𝝐𝑡 . In standard SGHMC, G is usually chosen as the identity matrix 𝐈𝐷 and 𝑪 is a constant matrix, which can be chosen as in [20], where other choices of the mass matrix than 𝑴 = 𝐈𝐷 are also examined. Rather than treating the matrices G and 𝑪 as constants, Gong et. al. [23] point out that making them functions of 𝒛 can lead to significant improvements in mixing as well as reduction of sample bias. They use two neural networks to learn functions for diagonal versions of the matrices 𝑮(𝒛) and 𝑪(𝒛) using as inputs the potential energy 𝑈(𝜽), the potential energy gradient 𝜕𝜽 𝑈(𝜽), and the momentum 𝒑. This choice of inputs actually ensures the translation-invariance of the sampling process to the posterior PDF by avoiding the direct use of 𝒛 = (𝜽, 𝒑) , which is a key factor contributing to the generalization ability. However, there are still some challenges in achieving meta-learning for Bayesian updating of structural dynamic models. It is found that, even for data from the same structure, both the arbitrariness of constant c ∗ in Equation (2) and the large variation in the scale of the posterior PDF will lead to a large variation in the network inputs, while the scale variation also demands a large variation in the network outputs. This phenomenon indicates a large difference in the input-output relationship between different tasks, which means that the trained neural networks are no longer suitable for

many new tasks. Therefore, a new network embedding design ensuring the scaleinvariance of the sampling process with respect to the posterior PDF is desired. In terms of network training, for Bayesian updating problems of structural dynamic models, the potential energy 𝑈(𝜽) changes more dramatically, necessitating more reliable training techniques. There are two main aspects that can be improved. Firstly, sample-based loss functions that disregard sample order are often selected to update network parameters through direct backpropagation. Without making full use of the property of Markov chains, both the naive loss function and its straightforward backpropagation path introduce some interference into the network training. Secondly, bounded priors, which constrain parameters to a reasonable range, are rarely employed in gradient-based MCMC methods due to the singularity of their gradients at the boundary. However, bounded priors are crucial for ensuring full exploration of the state space during training by avoiding the complex rare event regions. Our contributions to addressing these challenges are detailed in Section 3.

2.3 Blueprint of AM-SGHMC Before going into details, a higher-level blueprint is first given in this section for our proposed algorithm. In AM-SGHMC, two neural networks are proposed to substitute and optimize the matrices 𝑸(𝒛) and 𝑫(𝒛) for specific sampling problems. Considering the feasibility of network training, the parameterization of matrices 𝑸(𝒛) and 𝑫(𝒛) in SG-MCMC are designed as follows: 𝑸(𝒛) = [

𝟎 −𝑮(𝒛) ] 𝑮(𝒛) 𝟎 𝟎 ] 𝑪(𝒛)

(13)

𝑮(𝒛) = diag[𝒇𝜙𝑄 (𝒛)]

(14)

𝑫(𝒛) = [

𝟎 𝟎

(12)

where

𝑪(𝒛) = diag[𝒇𝜙𝐷 (𝒛)]

(15)

̂ (𝜽), 𝑝𝑖 , 𝐶𝑎𝑡𝑒𝑖 )) 𝒇𝜙𝑄 ,𝑖 (𝒛) = 𝜎𝑖 ∙ (𝑐1 + 𝑓𝜙𝑄 (𝑈

(16)

̂ (𝜽), 𝑝𝑖 , 𝜕𝜃 𝑈 ̂ ∗ (𝜽), 𝐶𝑎𝑡𝑒𝑖 ) 𝒇𝜙𝐷,𝑖 (𝒛) = 𝑐2 + 𝑓𝜙𝐷 (𝑈 𝑖

(17)

where each element is:

where 𝑓𝜙𝑄 > 0 and 𝑓𝜙𝐷 > 0 are the outputs of two embedded neural networks, and 𝑐1 and 𝑐2 are two small positive constants to prevent the vanishing of matrices 𝑸(𝒛) and 𝑫(𝒛) . In order to maintain scale-invariance with respect to the posterior PDF, which enhances meta-learning ability, adaptive scale estimators 𝜎𝑖 and adaptively ̂ (𝜽) and 𝜕𝜃 𝑈 ̂ ∗ (𝜽) are adopted, as detailed in Section normalized network inputs 𝑈 𝑖 3.1, based on their raw materials, 𝜃𝑖 , 𝑈(𝜽) and 𝜕𝜃𝑖 𝑈(𝜽) , which can be obtained during the sampling process. The 𝜃𝑖 and 𝑝𝑖 appearing in the network inputs are the 𝑖-th components of 𝜽 and 𝒑, respectively, corresponding to the 𝑖-th model parameter. Finally, the network input 𝐶𝑎𝑡𝑒𝑖 is an encoded category of the 𝑖-th parameter to allow different sampling strategies for different categories of parameters, which can also be viewed as a constant function of the state space 𝒛 = (𝜽, 𝒑). Thus, the inputs of the ̂ (𝜽), 𝑝𝑖 , 𝜕𝜃 𝑈 ̂ ∗ (𝜽) and 𝐶𝑎𝑡𝑒𝑖 , are all functions of state space 𝒛, which networks, 𝑈 𝑖 results in all elements in matrices 𝑸(𝒛) and 𝑫(𝒛) being functions of 𝒛 . More computational details and explanations about this adaptive input/output processing are shown in Section 3.1. Based on SDEs (3)-(5) with modified forward Euler discretization, the sampling process using discretized dynamics with step-size 𝜂 is: 𝒑𝑡+1 = (1 − 𝜂𝑪(𝒛𝑡 ))𝒑𝑡 − 𝜂𝑮(𝒛𝑡 )𝛁𝜽𝑡 𝑈(𝜽𝑡 ) ⃖ 𝑡 ) ∙ 𝛁𝜽 + ⃖𝑪(𝒛𝑡 ) ∙ 𝛁𝒑 ) + 𝝐𝑡 +𝜂 (𝑮(𝒛 𝑡 𝑡

(18)

⃖ ̂ 𝑡 ) ∙ 𝛁𝒑 𝜽𝑡+1 = 𝜽𝑡 + 𝜂𝑮(𝒛̂𝑡 )𝒑𝑡+1 − 𝜂𝑮(𝒛 𝑡+1

(19)

where 𝑡 is the discretized time; 𝒛𝑡 , 𝒛̂𝑡 and 𝝐𝑡 are defined as 𝒛𝑡 = (𝜽𝑡 , 𝒑𝑡 ), 𝒛̂𝑡 = (𝜽𝑡 , 𝒑𝑡+1 ) and 𝝐𝑡 ~𝒩(𝟎, 2𝜂𝑪(𝒛𝑡 )), respectively. Whenever the initial state 𝒛0 = (𝜽0 , 𝒑0 ) is determined, parameter samples from a Markov chain can be generated by performing Equations (18)-(19) iteratively. After simulating 𝐾 parallel Markov chains for 𝑇 steps, a set of samples 𝚯𝐾𝑇,𝜏 = 𝑇

[ ]

𝜏 {{𝜽𝑘𝑠𝜏 }𝑠=1 }𝐾 𝑘=1 with sampling interval 𝜏 can be obtained. The PDF of the samples

𝚯𝐾𝑇,𝜏

is

written

as

𝑞(𝜽|𝒟) ,

which

can

be

estimated

as

𝑞̅ (𝜽|𝒟) =

𝑇

[ ]

1 1

𝜏 𝐾 𝑘 ∑𝐾 𝑘=1 ∑𝑠=1 ϕ (𝜽; 𝜽𝑠𝜏 , 𝑐𝑜𝑝 𝚺𝚯𝐾 ) using the samples 𝚯 𝑇,𝜏 , where ϕ(∙; 𝝁, 𝚺) is a 𝐾 𝑇

[ ]

𝑇,𝜏

𝜏

multidimensional Gaussian PDF with mean 𝝁 and covariance matrix 𝚺, 𝚺𝚯𝐾 is the 𝑇,𝜏

covariance 1 1

matrix

of

samples

𝚯𝐾𝑇,𝜏 ,

and

𝑐𝑜𝑝

is

set

to

minimize

𝑇

[ ]

𝐾 𝜏 ̅(𝜽𝑘𝑠𝜏 |𝒟). Other more accurate kernel density functions, such as that 𝑇 ∑𝑘=1 ∑𝑠=1 log 𝑞

𝐾[ ] 𝜏

detailed in [10], can also be used. The loss function of AM-SGHMC algorithm can be approximated as: 𝐾

𝐿𝑜𝑠𝑠 = −ℒ𝑉𝐼 (𝚯𝐾𝑇,𝜏 ) =

𝑇 [ ] 𝜏

1 1 ∑ ∑[𝑈(𝜽𝑘𝑠𝜏 ) + log 𝑞̅(𝜽𝑘𝑠𝜏 |𝒟)] 𝐾 [𝑇 ] 𝜏 𝑘=1 𝑠=1

(20)

which is a Monte Carlo estimate of the negative Evidence Lower Bound (ELBO), ℒ𝑉𝐼 (𝑞) = −𝔼𝑞 [𝑈(𝜽)] + ℍ[𝑞] , consisting of an energy term −𝔼𝑞 [𝑈(𝜽)] and an information entropy term ℍ[𝑞]. While maximizing the variational lower-bound ELBO ℒ𝑉𝐼 (𝑞) , reducing the term 𝐿𝑜𝑠𝑠 = −ℒ𝑉𝐼 (𝑞) is like minimizing the KL-divergence KL[𝑞 ∥ 𝜋], which encourages the sampler to generate samples that converge faster to the stationary distribution 𝜋(𝜽) in a finite number of steps 𝑇. During network training,

only the gradients of ℍ[𝑞] are needed, so that the Stein gradient estimator detailed in [23] is used to directly estimate each gradient 𝛁𝜽𝑘𝑠𝜏 log 𝑞(𝜽𝑘𝑠𝜏 |𝒟) using the samples 𝚯𝐾𝑇,𝜏 , which is more effective than first approximating 𝑞̅ (𝜽|𝒟) and then applying the auto-differentiation technique. The specific form of the loss function and more computational details are shown in Section 3.2. For a specific sampling problem, by training with the loss function similar to Equation (20) through an efficient back-propagation path detailed in a later section, AM-SGHMC sampler can reliably achieve a good performance with considerable generalization ability for a set of homogeneous sampling problems. Its meta-learning ability is demonstrated in two examples of Bayesian updating of structural dynamic models in Section 5. An additional issue is that, how to rapidly assess whether the trained sampler can achieve satisfactory generalization performance, when faced with a new and significantly different task. Since the inputs to the networks consist solely of local component-wise information, the conditional posterior PDFs for each parameter, which more intuitively reflect the local component-wise information, are crucial in influencing the performance of the trained sampler. Due to the scale-invariance of the sampler with respect to the posterior PDF, one can briefly assess whether the shapes of the conditional posterior PDFs for the new task resemble those of parameters belonging to the same categories in the training task. To clearly observe the tail trends of the conditional distributions, we recommend visualizing the cross-section of the potential energy function, which is the negative logarithm of the conditional posterior PDF plus an unknown constant, as shown in Equation (2). If the shapes resemble each other, it can be initially concluded that the trained sampler can achieve satisfactory generalization performance for this new task.

3. MAIN CONTRIBUTIONS IN AM-SGHMC 3.1 Adaptive input/output processing with scale-invariance for the posterior PDF In the referenced NN-SGHMC algorithm, the neural networks are directly used to optimize the matrices 𝑸(𝒛) and 𝑫(𝒛), that is, each element in the matrices 𝑮(𝒛) and 𝑪(𝒛) is calculated as: 𝒇𝜙𝑄 ,𝑖 (𝒛) = 𝑐1 + 𝑓𝜙𝑄 (𝑈(𝜽), 𝑝𝑖 )

(21)

𝒇𝜙𝐷 ,𝑖 (𝒛) = 𝑐2 + 𝑓𝜙𝐷 (𝑈(𝜽), 𝑝𝑖 , 𝜕𝜃𝑖 𝑈(𝜽))

(22)

By training with a loss function similar to Equation (20), an NN-SGHMC sampler can achieve a satisfactory performance for specific sampling problems. However, given the arbitrariness of the potential energy reference point and the large variation in the scale of the posterior PDF, the generalization ability of the NN-SGHMC algorithm is limited, preventing it from achieving meta-learning for Bayesian updating of structural dynamic models. Therefore, an input/output processing procedure with scaleinvariance with respect to the posterior PDF is developed and detailed in this section. As shown in Equation (2), the constant c ∗ is of arbitrary value and so the potential energy input 𝑈(𝜽) can be very different even in the same model updating problem. To fix this issue, the normalization shown in Equation (23) is applied to the potential energy input 𝑈(𝜽). ̂ (𝜽) = 𝑈

𝑈(𝜽) − 𝜇𝑈 √2𝐷𝜎𝑈

(23)

̂ (𝜽) is the normalized potential energy, 𝐷 is the dimension of the parameter where 𝑈 space, and the lower-order moments 𝜇𝑈 and 𝜎𝑈 are the mean and standard deviation of the potential energy 𝑈(𝜽), respectively, which can be estimated adaptively during

burn-in similar to those adopted in Adam [29]. The numerator part eliminates the effect ̂ (𝜽) scale-invariant with respect to the posterior of the constant c, while making the 𝑈 ̂ (𝜽) as stable PDF. As for the denominator part, it can make the minimum value of 𝑈 as possible when the dimension 𝐷 of the parameter space changes from task to task. In addition to the normalization of the potential energy input, the following input/output processes are performed to achieve scale-invariance with respect to the posterior PDF, based on the adaptively estimated standard deviation 𝜎𝑖 , 𝑖 = 1, ⋯ , 𝐷 of each dimension of the parameter posterior PDF during burn-in. The potential energy gradient input 𝜕𝜃𝑖 𝑈(𝜽) in Equation (22) is changed into: ̂ ∗ (𝜽) = 𝜎𝑖 ∙ 𝜕𝜃 𝑈 ̂ (𝜽) 𝜕𝜃𝑖 𝑈 𝑖

(24)

And the output for elements of matrix 𝑸(𝒛) shown in Equation (21) is changed into: ̂ (𝜽), 𝑝𝑖 )) 𝒇𝜙𝑄 ,𝑖 (𝒛) = 𝜎𝑖 ∙ (𝑐1 + 𝑓𝜙𝑄 (𝑈

(25)

The above input/output processing approach is proposed based on the phenomenon that the shape of the conditional posterior PDF for parameters of the same category is relatively similar. Considering that the difference in PDF shape between different parameter categories may be relatively large, the parameter category is also encoded as one of the neural network inputs to accommodate different sampling strategies for different shapes. Specifically, the categories can be divided based on the functionality of the parameters, and subsequently, one-hot encoding, which is one of the most important encoding techniques for categorical data, is recommended. Of course, other ad-hoc encoding techniques can also be employed. Denoting the encoded category of the 𝑖-th parameter as 𝐶𝑎𝑡𝑒𝑖 , the Equations (25) and (22) can be finally updated into Equations (16)-(17) shown in Section 2.3. The process of adaptive estimates and a simple proof of the scale-invariance are

presented in Appendix A and Appendix B, respectively. With the computational details in Equations (23)-(24) as well as Appendix A, the newly developed adaptive meta-learning procedure, AM-SGHMC, is formed as shown in Figure 1, which aims to improve the generalization ability of the neural networkbased sampler for Bayesian updating problems of structural dynamic models.

Figure 1. Schematic of AM-SGHMC.

3.2 Loss function and its back-propagation path adapted to the Markov chain environment As for the back-propagation path during training, truncated back-propagate through time (BPTT) is applied in NN-SGHMC algorithm. That is, manually stop the gradient flow through the input of 𝑸 and 𝑫 matrices every 𝑇𝑇𝐵𝑃𝑇𝑇 steps to avoid computing higher order gradients, where the truncated step 𝑇𝑇𝐵𝑃𝑇𝑇 is typically 10 to 20. In addition, the sampling interval for training is typically 𝜏 = 1. The schematic of the forward calculation process involved in one backpropagation of NN-SGHMC algorithm is shown in Figure 2. During back-propagation, the gradient flow will accumulate from the loss function against the arrow direction to the neural network parameters 𝝓, guiding their update. Specifically, the loss function

will firstly guide the optimization of the samples, and then the optimization of each sample will guide the optimization of the previous sample and the corresponding matrices 𝑮(𝒛) and 𝑪(𝒛) based on the dynamic simulation process, Equations (18)(19), up to the initial sample. Finally, the optimization of all the matrices 𝑮(𝒛𝑘𝑡0+𝑠 ), 𝑪(𝒛𝑘𝑡0+𝑠 ) and 𝑮(𝒛̂𝑘𝑡0+𝑠 ), 𝑠 = 0, ⋯ , 𝑇𝑇𝐵𝑃𝑇𝑇 , will jointly guide the optimization of the network parameters 𝝓.

Figure 2. Schematic of the forward calculation process involved in one backpropagation of NN-SGHMC algorithm. Since the samples generated by MCMC methods are the states of a Markov chain, the generation of each sample can be fully determined by the previous adjacent sample only, which means that the information of all the other samples will be redundant when the previous adjacent sample is known. Accordingly, in AM-SGHMC algorithms, each sample 𝜽𝑘𝑡0+𝑠∗ is determined only by the matrices 𝑮(𝒛𝑘𝑡0 +𝑠∗−1 ) , 𝑪(𝒛𝑘𝑡0+𝑠∗−1 ) and 𝑮(𝒛̂𝑘𝑡0 +𝑠∗−1 ) output from the networks in the last step when the previous adjacent sample 𝜽𝑘𝑡0+𝑠∗−1 is known, i.e., each matrix 𝑮(𝒛) or 𝑪(𝒛) is responsible for only one subsequent sample. However, according to the back-propagation path in Figure 2, each matrix 𝑮(𝒛) or 𝑪(𝒛) is responsible for all subsequent samples within the

truncated step 𝑇𝑇𝐵𝑃𝑇𝑇 , which leads to large noise and affects the training effect. Thus, we firstly improve the back-propagation path by manually stop the gradient flow through the states at each step, which is similar to setting 𝑇𝑇𝐵𝑃𝑇𝑇 = 1, to ensure the stability of the training process. To distinguish, the number of steps involved in each back-propagation of AM-SGHMC is denoted as 𝑇𝑇 rather than 𝑇𝑇𝐵𝑃𝑇𝑇 . This improvement can be achieved directly for the energy term of loss function on the left of Figure 2, but for the information entropy term on the right, the estimate of the sample distribution 𝑞(𝜽|𝒟) needs to be considered. When estimating 𝑞̅(𝜽𝑘𝑡0+𝑠∗ |𝒟) to guide the sampling of 𝜽𝑘𝑡0 +𝑠∗ , considering that the samples after 𝑡0 + 𝑠 ∗ are affected by 𝜽𝑘𝑡0+𝑠∗ and have considerable uncertainty because they have not yet been generated when sampling 𝜽𝑘𝑡0 +𝑠∗ , we only utilize the previous samples, i.e., {{𝜽𝑘𝑡0+𝑠 }

𝑠∗

}𝐾 , rather than all the samples 𝑠=0 𝑘=1

{{𝜽𝑘𝑡0 +𝑠 }

𝑇𝑇

}𝐾 . 𝑠=0 𝑘=1

This change makes the order of samples reflected in the loss function. We denote the sample distribution estimated with samples {{𝜽𝑘𝑡0+𝑠 }

𝑠∗

}𝐾 𝑠=0 𝑘=1

as 𝑞̅𝑠∗ (𝜽|𝒟). Although

the repeated estimation of 𝑞̅𝑠∗ (𝜽|𝒟) with different 𝑠 ∗ takes more time, it is insignificant compared to the calculation of the energy term. Finally, considering that the estimated distribution is not accurate when the number of samples is small, the first 𝑀 estimates {𝑞̅𝑠 (𝜽𝑘𝑡0+𝑠 |𝒟)}

𝑀 𝑠=1

are not used in the

estimation of the information entropy term. So that the loss function for each backpropagation of AM-SGHMC algorithm can be written more concretely as: 𝐾

𝐿𝑜𝑠𝑠𝑡0 =

𝑇 [ 𝑇] 𝜏

𝐾

𝑇 [ 𝑇] 𝜏

1 1 1 1 ∑ ∑ 𝑈(𝜽𝑘𝑡0 +𝑠𝜏 ) + ∑ ∑ log 𝑞̅𝑠 (𝜽𝑘𝑡0 +𝑠𝜏 |𝒟) (26) 𝐾 [𝑇𝑇 ] 𝐾 [𝑇𝑇 ] − 𝑀 𝑘=1 𝑠=M+1 𝜏 𝑘=1 𝑠=1 𝜏

where 𝜏 = 1 is recommended.

In addition, once 𝑇𝑇 is too large, the framework for SG-MCMC already ensures that the long-term distribution of the samples is consistent with the target posterior PDF. This means that the loss function has already been minimized and is not suitable for further optimizing the sampling strategy networks. Thus, choosing 𝑇𝑇 in the range of 10 to 20, as in NN-SGHMC, remains a good choice. In our experiments below, for example, we set 𝑇𝑇 = 15 and 𝑀 = 3. To sum up, the schematic of the forward calculation process involved in one backpropagation of AM-SGHMC algorithm is shown in Figure 3, where the dashed lines indicate the stop gradient operations.

Figure 3. Schematic of the forward calculation process involved in one backpropagation of AM-SGHMC algorithm.

3.3 Treatment of prior boundaries enables adequate training In addition to the development of the AM-SGHMC algorithm itself, the treatment of PDFs defined by model updating problems is also crucial to ensure the well training of the neural network. A treatment of commonly used bounded priors is proposed in this section, which allows it to be used in MCMC methods with dynamics proposals. Bounded priors, such as uniform distributions and truncated Gaussian distributions,

can limit the parameters to a reasonable range, so as to reduce unnecessary exploration of the state space and avoid the failure of the MCMC method due to the complexity of the low probability regions. For the dynamics proposals utilizing the potential energy gradient, ideally the sample moving to the boundary will be bounced due to the infinite gradient at the boundary. However, due to the discretization of dynamics during calculation, the samples will hardly be exactly at the boundary, so that the infinite gradient at the boundary cannot be sensed. In general, the model parameter vector 𝔀 is equivalent to the state 𝜽 of the Markov chain, but we change this when there is a bounded prior. In order to enable the application of bounded priors in MCMC methods with dynamics proposals, a stochastic variable transformation is proposed, which can extend the parameter space near the boundary to infinity, thus forming a new state space suitable for dynamics proposals. Specifically, when encountering a bounded prior, for 𝑖-th parameter 𝓌𝑖 bounded by interval [𝑏𝑖,1 − 𝛿𝑖,1 , 𝑏𝑖,2 + 𝛿𝑖,2 ] where 𝑏𝑖,1 ≤ 𝑏𝑖,2 , denoting its corresponding state as 𝜃𝑖 ∈ ℝ, we define a monotonic increasing stochastic variable transformation from state space into parameter space as follows: 𝑓(𝜃𝑖 ; 𝑏𝑖,1 , 𝛿𝑖,1 ), 𝜃𝑖 < 𝑏𝑖,1 𝑏𝑖,1 ≤ 𝜃𝑖 ≤ 𝑏𝑖,2 𝓌𝑖 = {𝜃𝑖 , 𝑓(𝜃𝑖 ; 𝑏𝑖,2 , 𝛿𝑖,2 ), 𝜃𝑖 > 𝑏𝑖,2

(27)

where the function 𝑓(𝜃; 𝑏, 𝛿) is defined as: 𝑓(𝜃; 𝑏, 𝛿) = 2𝛿 ⋅ S𝑖𝑔𝑚𝑜𝑖𝑑 (

2(𝜃 − 𝑏) )+𝑏−𝛿 𝛿

(28)

An example of the proposed stochastic variable transformation with (𝑏1 , 𝛿1 ; 𝑏2 , 𝛿2 )𝑖 = (0.5,0.5; 1.5,0.5) is shown in Figure 4. As can be seen in Figure 4, the parts outside of 𝑏𝑖,1 and 𝑏𝑖,2 of the state space

are compressed into the 𝛿𝑖,1 and 𝛿𝑖,2 widths of the parameter space, respectively, while the part between 𝑏𝑖,1 and 𝑏𝑖,2 is not affected.

Figure 4. Example of stochastic variable transformation. After applying this transformation to all 𝐷 parameters, the potential energy 𝑈(𝜽) for the new state 𝜽 can be easily modified from the potential energy 𝑈(𝔀) of the parameter vector 𝔀: 𝐷

𝑈(𝜽) = − log 𝑝𝜣 (𝜽) + c = 𝑈(𝔀) − ∑ 𝑇(𝜃𝑖 ) + c ′

(29)

𝑖=1

where c′ is an arbitrary constant to facilitate calculation, − ∑𝐷 𝑖=1 𝑇(𝜃𝑖 ) is the modification term of the transformation, the potential energy of the parameter vector is given by Equation (2): 𝑈(𝔀) = − log(𝑝(𝓓|𝔀)) − log(𝑝(𝔀)) + c ∗

(30)

where the elements of the parameter vector 𝔀 can be derived from Equation (27). As for the modification term, the function 𝑇(𝜃𝑖 ) is: 𝑔(𝜃𝑖 ; 𝑏𝑖,1 , 𝛿𝑖,1 ), 𝜃𝑖 < 𝑏𝑖,1 𝑏𝑖,1 ≤ 𝜃𝑖 ≤ 𝑏𝑖,2 𝑇(𝜃𝑖 ) = {0, 𝑔(𝜃𝑖 ; 𝑏𝑖,2 , 𝛿𝑖,2 ), 𝜃𝑖 > 𝑏𝑖,2

(31)

where the function 𝑔(𝜃; 𝑏, 𝛿) is: 𝑔(𝜃; 𝑏, 𝛿) = 2 log 𝑆𝑖𝑔𝑚𝑜𝑖𝑑 (

2(𝜃 − 𝑏) 2(𝜃 − 𝑏) )− + log 4 𝛿 𝛿

The derivation is detailed in Appendix C.

(32)

Based on the new potential energy 𝑈(𝜽) calculated by Equations (27)-(32), the sampler using dynamics proposals can sample the whole real number space, and then transform into the bounded parameter space by Equation (27) to obtain the parameter samples exactly following the posterior PDF 𝑝(𝔀|𝓓). It can be seen that the proposed prior boundary treatment approach avoids the infinite gradient at the boundary and enables the application of bounded priors in MCMC methods with dynamics proposals, which is crucial for the free exploration in the early stage of network training in AM-SGHMC. There are several additional advantages of this approach: (1) The piecewise form of the transformation allows the boundaries on both ends to be treated separately. Furthermore, if there is no boundary on one end, the corresponding segment can be removed. (2) The identity transformation form of the middle segment ensures that the PDF shape of the state vector 𝜽 at the high probability region remains the same as that of the parameter vector 𝔀, which fits the needs of the scale-invariant input/output processing approach in section 3.1. (3) The involved nonlinear functions 𝑆𝑖𝑔𝑚𝑜𝑖𝑑(∙) and log𝑆𝑖𝑔𝑚𝑜𝑖𝑑(∙) are commonly used in neural networks, thus ensuring the stability and robustness of the computations. This boundary treatment approach is specifically designed for component-wise prior boundaries, which are the most intuitive and commonly used. As for the application of this approach, if there is no prior boundary for some parameters that may interrupt the early stage of the training process, we can add them based on their parameter definitions or common sense, without affecting regions with reasonable parameter values. Then, for each prior boundary, the start point 𝑏𝑖 of corresponding transformation segment should be set inside and near the boundary, to keep the PDF shape at the high probability region unchanged. And the width 𝛿𝑖 is the distance between 𝑏𝑖 and the boundary.

4. AM-SGHMC ALGORITHM According to the description in the previous section, the AM-SGHMC algorithm can be summarized as follows: 1. Initialize 𝒛0 = (𝜽0 , 𝒑0 ); 2. Repeat the following for 𝑡 = 0, ⋯ , 𝑁 − 1: In iteration 𝑡 , let the most recent sample be 𝒛𝑡 = (𝜽𝑡 , 𝒑𝑡 ) , then do the following to simulate a new sample 𝒛𝑡+1 = (𝜽𝑡+1 , 𝒑𝑡+1 ): a. Get sample information including 𝜃𝑡,𝑖 , 𝑝𝑡,𝑖 , 𝑈(𝜽𝑡 ), 𝜕𝜃𝑡,𝑖 𝑈(𝜽𝑡 ), 𝐶𝑎𝑡𝑒𝑖 where 𝑖 = 1, ⋯ , 𝐷 (𝑈(𝜽) is detailed in Section 3.3); b. If in "burn-in" phase, update the estimation of the lower-order moments 𝜇𝑈 and 𝜎𝑈 of potential energy 𝑈(𝜽) as well as standard deviation 𝜎𝑖 , 𝑖 = 1, ⋯ , 𝐷 of each dimension of the parameter posterior PDF (discussion of this adaptive estimation is presented in Appendix A); ̂ (𝜽𝑡 ) and 𝜕𝜃 𝑈 ̂ ∗ (𝜽𝑡 ) as in c. Normalize the inputs of neural networks 𝑈 𝑡,𝑖 Equations (23)-(24); d. Calculate elements of matrices 𝑸(𝒛𝑡 ), 𝑫(𝒛𝑡 ) and 𝑸(𝒛̂𝑡 ) using the sampling strategy learned by neural networks to perform discretized dynamics simulation of Itô diffusion process defined by 𝑸(𝒛) and 𝑫(𝒛) to obtain a new sample 𝒛𝑡+1 = (𝜽𝑡+1 , 𝒑𝑡+1 ) as in Equations (12)-(19); e. If training, recreate a copy of 𝒛𝑡+1 to stop the gradient flow; f. If training, every 𝑇𝑇 steps, calculate the loss function defined in Equation (26) and update the neural network by back-propagation.

5. ILLUSTRATIVE EXAMPLES 5.1 Multi-story shear-building model In order to verify the effectiveness and generalization ability of AM-SGHMC, we first consider some multi-story buildings excited by an earthquake. Three sets of noisy accelerometer data are simulated: Dataset 1 are the total acceleration data of a duration of 3 s collected from a 5-story building; Dataset 2 are the total acceleration data of a duration of 1 s collected from a 2-story building; Dataset 3 are the total acceleration data of a duration of 10 s collected from a 10-story building. The acceleration data are all collected with a sample interval of 0.01 s from the base, the first floor, and the roof and contaminated by a large amount of noise (i.e., 100% rms noise-to-signal ratio) as shown in Figure 5. In the data generating process, the ground motion is first generated using a zero-mean Gaussian distribution, and the structural parameters are randomly generated near their nominal values. Then, the structural response outputs are calculated and their average root-mean-square (rms) value is determined. Finally, the noise is generated using a zero-mean Gaussian distribution and added to the previously generated data, with the standard deviation being calculated as this rms value multiplied by the rms noise-to-signal ratio. Multi-story (𝑁-story) linear shear building models are utilized for Bayesian inference and the stiffness parameter 𝑘𝑖 and damping coefficient 𝑐𝑖 for each story, 𝑖 = 1, ⋯ , 𝑁, of the models are estimated.

(a) Acceleration Dataset 1 (5-story).

(b) Acceleration Dataset 2 (2-story).

(c) Acceleration Dataset 3 (10-story).

Figure 5. Acceleration datasets in multi-story shear-building example.

The likelihood function 𝑝(𝓓|𝔀) is set the same as [8]. The number of observed degrees of freedom (2, i.e., first floor and roof) and the length of the discrete time history data (300, 100, 1000 for Datasets 1-3) are denoted as 𝑁𝑜 and 𝑁𝑇 , respectively. For the 𝑛 th observed degree of freedom at time 𝑡𝑗 , the output predicted by the proposed structural model and the corresponding measured output are denoted as 𝑦𝑛 (𝑡𝑗 ; 𝔀) and 𝑦̂𝑛 (𝑡𝑗 ), respectively. The prediction and measurement errors 𝜀𝑛 (𝑡𝑗 ) = 𝑦̂𝑛 (𝑡𝑗 ) − 𝑦𝑛 (𝑡𝑗 ; 𝔀) for 𝑛 = 1,2, ⋯ , 𝑁𝑜 and 𝑗 = 1,2, ⋯ , 𝑁𝑇 , are modeled as independent and identically distributed Gaussian variables with mean zero and some unknown variance 𝜎 2 , based on the Principle of Maximum Entropy [30-31]. Altogether, we need to estimate 𝐷 = 2𝑁 + 1 model parameters with 𝜎 included, and thus the likelihood function 𝑝(𝓓|𝔀) for this problem is: 1

𝑁𝑜 𝑁𝑇

1 2 𝑝(𝓓|𝔀) = ∑ ∑[𝑦̂𝑛 (𝑡𝑗 ) − 𝑦𝑛 (𝑡𝑗 ; 𝔀)] ) 𝑁𝑜 𝑁𝑇 × exp (− 2 2𝜎 (2𝜋𝜎 2 ) 2 𝑛=1 𝑗=1

(33)

Based on the nominal values 𝑘0 = 2 × 107 Nm−1 , 𝑐0 = 6 × 104 Nm−1 s (not equal to the exact values), and a roughly estimated error magnitude 𝜎0 = 1.0ms−2 , the uncertain parameters 𝓌𝑖 , 𝑖 = 1, ⋯ ,2𝑁 + 1, are defined in dimensionless form: 𝓌𝑖 = 𝑘𝑖 /𝑘0 for 𝑖 = 1, ⋯ , 𝑁 ; 𝓌𝑁+𝑖 = 𝑐𝑖 /𝑐0 for 𝑖 = 1, ⋯ , 𝑁 ; and 𝓌2𝑁+1 = 𝜎/𝜎0 . We set the prior PDF 𝑝(𝔀) as truncated independent distributions. The parameters 𝓌𝑖 = 𝑘𝑖 /𝑘0 and 𝓌𝑁+𝑖 = 𝑐𝑖 /𝑐0 follow Gaussian distributions with means of 1 and coefficients of variation (c.o.v.) of 30%, truncated by intervals [0.499,1.501] and [−0.502,3.002], respectively, and 𝓌2𝑁+1 = 𝜎/𝜎0 follows a lognormal distribution with median 1 and a logarithmic standard deviation of 𝑠0 = 0.3 (the c.o.v. is about 30%), truncated by the interval [0.098,3.002]. Then the relationship between the parameter vector 𝔀 and the state 𝜽 of the

Markov chain is established according to Equations (27)-(28), with (𝑏1 , 𝛿1 ; 𝑏2 , 𝛿2 )𝑘𝑖 = (0.5, 0.001; 1.5, 0.001)

for all parameters

𝓌𝑖 = 𝑘𝑖 /𝑘0 ,

(𝑏1 , 𝛿1 ; 𝑏2 , 𝛿2 )𝑐𝑖 =

(−0.5, 0.002; 3, 0.002) for all parameters 𝓌𝑁+𝑖 = 𝑐𝑖 /𝑐0 and (𝑏1 , 𝛿1 ; 𝑏2 , 𝛿2 )𝜎 = (0.1, 0.002; 3, 0.002) for parameter 𝓌2𝑁+1 = 𝜎/𝜎0 .

k10

k3

k2

a

k1

Training case

k3

k2

y1

a

y10

k3 k2

y1

k1

Test case 1 non-generalization task

a

k1

Test case 2

k2

y1

a

y1

k1

Test case 3

generalization tasks

Figure 6. Experimental arrangement for multi-story shear-building example. For AM-SGHMC, the architecture of two embedded neural networks 𝑓𝜙𝑄 (∙) and 𝑓𝜙𝐷 (∙) and training setup are detailed in Appendix D. The experimental arrangement for this example is shown in Figure 6. The sampler is first trained by sampling on Dataset 1 (5-story), and then tested on Datasets 1-3 (5-story as non-generalization task, 2-story and 10-story as generalization tasks). During each test of AM-SGHMC, K = 32 parallel chains for T = 9000 steps are simulated, where the first 3000 steps are burn in. The adaptive estimates are updated only during steps 300 to 2800 with (𝛽1 , 𝛽2 )𝜽 = (0.99,0.995) and (𝛽1 , 𝛽2 )𝑈 = (0.99,0.998) for parameter samples 𝜽 and potential energy 𝑈(𝜽) . The sampling trajectories (after burn in) for three cases by AM-SGHMC for some pair-wise

parameters are shown in Figure 7 and compared with those by HMC method. For each test of HMC, similar with [8], the initial point is determined by 4000 steps of an efficient SPSA (simultaneous perturbation stochastic approximation) optimization algorithm, then K = 32 parallel chains for T = 9000 steps are simulated, where the first 500 steps are burn in. The colors in the figures are the two-dimensional marginal PDFs of the posterior PDFs estimated from the samples. The brighter the color, the higher the sample density is. As can be seen from the figure, the posterior sample distributions obtained by the two methods is similar, but the samples obtained by AM-SGHMC are more evenly distributed in the low probability part at the edge. As for the central high probability part obtained by HMC, the color of stiffness samples is less bright because many samples are wasted when exploring the low probability region, and the contour of damping samples is non-smooth because the initial exploration is relatively slow. The exact values of the dimensionless parameters used in the simulation are marked as pink stars in Figure 7. Comparing the sampling trajectories of the two methods, it can be seen that the deviation between the posterior sample distribution and the true value is mainly caused by the data, and has little to do with the difference of the two methods.

(a) Sampling trajectories of Dataset 1 (5-story).

(b) Sampling trajectories of Dataset 2 (2-story).

(c) Sampling trajectories of Dataset 3 (10-story).

Figure 7. Sampling trajectory plots for some pair-wise parameters by HMC (left column) and AM-SGHMC (right column) for three test cases. It is worth clarifying that, due to the small number of measurement points, the posterior PDF is only strongly constrained by the data information in certain directions. This results in the shape of the two-dimensional projection resembling the prior distribution, an effect that is more pronounced in high-dimensional examples. Therefore, here we attempt to visualize the following two aspects: firstly, the significant

difference in the scale of the posterior PDF across different directions; secondly, the nonlinear correlation between parameters within the high-dimensional space. For the first aspect, which is also a main difficulty that classical HMC method attempts to solve [32-35], Principal Component Analysis (PCA) is used to find orthogonal and linearly independent directions of the samples, and to order the scales of the corresponding posterior PDF along these directions. Figure 8 shows the twodimensional projections of the samples obtained by AM-SGHMC algorithm in the directions of the first and last principal components (PCs) of the three cases, and the projections of the samples obtained by the HMC algorithm in the same directions are also shown for comparison. To facilitate comparison, the sample intervals of the AMSGHMC algorithm are accordingly increased so that the sample numbers of the two methods in the three cases are similar.

(a) Projection from test case 1 (5-story).

(b) Projection from test case 2 (2-story).

(c) Projection from test case 3 (10-story). Figure 8. Projection of samples obtained by HMC (left column) and AM-SGHMC (right column) in the first and last PC directions for three test cases. The figures at the top of each case show the samples from all 32 chains. The colors in these figures correspond to the two-dimensional marginal PDF of the posterior PDF estimated from the samples. It can be seen that the overall distribution obtained by the two methods is consistent, and even though the projection will obscure the highdimensional shape and make the vertical scale larger, the horizontal and vertical scales are still an order of magnitude different. Each figure at the bottom of each case shows samples from only one chain. It can be seen that in test cases 1 and 2, AM-SGHMC algorithm is more effective for the exploration in the direction of the first PC; in test case 3, due to the large sampling interval required by HMC, the exploration efficiency of the two methods is similar. To illustrate how projection transformation obscures the geometric characteristic of high-dimensional samples, consider the ideal distribution of their two-dimensional projection in the figure: 𝑝 (𝜃𝑃𝑖 , 𝜃𝑃𝑗 |𝓓) = ∫ 𝑝(𝜽|𝓓) 𝑑𝜽̅̅̅̅̅ 𝑃𝑖,𝑗 = ∫ 𝑝 (𝜃𝑃𝑖 , 𝜃𝑃𝑗 |𝜽̅̅̅̅̅ 𝑃𝑖,𝑗 ; 𝓓) 𝑝 (𝜽̅̅̅̅̅ 𝑃𝑖,𝑗 |𝓓) 𝑑𝜽̅̅̅̅̅ 𝑃𝑖,𝑗 (34) where 𝜃𝑃𝑖 , 𝜃𝑃𝑗 , and 𝜽̅̅̅̅̅ 𝑃𝑖,𝑗 are the components of the parameter in the 𝑖-th, 𝑗-th and all other PC directions, 𝑃𝑖 , 𝑃𝑗 and ̅̅̅̅ 𝑃𝑖,𝑗 , respectively. When there is a high-dimensional

correlation, the randomness of the unknown correlation components causes the conditional PDF 𝑝 (𝜃𝑃𝑖 , 𝜃𝑃𝑗 |𝜽̅̅̅̅̅̅ 𝑃𝑖 𝑃𝑗 ; 𝓓) to change, making the two-dimensional projection PDF become the weighted average of a series of continuously varying conditional PDFs, and finally form a relatively smooth distribution in the figure. Since the PCs extracted by PCA are linearly independent, when the conditional PDF varies enough to obscure its own shape, the two-dimensional projection in the figure appears to be close to the Gaussian distribution, which is more common in higher dimensional problems. To visualize the high-dimensional nonlinear correlations between parameters in the posterior PDFs, the following three-dimensional relationships are shown. Under the conditions of the selected two-dimensional PC values 𝜃𝑃𝑖 and 𝜃𝑃𝑗 , consider the posterior mean of a third dimensional PC 𝜃𝑃𝑘 of the parameters, which can be approximated using Monte Carlo method with 𝑁𝛿 samples 𝜃𝑛,𝑃𝑘 , 𝑛 = 1, ⋯ , 𝑁𝛿 : 1 𝔼 (𝜃𝑃𝑘 |𝜃𝑃𝑖 , 𝜃𝑃𝑗 ) = ∫ 𝜃𝑃𝑘 𝑝 (𝜃𝑃𝑘 |𝜃𝑃𝑖 , 𝜃𝑃𝑗 ) 𝑑𝜃𝑃𝑘 ≈ 𝑁𝛿

𝑁𝛿

∑

𝜃𝑛,𝑃𝑘

(35)

𝑛=1 𝜃𝑛,𝑃𝑘 ~𝑝(𝜃𝑃𝑘 |𝜃𝑃𝑖 ,𝜃𝑃𝑗 )

However, based on the sampling results presented earlier, due to the absence of multiple 𝜃𝑛,𝑃𝑘 samples with exactly the same (𝜃𝑃𝑖 , 𝜃𝑃𝑗 ) values, a further approximation is made by utilizing 𝜃𝑛,𝑃𝑘 samples within the neighborhood of (𝜃𝑃𝑖 , 𝜃𝑃𝑗 ) and applying multiple iterations of averaging to mitigate noise. Specifically, samples within a threshold of 0.3 standard deviations are used, and 3 iterations of averaging are performed. This iterative operation can retain the overall trend of the conditional posterior mean 𝔼 (𝜃𝑃𝑘 |𝜃𝑃𝑖 , 𝜃𝑃𝑗 ) changing with conditions 𝜃𝑃𝑖 and 𝜃𝑃𝑗 .

Two obvious trends from test case 3 are selected and shown in Figure 9. The colors in the figure represent the two-dimensional marginal PDF of the posterior PDF, estimated from the samples, along the 𝑃𝑖 and 𝑃𝑗 directions.

(a) Nonlinear relationship close to the paraboloid.

(b) Nonlinear relationship close to the saddle surface. Figure 9. Top (left column) and oblique (right column) view of two obvious nonlinear relations obtained by filtering in test case 3. As can be seen from the figures, even if the high-dimensional samples are compressed into three-dimensional space by projection, the nonlinear relationship between the parameters can still be seen after being smoothed by iterative averaging. It is the high-dimensional complex manifold behind this nonlinear phenomenon that leads to the main difficulties faced by the sampling methods. Back to the sampling results, to further confirm the consistency of the results of the two methods, we use the naive version loss function in Equation (20). As mentioned

in Section 3.2, this naive version sample-based loss is independent of the sample order, so that it directly indicates the correctness of the local distribution of samples. The Naive Loss of the two methods on Datasets 1-3 are evaluated and recorded in Tables 13, respectively. It can be seen that the mean Naive Loss values of the two methods are always similar, which indicates the consistency of sample distributions obtained by the two methods. In order to compare the performance of the different samplers, Effective Sample Size (ESS) is employed, which is often used in the evaluation of samplers. ESS for a sequence of correlated samples can be viewed as the number of independent samples generated by the same PDF, where the two sets of samples can achieve the same accuracy in estimating the mean of the PDF. Our implementation of ESS estimation follows [23] and is briefly presented here. For a set of correlated samples 𝚯 𝑇 = {𝜽𝜏 }𝑇𝜏=1 generated by a Markov chain, its ESS is estimated by: ESS(𝚯 𝑇 ) =

𝑇 𝑇 [ ]−1 𝑠 𝜌 3 ∑ 1 + 2 𝑠=1 (1 − 𝑇) 𝜌𝑠 0

(36)

𝜌

where 𝜌𝑠 is the autocorrelation of 𝚯 𝑇 at lag 𝑠 estimated by: 0

𝑇

1 𝜌𝑠 = ∑ (𝜽𝜏 − 𝝁 ̂ )T (𝜽𝜏−𝑠 − 𝝁 ̂) 𝑇−𝑠

(37)

𝜏=𝑠+1

where 𝝁 ̂ is estimated by: 𝑇

1 𝝁 ̂ = ∑ 𝜽𝜏 𝑇

(38)

𝜏=1

The sum in Equation (36) is truncated whenever 𝑠 = 1000 or 𝑠 is even number and 𝜌𝑠−1 + 𝜌𝑠 < 0. The ESSs for both methods on each Dataset are evaluated. With each

corresponding time consumption recorded, we can define the sampling efficiency of the samplers as ESS per hour (ESS/h). The sampling efficiency of the two methods on Datasets 1-3 are also tabulated in Tables 1-3, respectively. As can be seen from the tables, the efficiency of AM-SGHMC is 3.2, 2.4 and 4.6 times that of HMC in 5-story, 2-story and 10-story tasks, respectively, and the higher the parameter dimension, the more the improvement is. These results show that AMSGHMC has the generalization ability and is more advantageous in higher-dimensional sampling problems. Table 1. Sampling performance comparison of Dataset 1 (5-story). Methods

Naive Loss

ESS

Time (h)

ESS/h

HMC

296.01

64.53

0.35

184.37

AM-SGHMC

294.93

198.79

0.33

602.39

Table 2. Sampling performance comparison of Dataset 2 (2-story). Methods

Naive Loss

ESS

Time (h)

ESS/h

HMC

83.06

158.67

0.20

793.35

AM-SGHMC

82.23

425.18

0.22

1932.6

Table 3. Sampling performance comparison of Dataset 3 (10-story). Methods

Naive Loss

ESS

Time (h)

ESS/h

HMC

1021.68

37.84

0.88

43.00

AM-SGHMC

1022.92

146.82

0.73

201.12

5.2 Multi-story braced-frame building model In this section, a more detailed building model is employed to further verify the ability of AM-SGHMC. We still consider the 𝑁 -story buildings excited by an earthquake. Different from the previous example, which considers only one degree of freedom (DOF) per floor, in this example, with reference to the benchmark structure [36-

39]

proposed by IASC (International Association for Structural Control) -ASCE

(American Society of Civil Engineers) Task Group on SHM, we consider the planar motion of 3 DOF for each floor, so that a 3𝑁-DOF model is obtained by assuming that the floor slabs of the braced frame “are rigid in and out of the plane”. The 4-story 2-bay by 2-bay steel braced-frame IASC-ASCE benchmark test structure is shown in Figure 10. The structural layout of the members on each floor of the benchmark structure is the same, and is adopted by the 𝑁-story buildings in this example.

Figure 10. The 4-story IASC-ASCE benchmark test structure [36]. Four sets of noisy accelerometer data are simulated: Dataset 1 are training data collected from a 4-story building contaminated by a typical amount of noise (i.e., 10% rms noise-to-signal ratio); Datasets 2-4 are test data collected from a 2-story building, a 4-story building, and a 6-story building, respectively, contaminated by a large amount of noise (i.e., 50% rms noise-to-signal ratio). The total acceleration data are all collected from the base, the first floor, and the roof along the horizontal direction of each outer wall, with a sample interval of 0.01s and a duration of 3 s, as shown in Figure 11.

(a) Acceleration Dataset 1 (4-story with typical noise).

(b) Acceleration Dataset 2 (2-story with large noise).

(c) Acceleration Dataset 3 (4-story with large noise).

(d) Acceleration Dataset 4 (6-story with large noise).

Figure 11. Acceleration datasets in multi-story braced-frame building example 𝑎3 (3)

𝐸𝑖 (4)

( )

𝐸𝑖

𝐸𝑖

𝑎4

𝑎2 (2) 𝐸𝑖 (1)

𝐸𝑖

𝑎1

Figure 12. Each floor is divided into 5 regions by dotted lines The 3𝑁-DOF braced-frame building models are utilized for Bayesian inference (𝑗)

and the Young's modulus parameters 𝐸𝑖

of the 5 regions, 𝑗 = 1, ⋯ ,5, as shown in

Figure 12, for each story, 𝑖 = 1, ⋯ , 𝑁, of the models are estimated. Specifically, each component takes the Young's modulus of the corresponding region, while the corner columns at the junction of two regions will use the average of the Young's modulus of the two regions. The likelihood function 𝑝(𝓓|𝔀) is set the same as the previous example, as shown in Equation (33), but with 𝑁𝑜 = 8 , 𝑁𝑇 = 300 and the 3𝑁-DOF braced-frame building model.

Based on the nominal value 𝐸0 = 2 × 1011 Pa (not equal to the exact value), and a roughly estimated error magnitude 𝜎0 = 0.8ms−2 , the uncertain parameters (𝑗)

𝓌𝑘 , 𝑘 = 1, ⋯ ,5𝑁 + 1, are defined in dimensionless form: 𝓌𝑖× +𝑗− = 𝐸𝑖 /𝐸0 for 𝑖 = 1, ⋯ , 𝑁; 𝑗 = 1, ⋯ ,5 ; and 𝓌 𝑁+1 = 𝜎/𝜎0 . We set the prior PDF 𝑝(𝔀) as truncated independent distributions. The dimensionless modulus of each component determined

by

parameters

(𝑗)

𝓌𝑘 = 𝐸𝑖 /𝐸0 , 𝑘 = 1, ⋯ ,5𝑁

follow

Gaussian

distributions with means of 1 and coefficients of variation (c.o.v.) of 10%, truncated by interval [0.1,2.0] , and 𝓌 𝑁+1 = 𝜎/𝜎0 follows a lognormal distribution with median 1 and a logarithmic standard deviation of 𝑠0 = 0.3 (the c.o.v. is about 30%), also truncated by interval [0.1,2.0]. Then the relationship between the parameter vector 𝔀 and the state 𝜽 of the Markov chain is established according to Equations (27)(28), with (𝑏1 , 𝛿1 ; 𝑏2 , 𝛿2 ) = (0.2,0.1; 1.9,0.1) for all parameters 𝓌𝑖 .

Training case low noise task

Test case 1

Test case 2

Test case 3

generalization tasks with large noise

Figure 13. Experimental arrangement for multi-story braced-frame building example. For AM-SGHMC, the architecture of two embedded neural networks 𝑓𝜙𝑄 (∙) and 𝑓𝜙𝐷 (∙) and training setup are also detailed in Appendix D. The experimental arrangement for this example is shown in Figure 13. The sampler is firstly trained by sampling on Dataset 1 (4-story with typical noise), and then tested on Datasets 2-4 (2-

story, 4-story and 6-story with large noise as generalization tasks). During each test of AM-SGHMC, K = 32 parallel chains for T = 9000 steps are simulated, where the first 3000 steps are burn in. The adaptive estimates are updated only during steps 500 to 2800 with (𝛽1 , 𝛽2 )𝜽 = (0.99,0.998) and (𝛽1 , 𝛽2 )𝑈 = (0.99,0.998) for parameter samples 𝜽 and potential energy 𝑈(𝜽). The samples (after burn in) for three cases by AM-SGHMC for some pair-wise parameters are shown in Figure 14 and compared with those by HMC method. The colors in the figures are the two-dimensional marginal PDFs of the posterior PDFs estimated from the samples. As can be seen from the figure, similar to the previous example, the posterior sample distributions obtained by the two methods is similar, but the samples obtained by AMSGHMC are more evenly distributed, resulting in smoother edges between different colors. As mentioned in the previous section, for such higher-dimensional problems, the shapes of the two-dimensional projections are very similar to those of the prior distributions.

(a) Samples of Dataset 2 (2-story).

(b) Samples of Dataset 3 (4-story).

(c) Samples of Dataset 4 (6-story).

Figure 14. Sample plots for some pair-wise parameters by HMC (left column) and AM-SGHMC (right column) for three test cases. Similar to the previous example, to prove and visualize the high-dimensional nonlinear correlations between parameters in the posterior PDFs, two obvious nonlinear relations obtained by conditional averaging from test case 2 are selected and shown in Figure 15. The colors in the figure represent the two-dimensional marginal PDF of the posterior PDF, estimated from the samples, along the two conditional PC directions.

(a) Nonlinear relationship close to the paraboloid.

(b) Nonlinear relationship close to the saddle surface. Figure 15. Oblique view of two obvious nonlinear relations obtained by filtering in test case 2. As in the previous example, to compare the performance of the different samplers, the Naive Loss, Effective Sample Size (ESS) and sampling efficiency (ESS/h) for both methods on each of the datasets are evaluated and tabulated in Tables 4-6, respectively. Table 4. Sampling performance comparison of Dataset 2 (2-story). Methods

Naive Loss

ESS

Time (h)

ESS/h

HMC

2635.97

7.96

3.20

2.49

AM-SGHMC

2636.83

31.32

1.89

16.59

Table 5. Sampling performance comparison of Dataset 3 (4-story). Methods

Naive Loss

ESS

Time (h)

ESS/h

HMC

2835.20

13.37

5.69

2.35

AM-SGHMC

2835.10

60.33

3.18

18.99

Table 6. Sampling performance comparison of Dataset 4 (6-story). Methods

Naive Loss

ESS

Time (h)

ESS/h

HMC

2898.67

14.24

8.29

1.71

AM-SGHMC

2898.09

69.89

4.59

15.23

As can be seen from the tables, the similar Naive Loss values indicate the similar sample distributions obtained by the two methods. As for the sampling efficiency, similar to the previous example, the sampling efficiency of AM-SGHMC is 6.66, 8.08 and 8.87 times that of HMC in 2-story, 4-story and 6-story tasks, respectively, and the higher the parameter dimension, the more the improvement is. These results show again the generalization ability of AM-SGHMC and its advantages in higher-dimensional sampling problems.

6. CONCLUDING REMARKS In recent decades, MCMC methods have been extensively used for Bayesian updating of structural dynamic models in structural health monitoring. Recently, neural network-enhanced MCMC algorithms were developed to improve performance for specific Bayesian updating tasks. However, a key challenge is the need for retraining the embedded neural networks for new tasks, which is time-consuming and thereby diminishes their competitiveness. In this paper, a novel meta-learning stochastic simulation approach, called the AMSGHMC, is developed to implement meta-learning for Bayesian updating of structural dynamic models where the MCMC sampler is based on trained adaptive neural networks that utilize the characteristics of the various posterior PDFs that are involved. Specifically, an input/output processing approach is proposed that takes the parameter types into consideration and guarantees scale invariance with respect to the posterior PDF of the model parameters. This approach can improve the generalization ability of

the sampler because the shape of the conditional posterior PDF for parameters of the same type will be relatively similar than those of different types. In addition, the loss function and its back-propagation path are improved so that they are more suitable for the training of networks in Markov chain environments. An approach is also proposed to enable the application of bounded priors in MCMC methods with dynamics proposals. Bounded priors, such as uniform distributions and truncated Gaussian distributions, can limit the parameters to a reasonable range, so as to avoid complex low probability regions. However, the infinite gradient at the boundary cannot be sensed by dynamics proposals. The proposed approach avoids this problem by performing carefully designed stochastic variable transformations that extend the parameter space near the boundary to infinity, which is crucial for free exploration in the early stage of network training. Two examples involving Bayesian updating of a multi-story shear-building model and a multi-story braced-frame building model under ground motions are used to demonstrate the effectiveness and generalization ability of AM-SGHMC. Specifically, in the shear-building model example, the adaptive meta-learning sampler of AMSGHMC is trained on a five-story building model updating task and then tested on a five-story (non-generalization), a two-story (generalization) and a ten-story (generalization) tasks. The test results show that, with similar accuracy, the sampling efficiency (Effective Sample Size per hour, ESS/h) of AM-SGHMC is 3.2 times, 2.4 times and 4.6 times, respectively, that of the classical Hamiltonian Monte Carlo method (HMC) for Bayesian updating of structural dynamic models. In the braced-frame building model example, the AM-SGHMC sampler is trained on a four-story building model updating task with typical amounts of noise, and then tested on a two-story, a four-story and a six-story generalization tasks with larger noise levels. The test results

show that, with similar accuracy, the sampling efficiency (ESS/h) of AM-SGHMC is 6.6 times, 8.0 times and 8.8 times, respectively, of that of the HMC method. In these examples, due to savings in training time, the computational effort for AM-SGHMC is less than that of the HMC method. Note that even though both examples involve linear structural dynamics, their relationships among the parameters indicated by the posterior PDFs are nonlinear, which can be seen intuitively from Figs. 9 and 15, indicating their essential system identification problems are actually also nonlinear. Therefore, there is no fundamental theoretical barrier for the application to nonlinear structural system identification problems, as long as there is a suitable model class that includes a likelihood function and a prior PDF of its parameters.

Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements The study was supported by National Natural Science Foundation of China under Grant No. 52078174. References 1 2

3

4

Beck J. L. (2010). Bayesian system identification based on probability logic. Struct Control Health Monit; 17(7): 825-847. Huang Y., Shao C., Wu B., Beck J. L. and Li H. (2019). State-of-the-art review on Bayesian inference in structural system identification and damage assessment. Adv Struct Eng; 22(6): 1329-1351. Zhao M., Huang Y., Zhou W. and Li H. (2021). Bayesian uncertainty quantification for guidedwave-based multidamage localization in plate-like structures using Gibbs sampling. Struct Health Monit; 20(6): 3092-3112. Yuen K. V., Ching J. and Phoon K. K. (2021). Bayesian Learning Methods for Geotechnical Data. ASCE-ASME J Risk Uncertain Eng Syst, Part A-Civ Eng; 7(1): 02020002.

5 6

7 8

9 10 11

12

13 14 15

16 17 18 19

20 21 22 23

Zhu Z., Au S. K., Li B. and Xie Y. L. (2021). Bayesian operational modal analysis with multiple setups and multiple (possibly close) modes. Mech Syst and Signal Proc; 150: 107261. Jia X., Sedehi O., Papadimitriou C., Katafygiotis L. S. and Moaveni B. (2022). Hierarchical Bayesian modeling framework for model updating and robust predictions in structural dynamics using modal features. Mech Syst and Signal Proc; 170: 108784. Li J., Huang Y. and Asadollahi P. (2021). Sparse Bayesian learning with model reduction for probabilistic structural damage detection with limited measurements. Eng Struct; 247: 113183. Cheung S. H. and Beck J. L. (2009). Bayesian model updating using Hybrid Monte Carlo Simulation with application to structural dynamics models with many uncertain parameters. Journal of Engineering Mechanics 135: 243–255. Robert, C. P., and Casella, G. (1999). Monte Carlo statistical methods, Springer, New York. Beck, J. L., and Au, S. K. (2002). “Bayesian updating of structural models and reliability using Markov chain Monte Carlo simulation.” J. Eng. Mech., 128(2), 380–391. Ching, J., Muto, M., and Beck, J. L. (2006). “Structural model updating and health monitoring with incomplete modal data using Gibbs Sampler.” Comput. Aided Civ. Infrastruct. Eng., 21(4), 242–257. Ching, J., and Chen, Y. J. (2007). “Transitional Markov chain Monte Carlo method for Bayesian model updating, model class selection, and model averaging.” J. Eng. Mech., 133(7), 816–832. Muto, M., and Beck, J. L. (2008). “Bayesian updating of hysteretic structural models using stochastic simulation.” J. Vib. Control, 14(1–2), 7–34. Ching, J., and Wang, J. S. (2016). Application of the transitional Markov chain Monte Carlo algorithm to probabilistic site characterization. Engineering Geology; 203:151-167. Wu S, Angelikopoulos P, Papadimitriou C, et al. (2018). Bayesian Annealed Sequential Importance Sampling: An Unbiased Version of Transitional Markov Chain Monte Carlo. ASCE-ASME J Risk Uncertain Eng Syst Part B Mech Eng; 4(1): 011008. Duane, S., Kennedy, A. D., Pendleton, B. J. and Roweth, D. (1987). Hybrid Monte Carlo. Physics Letters B, 195(2):216–222. Neal, R. M. et al. (2011). MCMC using Hamiltonian dynamics in Handbook of Markov Chain Monte Carlo, 2(11). Ahn, S., Korattikara, A. and Welling, M. (2012). Bayesian posterior sampling via stochastic gradient fisher scoring. arXiv:1206.6380. Patterson S. and Teh, Y. W. (2013). Stochastic gradient Riemannian Langevin dynamics on the probability simplex. In Advances in Neural Information Processing Systems, pp. 3102– 3110. Catanach T. A. and Beck J. L. (2017). Bayesian system identification using auxiliary stochastic dynamical systems. International Journal of Nonlinear Mechanics 94: 72–83. Levy, D., Hoffman, M. D. and Sohl-Dickstein, J. (2017). Generalizing Hamiltonian Monte Carlo with Neural Networks. arXiv:1711.09268. Song, J., Zhao, S., and Ermon, S. (2017). A-NICE-MC: Adversarial Training for MCMC. Advances in neural information processing systems, 30. Gong, W., Li, Y., and Hernández-Lobato, J. M. (2018). Meta-Learning for Stochastic Gradient MCMC. arXiv:1806.04522.

24 25 26 27

28 29 30 31 32 33 34 35 36 37

38 39

Schmidhuber, J. (1987). Evolutionary principles in self-referential learning, or on learning how to learn: the meta-meta-... hook. PhD thesis, Technische Universität München. Naik, D. K. and Mammone R. J. (1992). Meta-neural networks that learn by learning. In International Joint Conference on Neural Networks, v.1, pp. 437–442. IEEE. Zhai, G., Spencer, B. F., Yan, J. et al. (2025). Coupled data/physics-driven framework for accurate and efficient structural response simulation. Eng Struct; 327: 119636. Zhai, G., Xu, Y. and Spencer, B. F. (2025). Bidirectional graphics-based digital twin framework for quantifying seismic damage of structures using deep learning networks. Struct Health Monit; 24 (1): 86-110. Ma, Y.A., Chen, T. and Fox, E. (2015). A complete recipe for stochastic gradient MCMC. In Advances in Neural Information Processing Systems, pp. 2917–2925. Kingma, D. P. and Ba, J. (2014). Adam: A Method for Stochastic Optimization. Computer Science. Jaynes, E. T. (2003). Probability Theory: The Logic of Science, Cambridge University Press, 1329–1330. Jaynes, E. T. (1957). Information Theory and Statistical Mechanics, Phys. Rev. 106 (4) 620– 630. Duane, S., Kennedy, A. D., Pendleton, B. J., and Roweth, D. (1987). Hybrid Monte Carlo. Physics letters B, 195(2), 216-222. Neal, R. M. (2011). MCMC Using Hamiltonian Dynamics. In Handbook of Markov Chain Monte Carlo (pp. 113-162). Chapman and Hall/CRC. Girolami, M., and Calderhead, B. (2011). Riemann manifold Langevin and Hamiltonian Monte Carlo methods. J. Royal Statistical Society Series B: Statistical Methodology, 73(2), 123-214. Betancourt, M. (2017). A conceptual introduction to Hamiltonian Monte Carlo. arxiv:1701.02434. Beck, J. L., and Bernal, D. (2001). A benchmark problem for structural health monitoring. Experimental Techniques 25.3, 49-52. Johnson, E. A., Lam, H. F., Katafygiotis, L. S., and Beck, J. L. (2004). Phase I IASC-ASCE structural health monitoring benchmark problem using simulated data. Journal of Engineering Mechanics, 130(1), 3-15. Yuen, K. V., Au, S. K., and Beck, J. L. (2004). Two-stage structural health monitoring approach for phase I benchmark studies. Journal of Engineering Mechanics, 130(1), 16-33. Ching, J., and Beck, J. L. (2003). Two-step Bayesian structure health monitoring approach for IASC-ASCE phase II simulated and experimental benchmark studies. California Institute of Technology, Earthquake Engineering Research Laboratory.

Appendix A Adaptive estimate process for AM-SGHMC The adaptive estimate process mentioned in the main text is an extension of the lower-order moment estimate process adopted in Adam [29], one of the most popular optimizers. Adam implements optimization based on adaptive estimates of the 1st moment and raw 2nd moment of the first-order gradient, where the gradients 𝑔𝑡 , 𝑡 = 1, ⋯ , 𝑇 can only be obtained at the corresponding timestep 𝑡 based on the parameters updated in the previous step. The same situation occurs when we estimate the first two moments of the parameters 𝜽 and potential energy 𝑈(𝜽) during sampling, that is, 𝜽𝑡 and 𝑈(𝜽𝑡 ) are obtained at successive timesteps 𝑡 = 1, ⋯ , 𝑇. Therefore, in this section, the samples, such as 𝑔𝑡 , 𝜽𝑡 , and 𝑈(𝜽𝑡 ), obtained for estimating the first two moments are uniformly denoted as 𝑦𝑡 to focus on the adaptive estimate process. The lower-order moment estimate process adopted in Adam is extracted and given in Algorithm 1. In this algorithm, two exponential moving averages are used to update biased 1st and raw 2nd moments, 𝑚𝑡 and 𝑣𝑡 , with exponential decay rates 𝛽1 and 𝛽2 (recommended setting 𝛽1 = 0.9 , 𝛽2 = 0.999 ), respectively. Initialization bias correction terms are subsequently utilized to correct the bias caused by initializing the running average with zeros. Algorithm 1: First-two-moment estimate in Adam. All operations on vectors are elementwise. Require: 𝛽1 , 𝛽2 ∈ [0,1): Exponential decay rates for the moment estimates Require: 𝑦1 , ⋯ , 𝑦𝑇 : Samples obtained at subsequent timesteps 𝑡 = 1, ⋯ , 𝑇 𝑚0 ← 0 (Initialize 1st moment vector) 𝑣0 ← 0 (Initialize 2nd moment vector) 𝑡 ← 0 (Initialize timestep) for 𝑡 ′ = 1, ⋯ , 𝑇 do 𝑡 ←𝑡+1 𝑚𝑡 ← 𝛽1 𝑚𝑡−1 + (1 − 𝛽1 )𝑦𝑡 (Update biased 1st moment estimate) 𝑣𝑡 ← 𝛽2 𝑣𝑡−1 + (1 − 𝛽2 )𝑦𝑡2 (Update biased raw 2nd moment estimate) 𝑚 ̂ 𝑡 ← 𝑚𝑡 /(1 − 𝛽1𝑡 ) (Compute bias-corrected 1st moment estimate)

𝑣̂𝑡 ← 𝑣𝑡 /(1 − 𝛽2𝑡 ) (Compute bias-corrected raw 2nd moment estimate) end for return 𝑚 ̂ 𝑡 , 𝑣̂𝑡 at timesteps 𝑡 = 1, ⋯ , 𝑇 (Estimating the first two moments)

The application scenario of the estimation processes in our method differs from that in Adam in two main aspects. On the one hand, when sampling by MCMC method, 𝐾 parallel Markov chains are simulated at the same time, which means that at each 𝐾

timestep 𝑡, a batch of samples {𝑦𝑡,𝑘 }𝑘=1 of size 𝐾 will be obtained instead of one sample. On the other hand, because the means of the parameters 𝜽 and potential energy 𝑈(𝜽) are not zero, the centered second moment is a more accurate estimate of the variance than the raw second moment. Considering these differences, the iterative update formulas in Algorithm 1 need to be modified accordingly. Same as Algorithm 1, all operations on vectors in this section are elementwise. We start with the basic statistical formulas for the mean and the variance. When the batch size is 𝐾, the mean 𝜇𝑡 and the variance 𝜎𝑡2 at timestep 𝑡 can be calculated as: 𝑡

𝐾

1 𝜇𝑡 = ∑ ∑ 𝑦𝑖,𝑘 𝑡𝐾

(39)

𝑖=1 𝑘=1

𝑡

𝐾

1 2 𝜎𝑡2 = ∑ ∑(𝑦𝑖,𝑘 − 𝜇𝑡 ) 𝑡𝐾 − 1

(40)

𝑖=1 𝑘=1

According to these general term formulas, the recurrence formulas of these two series of numbers can be obtained: 𝐾

𝑡−1 1 1 𝜇𝑡 = 𝜇𝑡−1 + [ ∑ 𝑦𝑡,𝑘 ] 𝑡 𝑡 𝐾

(41)

𝑘=1

𝐾

𝑡𝐾 − 𝐾 − 1 2 𝐾 1 𝑡𝐾 − 𝐾 2 (𝜇 − 𝜇𝑡−1 )2 (42) 𝜎𝑡2 = 𝜎𝑡−1 + [ ∑(𝑦𝑡,𝑘 − 𝜇𝑡 ) ] + 𝑡𝐾 − 1 𝑡𝐾 − 1 𝐾 𝑡𝐾 − 1 𝑡 𝑘=1

Then by defining 𝛽1 and 𝛽2 as: 𝑡−1 𝑡

(43)

𝑡𝐾 − 𝐾 − 1 𝑡𝐾 − 1

(44)

𝛽1 = 𝛽2 =

Equations (41)-(42) can be finally rewritten as: 𝐾

1 𝜇𝑡 = 𝛽1 𝜇𝑡−1 + (1 − 𝛽1 ) [ ∑ 𝑦𝑡,𝑘 ] 𝐾

(45)

𝑘=1

𝐾

1 2 2 𝜎𝑡2 = 𝛽2 𝜎𝑡−1 + (1 − 𝛽2 ) [ ∑(𝑦𝑡,𝑘 − 𝜇𝑡 ) ] 𝐾 𝑘=1

1 + [𝛽2 + (1 − 𝛽2 ) ] (𝜇𝑡 − 𝜇𝑡−1 )2 𝐾

(46)

Based on the recurrence formulas (45)-(46), we modified Algorithm 1 into our extended version as shown in Algorithm 2. Algorithm 2: First-two-moment estimate with batch input and non-zero mean. All operations on vectors are elementwise. Require: 𝛽1 , 𝛽2 ∈ [0,1): Exponential decay rates for the moment estimates Require: 𝐾: Batch size of samples obtained at each timestep 𝐾

𝐾

Require: {𝑦1,𝑘 }𝑘=1 , ⋯ , {𝑦𝑇,𝑘 }𝑘=1 : Samples obtained at subsequent timesteps 𝑡 = 1, ⋯ , 𝑇 𝑚0 ← 0 (Initialize 1st moment vector) 𝑣0 ← 0 (Initialize 2nd moment vector) 𝑡 ← 0 (Initialize timestep) for 𝑡 ′ = 1, ⋯ , 𝑇 do 𝑡 ←𝑡+1 1

𝑦̅𝑡 ← ∑𝐾 𝑘=1 𝑦𝑡,𝑘 𝐾

𝑚𝑡 ← 𝛽1 𝑚𝑡−1 + (1 − 𝛽1 )𝑦̅𝑡 (Update biased 1st moment estimate) 𝑚 ̂ 𝑡 ← 𝑚𝑡 /(1 − 𝛽1𝑡 ) (Compute bias-corrected 1st moment estimate) if 𝑡 = 1 do 𝑚 ̂0 ← 𝑚 ̂1 end if ′ 𝑣𝑡−1 ← (𝑚 ̂𝑡 − 𝑚 ̂ 𝑡−1 )2 + 𝑣𝑡−1 2 1 ̅̅̅ 𝑦𝑡2 ← [(𝑚 ̂𝑡 − 𝑚 ̂ 𝑡−1 )2 + ∑𝐾 ̂𝑡) ] 𝑘=1(𝑦𝑡,𝑘 − 𝑚 𝐾

′ ̅̅̅2 (Update biased centered 2nd moment estimate) 𝑣𝑡 ← 𝛽2 𝑣𝑡−1 + (1 − 𝛽2 )𝑦 𝑡 𝑣̂𝑡 ← 𝑣𝑡 /(1 − 𝛽2𝑡 ) (Compute bias-corrected centered 2nd moment estimate)

end for return 𝑚 ̂ 𝑡 , 𝑣̂𝑡 at timesteps 𝑡 = 1, ⋯ , 𝑇 (Estimating the first two moments)

In AM-SGHMC, the proposed first-two-moment estimate process shown in Algorithm 2 is used in the adaptive estimation of the lower-order moments 𝜇𝑈 and 𝜎𝑈2 of the potential energy 𝑈(𝜽). As for the adaptive estimation of the variances 𝜎𝑖2 , 𝑖 = 1, ⋯ , 𝐷 of each component of the parameter samples during the middle part of the burn-in, the situation is somewhat different. In the testing process, based on the estimation in the training process, an initial rough estimate of the variances can usually be made based on the qualitative and quantitative changes in the new data, which can be used as a prior initial value vector of the adaptive estimation. Denoting the prior initial value vector as 𝑣0∗ , and set 𝑣0 ← 𝑣0∗ , the bias correction term 𝑣̂𝑡 ← 𝑣𝑡 /(1 − 𝛽2𝑡 ) in Algorithm 2 should be modified as: 𝑣𝑡 − 𝛽2𝑡 𝑣0∗ 𝑣̂𝑡 = 1 − 𝛽2𝑡

(47)

However, same as the initialization bias correction term utilized in Adam, the new bias correction term Equation (47) can completely correct the effect of the initial value, making the initial value meaningless. Therefore, the first-order Taylor expansion of Equation (47) at 𝛽2𝑡 = 0 is considered: 𝑣̂𝑡 = 𝑣𝑡 + 𝛽2𝑡 (𝑣𝑡 − 𝑣0∗ )

(48)

On the one hand, when a sufficient number of samples are obtained, that is, 𝑡 → +∞ and 𝛽2𝑡 → 0+ , Equation (48) approaches Equation (47), which is little affected by the initial value. On the other hand, during the initial timesteps, the estimates of

Equation (48) will be biased towards the prior initial value vector, making the estimates more robust. For the training process, the robustness of the adaptive estimates during the initial timesteps is more important due to the uncertainty of the sampler initialization. Thus, in addition to the variances, considering the initial value vector as 𝑚0∗ = 0, the bias correction term 𝑚 ̂ 𝑡 ← 𝑚𝑡 /(1 − 𝛽1𝑡 ) of the 1st moment vector in Algorithm 2 can also be similarly changed into: 𝑚 ̂ 𝑡 = 𝑚𝑡 ∗ (1 + 𝛽1𝑡 )

(49)

By introducing the prior initial value vector 𝑣0∗ and taking Equations (48)-(49) as the initialization bias correction term, the adaptive estimate process is further updated as Algorithm 3 and it is utilized in the adaptive estimation of the variances 𝜎𝑖2 , 𝑖 = 1, ⋯ , 𝐷 of each component of the parameter samples. Algorithm 3: Second moment estimate with prior initial value vector. All operations on vectors are elementwise. Require: 𝛽1 , 𝛽2 ∈ [0,1): Exponential decay rates for the moment estimates Require: 𝐾: Batch size of samples obtained at each timestep Require: 𝑣0∗ : Prior initial value vector of the 2nd moment vector 𝐾

𝐾

Require: {𝑦1,𝑘 }𝑘=1 , ⋯ , {𝑦𝑇,𝑘 }𝑘=1 : Samples obtained at subsequent timesteps 𝑡 = 1, ⋯ , 𝑇 𝑚0 ← 0 (Initialize 1st moment vector) 𝑣0 ← 𝑣0∗ (Initialize 2nd moment vector) 𝑡 ← 0 (Initialize timestep) for 𝑡 ′ = 1, ⋯ , 𝑇 do 𝑡 ←𝑡+1 1

𝑦̅𝑡 ← ∑𝐾 𝑘=1 𝑦𝑡,𝑘 𝐾

𝑚𝑡 ← 𝛽1 𝑚𝑡−1 + (1 − 𝛽1 )𝑦̅𝑡 (Update biased 1st moment estimate) if training process do 𝑚 ̂ 𝑡 ← 𝑚𝑡 ∗ (1 + 𝛽1𝑡 ) (Compute robust bias-corrected 1st moment estimate) else if testing process do 𝑚 ̂ 𝑡 ← 𝑚𝑡 /(1 − 𝛽1𝑡 ) (Compute bias-corrected 1st moment estimate) end if if 𝑡 = 1 do 𝑚 ̂0 ← 𝑚 ̂1

end if ′ 𝑣𝑡−1 ← (𝑚 ̂𝑡 − 𝑚 ̂ 𝑡−1 )2 + 𝑣𝑡−1 2 1 ̅̅̅ 𝑦𝑡2 ← [(𝑚 ̂𝑡 − 𝑚 ̂ 𝑡−1 )2 + ∑𝐾 ̂𝑡) ] 𝑘=1(𝑦𝑡,𝑘 − 𝑚 𝐾

′ ̅̅̅2 (Update biased centered 2nd moment estimate) 𝑣𝑡 ← 𝛽2 𝑣𝑡−1 + (1 − 𝛽2 )𝑦 𝑡 𝑣̂𝑡 ← 𝑣𝑡 + 𝛽2𝑡 (𝑣𝑡 − 𝑣0∗ ) (Compute robust bias-corrected centered 2nd moment estimate)

end for return 𝑣̂𝑡 at timesteps 𝑡 = 1, ⋯ , 𝑇 (Estimating the first two moments)

To more intuitively show the difference between Equation (47) and Equation (48) during the initial timesteps, the estimate for the first timestep 𝑡 = 1 is calculated below. According to Algorithm 3, when 𝑡 = 1 , we have 𝑣0′ = 𝑣0 = 𝑣0∗ determined by the 2

1 1 𝐾 prior initial value vector 𝑣0∗ and ̅̅̅ 𝑦12 = 𝐾 ∑𝐾 𝑘=1 (𝑦1,𝑘 − 𝐾 ∑𝑘 ′ =1 𝑦1,𝑘 ′ )

calculated

𝐾 ̅̅̅2 = based on the obtained samples {𝑦1,𝑘 }𝑘=1, then we have 𝑣1 = 𝛽2 𝑣0′ + (1 − 𝛽2 )𝑦 1

̅̅̅2. 𝛽2 𝑣0∗ + (1 − 𝛽2 )𝑦 1 The bias-corrected centered second moment estimate computed by Equations (47) and (48) are shown in equations (50) and (51), respectively: 𝑣̂1 = ̅̅̅ 𝑦12

(50)

̅̅̅2 𝑣̂1 = 𝛽22 𝑣0∗ + (1 − 𝛽22 )𝑦 1

(51)

It can be seen that Equation (47) is not affected by the prior initial value vector 𝑣0∗ at all, while Equation (48) gives a large weight 𝛽22 to the prior 𝑣0∗ at the first timestep.

Appendix B Simple proof of scale-invariance of the meta-learning procedure in AM-SGHMC The scale-invariance with respect to the posterior PDF mentioned in the main text means that when each component of the posterior PDF is linearly scaled by 𝜆𝑖 , 𝑖 = 1, ⋯ , 𝐷 and shifted by 𝑏𝑖 , 𝑖 = 1, ⋯ , 𝐷 (the shape of the posterior PDF in each dimension is thus unchanged), the displacement between adjacent samples generated

by the sampler is also correspondingly scaled, so that the sampling effect remains invariant. Here the model parameters before and after scaling are represented as random variables 𝜣+ and 𝜣′ , respectively, with a relation 𝛩𝑖′ = 𝜆𝑖 𝛩𝑖+ + 𝑏𝑖 , 𝑖 = 1, ⋯ , 𝐷. For any 𝜽+ in the sample space of 𝜣+ , the corresponding 𝜽′ is defined by 𝜃𝑖′ = 𝜆𝑖 𝜃𝑖+ + 𝑏𝑖 , 𝑖 = 1, ⋯ , 𝐷 . According to the multivariate change of variables theorem in probability theory, which is detailed in Appendix C, the relationship between their PDFs 𝑝𝜣+ (𝜽+ ) and 𝑝𝜣′ (𝜽′ ) is: 𝑝𝜣′ (𝜽′ ) =

1

𝑝𝜣 ∏𝐷 𝑖=1 𝜆𝑖

+

(𝜽+ )

(52) 1

Considering that the potential energy 𝑈(𝜽) = − log 𝑝𝜣 (𝜽) + c and ∏𝐷 𝜆 is a 𝑖=1 𝑖

constant, the relationship between their potential energy 𝑈 + (𝜽+ ) and 𝑈 ′ (𝜽′ ) is: 𝑈 ′ (𝜽′ ) = − log 𝑝𝜣′ (𝜽′ ) + c = − log 𝑝𝜣+ (𝜽+ ) − log

1

+c ∏𝐷 𝑖=1 𝜆𝑖

= − log 𝑝𝜣+ (𝜽+ ) + c ′ = 𝑈 + (𝜽+ ) + 𝑐 ′′

(53)

i.e., 𝑈 ′ (𝜽′ ) = 𝑈 + (𝜽+ ) + 𝑐′′ . Thus, 𝜇𝑈 ′ = 𝜇𝑈 + + 𝑐′′ , 𝜎𝑈 ′ = 𝜎𝑈 + and subsequently, the normalized potential energy inputs have the relationship: ̂ ′ (𝜽′ ) = 𝑈

=

=

𝑈 ′ (𝜽′ ) − 𝜇𝑈 ′ √2𝐷𝜎𝑈 ′ (𝑈 + (𝜽+ ) + 𝑐 ′′ ) − (𝜇𝑈 + + 𝑐 ′′ ) √2𝐷𝜎𝑈 + 𝑈 + (𝜽+ ) − 𝜇𝑈 ∗ √2𝐷𝜎𝑈 +

̂ + (𝜽+ ) =𝑈

(54)

i.e.,

̂ ′ (𝜽′ ) = 𝑈 ̂ + (𝜽+ ) . 𝑈

Thus,

̂ ′ (𝜽′ ) = 𝜕𝜃′ 𝑈 ̂ + (𝜽+ ) = 1 𝜕𝜃+ 𝑈 ̂ + (𝜽+ ) . 𝜕𝜃𝑖′ 𝑈 𝑖 𝜆 𝑖

𝑖

Considering that the relationship between standard deviations of 𝛩𝑖′ and 𝛩𝑖+ is 𝜎𝑖′ = 𝜆𝑖 𝜎𝑖+ due to the scaling, the normalized potential energy gradient inputs have the relationship: ′

̂ ∗ (𝜽′ ) = 𝜎𝑖′ ∗ 𝜕𝜃′ 𝑈 ̂ ′ (𝜽′ ) 𝜕𝜃𝑖′ 𝑈 𝑖 = 𝜆𝑖 𝜎𝑖+ ∗

1 ̂ + (𝜽+ ) 𝜕 +𝑈 𝜆𝑖 𝜃𝑖

̂ + (𝜽+ ) = 𝜎𝑖+ ∗ 𝜕𝜃+ 𝑈 𝑖

+

̂ ∗ (𝜽+ ) = 𝜕𝜃 + 𝑈

(55)

𝑖

̂ ∗′ (𝜽′ ) = 𝜕𝜃+ 𝑈 ̂ ∗+ (𝜽+ ). i.e., 𝜕𝜃𝑖′ 𝑈 𝑖

Considering that the auxiliary variables 𝒑~𝒩(𝟎, 𝐈) is invariant during scaling, i.e., the corresponding auxiliary variables before and after scaling will keep 𝒑′ = 𝒑+ , it has been shown that the inputs of both networks 𝑓𝜙𝑄 and 𝑓𝜙𝐷 remain unchanged after the scaling when the corresponding 𝒛′ = (𝜽′ , 𝒑′ ) is used. Then the outputs can be derived based on Equations (16)-(17): ̂ ′ (𝜽′ ), 𝑝𝑖′ , 𝐶𝑎𝑡𝑒𝑖 )) 𝒇′𝜙𝑄 ,𝑖 (𝒛′ ) = 𝜎𝑖′ ∗ (𝑐1 + 𝑓𝜙𝑄 (𝑈 ̂ + (𝜽+ ), 𝑝𝑖+ , 𝐶𝑎𝑡𝑒𝑖 )) = 𝜆𝑖 𝜎𝑖+ ∗ (𝑐1 + 𝑓𝜙𝑄 (𝑈 + = 𝜆𝑖 𝒇+ 𝜙𝑄 ,𝑖 (𝒛 )

(56) ′

̂ ′ (𝜽′ ), 𝑝𝑖′ , 𝜕𝜃′ 𝑈 ̂ ∗ (𝜽′ ), 𝐶𝑎𝑡𝑒𝑖 ) 𝒇′𝜙𝐷 ,𝑖 (𝒛′ ) = 𝑐2 + 𝑓𝜙𝐷 (𝑈 𝑖 ̂ + (𝜽+ ), 𝑝𝑖+ , 𝜕𝜃+ 𝑈 ̂ ∗+ (𝜽+ ), 𝐶𝑎𝑡𝑒𝑖 ) = 𝑐2 + 𝑓𝜙𝐷 (𝑈 𝑖

+ = 𝒇+ 𝜙𝐷 ,𝑖 (𝒛 )

(57)

′ + + ′ + + + + i.e., 𝒇′𝜙𝑄 ,𝑖 (𝒛′ ) = 𝜆𝑖 𝒇+ 𝜙𝑄 ,𝑖 (𝒛 ) and 𝒇𝜙𝐷 ,𝑖 (𝒛 ) = 𝒇𝜙𝐷 ,𝑖 (𝒛 ), where 𝒛 = (𝜽 , 𝒑 ).

Now, we focus on one sampling step before scaling in each dimension 𝑖 at the discretized time 𝑡. Denoting a sample generated from standard Gaussian distribution 𝒩(0, 1) as 𝜖0 , the 𝑖-th component of 𝝐𝑡 ~𝒩(𝟎, 2𝜂𝑪(𝒛𝑡 )) can be written as 𝜖𝑡,𝑖 = + + √2𝜂𝒇𝜙𝐷 ,𝑖 (𝒛𝑡 )𝜖0. Starting with state 𝒛+ 𝑡 = (𝜽𝑡 , 𝒑𝑡 ), according to Equations (18)-(19),

this sampling step can be written as: + (𝜽+ ) 𝑡 + 𝜕𝜃𝑡,𝑖

𝜕𝑈 + + + + + 𝑝𝑡+1,𝑖 = (1 − 𝜂𝒇+ 𝜙𝐷 ,𝑖 (𝒛𝑡 )) 𝑝𝑡,𝑖 − 𝜂𝒇𝜙𝑄 ,𝑖 (𝒛𝑡 )

+ + 𝜕𝒇+ 𝜕𝒇+ 𝜙𝑄 ,𝑖 (𝒛𝑡 ) 𝜙𝐷 ,𝑖 (𝒛𝑡 ) + +𝜂 ( + ) + √2𝜂𝒇+ 𝜙𝐷 ,𝑖 (𝒛𝑡 )𝜖0 + + 𝜕𝜃𝑡,𝑖 𝜕𝑝𝑡,𝑖

+ + + 𝜃𝑡+1,𝑖 = 𝜃𝑡,𝑖 + 𝜂𝒇+ ̂+ 𝑡 )𝑝𝑡+1,𝑖 − 𝜂 𝜙𝑄 ,𝑖 (𝒛

𝜕𝒇+ ̂+ 𝑡) 𝜙𝑄 ,𝑖 (𝒛 + 𝜕𝑝𝑡+1,𝑖

(58)

(59)

+ + + + + where 𝒛̂+ 𝑡 = (𝜽𝑡 , 𝒑𝑡+1 ), and the new state 𝒛𝑡+1 = (𝜽𝑡+1 , 𝒑𝑡+1 ) is generated.

As for the sampling step after scaling, starting with the corresponding state 𝒛′𝑡 = + ′ (𝜽′𝑡 , 𝒑′𝑡 ), i.e., 𝜃𝑡,𝑖 = 𝜆𝑖 𝜃𝑡,𝑖 + 𝑏𝑖 and 𝒑′𝑡 = 𝒑+ 𝑡 , it can be correspondingly written as:

′ ′ 𝑝𝑡+1,𝑖 = (1 − 𝜂𝒇′𝜙𝐷,𝑖 (𝒛′𝑡 )) 𝑝𝑡,𝑖 − 𝜂𝒇′𝜙𝑄 ,𝑖 (𝒛′𝑡 )

𝜕𝑈 ′ (𝜽′𝑡 ) ′ 𝜕𝜃𝑡,𝑖

𝜕𝒇′𝜙𝑄 ,𝑖 (𝒛′𝑡 ) 𝜕𝒇′𝜙𝐷 ,𝑖 (𝒛′𝑡 ) +𝜂 ( + ) + √2𝜂𝒇′𝜙𝐷 ,𝑖 (𝒛′𝑡 )𝜖0 ′ ′ 𝜕𝜃𝑡,𝑖 𝜕𝑝𝑡,𝑖 ′ ′ ′ 𝜃𝑡+1,𝑖 = 𝜃𝑡,𝑖 + 𝜂𝒇′𝜙𝑄 ,𝑖 (𝒛̂′𝑡 )𝑝𝑡+1,𝑖 −𝜂

𝜕𝒇′𝜙𝑄 ,𝑖 (𝒛̂′𝑡 ) ′ 𝜕𝑝𝑡+1,𝑖

where 𝒛̂′𝑡 = (𝜽′𝑡 , 𝒑′𝑡+1 ). For Equation (60), it can be further written as: + + ′ + + 𝑝𝑡+1,𝑖 = (1 − 𝜂𝒇+ 𝜙𝐷 ,𝑖 (𝒛𝑡 )) 𝑝𝑡,𝑖 − 𝜂𝜆𝑖 𝒇𝜙𝑄 ,𝑖 (𝒛𝑡 )

′′ 1 𝜕(𝑈 + (𝜽+ 𝑡)+𝑐 ) + 𝜆𝑖 𝜕𝜃𝑡,𝑖

+ + + 1 𝜕 (𝜆𝑖 𝒇𝜙𝑄 ,𝑖 (𝒛𝑡 )) 𝜕𝒇+ 𝜙𝐷 ,𝑖 (𝒛𝑡 ) + ) + √2𝜂𝒇+ +𝜂 ( + 𝜙𝐷 ,𝑖 (𝒛𝑡 )𝜖0 + + 𝜆𝑖 𝜕𝜃𝑡,𝑖 𝜕𝑝𝑡,𝑖

(60)

(61)

+ + + + = (1 − 𝜂𝒇+ 𝜙𝐷 ,𝑖 (𝒛𝑡 )) 𝑝𝑡,𝑖 − 𝜂𝒇𝜙𝑄 ,𝑖 (𝒛𝑡 )

+𝜂 (

+ 𝜕𝒇+ 𝜙𝑄 ,𝑖 (𝒛𝑡 ) + 𝜕𝜃𝑡,𝑖

+

+ 𝜕𝒇+ 𝜙𝐷 ,𝑖 (𝒛𝑡 ) + 𝜕𝑝𝑡,𝑖

𝜕(𝑈 + (𝜽+ 𝑡 )) + 𝜕𝜃𝑡,𝑖

+ ) + √2𝜂𝒇+ 𝜙𝐷 ,𝑖 (𝒛𝑡 )𝜖0

+ = 𝑝𝑡+1,𝑖 −−−−−−−−−−−−−−−−−−−

(62)

+ ′ i.e., 𝑝𝑡+1,𝑖 = 𝑝𝑡+1,𝑖 for 𝑖 = 1, ⋯ , 𝐷 , that is, 𝒑′𝑡+1 = 𝒑+ ̂ ′𝑡 = 𝑡+1 . That means 𝒛 + + (𝜽′𝑡 , 𝒑′𝑡+1 ) is also the corresponding state of 𝒛̂+ 𝑡 = (𝜽𝑡 , 𝒑𝑡+1 ).

Thus, for Equation (61), it can be further written as:

+ ′ ′ 𝜃𝑡+1,𝑖 − 𝜃𝑡,𝑖 = 𝜂𝜆𝑖 𝒇+ ̂+ 𝑡 )𝑝𝑡+1,𝑖 − 𝜂 𝜙𝑄 ,𝑖 (𝒛

𝜕 (𝜆𝑖 𝒇+ ̂+ 𝑡 )) 𝜙𝑄 ,𝑖 (𝒛

+ = 𝜆𝑖 (𝜂𝒇+ ̂+ 𝑡 )𝑝𝑡+1,𝑖 − 𝜂 𝜙𝑄 ,𝑖 (𝒛

+ 𝜕𝑝𝑡+1,𝑖

𝜕𝒇+ ̂+ 𝑡) 𝜙𝑄 ,𝑖 (𝒛 + 𝜕𝑝𝑡+1,𝑖

)

+ + = 𝜆𝑖 (𝜃𝑡+1,𝑖 − 𝜃𝑡,𝑖 )

(63)

+ + + ′ ′ ′ i.e., 𝜃𝑡+1,𝑖 − 𝜃𝑡,𝑖 = 𝜆𝑖 (𝜃𝑡+1,𝑖 − 𝜃𝑡,𝑖 ). By substituting 𝜃𝑡,𝑖 = 𝜆𝑖 𝜃𝑡,𝑖 + 𝑏𝑖 into it, we can + ′ get 𝜃𝑡+1,𝑖 = 𝜆𝑖 𝜃𝑡+1,𝑖 + 𝑏𝑖 . That means 𝒛′𝑡+1 = (𝜽′𝑡+1 , 𝒑′𝑡+1 ) is also the corresponding + + state of 𝒛+ 𝑡+1 = (𝜽𝑡+1 , 𝒑𝑡+1 ).

According to the derivations above, when the posterior PDF is linearly scaled, + + starting from the initial state 𝒛′0 = (𝜽′0 , 𝒑′0 ) corresponding to 𝒛+ 0 = (𝜽0 , 𝒑0 ) and

using a fixed random seed, parameter samples 𝜽′𝑡 corresponding to 𝜽+ 𝑡 can be generated by performing Equations (60)-(61) iteratively from a Markov chain. Obviously, for any 𝜽′0 in the sample space of 𝜣′ , there is always a 𝜽+ 0 in the 𝜃′ −𝑏𝑖

+ sample space of 𝜣+ calculated by 𝜃0,𝑖 = 0,𝑖𝜆

𝑖

, 𝑖 = 1, ⋯ , 𝐷, so that 𝜽′0 corresponds

to 𝜽+ 0 . Since the initial state of AM-SGHMC can be arbitrary, it is proved that the input/output processing proposed in Section 3.1 has the scale-invariance with respect to the posterior PDF.

Considering that the shape of the conditional posterior PDF for parameters of the same category is relatively similar, it is feasible for parameters of same category 𝐶𝑎𝑡𝑒𝑖 to share one pair of meta-strategy networks 𝑓𝜙𝑄 (∙,∙, 𝐶𝑎𝑡𝑒𝑖 ) and 𝑓𝜙𝐷 (∙,∙,∙, 𝐶𝑎𝑡𝑒𝑖 ).

Appendix C Derivation of the potential energy after transformation According to the multivariate change of variables theorem in probability theory, suppose that 𝑾 is a random variable taking values in 𝒮 ⊆ ℝ𝐷 , and that 𝑾 has a continuous distribution with PDF 𝑝𝑾 . Suppose another random variable 𝜣 = 𝑟(𝑾) where 𝑟 is a one-to-one and differentiable function from 𝒮 onto 𝒯 ⊆ ℝ𝐷 . Then the PDF 𝑝𝜣 of 𝜣 is given by: 𝑑𝔀 )| 𝑝𝜣 (𝜽) = 𝑝𝑾 (𝔀) |det ( 𝑑𝜽

(64)

𝑑𝔀

where ( 𝑑𝜽 ) is the first derivative of the inverse function 𝔀 = 𝑟 −1 (𝜽) defined as the 𝐷 × 𝐷 matrix of first partial derivatives: 𝑑𝔀 𝜕𝓌𝑖 ( ) = 𝑑𝜽 𝑖,𝑗 𝜕𝜃𝑗

(65)

𝑑𝔀

and det ( 𝑑𝜽 ) is the determinant of the first derivative matrix called Jacobian of the inverse function 𝑟 −1. Back to the transformation between states 𝜽 and parameters 𝔀 in the main text, where 𝑝𝑾 is known as the posterior PDF 𝑝𝑾 (𝔀) = 𝑝(𝔀|𝓓) , and the inverse function 𝑟 −1 is defined by Equations (27)-(28): 𝑓(𝜃𝑖 ; 𝑏𝑖,1 , 𝛿𝑖,1 ), 𝜃𝑖 < 𝑏𝑖,1 𝑏𝑖,1 ≤ 𝜃𝑖 ≤ 𝑏𝑖,2 𝑓(𝜃𝑖 ; 𝑏𝑖,2 , 𝛿𝑖,2 ), 𝜃𝑖 > 𝑏𝑖,2

𝓌𝑖 = 𝑟𝑖−1 (𝜃𝑖 ) = {𝜃𝑖 ,

(66)

where the function 𝑓(𝜃; 𝑏, 𝛿) of 𝜃 is defined as: 𝑓(𝜃; 𝑏, 𝛿) = 2𝛿 ⋅ S𝑖𝑔𝑚𝑜𝑖𝑑 (

2(𝜃 − 𝑏) )+𝑏−𝛿 𝛿

(67)

Since the inverse function 𝑟 −1 is elementwise and monotonically increasing in 𝑑𝔀

each component, the first derivative matrix ( 𝑑𝜽 ) becomes a positive definite diagonal matrix and the absolute value of Jacobian in Equation (64) can be simplified to: 𝐷

𝑑𝔀 𝑑𝓌𝑖 )| = ∏ |det ( 𝑑𝜽 𝑑𝜃𝑖

(68)

𝑖=1

and thus Equation (64) becomes: 𝐷

𝑝𝜣 (𝜽) = 𝑝𝑾 (𝔀) ∏ 𝑖=1

𝑑𝓌𝑖 𝑑𝜃𝑖

(69)

Considering the relationship between potential energy 𝑈(𝜽) and target PDF 𝑝𝜣 (𝜽) : 𝑈(𝜽) = − log 𝑝𝜣 (𝜽) + c , we take the negative logarithm of both sides of Equation (69) and get: 𝐷

𝑑𝓌𝑖 ) + c′ 𝑈(𝜽) = 𝑈(𝔀) − ∑ log ( 𝑑𝜃𝑖

(70)

𝑖=1

𝑑𝓌

where 𝑈(𝔀) can be calculated following Equation (30). Denoting log ( 𝑑𝜃 𝑖 ) as 𝑖

𝑇(𝜃𝑖 ), Equation (29) can be obtained. Then we focus on the calculation of 𝑇(𝜃𝑖 ). Recalling Equation (66), we get: 𝑑𝓌𝑖 𝑑𝑟𝑖−1 (𝜃𝑖 ) ) = log ( 𝑇(𝜃𝑖 ) = log ( ) 𝑑𝜃𝑖 𝑑𝜃𝑖 log[𝑓 ′ (𝜃𝑖 ; 𝑏𝑖,1 , 𝛿𝑖,1 )], 𝜃𝑖 < 𝑏𝑖,1 𝑏𝑖,1 ≤ 𝜃𝑖 ≤ 𝑏𝑖,2 = {0, ′ log[𝑓 (𝜃𝑖 ; 𝑏𝑖,2 , 𝛿𝑖,2 )], 𝜃𝑖 > 𝑏𝑖,2

(71)

Denoting log[𝑓 ′ (𝜃; 𝑏, 𝛿)] as 𝑔(𝜃; 𝑏, 𝛿), Equation (31) can be obtained. Finally, after calculating the first derivative 𝑓 ′ (𝜃; 𝑏, 𝛿) =

𝑑𝑓(𝜃; 𝑏,𝛿)

𝑓(𝜃; 𝑏, 𝛿) defined as Equation (67), Equation (32) can be obtained:

𝑑𝜃

of function

𝑔(𝜃; 𝑏, 𝛿) = log [

𝑑𝑓(𝜃; 𝑏, 𝛿) ] 𝑑𝜃

2(𝜃 − 𝑏) 2(𝜃 − 𝑏) = 2 log 𝑆𝑖𝑔𝑚𝑜𝑖𝑑 ( )− + log 4 𝛿 𝛿

(72)

Now the potential energy after transformation is fully derived.

Appendix D Network architecture and training setup Thanks to the scale-invariant input/output processing in Section 3.1, the metastrategy for the sampler for various Bayesian updating problems with the same type of structure can be learned well by simple neural networks embedded in AM-SGHMC. The architecture of the meta-strategy network is shown in Figure 16, which can be divided into three parts according to the dashed boxes. Parts I and III are additional input/output processing to make it suitable for ̂ (𝜽), 𝑝𝑖 , 𝜕𝜃 𝑈 ̂ ∗ (𝜽) and 𝐶𝑎𝑡𝑒𝑖 are converted into network training. In part I, inputs 𝑈 𝑖 network inputs 𝒊𝑈 , 𝒊𝑝𝑖 , 𝒊𝐺𝑖 and 𝒊𝐶𝑖 by functions: 2

̂ (𝜽) + 1)) + e − 1] − 1 𝒊𝑈 = log [(𝑅𝑒𝐿𝑈(𝑈 𝑝𝑖 ) − 1.5 10

(74)

̂ ∗ (𝜽) 𝜕𝜃𝑖 𝑈 ) − 1.5 30

(75)

𝒊𝑝𝑖 = 3 𝑆𝑖𝑔𝑚𝑜𝑖𝑑 ( 𝒊𝐺𝑖 = 3 𝑆𝑖𝑔𝑚𝑜𝑖𝑑 (

(73)

𝒊𝐶𝑖 = 𝑜𝑛𝑒_ℎ𝑜𝑡(𝐶𝑎𝑡𝑒𝑖 )

(76)

where 𝑅𝑒𝐿𝑈(∙) and 𝑆𝑖𝑔𝑚𝑜𝑖𝑑(∙) are two activation functions widely used in neural networks, and 𝑜𝑛𝑒_ℎ𝑜𝑡(∙) is one-hot encoding, one of the most important encoding techniques for categorical data. In the shear-building example, the model parameters are divided into three categories: stiffness for 𝑘𝑖 , damping for 𝑐𝑖 , and noise level for 𝜎. In the braced-frame building example, they are also divided into three categories:

(𝑗)

( )

outer modulus for 𝐸𝑖 , 𝑗 = 1,2,3,4, inner modulus for 𝐸𝑖 , and noise level for 𝜎. In Part III, the network outputs 𝒐𝑄 and 𝒐𝐷 are converted into function outputs by: ̂ (𝜽), 𝑝𝑖 , 𝐶𝑎𝑡𝑒𝑖 ) = 𝑀𝑄 𝑆𝑖𝑔𝑚𝑜𝑖𝑑(5𝒐𝑄 ) 𝑓𝜙𝑄 (𝑈

(77)

̂ (𝜽), 𝑝𝑖 , 𝜕𝜃 𝑈 ̂ ∗ (𝜽), 𝐶𝑎𝑡𝑒𝑖 ) = 𝑀𝐷 𝑆𝑖𝑔𝑚𝑜𝑖𝑑(5𝒐𝐷 ) 𝑓𝜙𝐷 (𝑈 𝑖

(78)

where the maximum values of functions 𝑓𝜙𝑄 and 𝑓𝜙𝐷 are controlled by the parameters 𝑀𝑄 and 𝑀𝐷 . After setting step-size 𝜂 , we recommend setting the parameters 𝑀𝑄 and 𝑀𝐷 to satisfy 𝜂𝑀𝑄 ≈ 3 and 𝜂𝑀𝐷 ≲ 1 , respectively. For example, we set 𝜂 = √0.001, 𝑀𝑄 = 100 and 𝑀𝐷 = 30.

Figure 16. Schematic of meta-strategy network architecture. Part II is two neural networks embedded in functions 𝑓𝜙𝑄 and 𝑓𝜙𝐷 , respectively. In the shear-building example, both of them are 3-hidden-layer multilayer perceptrons (MLPs) with 10 units for each hidden layer, and Leaky ReLU is utilized as the activation function for each unit of the hidden layers. In the braced-frame building example, due to the complexity of the posterior PDF, the network input-output (𝒊𝑄 -𝒐𝑄 and 𝒊𝐷 -𝒐𝐷 ) relationship will change rapidly during the adaptive estimates, so that the slow speed of the network training will slow down, or even disrupt, the process of adaptive estimates. In order to make the networks more flexible in the rapidly changing environment, we

add a shortcut connection beside the 3-hidden-layer MLP (NN), as shown in Figure 17. In the added shortcut connection, the linear transform part (Lin) and the radial basis function part (RBFs) can respond quickly to global and local changes in the 𝒊 -𝒐 relationship, respectively. They will be trained only during the adaptive estimate process and then fixed to remain unchanged. Note that the added short-cut part in Example II is also applicable for Example I.

Figure 17. Schematic of Part II networks for the two examples. For the training process, K 0 = 64 parallel chains are simulated for 100 epochs and each one consists of 10 sub-epochs. After each epoch, states 𝒛 = (𝜽, 𝒑) of the Markov chains are re-initialized using replay techniques with probability 0.2, that is, 20% of the states are re-initialized using the states simulated earlier. For each subepoch, the sampler is simulated 90 steps and then is updated by the Adam optimizer with learning rate of 0.01 and exponential decay rates (𝛽1 , 𝛽2 )Adam = (0.5,0.75). For each 𝑇𝑇 = 15 steps, 𝐾 = 10 chains are randomly selected to estimate the loss function 𝐿𝑜𝑠𝑠𝑡0 ({{𝜽𝑘𝑡0 +𝑠 }

𝑇𝑇

}𝐾 ) 𝑠=0 𝑘=1

as in Figure 3, and the gradient is backpropagated

and accumulated. The adaptive estimates are only updated during the last 6 sub-epochs of the first 50 epochs with (𝛽1 , 𝛽2 )𝜽 = (0.99,0.999) and (𝛽1 , 𝛽2 )𝑈 = (0.99,0.998) for parameter samples 𝜽 and potential energy 𝑈(𝜽).