# Elastic Bayesian Model Calibration

Devin Francom  Los Alamos National Laboratory  J. Derek Tucker  Sandia
National Laboratories  Gabriel Huerta  Sandia National Laboratories  Kurtis
Shuler  Sandia National Laboratories  Daniel Ries  Sandia National
Laboratories

###### Abstract

Functional data are ubiquitous in scientific modeling. For instance,
quantities of interest are modeled as functions of time, space, energy,
density, etc. Uncertainty quantification methods for computer models with
functional response have resulted in tools for emulation, sensitivity
analysis, and calibration that are widely used. However, many of these tools
do not perform well when the model’s parameters control both the amplitude
variation of the functional output and its alignment (or phase variation).
This paper introduces a framework for Bayesian model calibration when the
model responses are misaligned functional data. The approach generates two
types of data out of the misaligned functional responses: one that isolates
the amplitude variation and one that isolates the phase variation. These two
types of data are created for the computer simulation data (both of which may
be emulated) and the experimental data. The calibration approach uses both
types so that it seeks to match both the amplitude and phase of the
experimental data. The framework is careful to respect constraints that arise
especially when modeling phase variation, but also in a way that it can be
done with readily available calibration software. We demonstrate the
techniques on a simulated data example and on two dynamic material science
problems: a strength model calibration using flyer plate experiments and an
equation of state model calibration using experiments performed on the Sandia
National Laboratories’ Z-machine.

Keywords: amplitude/phase variability, Bayesian model calibration, functional
data analysis, strength material calibration.

##  1 Introduction

In domains of science and engineering where modeling is an important part of
investigation and discovery, quantifying uncertainty in model inferences and
predictions can be essential in order for model performance to be trusted.
When a model has uncertain parameters (or inputs), model calibration is the
act of tuning the parameters so that the model produces a desired response.
Most often, models are calibrated to experimental or observational
measurements, so that calibration seeks to make the model response reflect
reality. Model calibration (also known as inversion) is often a poorly
identified problem, where multiple combinations of inputs produce equally
valid solutions. Kennedy and O’Hagan (2001) proposed Bayesian model
calibration as a systematic approach to calibrating a model in the face of all
of the sources of uncertainty, so that calibration uncertainty is quantified.
These sources of uncertainty are parameter uncertainty, measurement
uncertainty, emulation or surrogate model uncertainty (an emulator is a fast
surrogate for a more expensive model), and model form error (i.e., model
misspecification, model discrepancy or model bias).

Often, computer model outputs are functional in nature, producing an output
measured over space and/or time. The majority of calibration work on these
types of outputs has been done on features extracted from the outputs such as
peak points or critical values (Walters et al., 2018). However, the process
for extracting these features can be tedious, prone to error, and is problem
specific. Other calibration approaches have been developed for use with
multivariate or functional response (Higdon et al., 2008; Bayarri et al.,
2007; Francom et al., 2019), though they are prone to problems when presented
with functional response data that are misaligned. Our goal is to extend
functional emulation and calibration methods for use when the response is
misaligned, and to do so without human-intensive feature engineering.

There has been considerable effort in statistics to develop methods that can
analyze functional data objects without loss of information. Such methodology
is known as functional data analysis and has a rich history. An excellent
introduction to this field is given in several books including Ramsay and
Silverman (2005), Horvath and Kokoszka (2012), and Srivastava and Klassen
(2016). An interesting aspect of most functional data is that the underlying
variability can be ascribed to two sources. These two sources are termed the
amplitude ($y$ or vertical) variability and the phase ($x$ or horizontal or
warping) variability. Capturing these two sources of variability is crucial
when modeling and monitoring functional data in a process control
architecture, and can greatly affect the construction of statistics (e.g.,
tolerance bounds, model parameters). In this work, we refer to methods that
handle both amplitude and phase variability in functional data as _elastic_.
This important concept is illustrated in Figure 1 through a simulated example.

In Figure 1 we have two functions that contain both a peak and a trough. Each
of the functions contain variability in the height of the peak and valley and
large variability in its placement. The relative heights of the peaks can be
attributed to the amplitude variability, while the different locations of the
peaks constitute the phase variability.

The phase variability can be accounted for by first aligning the functions. As
an example, the right panel shows time-aligned functions. The alignment
involves a transformation of the horizontal axis via the warping function
shown in the middle panel. The aligned function ($f_{1}(\gamma(t))$) captures
the amplitude variability while the warping functions ($\gamma$) captures the
phase variability. The Bayesian model calibration method introduced in this
paper considers the shape of the data by accounting for both directions of
variability.

Figure 1: Demonstration of amplitude and phase variability in functional data.
Left: Original functions. Middle: Warping Function (”phase variability”).
Right: Aligned functions (”amplitude variability”).

While standard calibration solutions are still applicable to misaligned
functions, especially when emulation is unnecessary, we demonstrate that using
functional metrics that consider the misaligned nature of the data can produce
more accurate calibration solutions and more reasonable model formulations.
Specifically, we use elastic functional data analysis methods, Srivastava and
Klassen (2016); Tucker et al. (2013); Marron et al. (2015), to construct a
metric to measure the distance between functions and to specify a likelihood.
We do this by decomposing our functional responses into aligned functions and
warping functions, as shown in Figure 1. For expensive models we can then
build an emulator or surrogate for the aligned functional responses, and,
under a suitable transformation, an emulator for the warping functions. These
two emulators can be used to calibrate in such a way that a proper distance is
used. The strengths of this approach are that (1) emulation is likely to be
more accurate when applied to aligned data instead of misaligned data (Francom
et al., 2022), (2) discrepancy modeling can be done in such a way that it is
isolated to the phase or the amplitude part of the model, and (3) the
alignment procedure ensures the isometry property holds and therefore, the
calibration is performed using a proper distance. Similarly, Kleiber et al.
(2014) proposed a functional calibration approach using a deformation that
relies on the $\mathbb{L}^{2}$ metric in the standard function space. However,
this metric is degenerate and exhibits a “pinching effect” which the elastic
methods do not have, as demonstrated in Srivastava and Klassen (2016). We
apply our approach to two dynamic material science applications described in
the following subsections.

###  1.1 Equation of state exploration using Z-machine experiments

This problem is motivated by data from Brown et al. (2014) which conducted
dynamic materials experiments on tantalum (Ta) generated with pulsed magnetic
fields using Sandia National Laboratories’ Z-machine, Savage et al. (2007). A
basic version of the experiments performed at the Z-machine can be seen in
Figure 2 (left). The Z-machine delivers electrical currents along an aluminum
(Al) panel creating massive pressure driving an impulse into a tantalum
sample. Stress waves flow through the tantalum, and the experiment results in
measurements of velocity on the far side of the sample as a function of time.
More details on the execution of the experiment can be found in Brown et al.
(2014) and Lemke et al. (2005). Figure 2 (right) shows the functional response
of these nine experiments.

Figure 2: Left: Example of experiments at the Z-machine. Right: Velocity
response functions from nine experiments.

The computer model input parameters consist of calibration parameters of key
interest and experiment-specific parameters that are not completely known.

###  1.2 Material strength exploration using gas gun experiments

Boteler and Dandekar (2006) performed a series of experiments where a flyer
was shot from a gas gun into a plate and the velocity of the opposite surface
of the plate was measured (as a function of time) as the resultant shock wave
moved through it. This experiment has similarities to the Z-machine example
above, except that the shock waves that propagate through the plate are driven
by the flyer’s impact instead of the more continuous electrical current drive.
Figure 3 shows three velocimetry curves measured during three of these flyer
plate gas gun experiments. Walters et al. (2018) used these experiments to
parameterize a strength model for aluminum. Also similar to the Z-machine
computer model, flyer plate impact simulations include both inputs to be
calibrated (strength model parameters) and experiment-specific parameters that
are uncertain (e.g., the exact velocity of the flyer).

Figure 3: From Walters et al. (2018), the response from three gas gun
experiments where both flyer and plate were aluminum alloys.

The structure of the rest of this paper is as follows. Section 2 gives a high
level overview of model calibration based on the Kennedy and O’Hagan (2001)
framework. Section 3.1 describes the types of variability for functional
output, functional alignment with a proper distance metric, and how to measure
the amount of variability in amplitude and phase space. Section 3.2 describes
the elastic Bayesian model calibration method as well as various modeling
choices in Section 3.3. Section 4 gives two simulation studies comparing the
proposed elastic functional calibration method with functional response
calibration that does not account for misalignment. Section 5 applies the
proposed method to the two material model calibration problems detailed above.
Lastly, Section 6 provides general conclusions and considerations on future
work.

##  2 Bayesian Model Calibration

