A sequential sampling strategy for extreme event statistics in nonlinear dynamical systems

1 Introduction

• Rare extreme events lie in the tails of probability distributions, so resolving them can require many experiments or ensemble realizations.The difficulty is especially acute for nonlinear systems with intermittent events and large intrinsic dimensionality.Source: §1 · p. 1
• The setup assumes random parameters with known distributions and a parameter dimensionality that is, or can be reduced to, moderate size O(10).The passage frames this as the setting for stochastic dynamical systems with transient extreme responses.Source: §1 · p. 1
• Existing approaches have limits: large-deviation methods can be costly and do not characterize the full distribution, while generic efficient statistical inference from dynamics is lacking.The cited event-decomposition approach requires explicit knowledge of the dynamics triggering extremes [15].Source: §1 · p. 2
• The proposed sequential method uses Gaussian process regression and a logarithmic pdf-distance criterion to select samples that reduce uncertainty, with emphasis on distribution tails.Its optimization uses only the inferred parameter-to-observation map; the parameter space must be moderate in size, on the order of O(10).Source: §1 · p. 2

2 Problem setup and method overview

• The objective is to estimate the pdf of a scalar observable induced by an unknown map, while minimizing evaluations because the map can represent an expensive simulation or experiment.The map takes uncertain parameters to a scalar quantity of interest, and the parameter distribution is known.Source: §2 · pp. 2–3
• Given existing samples, the next parameter is chosen to reduce the learned density's error against the target pdf, with special emphasis on large-magnitude map values.The target quantity's pdf is induced by the unknown parameter-to-observation map.Source: §2 · p. 3
• Gaussian process regression supplies a surrogate and posterior uncertainty bounds, and a logarithmic distance between bound-associated pdfs selects the next sample without using the true pdf.The optimization uses the inferred map rather than additional simulations.Source: §2 · p. 3
• Sampling should account for both parameter probability and map magnitude because high-probability regions need not have large map values.The proposed criterion targets regions where both probability and |T_hat| are significant.Source: Figure 1 · p. 3, §2 · p. 3

3 Method description

• GPR learns the parameter-to-observation map, and its posterior variance serves as an uncertainty proxy for selecting the next sample.The posterior distribution provides an explicit mean and kernel for the inferred map.Source: §3 · p. 3
• The method estimates the quantity-of-interest pdf from the posterior mean and selects a sample by minimizing distance between pdfs of confidence bounds.The bounds are formed using the posterior standard deviation; the augmented dataset approximates adding a candidate point without evaluating the expensive map.Source: §3 · p. 4
• The L1 distance between log-transformed pdfs emphasizes extreme and rare events, while the criterion avoids direct evaluation of the target pdf.The criterion instead uses confidence-bound pdfs and inexpensive GPR predictions.Source: §3 · p. 4
• For low-dimensional parameter spaces, particle-swarm optimization and explicit integration are practical; high-dimensional integration may require parameter reduction or importance sampling.The authors recommend an L2-based pre-training phase for main probability mass and allow fixed GPR hyperparameters after calibration.Source: §3 · pp. 4–5

3.1 Asymptotic behavior

• The proposed method converges to the true pdf as the number of samples tends to infinity.The asymptotic criterion result assumes sufficiently smooth GPR mean and standard deviation, large n, and a small ratio ∥∇σ∥/∥∇T∥.Source: §3.1 · p. 5
• Under the stated asymptotic conditions, the criterion samples uncertainty along log-density contours, reaching low-probability regions that low-order-moment criteria may undersample.The denominator involving the map pdf and integration over log fθ contours preserve sampling in rare-event regions.Source: §3.1 · p. 5, Figure 2 · p. 6

4 Applications

• The applications comprise an illustrative colored-noise oscillator and an offshore-platform wave-impact problem with millions of degrees of freedom and simulation runs lasting several hours.The latter motivates estimating extreme-event statistics with fewer expensive simulations.Source: §4 · p. 6

4.1 Nonlinear oscillator driven by correlated stochastic noise

• After only 100 samples, the estimated pdf captures the exact pdf’s qualitative features with very good quantitative agreement; the reference pdf was obtained by sampling the true map at 64,000 points on a 40 × 40 × 40 grid.The experiment starts with 6 Latin-Hypercube samples, performs 12 iterations using the dL2 error metric, then switches to the dL1 error of the logarithms of the pdfs.Source: §4.1 · p. 6
• Nearly an order of magnitude lower L1 error in the logarithmic pdf than Latin-Hypercube sampling is observed across the proposed method’s tested settings.The comparison is against the exact pdf and includes variations in initial dataset size and core iterations.Source: §4.1 · p. 7, Figure 5 · p. 8

4.2 Hydrodynamic forces and moments on an offshore platform

• Offshore wave inputs are parameterized by primary wave-group length L and height A, yielding a two-variable input space; platform loads are computed with three-dimensional SPH simulations.The simulations model irregular waves using the JONSWAP spectrum.Source: §4.2 · pp. 7–8, Figure 6 · p. 8
• After 15 iterations (i.e. 14 total sample points), together with 4 initial-configuration samples, the method approximates the pdf to good agreement with the exact pdf computed from a densely sampled grid.The GP regression operates on the logarithm of the observable because the underlying function is always positive.Source: §4.2 · p. 9

5 Conclusions

