An Elastic Approach

February 24, 2024

- Department of Energy National Lab
- 14,000 Staff across 7 Locations
- Two Statistics Departments
- 30 Full Time Staff, 4 Post-Docs, 10 year round interns

- Collaborators: Devin Francom (LANL), Gabriel Huerta, Kurtis Shuler, Gavin Collins

- Introduction
- Functional Data Analysis
- Elastic Metric

- Elastic Bayesian Model Calibration
- Results
- Simulation
- Sandia Z-Machine Equation of State

- Question: How can we model functions
- Can we use them as predictors in a regression model?
- Can we calibrate a computer model?

- It is the same goal (question) of any area of statistical research
- One problem occurs when performing these types of analysis is that functional data can contain variability in time (x-direction) and amplitude (y-direction)
- How do we characterize and utilize this variability in the models that are constructed from functional data?

In FDA, one develops the theory on function spaces and not finite vectors, and discretizes the function only at the final step – compute implementation

Ulf Grenander: “Discretize as late as possible” (1924-2016)

Even after discretization, we retain the ability to interpolate resample as needed!

- Cross-sectional statistics ignore horizontal component, often destroys structures

- Let \(f\) be a real valued-function with the domain \([0,1]\), can be extended to any domain
- Only functions that are absolutely continuous on \([0,1]\)

- Let \(\Gamma\) be the group of all warping functions \[ \Gamma = \{\gamma:[0,1]\rightarrow[0,1]|\gamma(0)=0,\gamma(1)=1,~\gamma \mbox{ is a diffeo}\} \]
- It acts on the function space by composition \[(f,\gamma) = f\circ\gamma\]
- It is common to use the following objective function for alignment \[\min_{\gamma\in\Gamma}\|f_1\circ\gamma-f_2\|\]
- Note: It is not a distance function since it is not symmetric.

Define the Square Root Velocity Function \[q:[0,1]\rightarrow\mathbb{R}^1,~q(t)=sign(\dot{f}(t))\sqrt{|\dot{f}(t)|}\] Fisher Rao Distance is \(\mathbb{L}^2\) in SRVF space \[d_a(f_1,f_2) = \inf_\gamma\|(q_1\circ\gamma)\sqrt{\dot{\gamma}}-q_2\|\] Distance is a proper distance

Can compute distance on warping functions (how much alignment) \[d_p(\gamma) = \arccos\left(\int_0^1\sqrt{\dot{\gamma}}\,dt\right)\]

\(\Gamma\) is a nonlinear manifold and it is infinite dimensional

Represent an element of \(\gamma\in\Gamma\) by the square-root of its derivative: \(\psi = \sqrt{\dot{\gamma}}\)

Important advantage of this transformation is the set of all such \(\psi\)’s is a Hilbert Sphere \(\mathbb{S}_\infty\)

- We wish to calibrate a computer model with parameters \(\theta\) to an experiment
- Can compute computer model (simulations) over wide range of \(\theta\)
- The data is functional in nature and has phase and amplitude variability
- Utilize elastic metrics in a Bayesian Model Calibration Framework

Let \(z(t,\boldsymbol{x_i})\) denote an experimental measurement from the \(i^{th}\) experiment

Similarly, let \(y(t,\boldsymbol{x_i}, \boldsymbol{u})\) denote a simulation of the \(i^{th}\) experiment with input parameters \(\boldsymbol{u}\)

An approach to Bayesian model calibration with functional responses is to specify \[z(t,\boldsymbol{x_i}) = y(t,\boldsymbol{x_i},\boldsymbol{\theta}) + \delta(t, \boldsymbol{x_i}) + \epsilon_i(t,\boldsymbol{x_i}), ~~\epsilon(t,\boldsymbol{x_i})\sim \mathcal{N}(0,\sigma_\epsilon^2)\]

Where \(\delta\) is a model discrepancy term and \(\epsilon\) represents all other error

This model will suffer from the aforementioned problems with phase variability

- Decompose observation into aligned functions and warping functions \[z(t,\boldsymbol{x_i}) = \tilde z(t, \boldsymbol{x_i}) \circ\gamma_{z\rightarrow z}(t,\boldsymbol{x_i})\]
- and decompose the simulations \[y(t,\boldsymbol{x_j}, \boldsymbol{u_j}) = \tilde y(t,\boldsymbol{x_j}, \boldsymbol{u_j}) \circ \gamma_{y\rightarrow z}(t, \boldsymbol{x_j}, \boldsymbol{u_j})\]
- To facilitate modeling, we transform the warping functions into shooting vector space with \[\boldsymbol{v}_{z\rightarrow z}(\boldsymbol{x_i}) = \exp_{\psi}^{-1}\left(\sqrt{\dot\gamma_{z\rightarrow z}(\boldsymbol{x_i})}\right)\] \[\boldsymbol{v}_{y\rightarrow z}(\boldsymbol{x_j},\boldsymbol{u_j}) = \exp_{\psi}^{-1}\left(\sqrt{\dot\gamma_{y\rightarrow z}(\boldsymbol{x_j},\boldsymbol{u_j})}\right)\]