###  2.1 Univariate Response

The traditional approach to Bayesian model calibration, introduced by Kennedy
and O’Hagan (2001) and expanded in Higdon et al. (2004), seeks to calibrate
parameters of a model using observations. Let $y(\bm{x},\bm{u})$ denote the
model where $\bm{x}$ denotes conditions that are certain, often fixed
conditions of an experiment, and $\bm{u}$ denotes uncertain parameters in need
of calibration. Let $z(\bm{x})$ denote an observation. In these approaches,
$\bm{x}$ could include functional variables like space or time and/or
conditions at which an experiment was performed. Let $n$ be the number of
observations or experimental conditions measured, so that our observed data
are $z(\bm{x}_{1}),\dots,z(\bm{x}_{n})$. Then the calibration model is

$\displaystyle z(\bm{x}_{i})$
$\displaystyle=y(\bm{x}_{i},\bm{\theta})+\delta(\bm{x}_{i})+\epsilon(\bm{x}_{i}),~{}~{}\epsilon(\bm{x}_{i})\sim\mathcal{N}(0,\sigma_{\epsilon}^{2})$
(1)

where $\bm{\theta}$ denotes the best set of calibration parameters, $\delta$
denotes (latent) error in the form of the model (often called model
discrepancy), and $\epsilon$ denotes measurement or observation error in $z$.
The Gaussian likelihood specified here assumes that measurement errors are
independent and identically distributed, but other measurement error models
can be used. The unknowns in this model are the calibration parameters
$\bm{\theta}$, the measurement error variance $\sigma^{2}_{\epsilon}$, and the
form of the discrepancy function. In order for this model to be identifiable,
priors for each of the unknowns need to be chosen carefully. If $n$ is small,
the prior for $\sigma^{2}_{\epsilon}$ will be influential. Additionally, there
is a natural trade off between the calibration parameters and the model
discrepancy, so that at least one of these needs to be well constrained by the
prior. For instance, we could use a Gaussian process prior for the discrepancy
function $\delta$ in such a way that we prefer small-valued smooth functions,
or we could try to specify priors that prefer positive values or monotone
functions (Higdon et al., 2004; Brynjarsdóttir and O’Hagan, 2014).

In many realistic scenarios the model $y$ is expensive to evaluate, requiring
powerful computing and non-trivial amounts of time. This slows down the
evaluation of the likelihood and makes inference impractical. In these
scenarios, we build a surrogate model (or emulator) to use in place of $y$,
which is trained using a (relatively) small number of model evaluations,
$\\{y(\bm{x}_{j},\bm{u}_{j})\\}_{j=1}^{N_{sim}}$ (Sacks et al., 1989). The
full Bayesian approach to inference then seeks to learn all the unknowns
(calibration parameters, model discrepancy, measurement error, and emulator
parameters) conditional on all of the data ($n$ observations and $N_{sim}$
model runs). Higdon et al. (2004) use this approach with a Gaussian process
emulator, while Kennedy and O’Hagan (2001) fix some of the Gaussian process
emulator and/or discrepancy parameters in advance for computational and
identifiability reasons. More generally, Liu et al. (2009) describe cases
where modularization of the model will result in better inference, and many
practitioners use modularized approaches either for philosophical or
computational reasons.

###  2.2 Functional Response

Various authors have extended the Bayesian model calibration approach above to
be more explicitly suited to functional response. This can vastly increase the
training data size, but Gaussian processes with Kronecker covariance
structures can help produce scalable models (Williams et al., 2006). A
somewhat different approach is used by Gu and Berger (2016), involving fitting
a Gaussian process for each output while sharing some parameters across
models, and, with a suitable discrepancy model, the calibration can be
obtained.

The most common approach is to project the functional response onto basis
functions and build functional models for calibration in the reduced dimension
basis coefficient space (Higdon et al., 2008; Bayarri et al., 2007; Francom et
al., 2019). Let $t$ denote the functional variable (e.g., time), and let
$z(t,\bm{x}_{i})$ denote an experimental measurement from the $i^{th}$
experiment at functional variable $t$. Similarly, let $y(t,\bm{x}_{i},\bm{u})$
denote a simulation of the $i^{th}$ experiment at the functional variable $t$
with input parameters $\bm{u}$. Then an approach to Bayesian model calibration
with functional response specifies

$\displaystyle z(t,\bm{x}_{i})$
$\displaystyle=y(t,\bm{x}_{i},\bm{\theta})+\delta(t,\bm{x}_{i})+\epsilon_{i}(t,\bm{x}_{i}),~{}~{}\epsilon(t,\bm{x}_{i})\sim\mathcal{N}(0,\sigma_{\epsilon}^{2}).$
(2)

The typical approach to inference is to discretize $t$ onto a grid
$t_{1},\dots,t_{N_{T}}$, which generates $N_{T}$ dimensional vectors
$\bm{z}(\bm{x}_{i})$, $\bm{y}(\bm{x}_{i},\bm{\theta})$, and
$\bm{\delta}(\bm{x}_{i})$ of the respective functions evaluated on the
discretized grid, which simplifies the model to a multivariate representation,

$\displaystyle\bm{z}(\bm{x}_{i})$
$\displaystyle=\bm{y}(\bm{x}_{i},\bm{\theta})+\bm{\delta}(\bm{x}_{i})+\bm{\epsilon}(\bm{x}_{i}),~{}~{}\bm{\epsilon}(\bm{x}_{i})\sim
N(\bm{0},\sigma_{\epsilon}^{2}\bm{I}).$ (3)

As in the univariate case, when $y$ is expensive to evaluate, we will require
an emulator or surrogate model in order to evaluate the likelihood function
quickly. Higdon et al. (2008) projects the model runs onto functional
principal components and projects the discrepancy onto a separate flexible
basis, and inference is carried out in the resulting low dimensional space.
Bayarri et al. (2007) project the model runs and discrepancy onto a wavelet
basis and carry out the inference in the low dimensional space. Francom et al.
(2019) projects the model runs onto functional principal components, allows
the discrepancy to be full-dimensional, and fits the emulator in a modular
fashion.

All of these calibration approaches can be applied to misaligned functional
data, but all will have difficulty with emulation and the specification of
model discrepancy, as they will only deal explicitly in amplitude variation.
Additionally, distances computed between the different functions will not be
proper distances since they would not account for the phase variability.

##  3 Elastic Bayesian Model Calibration

In this section, we will review elastic functional data analysis, introduce
our elastic approach to calibration, and suggest various modeling choices and
practical consierations.

###  3.1 Elastic Functional Data Analysis

####  3.1.1 Types of variability and model discrepancy in functional outputs

When model outputs are functional, we must consider two types of variability
in the functional outputs: amplitude and phase variability. Amplitude
variability is variability in the output for a fixed time (t), or, more
simply, y-axis variation. Phase variability is variability in time, or, more
simply, x-axis variability.

For computer model applications, input variables can induce both phase and
amplitude variability in the computer model outputs, resulting in
fundamentally different shapes over the range of plausible inputs.
Additionally, model discrepancy can result in an imperfect match between the
computer model prediction and observed data at the correct value of the model
inputs. With misaligned functional data, we must consider model discrepancy in
both phase and amplitude variability to accurately represent discrepancy-
induced shape distortions in the functional predictions.

When model discrepancy is not solely driven by amplitude variability,
pointwise (non-elastic) calibration (as in Section 2.2) will produce biased
estimates of the calibration parameters. Another large problem with the use of
the $\mathbb{L}^{2}$ metric is known as the pinching problem, Ramsay and Li
(1998). Specifically, if we have two functions, $f_{1}$ and $f_{2}$ and the
range$(f_{1})$ is entirely above the range$(f_{2})$, the $\mathbb{L}^{2}$
metric becomes degenerate and pinches the warped function. To address this
problem, Srivastava et al. (2011) introduced a mapping for functional data
called the square-root velocity function or SRVF that improves functional
alignment, and provides fundamental mathematical equalities that lead to the
formal development of this topic. Moreover, the metric used in the alignment
is a proper distance and avoids the pinching effects of the standard
$\mathbb{L}^{2}$ metric in function space without the use of a penalty. We
propose functional calibration metrics that account for both phase and
amplitude variation while properly measuring the distance between functions.

####  3.1.2 Functional Alignment