• The algorithm selects a sequence of design points that improves estimates of a scalar quantity’s probability distribution, with the selection criterion emphasizing tail statistics.The paper describes the method as practical for very high-dimensional systems when parameter spaces are moderate in dimensionality.Source: §5 · p. 9
• The authors prove asymptotic convergence of the algorithm and analyze its asymptotic behavior.Source: §5 · p. 9
• The method is demonstrated on two problems, including a demanding system with millions of degrees of freedom.Source: §5 · p. 9

A Overview of Gaussian process regression

• GPR places a Gaussian-process prior on the parameter-to-observation map and uses sampled input-output pairs to form a posterior mean and covariance.The posterior variance provides an uncertainty estimate for guiding optimization and exploration.Source: A Overview of Gaussian process… · p. 11
• The GPR mean interpolates the exact map at design points, and adding data points produces non-increasing local posterior variance.The estimated variance vanishes at the design points.Source: A Overview of Gaussian process… · p. 11
• For uniform tensor grids, the map convergence rate is at least O(n^-1/2), while the precise rate depends on the criterion; the convergence result assumes a bounded Lipschitz domain satisfying an interior cone condition.The convergence analysis also uses the relation between convergence of the mean estimate and convergence of estimated variance to zero.Source: A Overview of Gaussian process… · p. 11

• A full Cholesky decomposition costs O(n^3), while updating it for one new point with fixed covariance parameters costs O(n^2).The update avoids recomputing the factorization from scratch at each iteration.Source: B Iterative updates to a Gaussian… · p. 12
• Adding a new design point extends the covariance matrix with its covariances against the existing dataset and its own variance.The passage denotes these terms c21 and c22, respectively.Source: B Iterative updates to a Gaussian… · p. 12
• The update obtains the new Cholesky-factor blocks by solving a linear system using the previously known factor.The solve costs O(n^2).Source: B Iterative updates to a Gaussian… · p. 12

• Each iteration selects a test parameter by minimizing the distance criterion Q, evaluates the response there, and appends the new pair to the dataset.Source: C Algorithm pseudocode · p. 12
• The loop stops when the distance between the estimated pdf bounds meets the user-defined error tolerance.The distance metric is user-defined.Source: C Algorithm pseudocode · p. 12
• The procedure produces a series of pdf estimates for the quantity of interest that converges to the true pdf under appropriate conditions.Source: C Algorithm pseudocode · p. 12

• With no resampling of existing points, L2 convergence of the maps implies convergence in probability of qn and convergence in distribution of Fn.The no-resampling condition yields hΘ → 0, which supports map convergence.Source: D Convergence and asymptotic behavior… · p. 13
• For large n, Q measures and minimizes surrogate-map variance while weighting regions with large |Tn| and small fθ more heavily.The asymptotic statement assumes sufficiently smooth GPR functions, small σ, and a small ratio ∥∇σ∥/∥∇T∥.Source: D Convergence and asymptotic behavior… · p. 13
• The pdf denominator in the log-density criterion gives low-probability regions greater sampling weight, helping capture tail behavior.Source: D Convergence and asymptotic behavior… · p. 14

E Parameterization of a random process using the Karhunen-Loève expansion

• The Karhunen–Loève expansion provides an optimal linear parameterization of a random process in the mean-square sense.Source: E Parameterization of a random process… · p. 15
• Depending on the stochastic process structure, a nonlinear parameterization may capture the desired physics with fewer random parameters.The KL expansion is not the only strategy and may be ill-suited for some systems.Source: E Parameterization of a random process… · p. 15
• The expansion represents a process using covariance eigenfunctions and orthogonal random coefficients; for a Gaussian process, the coefficients are independent and the series can be truncated to a desired accuracy.Source: E Parameterization of a random process… · p. 16

F.1 Numerical experiments

• DualSPHysics supplies the numerical simulations through a meshless Lagrangian smoothed particle hydrodynamics method.The code is open source and has been validated on offshore engineering cases involving forces on structures and wave propagation.Source: §F.1 · p. 16
• The tank has a 100 m horizontal section, 62 m water depth, and a 17.5° sloping beach, with a 45 m-wide, 25 m-high gravity platform.The platform also has three columns that narrow from 10 m base diameter to 4 m width at 42.5 m height.Source: §F.1 · pp. 16–17, Figure 8 · p. 16
• Irregular seas are modeled using a JONSWAP spectral density; its spectral-width parameter δ is 0.07 for k ≤ k0 and 0.09 for k > k0, with fp denoting peak frequency.Source: §F.1 · p. 17

F.2 Additional cases

• The two sea cases use α = 0.060, γ = 3.0, Hs = 13.2 m for JONSWAP 1 and α = 0.100, γ = 2.0, Hs = 15.7 m for JONSWAP 2.Both use a peak wave period of 10 s and are described as representative of North Sea conditions.Source: §F.2 · p. 17
• The force and moment response maps relate wave-group lengthscale θ1 = L and amplitude θ2 = A to the observations.The maps are shown for force Fx and moment My.Source: Figure 9 · p. 17, Figure 10 · p. 18
• Figure 11 gives the target force and momentum pdfs for both spectra, computed from the densely sampled parameter-to-observation map.The response map is fixed regardless of the wave-field spectrum.Source: §F.2 · pp. 17–18, Figure 11 · p. 18
• Algorithm progressions are shown for force and momentum under JONSWAP 1 and JONSWAP 2 in Figures 12–15.The accompanying text specifically identifies the JONSWAP 1 momentum results and points to Figures 14 and 15 for JONSWAP 2.Source: §F.2 · p. 18, Figure 12 · p. 19, Figure 13 · p. 19