- Calibrate the aligned data and shooting vectors using the following model \[\tilde{z}(t,\boldsymbol{x_i}) = \tilde{y}(t, \boldsymbol{x_i}, \boldsymbol{\theta}) + \delta_{\tilde{y}}(t, \boldsymbol{x_i}) + \epsilon_{\tilde{z}}(t, \boldsymbol{x_i}), ~~\epsilon_{\tilde{z}}(t, \boldsymbol{x_i})\sim \mathcal{N}(0,\sigma_{\tilde{z}}^2)\] \[\boldsymbol{v}_{z\rightarrow z}(\boldsymbol{x_i}) = \boldsymbol{v}_{y\rightarrow z}(\boldsymbol{x_i}, \boldsymbol{\theta}) + \boldsymbol{\delta_v(\boldsymbol{x_i})} + \boldsymbol{\epsilon_v(\boldsymbol{x_i})}, ~~\boldsymbol{\epsilon_v(\boldsymbol{x_i})} \sim \mathcal{N}(0,\sigma_v^2\boldsymbol{I})\]
- Note: The shooting vector will be identity if the data is aligned to the observation (experiment)
- Then if \(\theta\) is calibrated correctly the shooting vectors will be identity

For each experiment the likelihood is a Gaussian likelihood

- We fit an emulator (Gaussian Process, BASS, MARS) to the simulated data
- Uniform priors on \(\theta\)
- Sample posterior using delayed rejection adaptive Metropolis Hastings
- Implemented using Impala (LANL) or Dakota (SNL) calibration framework

- Simulation study where each function is parameterized Gaussian pdf
- A set of 100 functions were simulated with \(\theta_0,\theta_1\) being drawn from a \(U[0,1]\)
- Third nuisance parameter \(\theta_2\) also drawn from \(U[0,1]\)

\[f_i(t) = \frac{\theta_1}{0.05\sqrt{2\pi}}\exp\left(-\frac{1}{2}\left(\frac{t-(\sin(2\pi \theta_0^2)/4-\theta_0/10+0.5)}{0.05}\right)^2\right)\]

- Trained BASS Emulator on Aligned Functions and Shooting Vectors (using elastic fPCA)
- Calibrated using framework with tempering and adaptive MCMC
- Blue shows draws from posterior distribution at 95% credible interval

- Repeat same simulation model with with the experiment having a timing discrepancy
- Discrepancy modeling with basis functions in shooting vector space

- Trained BASS Emulator on Aligned Functions and Shooting Vectors (using elastic fPCA)
- Calibrated using framework with tempering and adaptive MCMC
- Blue shows draws from posterior distribution at 95% credible interval

Calibration equation of state of platinum generated with pulse magnetic fields

Estimate the parameters describing the compressibility (relationship between pressure and density) to understand extreme pressures

Conducted using Sandia Z-machine, a pulse power drive reactor

We have a good fit between experiment and simulation

Better residuals over standard calibration method

Comparison to numerical method (ILA), we have good agreement

Small percentage difference between ILA and Calibrated Model

With the elastic method we have tighter posteriors over previous methods

Additionally, scientists at Z have confirmed the parameters conform more to physical understanding of the material

Functional metrics provide a global measure of the difference of a function in terms of amplitude and phase

Integrated elastic functional metrics into Bayesian Model Calibration framework utilizing aligned data and shooting vector representation

Demonstrated ability on simulated and platinum equation of state calibration problems

Future Work

- Additional testing on real world examples

D. Francom, J. D. Tucker, J. G. Huerta, K. Shuler, and D. Ries, “Elastic Bayesian Model Calibration”, arXiv:2305.08834 [stat.ME], in revision, 2023.

J. Brown, J. P. Davis, J. D. Tucker, J. G. Huerta, and K. Shuler, “Quantifying uncertainty in analysis of shockless dynamic compression experiments on platinum, Part 2: Bayesian model calibration”, Journal of Applied Physics, vol. 134, no. 23, 2023.

Questions?

jdtuck@sandia.gov

http://research.tetonedge.net