To explain metrics for comparing functional data in a calibration setting, we
simplify to the comparison of two functions of $t$, $z(t,\bm{x}_{i})$ and
$y(t,\bm{x}_{i},\bm{\theta})$. Varying $\bm{\theta}$ will change the shape of
$y(t,\bm{x}_{i},\bm{\theta})$ so we seek to find $\bm{\theta}$ such that
$z(t,\bm{x}_{i})$ and $y(t,\bm{x}_{i},\bm{\theta})$ are optimally matched,
where the optimality criterion considers the distance between the functions in
both amplitude and phase. For notational convenience, and to emphasize that we
are comparing these functions for a fixed $\bm{x}_{i}$ and $\bm{\theta}$, we
rewrite $z(t,\bm{x}_{i})$ and $y(t,\bm{x}_{i},\bm{\theta})$ as $z(t)$ and
$y(t)$ for the rest of this subsection. To measure the distance between $z(t)$
and $y(t)$, we use elastic functional data analysis (EFDA) (Srivastava and
Klassen, 2016). The main premise behind EFDA is to construct a proper distance
metric between the computational prediction $y(t)$ and the experimental data
$z(t)$. To construct this metric, a continuous mapping $\gamma_{y\rightarrow
z}(t):[0,1]\rightarrow[0,1]$ between $y(t)$ and the experimental data $z(t)$
is constructed such that $\gamma$ is a diffeomorphism. The function
$\gamma_{y\rightarrow z}(t)$ is referred to as a warping function, as it
measures phase distortions in $y(t)$ such that $y\circ\gamma_{y\rightarrow
z}(t)=y(\gamma_{y\rightarrow z}(t))$ aligns with $z(t)$. Srivastava et al.
(2011); Tucker et al. (2013) show that by applying a specific transformation
to the original functions $z(t)$ and $y(t)$, there exist simple expressions
for the amplitude and phase distance between functions. We now describe how to
construct this transformation and how to define the distance metrics on the
transformed data.

The functions $z(t)$ and $y(t)$ are transformed to their square-root velocity
functions (SRVFs). That is, we define the SRVF of $f(t)$ as:

$\displaystyle q(t)$
$\displaystyle=\operatorname{sgn}(\dot{f}(t))\sqrt{|\dot{f}(t)|}$ (4)

where $\dot{f}$ denotes the time derivative of $f$. The SRVF is a bijective
mapping up to a translation; that is, $f(t)$ can be uniquely determined from
$q(t)$ and a single point on the curve $f(t)$.

Let $q_{z}(t)$ and $q_{y}(t)$ denote the SRVFs of $z(t)$ and $y(t)$. The
warping function that aligns $y$ to $z$, denoted $\gamma_{y\rightarrow z}$,
can be estimated solving the following optimization problem via Dynamic
Programming (Tucker et al., 2013)

$\displaystyle\gamma_{y\rightarrow
z}=\operatorname*{arg\,inf}_{\gamma\in\Gamma}||q_{z}-(q_{y}\circ\gamma)\sqrt{\dot{\gamma}}||$
(5)

or using a Bayesian approach (Cheng et al., 2016; Lu et al., 2017).

An advantage of this approach to warping function estimation is that the
analyst does not have to specify landmarks for function alignment; the
estimation of $\gamma$ is achieved by using the group structure of $\Gamma$.
When we use the optimization in Equation 5 in later sections, we will refer to
this as a decomposition (specifically the warping decomposition) because it
decomposes a misaliged function into an aligned function and a warping
function.

Given an alignment function $\gamma(t)$ and aligned SRVF,
$(q_{y}\circ\gamma\sqrt{\dot{\gamma}})(t)$, we can construct measures of phase
and amplitude variability. Specifically, amplitude variability is measured as:

$\displaystyle d_{a}(q_{z},q_{y})$
$\displaystyle=||q_{z}-(q_{y}\circ\gamma_{y\rightarrow
z})\sqrt{\dot{\gamma}_{y\rightarrow z}}||^{2},$ (6)

Srivastava et al. (2011) show that this $\mathbb{L}^{2}$ distance on the
transformed and aligned SRVF functions is a proper distance metric for
amplitude.

####  3.1.3 Measuring phase distance

Defining a measure of phase variability is more difficult than for amplitude
variability, because the space of warping functions, $\Gamma$, is an infinite-
dimensional nonlinear manifold, and therefore cannot be treated as a standard
Hilbert space. Since we would like to exploit Riemannian-geometric structure
when making inferences about the warping functions, we again apply a specific
transformation to the warping functions such that we can use a standard
distance metric (the $\mathbb{L}^{2}$ norm) to measure distance. Specifically,
we represent an element $\gamma\in\Gamma$ by the square-root of its derivative
$\psi=\sqrt{\dot{\gamma}}$. Note that this is the same as the SRVF defined
earlier, and takes this form since $\dot{\gamma}>0$. The identity
$\gamma_{z\rightarrow z}$ maps to a constant function with value
$\psi_{z\rightarrow z}(t)=1$, which corresponds to no warping. Since
$\gamma(0)=0$, the mapping from $\gamma$ to $\psi$ is a bijection and one can
reconstruct $\gamma$ from $\psi$ using $\gamma(t)=\int_{0}^{t}\psi(s)^{2}ds$.
An important advantage of this transformation is that since
$\|\psi\|^{2}=\int_{0}^{1}\psi(t)^{2}dt=\int_{0}^{1}\dot{\gamma}(t)dt=\gamma(1)-\gamma(0)=1$,
the set of all such $\psi$’s is the positive orthant of the Hilbert sphere
$\Psi=\mathbb{S}_{\infty}^{+}$ (i.e., a unit sphere in the Hilbert space
$\mathbb{L}^{2}$). In other words, the square-root representation simplifies
the complicated geometry of $\Gamma$ to a unit sphere. The distance between
any two warping functions, i.e., the phase distance, is exactly the arc-length
between their corresponding SRVFs on the unit sphere $\mathbb{S}_{\infty}$:

$d_{p}(\gamma_{1},\gamma_{2})=d_{\psi}(\psi_{1},\psi_{2})\equiv\cos^{-1}\left(\int_{0}^{1}\psi_{1}(t)\psi_{2}(t)dt\right)\
.$

While the geometry of $\Psi\subset\mathbb{S}_{\infty}$ is more tractable, it
is still a nonlinear manifold and computing distances remains difficult.
Instead, we use a tangent (vector) space at a certain fixed point for further
analysis. The tangent space at any point $\psi\in\Psi$ is given by:
$T_{\psi}(\Psi)=\\{v\in\mathbb{L}^{2}|\int_{0}^{1}v(t)\psi(t)dt=0\\}$. To map
between the representation space $\Psi$ and tangent spaces, one requires the
exponential and inverse-exponential mappings. The exponential map at a point
$\psi\in\Psi$ denoted by $\exp_{\psi}:T_{\psi}(\Psi)\mapsto\Psi$, is defined
as

$\displaystyle\exp_{\psi}(v)$
$\displaystyle=\cos(\|v\|)\psi+\sin(\|v\|)\frac{v}{\|v\|},$ (7)

where $v\in T_{\psi}(\Psi)$. Thus, $\exp_{\psi}(v)$ maps points from the
tangent space at $\psi$ to the representation space $\Psi$. Similarly, the
inverse-exponential map, denoted by $\exp_{\psi}^{-1}:\Psi\mapsto
T_{\psi}(\Psi)$, is defined as

$\displaystyle\exp_{\psi}^{-1}(\psi_{1})$
$\displaystyle=\frac{\kappa}{\sin(\kappa)}(\psi_{1}-\cos(\kappa)\psi),$ (8)

where $\kappa=d_{p}(\gamma_{1},\gamma)$. This mapping takes points from the
representation space to the tangent space at $\psi$.

The tangent space representation $v$ is sometimes referred to as a _shooting
vector_. A sensible point for calibration on $\Psi$ is to define the tangent
space at the identity element $\psi_{z\rightarrow z}$. Since we are aligning
the simulations relative to the experimental data $z$, phase variability is
defined relative to $q_{z}$. Therefore, the warping function for the
experimental data is $\gamma_{z\rightarrow z}$ (no warping is required) and
the exponential map of this function is 0 by definition. Therefore, deviations
in $v$ from 0 represent deviation in phase from the experimental data.

###  3.2 Using Elastic FDA within Bayesian Model Calibration

With the metrics for distance defined in the previous section, we are ready to
detail how to do Bayesian model calibration with misaligned functional
responses, which we call elastic Bayesian model calibration.

We decompose the observations from experiment $i$ into aligned functions and
warping functions (each specific to experiment $i$) so that

$\displaystyle z(t,\bm{x}_{i})$
$\displaystyle=\tilde{z}(t,\bm{x}_{i})\circ_{t}\gamma_{z\rightarrow
z}(t,\bm{x}_{i})$ (9)

where $\circ_{t}$ emphasizes that the composition is only in $t$, such that
$f(a,b)\circ_{a}g(a,c)=f(g(a,c),b)$. This is, of course, the identity warping,
so nothing is happening in this step, but we are explicit about this to allow
for generalizations later. We similarly decompose each simulation (e.g.,
simulation $j$) with

$\displaystyle y(t,\bm{x}_{j},\bm{u}_{j})$
$\displaystyle=\tilde{y}(t,\bm{x}_{j},\bm{u}_{j})\circ_{t}\gamma_{y\rightarrow
z}(t,\bm{x}_{j},\bm{u}_{j}).$ (10)

using the warping decomposition of Equation 5. To facilitate modeling with
proper distance metrics, we transform the warping functions into shooting
vector space with

$\displaystyle\bm{v}_{z\rightarrow z}(\bm{x}_{i})$
$\displaystyle=\exp_{\psi}^{-1}\left(\sqrt{\dot{\gamma}_{z\rightarrow
z}(\bm{x}_{i})}\right)$ (11) $\displaystyle\bm{v}_{y\rightarrow
z}(\bm{x}_{j},\bm{u}_{j})$
$\displaystyle=\exp_{\psi}^{-1}\left(\sqrt{\dot{\gamma}_{y\rightarrow
z}(\bm{x}_{j},\bm{u}_{j})}\right).$ (12)

We can then use these aligned model runs and associated shooting vectors for
emulator building, if desired.

####  3.2.1 Calibration when Emulation is Unnecessary

Our calibration model using the aligned data and shooting vectors then has
likelihood

$\displaystyle\tilde{z}(t,\bm{x}_{i})$
$\displaystyle=\tilde{y}(t,\bm{x}_{i},\bm{\theta})+\delta_{\tilde{y}}(t,\bm{x}_{i})+\epsilon_{\tilde{z}}(t,\bm{x}_{i}),~{}~{}\epsilon_{\tilde{z}}(t,\bm{x}_{i})\sim\mathcal{N}(0,\sigma_{\tilde{z}}^{2})$
(13) $\displaystyle\bm{v}_{z\rightarrow z}(\bm{x}_{i})$
$\displaystyle=\bm{v}_{y\rightarrow
z}(\bm{x}_{i},\bm{\theta})+\bm{\delta}_{v}(\bm{x}_{i})+\bm{\epsilon}_{v}(\bm{x}_{i}),~{}~{}\bm{\epsilon}_{v}(\bm{x}_{i})\sim\mathcal{N}(0,\sigma_{v}^{2}\bm{I}).$
(14)

In order to do inference with this model, we discretize $t$ so that equation
13 can be re-written in a vector form,
$\tilde{\bm{z}}(\bm{x}_{i})=\tilde{\bm{y}}(\bm{x}_{i},\bm{\theta})+\bm{\delta}_{\tilde{y}}(\bm{x}_{i})+\bm{\epsilon}_{\tilde{z}}(\bm{x}_{i})$,
and the likelihood can then be written as

$\displaystyle
f\left(\tilde{\bm{z}}(\bm{x}_{1}),\dots,\tilde{\bm{z}}(\bm{x}_{n}),\bm{v}_{z\rightarrow
z}(\bm{x}_{1}),\dots,\bm{v}_{z\rightarrow
z}(\bm{x}_{n})~{}|~{}\bm{\theta},\sigma^{2}_{\tilde{z}},\sigma^{2}_{v},\bm{\beta}_{\tilde{y}},\bm{\beta}_{v}\right)=$
$\displaystyle\prod_{i=1}^{n}\left[\mathcal{N}\left(\tilde{\bm{z}}(\bm{x}_{i})~{}|~{}\tilde{\bm{y}}(\bm{x}_{i},\bm{\theta})+\bm{\delta}_{\tilde{y}}(\bm{x}_{i}),~{}\sigma_{\tilde{z}}^{2}\bm{I}\right)~{}\mathcal{N}\left(\bm{v}_{z\rightarrow
z}(\bm{x}_{i})~{}|~{}\bm{v}_{y\rightarrow
z}(\bm{x}_{i},\bm{\theta})+\bm{\delta}_{v}(\bm{x}_{i}),~{}\sigma_{v}^{2}\bm{I}\right)\right]$

where $\bm{\beta}_{\tilde{y}}$ and $\bm{\beta}_{v}$ parameterize the
discrepancy. With priors specified for $\bm{\theta}$,
$\sigma^{2}_{\tilde{z}}$, $\sigma^{2}_{v}$, $\bm{\beta}_{\tilde{y}}$, and
$\bm{\beta}_{v}$, the posterior

$\pi\left(\bm{\theta},\sigma^{2}_{\tilde{z}},\sigma^{2}_{v},\bm{\beta}_{\tilde{y}},\bm{\beta}_{v}~{}|~{}\tilde{\bm{z}}(\bm{x}_{1}),\dots,\tilde{\bm{z}}(\bm{x}_{n}),\bm{v}_{z\rightarrow
z}(\bm{x}_{1}),\dots,\bm{v}_{z\rightarrow z}(\bm{x}_{n})\right)$

is proportional to likelihood multiplied by priors and can be sampled with
MCMC. The main benefit of using this approach when an emulator is unnecessary
is that the discrepancy model can be specified in a more reasonable way, and
will be used with a proper distance metric in the likelihood. Note that this
approach requires the warping decomposition be called for each likelihood
evaluation. Even though this is just the decomposition of a single curve and
fairly fast, this expense can make emulation more desirable.

####  3.2.2 Calibration when Emulation is Necessary

Assume that the set of $N_{sim}$ model runs is decomposed into aligned
functions and warping functions. We can then use separate or joint emulators
with training inputs $\left\\{\bm{x}_{j},\bm{u}_{j}\right\\}_{j=1}^{N_{sim}}$
and outputs $\left\\{\tilde{\bm{y}}(\bm{x}_{j},\bm{u}_{j}),v_{y\rightarrow
z}(\bm{x}_{j},\bm{u}_{j})\right\\}_{j=1}^{N_{sim}}$. A full Bayesian approach
to emulation and calibration results in a joint likelihood of the observation
data and model runs while a modular Bayesian approach fits the emulator first
and then uses the emulator in the calibration model. Because these models
result in likelihoods and posteriors that are not very different from the
likelihood without an emulator, we omit the explicit likelihood here. For
Gaussian process emulators, the likelihood under the full Bayesian model can
be seen in Higdon et al. (2004, 2008).

###  3.3 Modeling Choices and Practical Considerations

To this point, we have discussed how using standard functional response
calibration tools with misaligned functional responses is possible but prone
to problems. We also mentioned that emulation techniques for functional
response do not perform well when applied to misaligned functional responses,
as demonstrated in Francom et al. (2022). We then showed how elastic
functional data analysis methods use proper distance metrics by aligning
functions in SRVF space and using (1) distance between aligned functions and
(2) distance between shooting vectors (transformed warping functions). Then we
introduced how to frame functional response computer model calibration such
that it uses these proper elastic metrics. Simply put, we use the warping
decomposition to separate our misaligned functional responses into aligned
functional responses and shooting vectors, and we use these two new datasets
instead of the original data when calibrating.

Below we discuss additional modeling choices and practical considerations for
using this methodology.

####  3.3.1 Warping decomposition

Uncertainty: Because we are using an optimization technique to obtain the
warping functions, we are fixing a part of the model that could be considered
uncertain. This is a modeling choice, similar to the choice of Kennedy and
O’Hagan (2001) to fix some emulator parameters at maximum likelihood values or
to the choice of Higdon et al. (2008) to not allow for uncertainty in the
basis representation of the functional response (i.e., the functional
principal components are fixed). If the warping functions were to be inferred
jointly with all of the other unknowns in a Bayesian framework, this would
lead to a much greater computational burden and would require much more
specialized calibration software. A possible shortcut could be to use modular
Bayes techniques to propagate uncertainty from the warping decomposition to
calibration uncertainty while cutting the feedback from the calibration to the
warping decomposition (Liu et al., 2009; Plummer, 2015).

Regularization: Another practical consideration is for the family of warping
functions allowed in the warping decomposition, $\Gamma$. If these functions
are not smooth or regularized, their variation can be difficult to predict
using the parameters, resulting in emulators with large residual variance.
However, if they are over regularized, they may not align the functional
responses enough. In case one wants to control the amount of warping or
elasticity this can be done as described in Wu and Srivastava (2011) using a
penalty on equation 5.

Choice of alignment reference: This decomposition can occur by either aligning
the model to the experiment, as described above, or by aligning to some other
common element (e.g., one of the model runs). Hence, the quantities
$\gamma_{y\rightarrow z}$, $\gamma_{z\rightarrow z}$, $\bm{v}_{y\rightarrow
z}$, and $\bm{v}_{z\rightarrow z}$ used above can be replaced with
$\gamma_{y\rightarrow y^{*}}$, $\gamma_{z\rightarrow y^{*}}$,
$\bm{v}_{y\rightarrow y^{*}}$ and $\bm{v}_{z\rightarrow y^{*}}$ for some
common reference $y^{*}(t)$ like the first model run. There are a few reasons
why this could be useful. First, the model runs have no discrepancy or noise,
which means that warping to them can be more stable. Second, exploration of
the computer model is often of interest even if there is no calibration data,
in which case emulation is frequently used and can be performed more
accurately when accounting for misalignment.

Modeling transformed aligned curves: Recall that we use the SRVF
transformation of the original curves to obtain the warping decomposition.
However, we opt not to build the model for the aligned data in the SRVF space.
Modeling in derivative space means that when models are transformed back to
the native space, errors are integrated and unrealistic heteroskedasticity
arises. Note that this problem is muted when transforming from shooting
vectors to warping functions because of the highly-constrained shape of the
warping functions.

####  3.3.2 Emulation

Francom et al. (2022) showed that taking alignment into account can
significantly improve emulator accuracy and efficiency, though their approach
relied on a model for landmarks (rather than shooting vectors) to build
warping functions. They found that, in many cases, training separate emulators
for aligned functions and warping functions is more desirable than training a
single emulator for both. This is because the warping function model will
often have larger unexplained variation (i.e., larger residual variance than
the aligned data model, even under appropriate standardization) because of the
latent nature of the warping functions. Subtle variations in the warping
function regularization level can result in fairly significant variations in
the shooting vectors, and while this variation is less pronounced when
transforming back into warping function space, the larger variation can
“corrupt” (Liu et al., 2009) the modeling of the aligned data unless careful
precautions are taken. The easiest precaution is to merely create the two
emulators independently. This works well in practice.

####  3.3.3 Discrepancy models

The decomposition of misaligned data into aligned data and shooting vectors
not only facilitates better emulation and calibration likelihood (distance
metrics), it also facilitates more realistic discrepancy modeling. For
instance, if the discrepancy is a time shift, that can be expressed through
the shooting vector discrepancy model $\bm{\delta}_{v}$. If the discrepancy is
a change to one feature of the curves, that can be handled directly through
the aligned data discrepancy model $\bm{\delta}_{\tilde{y}}$. The
specification of $\bm{\delta}_{\tilde{y}}$ can be reasoned about in a natural
way, for instance a set of basis functions could be used as in Higdon et al.
(2008). However, the specification of $\bm{\delta}_{v}$ is more nuanced
because it is in a transformed space. For instance, adding a constant in that
space (e.g., $\bm{\delta}_{v}(\bm{x}_{i})=\bm{1}$) has no effect because the
shooting vectors are scaled by their norm in the exponential map that
transforms them back to a SRVF on the unit sphere. Hence, a constant time
shift discrepancy is not achieved by adding a constant to the shooting
vectors. We have found that a piecewise linear basis for shooting vector
discrepancy is effective at capturing time shifts, as we will demonstrate in
the flyer plate analysis below.

####  3.3.4 Residual error models

We will typically assume independence between $\bm{\epsilon}_{\tilde{z}}$ and
$\bm{\epsilon}_{v}$, but this can be relaxed. Of greater interest is the
assumption of independence within $\bm{\epsilon}_{\tilde{z}}$ and
$\bm{\epsilon}_{v}$. As mentioned above, adding a correlation structure to the
residual model still results in appropriate distance metrics when the
Mahalanobis distance is used in the likelihood calculation. This can be done
explicitly to avoid the problems identified in Brown and Hund (2018).

####  3.3.5 Using Standard Functional Calibration Tools for Elastic
Calibration

Our formulation results in a typical Gaussian-likelihood pointwise-
calibration, though using the aligned responses and the shooting vectors
instead of the original misaligned functional responses (the likelihood
function for the elastic calibration is formed by combining vector versions of
equations 13 and 14). This means that the same tools that were discussed in
Section 2.2 can be used for elastic calibration, making the application of
these methods very practical. For instance, SEPIA (Gattiker et al., 2023)
implements the model of Higdon et al. (2008), and Francom et al. (2023)
implements a similar model but allows for modularization of the emulator as in
Francom et al. (2019) as well as MCMC that is robust to posterior
multimodality (via tempering). Both of these tools can be used to do elastic
calibration, and the warping decomposition can be achieved using fdasrsf
(Tucker, 2023).

Further, there are several choices available to emulate computer models with
functional response, some variations of which are described in Francom et al.
(2022); Hutchings et al. (2023). In the examples that follow, we use a
Bayesian multivariate adaptive regression spline (BMARS) (Francom et al.,
2018) emulator for both the aligned simulated curves and their corresponding
shooting vectors, paired with the calibration approach in Francom et al.
(2023). We chose BMARS due to its performance and the ability to scale to a
large number of training points which we have in our equation of state and
material strength data examples, though the methods are agnostic to emulator
choice.

##  4 Simulated Examples

We create two simulated examples where both will have three parameters to
calibrate and misaligned functional response simulated using a Gaussian
density function. One of the examples is formulated without model discrepancy
and the other with discrepancy.

###  4.1 Example 1

The first simulated data set is constructed from the following model where
each function is a parameterized Gaussian pdf:

$y(t,\bm{u})=\frac{u_{1}}{0.05\sqrt{2\pi}}\exp\left(-\frac{1}{2}\left(\frac{t-(\sin(2\pi
u_{0}^{2})/4-u_{0}/10+0.5)}{0.05}\right)^{2}\right)+0u_{2}.$

The calibration parameters are $\bm{u}=[u_{0},u_{1},u_{2}]$, though the model
only uses the first two of these parameters. There are no experimental
parameters ($\bm{x}$), so we omit these from the notation used in the previous
sections. We generated a set of 100 functions (model runs) using
$\bm{u}_{1},\dots,\bm{u}_{100}$ where each $\bm{u}_{j}$ was sampled uniformly
within the unit cube. The ‘experimental data’ $z(t)$ was generated using the
parameter values $\bm{u}^{*}=[0.1028,0.5930,0]$ without any model discrepancy,
so that $z(t)=y(t,\bm{u}^{*})$. Figure 4 presents the simulated model runs in
grey and the experimental data shown in black.

Figure 4: Simulated curves $y(t,\bm{u}_{1}),\dots,y(t,\bm{u}_{100})$ and
experiment data $z(t)$.

We separated the phase and amplitude by aligning the simulated curves to the
experimental data curve utilizing the elastic methodology from Section 3.1.
Figure 5(a) presents the aligned curves to the experimental data, while Figure
5(b) presents the corresponding warping functions of both the simulated curves
($\gamma_{y\rightarrow z}(t,\bm{u}_{j})$ for $j=1,\dots,100$) and the
experimental data ($\gamma_{z\rightarrow z}(t)$, which is the identity).
Figure 5(c) presents the corresponding shooting vectors for the simulation
($\bm{v}_{y\rightarrow z}(\bm{u}_{j})$) and the experimental data
($\bm{v}_{z\rightarrow z}$). As noted previously, if the simulated curves are
aligned to the experiment data, the corresponding warping function will be the
identity function and the shooting vectors will all be equal to zero.

(a) Aligned simulated functions.

(b) Warping functions.

(c) Shooting vectors.

Figure 5: Alignment of simulated curves to experimental data.

We then built separate BMARS emulators for the aligned curves and the shooting
vectors and performed a (modular) elastic Bayesian model calibration as
described in Section 3.2 to infer the posterior.

(a) Calibrated and aligned curves.

(b) Calibrated warping functions.

(c) Calibrated shooting vectors.

Figure 6: Predictive posterior samples after calibration of the simulated
data.

Figure 6 presents some posterior predictive samples after calibration of the
aligned curves, the corresponding warping functions and shooting vectors (in
blue). Figure 7 presents some posterior predictive samples in the original
data space (in blue) with the simulated curves/model runs (in grey) and the
experimental data (in black). The predictive samples cover the experimental
data with small uncertainty, showing good performance of the emulator at
$\bm{\theta}$. A pairwise plot summarizing the posterior samples of
$\bm{\theta}$ (marginal and bivariate) is shown in Figure 8. The true value of
$\bm{\theta}$ is represented by the mark ($x$) in the lower triangular and the
vertical line on the diagonal. The true parameter values are covered by high
density regions in the posterior, indicating they are recovered well by the
model. Additionally, the distribution of the nuisance parameter resembles a
uniform distribution (prior) and does not affect the calibration.

Figure 7: Predictive posterior samples after calibration of simulated data in
the original data space. Figure 8: Pair plot of parameters after calibration
of simulated data with truth shown by the x and vertical line.

To compare our elastic method to the standard method, we fitted a BASS
emulator on the original simulated computer model and directly performed a
modular Bayesian model calibration. Figure 9 presents the posterior predictive
samples after calibration of the simulated data in original data space. We see
that the predictions do not resemble the experimental data well and have more
than one peak, given that phase variability is not taken into account. The
underlying computed mean from the emulator is not representative of the
underlying data. Additionally, Figure 10 presents the pairwise plot of the
calibrated parameters with the true value shown by the mark. The resulting
posterior distribution covers the the true value of the parameters, but with
much larger uncertainty. In comparison to Figure 8, we can immediately notice
the benefit of adopting BMC with an elastic framework. This is more clearly
shown in Figure 11 where the marginal posteriors and the bivariate posterior
contours (at 95%) are being compared for the traditional calibration (not
aligned) and the elastic functional data calibration (aligned).

Figure 9: Predictive posterior samples after calibration of the simulated data
in original data space when the emulator was fitted on the original data
space.  Figure 10: Pairs plot of the posterior samples of the parameters after
calibration of the simulated data where the true value is marked by the x and
vertical line. The emulator was fitted directly on the original data. Figure
11: Pairs plot of the marginal posteriors and bivariate posterior contours
(95%) for the calibrations with alignment and no alignment.

###  4.2 Example 2

To explore how the methodology handles model discrepancy, another simulated
data set was generated based on the model,

$y(t,\bm{u})=\frac{u_{1}}{0.05\sqrt{2\pi}}\exp\left(-\frac{1}{2}\left(\frac{t-(u_{0}+0.1)}{0.05}\right)^{2}\right).$

for which 300 functions were generated with $u_{0}$ and $u_{1}$ being drawn
from the $U[0,0.3]$ distribution. These functions are considered as the
“computer model output.” An additional nuisance parameter, $u_{2}$, was drawn
from the standard uniform distribution $U[0,1]$. The experimental data was
generated according to

$z(t)=\frac{u^{*}_{1}}{0.05\sqrt{2\pi}}\exp\left(-\frac{1}{2}\left(\frac{t-(u^{*}_{0}+0.3)}{0.05}\right)^{2}\right).$

where $\bm{u}^{*}=[u^{*}_{0},u^{*}_{1}]=[0.3,0.2]$, where the shift in time of
$0.2$ indicates model discrepancy.

Similar to the previous example, the 300 generated functions were aligned to
the experimental data and emulators were built based on the aligned functions
and shooting vectors, which were used for calibration. A more traditional
approach to this problem would be to pick out important features of the
misaligned curves and do the emulation and calibration (including discrepancy)
using the features. As an example, we compare against a calibration that
reduces each curve to its maximum value and the timing of the maximum value.
We then build an emulator to jointly model these features, and we specify a
constant shift discrepancy model.

Due to the inclusion of the discrepancy model, posterior predictive samples
from both the elastic and standard method have good agreement and coverage of
the experimental curve. As shown in the right panel of Figure 12, the
posteriors of $\bm{\theta}$ for the standard method have a credible interval
that cover the true parameter values. However, the left panel shows that the
posteriors of $\theta$ from the elastic functional method have also good
coverage and provide tighter credible intervals. This is more clearly shown in
Figure 13 where the posteriors from both approaches are compared via marginal
distributions and bivariate contours.

(a) Elastic Method

(b) Standard Method

Figure 12: Predictive posterior samples after calibration of the simulated
data with discrepancy in the original data space. Figure 13: Pairs plot of the
marginal posteriors and bivariate marginal contours (95%) for the calibration
that used only features and the one that used the proposed methodology.

##  5 Dynamic Material Model Calibration

We explore the application of the elastic functional Bayesian calibration
method to two real world applications, a tantalum equation of state
calibration and a material strength calibration for aluminum alloy.

###  5.1 Tantalum Z-Machine

In this application, we seek to calibrate the equation of state (EoS) of
tantalum (Ta) generated from pulsed magnetic fields (Brown et al., 2014). We
are aiming to estimate parameters describing the compressibility (relationship
between pressure and density) to understand better how materials compress to
extreme pressures. Ta is an ideal material for this study as it is able to
remain in its initial crystal structure to pressures up to 10 million times
standard atmospheric pressure (Soderlind and Moriarty, 1998). A description of
the Ta experiments is shown in Figure 2. The experiments were conducted using
Sandia National Laboratories’ Z-machine, which is a pulsed power drive that
can deliver massive electric currents over short time scales. These currents
were forced to flow along an aluminum (Al) panel, producing a large magnetic
pressure which drives a time-dependent stress wave (impulse) into the system.
Ta samples and transparent lithium fluoride (LiF) windows were glued to the
panel such that the stress wave propagates sequentially through each of these
materials.

The material properties are modeled using a physically motivated form given by
Vinet et al. (1989). This form describes the pressure-density $(P-\rho)$
response as

$P(\rho)=3B_{0}\left(\frac{1-\eta}{\eta^{2}}\right)\exp\left(\frac{3}{2}(B^{\prime}_{0}-1)(1-\eta)\right),$

where $\eta=\sqrt[3]{\rho_{0}/\rho}$, $\rho_{0}$ is the initial density and
$B_{0}$ and $B^{\prime}_{0}$ are the bulk modulus and its pressure derivative
at ambient conditions. From the computer experiment perspective, we will work
with 6 inputs (3 EoS parameters ($\rho,B_{0},B^{\prime}_{0}$) and 3
experiment-specific parameters) and output velocity curves on a grid of 100
equidistant time points. As described in (Brown et al., 2014), the main goal
is to provide inferences on the 3 EoS parameters with uncertainty
quantification and to propagate these inferences to the Vinet model.

As it was done in the simulation examples, the computer model output was
aligned to the experimental data, and an emulator was fitted to the aligned
computer model and the corresponding shooting vectors. We then performed a
modular elastic Bayesian model calibration utilizing the formulation above.
There are a total of 9 experiments and the 3 EoS parameters were calibrated
across the experiments and 4 experiment specific parameters.

Figure 14 presents the calibration results for the Ta experiments. The black
curves shown in the figure correspond to the experimental velocity curves for
the 9 experiments. The shaded colored regions are the 95% prediction intervals
that result from the elastic functional Bayesian calibration. The prediction
intervals exhibit good agreement with each experimental curve. The residuals
defined as the difference between the experimental data and the calibrated
predicted mean, are shown in the right panel of Figure 14. Each of the
residuals are color coded to the corresponding experiment from the left panel.
The full functional approach has tighter coverage of the experimental curves
and smaller residual values compared to the resulting predictions and
residuals that stem from the scaled ESS approach of Brown and Hund (2018).

Figure 14: (Left) Experiment velocities of Ta shown in black compared with 95%
prediction intervals from the elastic functional Bayesian calibration.
(Right). Corresponding color-coded residuals (difference between experiment
and calibration mean) for each experiment

Furthermore, Figure 15 presents a pairwise plot of the posterior distribution
samples of the calibrated EoS parameters for Ta with the elastic approach. The
posteriors are well concentrated within the parameter space and the medians
compare to those found in Brown and Hund (2018), with a reduced uncertainty
around the same parameter values. It should be remarked that the elastic FDA
calibration approach does not involve any likelihood scaling. To assess the
calibration further, we compared the posterior estimates of the EoS to those
reported in the literature for the Vinet model. Figure 16 presents a plot of
the pressure versus density Vinet curve, with a 95% credible interval
generated from the elastic calibration shown in red. The blue curve
corresponds to the loading path determined analytically by using the standard
techniques in the dynamic materials community as described in (Brown et al.,
2014). The dashed curve is the state of the art theoretical calculation given
by Greeff et al. (2009), which has been used to simulate these types of Ta
experiments. The posterior estimate of the curve from the elastic calibration
lies just above the theoretical calculation and below the analytical result,
which corresponds with what is expected when using the analytic analysis
(Kraus et al., 2016). The estimates of the physical parameters resulting from
the calibration tend to agree well with those from previous work and provide
uncertainty quantification.

Figure 15: Pair plot of the posterior densities of EoS parameters for Tantalum
Figure 16: Calibrated Ta material response to the Vinet model compared with
the analytic results in Brown et al. (2014) and the theoretical calculations
in Greeff et al. (2009).

###  5.2 Flyer Plate Impact

Material strength characterizes how a material temporarily or permanently
deforms as it experiences pressure. This is of interest in various areas of
science and engineering, with applications in the aerospace, medical, and
automotive industries (Gray et al., 2005). When experimentation is difficult,
material strength models are an important tool for predicting how a material
will react to pressure. One such model is the Johnson-Cook strength model
(Johnson and Cook, 1983), for which the model behavior is dictated by a
relatively small collection of material-specific physical constants. To set
these constants, scientists and engineers rely on experimentation, hence
calibrating the material strength model parameters to experimental
measurements.

In this application, we consider the calibration of the Johnson-Cook material
strength model for an aluminum alloy using a set of plate impact experiments.
Walters et al. (2018) performed a Bayesian model calibration using plate
impact experiments by reducing the velocimetry curves to a small set of
features considered to be important to strength. We show how elastic Bayesian
model calibration relies on the entire curve without using human-intensive
feature engineering. Plate impact experiments achieve high pressure on a
material sample by shooting a flyer at high velocity into the sample (plate).
Lasers measure how the free surface of the sample moves as the shock wave
moves through it. The result is a trace of the velocity of the free surface of
the sample over time, called a velocimetry curve. The first column of Figure
17 shows measured velocimetry curves along with 1000 simulated velocimetry
curves using different settings of the Johnson-Cook model and for three
experiments. To obtain these simulations, the Johnson-Cook model is used
within a larger hydrodynamics code that is expensive to evaluate. Additionally
Figure 17 shows the aligned and warping functions along with the shooting
vectors corresponding to the measured and simulated curves. The misalignment
of the experimental curves to the simulated is clear.

Figure 17: Measured and simulated velocimetry curves. Aligned, warping
functions and shooting vectors Figure 18: Basis functions used to capture a
time shift discrepancy in shooting vector space for the flyer plate example.

(a) Original Data and 1000 simulated velocimetry curves

(b) Calibrated Posterior Predictive Samples for the aligned curves

(c) Calibrated Posterior Predictive Samples for the warping functions.

(d) Calibrated Posterior Predictive Samples for the shooting vectors.

Figure 19: (a) The original measured velocimetry curve with 1000 computer
model runs. (b) Calibrated posterior predictive samples of velocimetry curve;
experimental data and model runs.

(a) Posterior discrepancy in shooting vector space.

(b) Posterior discrepancy in warping function space.

From the computer experiment perspective, we will work with 11 inputs: 5
Johnson-Cook parameters plus 6 experiment-specific parameters (two for each
experiment) and the output velocimetry curves on a grid of 200 time points.
The 1,000 velocimetry simulations are aligned to the experimental data, used
to build an emulator and perform the elastic calibration. To allow for a time
shift discrepancy, especially relevant to experiments 2 and 3, we use three
piecewise constant and one piecwise linear basis functions in shooting vector
space as shown in Figure 18. More specifically, we parameterize
$\bm{\delta}_{v}(\bm{x}_{i})=\bm{D}\bm{d}(\bm{x}_{i})$ where $\bm{D}$ is the
$200\times 4$ matrix of basis functions and $\bm{d}(\bm{x}_{i})$ is a vector
of $4$ basis coefficients that can vary by experiment. In the approach of
Higdon et al. (2008), $d_{k}(\cdot)$ would be given a Gaussian process prior
so that discrepancy could be predicted for a new experiment (with experimental
settings $\bm{x}^{*}$). In this case, we are not interested in predicting the
discrepancy at new experimental settings, so we specify a standard normal
prior for each $d_{k}(\bm{x}_{i})$. This indicates that we want to allow for a
time shift, but would prefer to have none ($d_{k}(\bm{x}_{i})=0$) if possible.
The prior variance of $d_{k}(\bm{x}_{i})$ can be modified in order to favor
more or less discrepancy. The form of the basis functions in Figure 18 is
carefully chosen to be constant in domains where the aligned data are roughly
constant and a linear function in the region where the aligned data are most
active.

Figure 19(a) presents the experimental data, the computer simulated models
shown in grey and the posterior predictive samples of the velocimetry curve
after calibration shown in blue. In a similar way, Figures 19(b), 19(c) and
19(d) show the posterior predictive samples of the aligned curves, warping
functions and shooting vectors respectively. In all of these cases, the
predictive samples have good coverage of the experimental data. Furthermore,
Figure 21 presents a pairs plot of the posterior distribution of the 11
calibrated parameters with diagonal elements representing marginal
distributions while the lower diagonal are bivariate contours and pair sample
plots for the upper diagonal. Compared to Walters et al. (2018), we are able
to use the entirety of the functional data while they used a few hand selected
features.

Figure 20(a) shows the posterior distribution of shooting vector
discrepancies, $\bm{\delta}_{v}(\bm{x}_{1})$, $\bm{\delta}_{v}(\bm{x}_{2})$,
and $\bm{\delta}_{v}(\bm{x}_{3})$. Figure 20(b) shows the difference between
two sets of warping functions: those that include the discrepancy in Figure
20(a) and those that do not, indicating the type of shift that these
discrepancies induce.

Figure 21: Pairwise plot and marginal densities for the posterior distribution
of the 11 calibrated parameters for the flyer plate experiment.

##  6 Conclusion

Model calibration for functional responses is typically more difficult than in
more traditional settings. Adjusting input parameters so model output matches
experimental output despite discrepancy is seldom trivial, but can be even
more challenging for functional data when altering input parameters results in
amplitude and phase variability in the responses. Traditional methods ignore
these aspects which are unique to models with functional responses, putting
them at risk for higher levels of bias and lower levels of efficiency. In this
paper we develop methods to handle amplitude and phase variability in a
systematic way, resulting in more efficient and more accurate estimates of the
model parameters. The improvements in these estimators is achieved through
better handling of the functional responses as opposed to collecting more
data.

The elastic Bayesian model calibration procedure presented here uses
information from both the amplitude and phase space to calibrate the
parameters, in contrast with more common calibration methods which would only
look at the amplitude space. The benefits of this approach are demonstrated in
simulation studies where amplitude and phase variability are present. The
simulation studies show the elastic Bayesian model calibration approach
results in superior estimation of the model parameters and the predicted
functional responses. Information about the amount of warping necessary to
align the functions provides an additional indirect benefit, as it provides
valuable insight into the model’s input parameters. Standard methods which do
not handle the phase and amplitude variability separately, may result in
predicted functional responses which do not fit the data well and suggest much
larger uncertainty in calibration parameters. These results make a strong case
for the elastic modeling approach, as most functional data have some
misalignment and phase discrepancy.

Applying the method to the Z-machine data and flyer-plate data generated
similar results to previous studies without the need of likelihood scaling as
in Brown and Hund (2018) or feature selection as in Walters et al. (2018). The
elastic Bayesian calibration model yields credible intervals with good
agreement of the experimental data using a principled approach. From a
theoretical perspective, the elastic Bayesian calibration approach is
satisfying in that it treats the model’s output as functional throughout the
calibration procedure.

Future work may investigate some of the additional modeling assumptions made
in the elastic calibration approach. For example, the modeling assumption that
$\tilde{y}$ and $v$ are independent needs further assessment. Additionally,
uncertainty introduced from the warping procedure is not considered in the
analysis presented here. A more complete treatment of this uncertainty may be
integrated throughout the calibration procedure. However, the additional
uncertainty from warping is small in the noiseless cases presented in this
paper.

## Acknowledgment

This paper describes objective technical results and analysis. Any subjective
views or opinions that might be expressed in the paper do not necessarily
represent the views of the U.S. Department of Energy or the United States
Government. This work was supported by the Laboratory Directed Research and
Development program at Sandia National Laboratories; a multi-mission
laboratory managed and operated by National Technology and Engineering
Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell
International, Inc., for the U.S. Department of Energy’s National Nuclear
Security Administration under contract DE-NA0003525.

## References

  * Bayarri et al. [2007] M. J. Bayarri, J. O. Berger, R. Paulo, J. Sacks, J. A. Cafeo, J. Cavendish, C. Lin, and J. Tu.  A framework for validation of computer models.  _Technometrics_ , 49(2), 2007.
  * Boteler and Dandekar [2006] J. M. Boteler and D. P. Dandekar.  Dynamic response of two strain-hardened aluminum alloys.  _Journal of applied physics_ , 100(5):054902, 2006\.
  * Brown and Hund [2018] J. Brown and L. Hund.  Estimating material properties under extreme conditions by using bayesian model calibration with functional outputs.  _Journal of the Royal Statistical Society, Series C_ , 67(4):1023–1045, 2018.
  * Brown et al. [2014] J. L. Brown, C. S. Alexander, J. R. Asay, T. J. Vogler, D. H. Dolan, and J. L. Belof.  Flow strength of tantalum under ramp compression to 250 gpa.  _Journal of Applied Physics_ , 115(4):043530, 2014\.  doi: 10.1063/1.4863463.
  * Brynjarsdóttir and O’Hagan [2014] J. Brynjarsdóttir and A. O’Hagan.  Learning about physical parameters: The importance of model discrepancy.  _Inverse Problems_ , 30(11):114007, 2014.
  * Cheng et al. [2016] W. Cheng, I. L. Dryden, X. Huang, et al.  Bayesian registration of functions and curves.  _Bayesian Analysis_ , 11(2):447–475, 2016.
  * Francom et al. [2018] D. Francom, B. Sanso, A. Kupresanin, and G. Johannesson.  Sensitivity analysis and emulation for functional data using bayesian adaptive splines.  _Statistica Sinica_ , 28(2):791–816, 2018.
  * Francom et al. [2019] D. Francom, B. Sanso, V. Bulaevskaya, D. Lucas, and M. Simpson.  Inferring atmospheric release characteristics in a large computer experiment using bayesian adaptive splines.  _Journal of the American Statistical Association_ , 114(528):1450–1465, 2019.
  * Francom et al. [2022] D. Francom, B. Sanso, and A. Kupresanin.  Landmark-warped emulators for models with misaligned functional response.  _SIAM/ASA Journal on Uncertainty Quantification_ , 10(1):125–150, 2022.
  * Francom et al. [2023] D. Francom, P. Trubey, et al.  LANL/impala, 2023.  URL https://github.com/lanl/impala.
  * Gattiker et al. [2023] J. Gattiker, N. Klein, E. Lawrence, and G. Hutchings.  LANL/SEPIA, 2023.  URL https://doi.org/10.5281/zenodo.4048801.
  * Gray et al. [2005] G. T. R. Gray, P. J. Maudlin, L. M. Hull, Q. K. Zuo, and S. Chen.  Predicting material strength, damage, and fracture the synergy between experiment and modeling.  _Journal of Failure Analysis and Prevention_ , 5(3):7–17, 2005.
  * Greeff et al. [2009] C. W. Greeff, S. P. Rudin, S. D. Crockett, and J. M. Wills.  The cold equation of state of tantalum.  _Melville: American Institute of Physics_ , pages 681–684, 2009.
  * Gu and Berger [2016] M. Gu and J. O. Berger.  Parallel partial gaussian process emulation for computer models with massive output.  _The Annals of Applied Statistics_ , 10(3):1317–1347, 2016.
  * Higdon et al. [2004] D. Higdon, M. Kennedy, J. C. Cavendish, J. A. Cafeo, and R. D. Ryne.  Combining field data and computer simulations for calibration and prediction.  _SIAM Journal on Scientific Computing_ , 26(2):448–466, 2004.
  * Higdon et al. [2008] D. Higdon, J. Gattiker, B. Williams, and M. Rightley.  Computer model calibration using high-dimensional output.  _Journal of the American Statistical Association_ , 103(482), 2008.
  * Horvath and Kokoszka [2012] L. Horvath and P. Kokoszka.  _Inference for Functional Data with Applications_.  Springer, 2012.
  * Hutchings et al. [2023] G. Hutchings, B. Sansó, J. Gattiker, D. Francom, and D. Pasqualini.  Comparing emulation methods for a high-resolution storm surge model.  _Environmetrics_ , 34(3):e2796, 2023.
  * Johnson and Cook [1983] G. R. Johnson and W. H. Cook.  A constitutive model and data for metals subjected to large strains, high strain rates, and high temperatures.  _Proceedings 7th International Symposium on Ballistics_ , pages 541–547, 1983.
  * Kennedy and O’Hagan [2001] M. C. Kennedy and A. O’Hagan.  Bayesian calibration of computer models.  _Journal of the Royal Statistical Society. Series B, Statistical Methodology_ , pages 425–464, 2001.
  * Kleiber et al. [2014] W. Kleiber, S. R. Sain, and M. J. Wiltberger.  Model calibration via deformation.  _SIAM/ASA Journal on Uncertainty Quantification_ , 2(1):545–563, 2014.
  * Kraus et al. [2016] R. G. Kraus, J. P. Davis, C. T. Seagle, D. E. Fratanduono, D. C. Swift, J. L. Brown, and J. H. Eggert.  Dynamic compression of copper to over 450 gpa: a high-pressure standard.  _Phys. Rev. B_ , 93(13):134105, 2016.
  * Lemke et al. [2005] R. W. Lemke, M. D. Knudson, D. E. Bliss, K. Cochrane, J. Davis, A. A. Giunta, H. C. Harjes, and S. A. Slutz.  Magnetically accelerated, ultrahigh velocity flyer plates for shock wave experiments.  _Journal of Applied Physics_ , 98(7):073530, 2005\.  doi: 10.1063/1.2084316.
  * Liu et al. [2009] F. Liu, M. J. Bayarri, J. O. Berger, et al.  Modularization in Bayesian analysis, with emphasis on analysis of computer models.  _Bayesian Analysis_ , 4(1):119–150, 2009.
  * Lu et al. [2017] Y. Lu, R. Herbei, and S. Kurtek.  Bayesian registration of functions with a gaussian process prior.  _Journal of Computational and Graphical Statistics_ , 26(4):894–904, 2017.
  * Marron et al. [2015] J. S. Marron, J.O. Ramsay, L.M. Sangalli, and A. Srivastava.  Functional data analysis of amplitude and phase variation.  _Statistical Science_ , 30(4):468–484, 2015.
  * Plummer [2015] M. Plummer.  Cuts in bayesian graphical models.  _Statistics and Computing_ , 25(1):37–43, 2015\.
  * Ramsay and Li [1998] J. O. Ramsay and X. Li.  Curve registration.  _Journal of the Royal Statistical Society, Ser. B_ , 60(2):351–363, 1998.
  * Ramsay and Silverman [2005] J. O. Ramsay and B. W. Silverman.  _Functional Data Analysis_.  Springer, 2005.
  * Sacks et al. [1989] J. Sacks, W. J. Welch, T. J. Mitchell, and H. P. Wynn.  Design and analysis of computer experiments.  _Statistical Science_ , pages 409–423, 1989.
  * Savage et al. [2007] M. E. Savage, L. F. Bennett, D. E. Bliss, W. T. Clark, R. S. Coats, J. M. Elizondo, K. R. LeChien, H. C. Harjes, J. M. Lehr, J. E. Maenchen, et al.  An overview of pulse compression and power flow in the upgraded z pulsed power driver.  In _2007 16th IEEE International Pulsed Power Conference_ , volume 2, pages 979–984. IEEE, 2007.  doi: 10.1103/PhysRevSTAB.13.010402.
  * Soderlind and Moriarty [1998] P. Soderlind and J. A. Moriarty.  First-principles theory of ta up to 10 mbar pressure: structural and mechanical properties.  _Phys. Rev. B,_ , 57:10340–10350, 1998.
  * Srivastava and Klassen [2016] A. Srivastava and E. Klassen.  _Functional and Shape Data Analysis_.  Springer, 2016.
  * Srivastava et al. [2011] A. Srivastava, W. Wu, S. Kurtek, E. Klassen, and J. S. Marron.  Registration of functional data using Fisher-Rao metric.  _arXiv:1103.3817v2 [math.ST]_ , 2011.  URL http://arxiv.org/abs/1103.3817v2.
  * Tucker [2023] J. D. Tucker.  fdasrsf, 2023.  URL https://github.com/jdtuck/fdasrsf_python.
  * Tucker et al. [2013] J. D. Tucker, W. Wu, and A. Srivastava.  Generative models for functional data using phase and amplitude separation.  _Computational Statistics and Data Analysis_ , 61:50–66, 2013.
  * Vinet et al. [1989] P. Vinet, J. Rose, J. Ferrante, and J. Smith.  Universal features of the equation of state of solids.  _J. Phys. Condnsd Matt._ , 1:1941–1963, 1989.
  * Walters et al. [2018] D. J. Walters, A. Biswas, E. C. Lawrence, D. C. Francom, D. J. Luscher, D. A. Fredenburg, K. R. Moran, C. M. Sweeney, R. L. Sandberg, J. P. Ahrens, et al.  Bayesian calibration of strength parameters using hydrocode simulations of symmetric impact shock experiments of al-5083.  _Journal of Applied Physics_ , 124(20):205105, 2018.
  * Williams et al. [2006] B. Williams, D. Higdon, J. Gattiker, L. Moore, M. McKay, S. Keller-McNulty, et al.  Combining experimental data and computer simulations, with an application to flyer plate experiments.  _Bayesian Analysis_ , 1(4):765–792, 2006.
  * Wu and Srivastava [2011] W. Wu and A. Srivastava.  An information-geometric framework for statistical inferences in the neural spike train space.  _Journal of Computational Neuroscience_ , 31:725–748, 2011\.