# From Displacements to Distributions: A Machine-Learning Enabled Framework
for Quantifying Uncertainties in Parameters of Computational Models

Taylor Roper  University of Colorado Denver
Department of Mathematical and Statistical Sciences
1201 Larimer St
Denver, CO, 80204, USA<EMAIL_ADDRESS>Harri Hakula  Aalto
University
Department of Mathematics and System Analysis
P.O. Box 11100
FI–00076 Aalto, Finland<EMAIL_ADDRESS>Troy Butler  University of
Colorado Denver
Department of Mathematical and Statistical Sciences
1201 Larimer St
Denver, CO, 80204, USA<EMAIL_ADDRESS>
###### Abstract

This work presents novel extensions for combining two frameworks for
quantifying both aleatoric (i.e., irreducible) and epistemic (i.e., reducible)
sources of uncertainties in the modeling of engineered systems. The data-
consistent (DC) framework poses an inverse problem and solution for
quantifying aleatoric uncertainties in terms of pullback and push-forward
measures for a given Quantity of Interest (QoI) map. Unfortunately, a pre-
specified QoI map is not always available a priori to the collection of data
associated with system outputs. The data themselves are often polluted with
measurement errors (i.e., epistemic uncertainties), which complicates the
process of specifying a useful QoI. The Learning Uncertain Quantities (LUQ)
framework defines a formal three-step machine-learning enabled process for
transforming noisy datasets into samples of a learned QoI map to enable DC-
based inversion. We develop a robust filtering step in LUQ that can learn the
most useful quantitative information present in spatio-temporal datasets. The
learned QoI map transforms simulated and observed datasets into distributions
to perform DC-based inversion. We also develop a DC-based inversion scheme
that iterates over time as new spatial datasets are obtained and utilizes
quantitative diagnostics to identify both the quality and impact of inversion
at each iteration. Reproducing Kernel Hilbert Space theory is leveraged to
mathematically analyze the learned QoI map and develop a quantitative
sufficiency test for evaluating the filtered data. An illustrative example is
utilized throughout while the final two examples involve the manufacturing of
shells of revolution to demonstrate various aspects of the presented
frameworks.

###### keywords:

inverse problems , uncertainty quantification , quantity of interest , data-
consistent , push-forward measure , pullback measure

###### MSC:

28A50 , 60-04 , 60-08

##  1 Introduction

Thin structures and shells in particular are common in engineering [1].
Despite their ubiquity, they remain challenging from the point-of-view of
computational mechanics. One of the reasons for this is their sensitivity to
perturbations in materials, kinematic constraints, and shell geometry itself
due to manufacturing imperfections.

In the context of uncertainty analysis for structural design, [2] states that
“[a]ll variables in engineering structures are stochastic to a certain degree”
and that model uncertainty is “due to imperfections and idealizations made in
physical model formulations [$\ldots$] as well as in the choices of
probability distribution types for the representation of uncertainties.” The
authors go on to present a method for uncertainty quantification (UQ) as a
probabilistic quantity where one makes “many” observations and corresponding
predictions of a particular scalar quantity of interest, forming the ratio of
observations to predictions as a measure of discrepancy, and utilizing both
the mean and variance of this ratio to perform the UQ analysis. While we
consider a distinct measure-theoretic UQ framework, the perspective provided
in [2] nonetheless proves prescient in this work. At a high-level, we present
novel developments within a machine learning-enabled data-to-distribution
pipeline to transform (noisy) observed and predicted spatio-temporal data
clouds into samples of (vector-valued) quantities of interest (QoI) that
enables the construction of a data-consistent update to initially assumed
distributions on model parameters via the ratio of (non-parametric) observed
and predicted distributions of the QoI.

In general, any UQ analysis for a computational model of engineered or
physical system governed by principles of mechanics seeks to use principles of
probability theory to identify, classify, and quantify sources of uncertainty
between simulated predictions and observational data, e.g., see [3] for a
comprehensive framework analyzing the impacts of aleatoric (i.e., irreducible)
and epistemic (i.e., reducible) sources of uncertainty on model predictions.
To provide some context, suppose a manufacturing process is implemented to
produce an “ideal” design of some engineered system for which uncertainty is
inherent and irreducible in system inputs, e.g., structural design [2],
material science [4], or thin elastic membranes [5]. Some common sources of
input uncertainty in such systems include, but are not limited to, material
impurities, e.g., see [6, 7, 8], and mechanical tolerances in the machinery
used to construct system components and sub-systems, e.g., see [9, 10, 11].
Complicating matters further is that the system inputs impacted by these
uncertainties often define physical characteristics that have to be
parameterized within the model and are hidden from direct observation.

(a)

(b)

Figure 1: Trommel configuration, hole coverage 35%. (a) Three different
sections are indicated with different graylevels. (b) The sensor locations are
indicated with red colour. The boundaries at $x=\pm\pi$ are clamped and $y=0$
and $y=2\pi$ are periodic.

As an example, consider Young’s modulus, which is a common parameter appearing
in models of elastic materials that is a measure of material stiffness whose
value is often inferred by exerting a known amount of force to the material
and measuring its deformation, e.g., see [12]. Similar to such experimental
settings, the impact of perturbing model parameters is usually observed
indirectly via the simulation of model state data. It follows that the UQ
analysis requires the formulation and solution of a stochastic inverse problem
utilizing discrepancies between (noisy) observational and simulated data. One
example considered in this work involves the construction and monitoring of
trommel screens, which are mechanical screening machines often used to
separate materials of distinct sizes in mineral and solid-waste applications
[13]. Due to harsh operating conditions, the structural health of a population
of screens (defined by the distribution of Young’s modulus) should be
periodically analyzed. In this work, we consider trommel screens with three
sections that possess distinctive regular penetration patterns but constant
35% hole coverage, see Figure 1a for the computational mesh utilized in
simulations with shading indicating the distinct sections. Each section
possess a unique Young’s modulus, which must be inferred from a combination of
observed and simulated data at sensor locations indicated by the red circles
shown in Figure 1b.

###  1.1 Uncertainty Quantification and Stochastic Inverse Problems

The type of stochastic inverse problem that is formulated and the UQ methods
developed for solving such a problem are dictated by the assumptions made
about the type of uncertainty that is to be quantified. The last few decades
have led to significant advancements in UQ research and thousands of papers
published across many computational science journals as well as several
dedicated UQ journals. It is therefore impossible to provide a complete
literature review of UQ and the types of inverse problems at the heart of many
UQ methodologies. Below, we provide a brief literature review that is
intentionally chosen to (i) focus on the two types of uncertainties (epistemic
and aleatoric) considered (often separately) in the literature that are at the
heart of this work, and (ii) to help situate the specific contributions of
this manuscript (summarized in Section 1.2) within the vast UQ literature on
these topics.

In cases where the inverse problem arises from an initial assumption of
epistemic uncertainty in data and parameters, e.g., assuming noisy data are
observed for a single instance of a system associated with true, but unknown,
parameter values, the Bayesian framework is often deployed, e.g., see [14, 15,
16, 17, 18, 19, 20]. The solution to the resulting inverse problem within the
Bayesian framework is known as a posterior, which is a conditional density
defined by the product of a prior density on parameters and a data-likelihood
function that is usually constructed from the differences in simulated and
observed QoI data. The posterior is consequently interpreted as defining the
relative likelihoods that any particular estimate for the parameters could
have produced all of the observed (noisy) data. A commonly used parameter
estimate utilized in the literature is the maximum a posteriori (MAP)
estimate, e.g., see [21, 22, 23]. Under a typical Bayesian setup and
assumptions (such as an additive noise model on the data that follows a given
distribution, usually assumed to be Gaussian), the Bernstein-von Mises theorem
[24] guarantees that the posterior becomes more “spiked” around the true
parameter value and subsequently that the resulting uncertainty in a parameter
estimate, such as the MAP point, is reduced.

In cases where the inverse problem arises from an initial assumption of
aleatoric uncertainty in data and parameters, e.g., assuming natural
variability in system outputs and inputs as discussed above in the context of
an engineered system, the inferrential target is fundamentally a distribution
on parameters that quantifies this intrinsic and irreducible variability.
Alternative Bayesian formulations exist to address such problems such as
hierarchical Bayesian methods [25] that begin with specifying prior
distributions from a parametric family (e.g., normal distributions). Some non-
parametric, hierarchical Bayesian methods for modeling aleatoric uncertainty
also exist, e.g., see [26, 27]. These are an active research area within the
machine/statistical learning community, but they are currently more expensive
and complex to implement than what we utilize in this work.

In this manuscript, we utilize a measure-theoretic data-consistent (DC)
framework. In the DC framework, the inverse problem and its solution are
defined in terms of pullback and push-forward measures [28, 29, 30]. The DC
solution has what is referred to as the data-consistency property in that its
push-forward through the QoI map matches the observed probability measure. The
density-based approximation of the DC solution, as derived in [29] via the
Disintegration Theorem [31], has seen the most development, analysis, and
application in recent years, e.g., see [32, 33, 34, 35, 36, 4, 37]. It is
worth noting that the density form of the solution perhaps first appeared in
[38] where it was derived through heuristic arguments based on logarithmic
pooling and referred to as “Bayesian melding.” The fundamental distinction in
assumptions, form, and properties of the solution from the typical Bayesian
framework led to a distinction of the terminology used in the DC framework in
[39] (which is a follow-up to [29]). In [39] and many of the works that
chronologically follow it, an initial and predicted density are used to
describe the initial quantification of uncertainties on parameters and QoI,
respectively, independent of any observed data. The observed density describes
the quantification of uncertainty for the observed QoI data. An update to the
initial density is then obtained via the product of the initial density with
the ratio of observed to predicted densities evaluated on the outputs of the
QoI map. The updated density serves as the DC solution. Recently, [40]
extended the DC framework to problems that focus on epistemic uncertainties
and parameter estimation while [41, 42] extended the DC framework to problems
that simultaneously involve both aleatoric and epistemic uncertainties. The
work of [42] is most closely related to the current work, which is discussed
more below in the context of the primary contributions of this manuscript.

###  1.2 Mathematical and Algorithmic Contributions

The terms data and information are sometimes used interchangeably, but they
are in fact distinct. The distinction is critical in this work, and we refer
to the commonly accepted definitions that data is a collection of facts
whereas information is the organization and interpretation of those facts. At
a high-level, the interconnected goals of this work are as follows: (i)
develop a robust framework that can learn the most useful quantitative
information present in spatio-temporal data represented as a QoI map, (ii)
utilize the learned QoI map to transform simulated and observed data into the
predicted and observed distributions for DC-based inversion, and (iii)
mathematically represent and analyze this learned QoI map to develop
quantitative diagnostics that evaluate critical aspects of the framework.

When data are strictly temporal in nature, goals (i) and (ii) are achieved
with the Learning Uncertain Quantities (or LUQ111Pronounced as “luck.” for
short) framework introduced in [42], which we briefly summarize here to help
situate the contributions of this work. The basic motivation behind LUQ is
that individual measurement devices can easily produce $\mathcal{O}(10^{4})$
(or more) time steps for which data are collected while the dynamical behavior
is in fact characterized by some “low-dimensional” features. In other words,
organizing and interpreting the data in terms of these low-dimensional
features defines the useful information while the quantitative
characterization of this information defines the QoI map. This addresses goal
(i). Applying this QoI map to the data defines a transformation of samples of
model parameters into numerical samples of the QoI for which (predicted and
observed) distributions are estimable. This addresses goal (ii). LUQ proceeds
in three steps: (Step 1) filtering, (Step 2) clustering and classification,
and (Step 3) feature extraction. The filtering step utilizes predicted and
observed data (both of which may include noise) to approximate the underlying
temporal signal, which is then sampled to produce filtered data, i.e., the
output of this step is a temporal dataset where epistemic uncertainty is
reduced. The clustering and classification step is an optional step. It is
useful in cases where a dynamical system presents strikingly different
qualitative behaviors as model parameters are varied (e.g., due to
bifurcations) or over different windows of time (e.g., due to transient or
equilibrium behavior). In such cases, unsupervised learning approaches are
first applied to classify (i.e., label) the different types of dynamics that
are present in the predicted filtered data. This is followed by the training
and optimal selection of classifiers for this labeled data so that
subsequently filtered observed temporal data are appropriately labeled by
their prevalent dynamics. This leads to the final step of feature extraction
that applies kernel-based principal component analysis (kPCA) [43, 44, 45, 46]
to define, implicitly, the QoI maps that best describe the low-dimensional
nature of the dynamics.

The main contributions of this work are both mathematical and algorithmic in
nature with algorithmic developments encoded within the open-source software
package LUQ [47]. The LUQ package [47] underwent a complete refactoring to
provide (a) an agnostic interface in terms of spatial or temporal data types
and (b) advanced filtering options that apply to higher-dimensional data sets.
The refactored package also allows multiple points-of-entry within the
workflow including offline construction of LUQ objects from pre-computed
predicted data sets that can be subsequently applied to observed data sets as
they are made available. At a conceptual level, the new LUQ framework and
software significantly extend prior work to apply to both spatial and spatio-
temporal data. At a more specific algorithmic level, the filtering step in the
prior LUQ framework only applied to temporal data associated with observing
the evolution of a single state variable in a dynamical system since it
utilized adaptively constructed splines to approximate the underlying temporal
signal. As discussed in [42], the choice of splines in the filtering step is
motivated by a type of universal approximation theory for splines (see [48])
where it is possible to approximate the underlying dynamical response
(assuming a particular regularity or finite number of discontinuities) to
arbitrary pointwise accuracy if both a sufficiently high frequency for
collecting data and number of knots are used. However, the use of splines
prevents any practical application to spatial data in 2- or higher dimensions.
Motivated by the universal approximation theorem for neural nets (NNs) [49],
we provide a deep learning-inspired approach to filtering that utilizes radial
basis functions (RBFs) in the hidden layer [50] which allows for direct
implementation in arbitrary (but finite) dimension.

The contributions listed above address goals (i) and (ii). To address goal
(iii), we provide an explicit derivation and mathematical justification of the
form of learned QoI maps obtained in the feature extraction step. This
utilizes theoretical results from reproducing kernel Hilbert spaces (RKHS)
[51]. This provides significant insight into the mathematical structure of the
QoI map, which is then utilized to develop a quantitative scheme for
numerically determining the sufficiency of filtered data produced in the first
step of the LUQ framework.

###  1.3 Numerical Examples

The main contributions of this work are contextualized in several numerical
examples. First, to help illustrate various DC and LUQ concepts as they are
presented throughout the text, we utilize a simple 2-dimensional wave model
for a fluid where waves are excited by randomly located droplets. The shape of
the droplets and the choice of boundary conditions results in complicated wave
shapes in space and time. This serves to highlight the capabilities of the new
filtering within LUQ as we seek QoI and DC solutions that attempt to recover
the data-generating distribution of droplet locations.

The two main numerical examples are presented in Section 7. These concern
shells of revolution, have several engineering applications, and also serve to
highlight various new and improved capabilities of the LUQ software package.
Shells of revolution are 3D structures, but under assumptions on symmetries
and properties of loading, the problems can be reduced dimensionally to 2D and
even 1D problems. Consider a cylinder with constant radius $R=1$ and thickness
$d=1/100$, say. Many properties of the solutions of the shell problems, for
instance, the resulting displacement fields can be described as functions of
the dimensionless thickness $t=d/R$. Under the assumptions of linear
elasticity, the Young’s modulus $E$, if isotropic over the whole structure,
simply acts as the scaling parameter to the displacement field.

The first example in Section 7 involves the manufacturing of approximately
cylindrical shells subjected to persistent loading in the transverse direction
as would occur, for instance, when such shells are utilized in the
construction of wind turbine towers to be deployed in areas with strong
directional prevailing winds. In this example, it is assumed that
imperfections in materials impacts the Young’s modulus and imperfections in
the machining impacts the shape characteristics of each shell. This example
demonstrates how we can utilize noisy deformation data to update distributions
on Young’s modulus and shape characteristics, which is useful in monitoring
the overall structural integrity of the population of these shells.

The second example of Section 7 involves the construction and monitoring of
trommel screens as previously mentioned. This example demonstrates how we can
utilize spatial data obtained by different experiments at distinct times with
distinct loads to sequentially learn QoI (of potentially distinct dimension)
to learn the distributions of the Young’s modulus associated with each of the
distinctive spatial regions of the trommels.

###  1.4 Organization

The rest of this paper is organized as follows. Section 2 provides a brief
overview of the mathematics behind DC-based inversion including a critical
quantitative diagnostic for verifying a key predictability assumption. Section
3 summarizes the prior work in LUQ for temporal data. Section 4 describes the
extension of the LUQ framework to datasets with spatial attributes that
require filtering. Spatio-temporal datasets are then discussed in Section 5 in
the context of an iterative approach to DCI that utilizes the new spatial-
filtering options within LUQ at each temporal iteration. The mathematical
analysis of learned QoI is presented in Section 6, which is subsequently
utilized to develop a sufficiency diagnostic for filtered data. The main
numerical examples follow in Section 7. Concluding remarks and future research
directions are provided in Section 8. The appendices contain information for
accessing all software and datasets utilized in this work as well as
mathematical information related to the main numerical examples.

##  2 Data-Consistent Inversion: Overview

We begin by introducing the minimal notation and terminology required for
defining the stochastic inverse problem and its solution in the DC framework
considered in this work. Let $\Lambda$ denote an uncertain model parameter
taking values in the parameter space denoted by
$\Lambda\subseteq\mathbb{R}^{p}$. Let
$Q:\Lambda\to\mathcal{Q}\subseteq\mathbb{R}^{q}$ denote a quantities of
interest (QoI) map from the parameter space to the space of observable QoI
data associated with model outputs, which is denoted by
$\mathcal{Q}:=Q\left(\Lambda\right)$. Given an observed probability measure
$\mathbb{P}_{\text{obs}}$ on $\mathcal{Q}$, the goal of Data-Consistent
Inversion (DCI) is to compute a probability measure $\mathbb{P}_{\Lambda}$ on
$\Lambda$ such that the push-forward of this measure through $Q$ matches
$\mathbb{P}_{\text{obs}}$. More formally, if $\mathcal{B}_{\Lambda}$ and
$\mathcal{B}_{\mathcal{Q}}$ denote the Borel $\sigma$-algebras on $\Lambda$
and $\mathcal{Q}$, respectively, and $\mathbb{P}_{\text{obs}}$ is defined on
$\left(\mathcal{Q},\mathcal{B}_{\mathcal{Q}}\right)$, then the DCI solution is
given by a probability measure $\mathbb{P}_{\Lambda}$ on
$\left(\Lambda,\mathcal{B}_{\Lambda}\right)$ such that

$\mathbb{P}_{\Lambda}\left(Q^{-1}\left(E\right)\right)=\mathbb{P}_{\text{obs}}\left(E\right)\>\>\forall
E\in\mathcal{B}_{\mathcal{Q}},$ (1)

where $Q^{-1}$ is the measure-theoretic shorthand for the pre-image of $Q$,
i.e., the map $Q$ is not assumed to be invertible. This introduces an
immediate complication since it is possible (and in fact typical) for many DCI
solutions to exist. To address this issue, we first note it is self-evident
from (1) that DCI solutions are uniquely defined on the induced “contour”
$\sigma$-algebra $\mathcal{C}_{\Lambda}$ defined as
$\mathcal{C}_{\Lambda}=\left\\{Q^{-1}\left(E\right)\>|\>E\in\mathcal{B}_{\mathcal{Q}}\right\\}$.
For any measurable QoI map,
$\mathcal{C}_{\Lambda}\subseteq\mathcal{B}_{\Lambda}$. However, the
containment is usually proper since $\mathcal{B}_{\Lambda}$ often contains
events not present in $\mathcal{C}_{\Lambda}$. Since any event in
$\mathcal{B}_{\Lambda}$ can be written as an intersection with the smallest
contour event containing it222This can be defined in the usual measure-
theoretic way of defining smallest sets in a $\sigma$-algebra with a specific
property by taking the intersection of all such sets from the $\sigma$-algebra
with the specific property., specifying different solutions to the DCI problem
depends upon specifying different probabilities of
$\mathcal{B}_{\Lambda}$-measurable subsets of events in
$\mathcal{C}_{\Lambda}$.

###  2.1 A density-based solution

In [29], a density-based approach to solving DCI is introduced. We summarize
this approach using the terminology from [39] as mentioned in the
introduction. This approach utilizes a disintegration theorem [31] and an
initial density on $\Lambda$, denoted by $\pi_{\text{init}}$, to construct a
specific solution from the space of all DCI solutions as

$\pi_{\text{update}}(\lambda):=\pi_{\text{init}}(\lambda)\frac{\pi_{\text{obs}}(Q(\lambda))}{\pi_{\text{pred}}(Q(\lambda))}.$
(2)

Here, $\pi_{\text{pred}}$ denotes the _predicted density_ defined as the push-
forward of $\pi_{\text{init}}$ through the map $Q$ and $\pi_{\text{obs}}$
denotes the density associated with $\mathbb{P}_{\text{obs}}$. Note that in
more general cases, specifying dominating measures on
$(\Lambda,\mathcal{B}_{\Lambda})$ and
$(\mathcal{Q},\mathcal{B}_{\mathcal{Q}})$ admits the same form of solution but
with the densities interpreted as Radon-Nikodym derivatives. The solution,
$\pi_{\text{update}}$, is referred to as an updated density (or simply as the
update) because it serves to update $\pi_{\text{init}}$ only in directions
informed by the QoI map. This is easily verified by specifying a fixed QoI
value, $\bm{q}\in\mathcal{Q}$, and restricting $\lambda$ to the “contour”
defined by $Q^{-1}(\bm{q})$. For all such $\lambda$, $Q(\lambda)=\bm{q}$ and
the ratio of observed to predicted densities in (2) is a constant. It follows
that the initial likelihoods of the contours are updated while the initial
conditional densities along the contours remain unchanged in the update. To
emphasize the role of the ratio of observed to predicted densities in updating
the initial, we rewrite (2) as

$\pi_{\text{update}}(\lambda)=r(Q(\lambda))\,\pi_{\text{init}}(\lambda),\
r(Q(\lambda)):=\frac{\pi_{\text{obs}}(Q(\lambda))}{\pi_{\text{pred}}(Q(\lambda))}.$
(3)

As shown in [29], $\pi_{\text{update}}$ is unique (up to the choice of
$\pi_{\text{init}}$), and is also stable in the total-variation metric with
respect to perturbations in $\pi_{\text{init}}$, $\pi_{\text{obs}}$, or
$\pi_{\text{pred}}$. The last stability result is critical since even if
$\pi_{\text{init}}$ and $\pi_{\text{obs}}$ are specified exactly, we must
often estimate $\pi_{\text{pred}}$ in constructing an estimate of
$\pi_{\text{update}}$. In cases where the QoI map is itself approximated as
part of the process, [39, 36] provides convergence analysis for the sequence
of approximate updates.

###  2.2 Numerical approximation, predictability assumption, and a
quantitative diagnostic

It is often the case that $Q$ is nonlinear or there is a change in dimensions
between $\Lambda$ and $\mathcal{Q}$. In such cases, we expect
$\pi_{\text{pred}}$ and thus $\pi_{\text{update}}$ to both be non-parametric
even if $\pi_{\text{init}}$ and $\pi_{\text{obs}}$ are given exactly from some
well-studied parametric families of distributions (e.g., as Gaussian
distributions). In this work, for the sake of both simplicity and
reproducibility of results, we estimate $\pi_{\text{pred}}$ with standard
kernel density estimation (KDE) [52] techniques. This requires us to address a
key assumption involving $\pi_{\text{pred}}$ that allows for practical
approximation and generation of random samples with standard techniques333By
standard techniques, we refer to “off the shelf” approaches such as kernel
density estimation and rejection sampling as utilized in [29, 39].. In [29],
this is referred to as the predictability assumption given by

######  Assumption 1.

There exists a constant $C>0$ such that $\pi_{\text{obs}}(\bm{q})\leq
C\pi_{\text{pred}}(\bm{q})\>\>\forall\bm{q}\in\mathcal{Q}$.

Under this predictability assumption, $\pi_{\text{update}}$ defined in (2) or
(3) is a density. It follows that

$\mathbb{E}_{\text{init}}(r(Q(\lambda))=\int_{\Lambda}r(Q(\lambda))\,d\mathbb{P}_{\text{init}}=\int_{\Lambda}\pi_{\text{init}}(\lambda)r(Q(\lambda))\,d\mu_{\Lambda}=\int_{\Lambda}\pi_{\text{update}}(\lambda)\,d\mu_{\Lambda}=1,$
(4)

where $\mathbb{P}_{\text{init}}$ is the initial probability measure associated
with the density $\pi_{\text{init}}$ and $\mu_{\Lambda}$ is the Lebesgue
measure444In the case of more general Radon-Nikodym derivatives,
$\mu_{\Lambda}$ is a specified dominating measure on $\Lambda$. on $\Lambda$.
We turn this into a quantitative diagnostic for verifying the predictability
assumption as discussed in [29]. Specifically, given a set of independent
identically distributed (iid) parameter samples drawn from the initial
distribution, the corresponding sample average of the ratio $r(Q(\lambda))$
should be approximately unity. To make this a computationally cheap diagnostic
to compute in this work, we re-use the same parameter samples involved in
estimating $\pi_{\text{pred}}$. More precisely, we estimate
$\pi_{\text{pred}}$ from a KDE applied to a set of QoI samples obtained by
evaluating the QoI map on an iid set of parameters generated from the initial
distribution. The quantitative diagnostic defined by the sample average of
$r(Q(\lambda))$ is then estimated without any further model simulations. It is
the experience of the authors that values between $[0.9,1.1]$ are expected
when the predictability assumption is satisfied and sufficiently accurate
estimates of the predicted and observed densities are utilized.

###  2.3 Illustrative Example: Randomly Generated Waves (Part I)

As an illustrative example, we consider the random generation of 2-dimensional
water waves in a square domain where the goal is to determine the distribution
of initial droplets based on observations of an a priori defined QoI map. For
the mathematical model, we consider the prototypical wave equation
$u_{tt}=\nabla^{2}u$, where $u=u(x,y,t)$ denotes the heights of the waves on a
physical (dimensionless) domain defined by $(0,5)^{2}\subset\mathbb{R}^{2}$
with homogeneous Dirichlet boundary conditions. The model is solved with a
standard centered finite difference scheme on a $101\times 101$ regular
uniformly-spaced mesh ($99\times 99$ mesh when excluding the boundary) with
$0.005$ sized time-steps555See A on instructions for obtaining the code
utilized to generate the results in this example.. The initial condition
models a droplet of water at some uncertain location $(a,b)\in(0,5)^{2}$ given
by

$u(x,y,0)=0.2\exp\left(-10\left((x-a)^{2}+(y-b)^{2}\right)\right).$

Figure 2: Contour plot of the QoI map $Q$ over $[0,5]^{2}$ for illustrative
example (part I) with the specified QoI map
$Q(a,b)=u\left(4.0,1.0,2.5;a,b\right)$ which represents the wave height at
spatial location $(4.0,1.0)$ 2.5 time units after initial water droplet at
$(a,b)$ occurs.

In the more abstract notation given above,
$\lambda=(\lambda_{1},\lambda_{2})=(a,b)$ defines the parameters for this
problem with $\Lambda=(0,5)^{2}$ denoting all physically possible values for
these parameters. To emphasize the role of $(a,b)$ on the solution, we rewrite
it as $u(x,y,t;a,b)$, and we define the QoI map as
$Q\left(a,b\right)=u\left(4.0,1.0,2.5;a,b\right)$. In other words, we assume
that a single measurement of a wave height is obtained at $t=2.5$ at location
$(x,y)=(4.0,1.0)$. Note first that $Q$ maps from $\mathbb{R}^{2}$ into
$\mathbb{R}$ so that the contour structure of $Q$ significantly influences the
DCI solution. See Figure 2 for a contour plot of $Q$ over $\Lambda$.

To construct the observed distribution, we utilize a _data-generating (DG)
distribution_ defined as the exact distribution of $(a,b)$ values that leads
to the observations. Here, $a$ and $b$ are assumed independent with the DG
distribution for $a$ (the horizontal location of a droplet) taken as a
$\text{Beta}(2,5)$ distribution that is scaled and shifted to the interval
$[1,2]$ while the DG distribution for $b$ (the vertical location of a droplet)
is taken as a $N(2.5,0.5)$ distribution. We randomly generate 200 iid $(a,b)$
values from the joint DG distribution and numerically solve the model to
simulate the associated QoI values defining the observations.

To quantify the differences and similarities between densities, we make use of
the total variation (TV) metric throughout this work. In particular, we use
the standard form of the TV metric for densities defined as one-half the
$L^{1}$-distance between densities [53]. This produces values between $[0,1]$
with a value of $1$ obtained when the two densities have non-overlapping
support. To provide context on typical TV values obtained for KDE estimates of
the DG distribution, we generate ${1}\mathrm{E}{3}$ sets of 200 iid $(a,b)$
samples, compute the corresponding KDEs for each set of 200 samples, and then
compute the TV distance between these KDEs and the exact DG distribution.
These results are summarized in the histograms of Figure 3, which illustrates
that a “typical” TV distance for the joint KDE estimate is between
approximately $0.08$ and $0.14$ from the exact joint density while the
marginal estimates appear concentrated around a value of $0.05$. The vertical
lines in these plots correspond to a particular set of 200 iid $(a,b)$ samples
drawn from the DG distribution that we utilize throughout this work. This
particular set of samples from the DG distribution produces KDE estimates of
the joint, $a$-marginal, and $b$-marginal densities that are TV distances from
their exact counterparts of approximately $0.12$, $0.03$, and $0.05$,
respectively.

Figure 3: Histogram of TV values between KDEs and exact joint DG densities
using 1000 random sets of 200 random iid samples. On the left is the resulting
histogram. The middle and right plots represent the resulting histograms for
the horizontal and vertical marginals, respectively. The vertical lines in
each plot represent the TV metric between the KDEs of 200 i.i.d. samples used
for the illustrative example throughout the paper.

To construct the DCI solution, we take the uniform distribution on $[0,5]^{2}$
as the initial distribution for $(a,b)$, randomly generate ${1}\mathrm{E}{4}$
iid samples from this initial distribution, simulate the associated QoI
values, and then apply a standard KDE to estimate the predicted density. We
then calculate the associated $r(Q(\lambda))$ values for these predicted QoI
values based on the estimated observed and predicted densities. The sample
mean of these $r$-values (i.e., the quantitative diagnostic as defined by (4))
is, to three-decimals, 0.961, which indicates the predictability assumption is
verified and that the predicted and observed densities are reasonably
approximated.

We use a mixture of qualitative and quantitative methods for analyzing and
comparing the various densities. First, we utilize the $r$-values to visualize
the updated density and its marginals in Figure 4. The left plot shows the
${1}\mathrm{E}{4}$ initial samples colored by their associated $r$-values,
which illustrates the impact of the contour structure of the QoI map as well
as the uniform initial density. Clearly, the updated density is no longer
uniform although this initial distribution does essentially assign uniform
likelihoods across the contours of the QoI map. There is a clear concentration
of probability in a contour event containing most of the samples generated
from the DG distribution. This is reflected in the TV metric computations as
well where the TV metric between the exact uniform initial joint density and
the KDE of the DG joint density is 0.808 while the TV metric is 0.718 between
the KDEs of the updated and DG joint densities. The updated marginals are
computed as weighted KDE estimates with the $r$-values serving as the weights
and are shown in the middle and right plots (for $a$ and $b$, respectively).
These are of limited utility and do not appear significantly updated from the
uniform initial distributions due to the contour structure of the QoI map.
This is reflected in the TV metric computations where we see the TV metrics
between the exact uniform marginals and KDE of the DG marginals are 0.706 and
0.544 for $a$ and $b$ (resp.) while the TV metrics between the KDEs of the
updated marginals and DG marginals are 0.651 and 0.543 for $a$ and $b$
(resp.).

(a) Marginal for $a$.

(b) Marginal for $b$.

Figure 4: (Left) Color plot of the DCI solution using $10E3$ uniform samples
colored by the updated density for the illustrative example part I with data-
generating samples overlaid in black. (Middle and Right) Updated density
marginals of $a$ and $b$ for the illustrative wave equation example part I
along with the original data-generating marginals. The TV metric between the
initial predicted distribution and the data-generating distribution is 0.808
while the TV between the updated solution and the data-generating distribution
is 0.718 showing minor improvements in learning the data-generating
distribution.

The above qualitative and quantitative discrepancies between the updated and
DG distributions highlights the deficiency in both a priori specifying and
utilizing a single QoI. Unfortunately, given the complexity in the contour
structure of any particular measurement over the parameter space (e.g., as
shown in Figure 2), it is not at all obvious how one would a priori specify
multiple QoI that lead to an updated density that is a reasonable
approximation of the true data-generating distribution. This issue is
addressed in the next section.

##  3 Learning Uncertain Quantities (LUQ) for Temporal Data

Let $\mathcal{M}(u;\lambda)=0$ define a model that admits a solution
$u(\lambda)$ of state variables dependent upon model parameters $\lambda$
representing various uncertain inputs of the model taking values in parameter
space $\Lambda\subset\mathbb{R}^{n}$. Note that in this work $u$ also depends
on spatial or temporal variables, but we omit explicit dependence on those
variables here both for simplicity in notation and to emphasize the dependence
on the parameters. While the solution operator of the model defines a mapping
from parameters, $\lambda$, to state variables, $u$, in practice it is
impossible to observe the full solution across all of space or time. Instead,
we collect a finite amount of spatial or temporal data on the solution $u$
that is then analyzed and transformed to the quantitative information useful
for defining and solving a stochastic inverse problem. This is the goal of the
Learning Uncertain Quantities (LUQ) framework.

###  3.1 The LUQ framework

Originally introduced in [42], the LUQ framework defines a three-step data-to-
distribution pipeline to transform noisy temporal data clouds for dynamical
systems into distributions of learned QoI for the purposes of DCI. An open-
source Python module LUQ is also available [47]. We summarize the three steps
below and refer the interested reader to [42] for more details.

####  3.1.1 Step 1: Filter Noisy Data

The filtering step occurs in two stages. The first step is to approximate the
underlying dynamical signal responsible for the noisy data. Then, from this
approximate signal, “noise-free” (i.e., filtered) data are produced.

In [42], splines are utilized for filtering temporal measurements based on a
universal approximation theorem and their ubiquity in inverse problems (see
[48] for a thorough treatise of the subject). For each set of temporal
measurements obtained at a fixed point in space, piecewise linear splines are
optimized to approximate the underlying dynamical response associated with the
noisy temporal data. While the first and last knot are fixed to the beginning
and end of the time-interval, the interior knot placements are optimized along
with the values of the spline at each knot. The number of knots used is
incrementally increased from a user-specified minimum until either the
relative improvement in the sum square error falls below a threshold or a
maximum number of knots, specified by the user, is reached. From the resulting
piecewise linear splines, a specified number of usually equispaced data points
are then calculated defining the filtered time series data.

####  3.1.2 Step 2: Learn and Classify Dynamics

The second step of LUQ, which is optional, is to perform unsupervised learning
in order to label (i.e., define clusters) followed by supervised learning to
classify the various types of dynamics present in the filtered data. This is
an important step when different types of signals are present in the observed
dynamics (e.g., due to bifurcations), which may be quantified by entirely
distinct types and numbers of QoI.

The LUQ module fully supports some of the most popular clustering schemes such
as $k$-means [54], Gaussian mixture modeling [55], and spectral clustering
routines [56, 57]. All such methods begin with a specification of the number
of labels (i.e., clusters) to learn from the predicted data. While these
schemes often implicitly contain a classifier model that can be applied to
other data (e.g., the observed data), Support Vector Machines (SVMs) are a
more robust class of supervised learning methods that can be used for, among
other things, classification that has advantages in terms of operations in
high dimensions and memory efficiency. The LUQ module leverages the widely
used library `LIBSVM` [58] for SVMs to construct classifiers.

Denote by $K$ the number of labels to learn from the predicted data that are
then utilized to construct an SVM for classification. This implicitly defines
a partitioning of $\Lambda$ into $K$ subsets666In general, the $K$ subsets may
be both disconnected and non-convex., which we denote by $\Lambda_{k}$ for
$k\in\left\\{1,2,\ldots,K\right\\}$. It is important to note that it is not
technically necessary to explicitly identify $\Lambda_{k}$ for any $k$, but it
is computationally trivial to identify which initial subset of samples drawn
from the initial distribution belong to each $\Lambda_{k}$ by simply sorting
the labels of the associated sample set of prediction data. We make use of the
notation $\Lambda_{k}$ to describe the structure of the DCI solution at the
end of Step 3.

####  3.1.3 Step 3: Learn QoI via Feature Extraction and Compute the DCI
solution

The LUQ module makes use of kernel-based principal component analysis (kPCA)
[43, 44, 45, 46], which is one of the most popular approaches to performing
feature extraction. Specifically, kPCA is performed with multiple kernels on
each classification of predicted filtered time series data to learn the most
beneficial QoIs for each type of dynamics. We provide an in-depth mathematical
derivation of the QoI maps in terms of the chosen kernel in Section 6, which
is subsequently utilized in a proposed quantitative method for evaluating the
sufficiency of the filtered data for learning the QoI.

While the LUQ module technically ends after feature extraction is performed to
learn the QoI, it is now possible to begin computing the DCI solution based on
this learned QoI. To that end, let $Q^{(k)}(\lambda)$ denote the learned QoI
map defined on $\Lambda_{k}$ and $\pi_{\text{pred}}^{(k)}$ denote a predicted
density associated with a re-normalization of $\pi_{\text{init}}$ restricted
to $\Lambda_{k}$. Similarly, let $\pi_{\text{obs}}^{(k)}$ denote the observed
density associated with the $k$th label of data and $w_{k}$ denote the
observed weights associated with this label defined by the ratio of observed
data with this label to all observed data of any label, i.e., $0\leq w_{k}\leq
1$, $\sum_{k=1}^{K}w_{k}=1$, and

$\sum_{k=1}^{K}w_{k}\pi_{\text{obs}}^{(k)}(Q^{(k)}(\lambda))\mathbb{I}_{\Lambda_{k}}(\lambda)$
(5)

defines the observed density over $\Lambda$ where $\mathbb{I}_{\Lambda_{k}}$
denote the indicator function over $\Lambda_{k}$. With this notation, the DCI
solution is given by

$\pi_{\text{update}}(\lambda)=\sum_{k=1}^{K}w_{k}\left(\pi_{\text{init}}(\lambda)\frac{\pi_{\text{obs}}^{(k)}(Q^{(k))}(\lambda)}{\pi_{\text{pred}}^{(k)}(Q^{(k)}(\lambda))}\right)\mathbb{I}_{\Lambda_{k}}(\lambda).$
(6)

While using the same number of principal components (PCs) as the number of
desired QoIs will typically explain the majority of variance in the data, the
resulting QoIs may be heavily correlated leading to a poor construction of the
updated density. To decide whether a set of desired QoI are to be used for
DCI, we check the diagnostic discussed in 2.2 for each label (i.e., we compute
the diagnostic for each ratio of $\pi_{\text{obs}}^{(k)}$ to
$\pi_{\text{pred}}^{(k)}$ shown in (6)) and systematically reduce the number
of QoI extracted for each type of dynamic if the diagnostic deviates too far
from unity.

###  3.2 Illustrative Example: Randomly Generated Waves (Part II)

Returning to the illustrative example from Section 2.3 with the same DG and
initial distributions as well as the same random samples from these
distributions, we now consider temporal data recorded every 0.5 time units up
to $t=7.0$ at the spatial location $(4.0,1.0)$ producing 14 measurements of
wave heights for every experiment777See A on instructions for obtaining the
code utilized to generate the results in this example.. In the interest of
simplicity and to illustrate limitations in the analysis of only temporal data
in the previous LUQ framework, we avoid any introduction of noisy data. This
allows us to skip Step 1. Moreover, clustering is not found to be useful in
this case, so we also skip Step 2. Ideally, we seek a QoI map of the same
dimension as the parameter space to attempt to recover the DG distribution.
Hence, the goal is to learn a 2-dimensional QoI map for the DCI solution. To
this end, we apply Step 3 and perform kernel PCA on the 14-dimensional data to
learn the 2-dimensional QoI map for DCI. The sample average of the resulting
$r$-values (i.e., the quantitative diagnostic as defined by (4)) for this QoI
map is approximately 1.031 indicating that the predictability assumption is
satisfied and the observed and predicted densities are well approximated.

The QoI map learned is still producing a coarse contour structure for the
updated density as evidenced by Figure 5. While there appears to be improved
concentration of probability near the samples drawn from the DG distribution,
there is no improvement in terms of the TV metric. Here, the TV metric between
the updated joint density and the exact DG density is approximately 0.731. By
comparison, the TV metric between the updated joint density and the exact DG
density is 0.718 using the QoI from Section 2.3. It is self-evident from
Figure 5 that the marginals are not improved from those seen in Figure 4. When
using strictly temporal data from different spatial locations, similar results
occur. This highlights the need to extend the LUQ framework to admit spatio-
temporal data rather than simply temporal data.

(a) Joint density.

(b) Marginal for $a$.

(c) Marginal for $b$.

Figure 5: (Left) Color plot of the DCI solution using $10E3$ uniform samples
colored by the updated density for the illustrative example part II with data-
generating samples overlaid in black. (Middle and Right) Updated density
marginals of $a$ and $b$ for the illustrative wave equation example part II
using purely temporal data along with the original data-generating marginals.
The TV metric between the updated solution and the data-generating
distribution is 0.731. For part I of this example, the resulting TV metric was
0.718. While there is an increase in TV metrics, the increase is within
expected error tolerances with the quadrature used to compute the TV metric
integral suggesting that there are minor changes in the results from parts I
and II.

##  4 Extending LUQ for Spatial and Spatio-Temporal Data

Utilizing linear splines or even higher order splines to filter spatial or
spatio-temporal data introduces technical difficulties due to the increased
dimension that the spatial domain can introduce as well as the complexity in
defining a mesh of potential spline knots. At the same time, the filtering
step is essential in LUQ since we cannot reasonable expect a kPCA to learn the
QoI maps that best explain perturbations in parameters that are useful for DCI
from data with significant noise levels as this results in much of the
variation in the data being explained by the noise itself. We therefore extend
both the LUQ framework and module LUQ to learn a mixture of polynomial and
radial basis functions from noisy data for the purpose of filtering data in
other dimensions. This filtering of data utilizes a variety of user-specified
functions and is described in general terms through the use of a Neural
Network (NN) [59].

###  4.1 A deep-learning based perspective for filtering

Much like with the use of splines for filtering temporal data, the motivation
for using NNs comes from a Universal Approximation Theorem (UAT) for NNs that
states that there exists a NN with a single hidden layer that can approximate
a continuous function defined on a compact subset of $\mathbb{R}^{m}$ to
arbitrary accuracy in the sup-norm metric [49]. Assuming a scalar valued
function is to be approximated on an $m$-dimensional input space using a NN
with a single hidden layer of width $n$, then such a NN’s output, denoted by
$N(x)$ is given by a function of the form

$N(x)=\mathcal{A}\left(\alpha_{0}+\sum_{j=1}^{n}\alpha_{j}\phi_{j}\left(w_{j,0}+\sum_{i=1}^{m}w_{j,i}x_{i}\right)\right).$
(7)

Here, $\mathcal{A}$ is the activation function of the neuron on the output
layer while $\phi_{j}$ is the activation function associated with the $j$th
neuron in the hidden layer. The $\alpha_{j}$, for $0\leq j\leq n$, define the
hyperparameters between the hidden and output layer, while the $w_{j,i}$, for
$0\leq i\leq m$ for each $0\leq j\leq n$, define the hyperparameters between
the input and hidden layer. See Figure 6 for a schematic. All of these
hyperparameters are learned via backpropagation [59]. When $\mathcal{A}$ is
chosen to be the identity, as we initially utilize for simplicity, the
equation for $N(x)$ simplifies to

$N(x)=\alpha_{0}+\sum_{j=1}^{n}\alpha_{j}\phi_{j}\left(w_{j,0}+\sum_{i=1}^{m}w_{j,i}x_{i}\right).$
(8)

Figure 6: Neural network schematic illustrating location of hyperparameters
$w_{j,i}$ and $\alpha_{j}$ in (7) and (8).

While the UAT establishes the existence of a NN achieving a specified level of
accuracy (assuming significant data are provided to train the NN), it does not
provide any insight as to the necessary width (i.e., the number of neurons)
required to achieve such accuracy. Some insight is still generally required to
setup a range of potentially useful NNs. Below, we describe how such insight
may be obtained in the context of utilizing NNs in the filtering step for LUQ.

First, consider a NN with $n$ neurons where each neuron is activated via
rectified linear unit (ReLU) activation functions [60]. In the case of
temporal data, the input to the $j$th neuron in this layer is defined by
$w_{j,1}t+w_{j,0}$. The activation function of the $j$th neuron is then given
by $\phi_{j}(w_{j,1}t+w_{j,0})=\max\left\\{0,w_{j,1}t+w_{j,0}\right\\}$. The
output of the NN (assuming the identity activation output layer) is then
equivalent to utilizing linear splines where the values of
$\\{w_{j,0},w_{j,1}\\}_{j=1}^{n}$ are associated with node placement while the
values of $\alpha_{j}$ for $0\leq j\leq n$ are associated with the value of
the spline through the nodes. In this case, identifying a range of neurons to
use in the hidden layers of potential NNs is equivalent to choosing a range of
knot values to use in linear spline approximations. A pre-analysis of noisy
data (e.g., using visualization or other techniques) or domain-specific
knowledge of solutions to models and their general features, then allows the
user to specify a range of neurons necessary for approximating the underlying
temporal signal.

We adapt this conceptual perspective to the use of radial basis functions
(RBFs). Specifically, we consider the ubiquitous Gaussian-type functions of
the form $\phi_{j}(x)=e^{-\lVert x-r_{j}\rVert^{2}/\sigma_{j}}$ where $x$ is
the $m$-dimensional location variable. While this can be rewritten using the
$w_{j,i}$ notation for hyperparameters, it is conceptually simpler to explain
the interpretation of the hyperparameters in the notation of $r_{j}$ and
$\sigma_{j}$. The hyperparameters $r_{j}$ are interpreted as the
$m$-dimensional vectors defining the centers of the activation functions and
the $\sigma_{j}$ denotes a scaling factor associated with the width of the
activation. If $\sigma_{j}$ is “small,” then the associated $j$th neuron is
generally associated with a localized feature present in the function near the
$r_{j}$ location whereas if $\sigma_{j}$ is “large,” then the associated $j$th
neuron is generally associated with a non-localized trend (such as an average
value) present in the function “around” the $r_{j}$ location. This perspective
and interpretation is valuable in user-provided guidance in specifying the
depth of the NN required to approximate data with a NN instead of linear
splines. Moreover, for moderate levels of complexity in the function response
(i.e., for functions with a few localized features such as those resulting
from many smoothly modeled phenomena), we can use $n$ values that permit
direct optimization of the loss function associated with training a NN, which
is often equivalent to minimizing a sum square error. We discuss this more
below.

###  4.2 RBF-based filtering in LUQ

One of the redesigns of the LUQ package produced as part of this work is the
inclusion of an RBF fitting algorithm using a trust-region reflective least
squares algorithm [61] available within the scipy Python package for the
filtering step. This allowed us to circumnavigate the need to tune training
parameters in NNs and to speed up the optimization by not using the more
robust optimization algorithms standard in large-scale machine-learning
libraries. We note, however, that the RBF fitting technique we deploy is
mathematically equivalent to fitting a specific NN.

Before fitting the data, we perform some standard pre-processing of the data
by first centering the data according to its mean. The RBF fitting algorithm
we implement allows for centering the data over a linear or quadratic trend if
desired. With the data centered, we then fit a user-specified number $n$ of
Gaussian-type RBFs and an optional polynomial $p(x)$ of user-specified degree
as

$\sum_{i=1}^{n}w_{i}e^{-\left\lVert x-r_{i}\right\rVert^{2}/\sigma_{i}}+p(x)$

where for each $i=1,\dots,n$, the hyperparameters that are fitted are the
weights $w_{i}$, centers $r_{i}$, scales $\sigma_{i}$, and any coefficients
that define the optionally specified polynomial $p$. It is worth noting that
the options to center the data over a linear or quadratic trend as well as the
option to include a polynomial to learn as part of the RBF fitting serves to
generally reduce the number $n$ of RBFs required to approximate the underlying
function since more global trends in the data are approximated by these
polynomial objects, which are typically more well-suited to such tasks than
RBFs.

###  4.3 Filtered data from disparate datasets

It is worth noting an additional feature of the LUQ framework that is fully
incorporated within the new LUQ module developed as part of this work. The
predicted data is often simulated and therefore straightforward to define on a
spatial grid utilized within the simulation model. The observed data may be
more irregularly spaced and not aligned to any grid. The LUQ module allows for
such disparities between datasets by allowing the user to specify precisely
the coordinates where (predicted or observed) data are recorded and
identifying the unifying set of filtered coordinates for which the filtered
data are obtained. This flexibility allows for functionality of the framework
even in situations where the predicted and observed datasets are of different
dimensions. Iterating through the filtering method within LUQ for each
individual sample associated with a specific predicted or observed stream of
data further allows the user to specify data coordinates for each sample,
which is useful in situations where data coordinates or number of data differ
from one iteration to the next (e.g., due to adaptive mesh
refinement/coarsening in the predicted case or accounting for sensor outages
in an observation network). The filtering step is particularly crucial in such
situations to ensure a common dimensionality of the data vectors analyzed in
either the clustering or feature extraction steps.

###  4.4 Illustrative Example: Randomly Generated Waves (Part III)

Figure 7: Random sample taken from the 200 i.i.d. observed samples consisting
of exact wave height measurements over the spatial domain $[0,5]$ along with
the fitted surface computed by fitting a number of Gaussians up to 7 total
until a relative tolerance is reached from the mean-centered noisy wave height
data. The fitting learns the locations, scales, and weights of each Gaussian.
The sample used in this plot was fitted using 6 Gaussians.

Revisiting the illustrative example, we now assume we have sensors at 81
spatial locations given by $\left(0.5i,0.5j\right)$ for $i,j=1,\dots,9$ that
record the wave height at $t=2.5$. For the sake of illustrating the filtering
step for spatial data, we add “white” noise simulated by a normal distribution
centered at $0.0$ with standard deviation of ${2.5}\mathrm{E}{-3}$, which is a
significant amount of noise given that the wave height measurements are on the
order of ${1}\mathrm{E}{-2}$). We specify a range of 1-7 RBFs to approximate
the wave’s deviation from its approximate mean value computed from the noisy
data. Beginning with a single RBF, the number is incrementally refined until
either the relative error between fittings reaches an error tolerance of 0.01
or the a priori specified maximum number of 7 RBFs is reached. A
representative plot is shown in Figure 7 for an RBF-based surface that is
learned for a single randomly generated sample of observation data. Note that
the predicted data exist at each time step on the computational mesh of
$99\times 99$ spatial coordinates. However, we choose to evaluate the RBF-
based surface on a different mesh defining the filtered data coordinates given
by the mesh consisting of 2401 spatial locations given by
$\left(0.1i,0.1j\right)$ for $i,j=1,\dots,49$. This is not an arbitrary choice
for the filtered coordinates as we show in Section 6 where we provide both a
mathematical and statistical justification for choosing these filtered
coordinates. For now, we simply note that for each sample from the DG
distribution there are 81 noisy data points, for each sample from the initial
distribution there are 9801 data points, but for each sample, regardless of
the distribution it comes from, there are 2401 filtered data points.

The learned QoI map is obtained via kPCA on the predicted spatial data, and
the remaining steps are analogous to the prior illustrative example parts. The
resulting updated density along with the marginals of the updated density
compared to the data-generating densities are shown in Figure 8. The resulting
TV between the updated density and the exact DG density is about 0.390 which
is a significant improvement from the values of 0.718 and 0.731 obtained in
Part I using a single measurement at the spatial location $(4.0,1.0)$ and Part
II using temporal data at the same spatial location, respectively. These
qualitative and quantitative results demonstrate a significant improvement
over the results of Parts I and II of this illustrative example in terms of
the updated density providing a better approximation to the DG density, though
the value obtained here is still not within the expected range of TV values
between KDEs of the DG density and the exact DG density which is within about
0.01 and 0.13.

(a) Joint density.

(b) Marginal for $a$.

(c) Marginal for $b$.

Figure 8: (Left) Color plot of the DCI solution using $10E3$ uniform samples
colored by the updated density for the illustrative example part III with
data-generating samples overlaid in black. (Middle and Right) Updated density
marginals of $a$ and $b$ for the illustrative wave equation example part III
using noisy spatial data along with the original data-generating marginals.
The TV metric between the updated solution and the data-generating
distribution is about 0.390 showing a major improvement from parts I and II of
this example that produced TV metrics of about 0.718 and 0.731, respectively.
The TV metric between the KDE of the data-generating distribution and the
exact data-generating distribution is about 0.119 suggesting that the DCI
solution has room for improvement which is achieved in part IV by utilizing a
combination of spatial and temporal data.

##  5 Temporal Iterations of Spatially Filtered Data to Update Distributions

###  5.1 Iteratively obtained QoI and DCI

Once the spatial data are filtered at a particular time as discussed in the
previous section, we learn the QoI maps to use for DCI. If the data are purely
temporal or spatial, then QoI maps are learned and utilized as described in
Steps 2 and 3 of Section 3. For spatio-temporal data, we propose an iterative
approach to calculating a DCI solution. At a conceptual level, after the first
updated density is computed in the first iteration, the initial density at
each subsequent iteration is defined as the updated density of the prior
iteration.

Mathematically, given the $i$th iteration that produces an updated density
denoted by $\pi_{\text{update}}^{(i)}$, the $i+1$ iteration for $i\geq 1$
starts with an initial density defined by
$\pi_{\text{init}}^{(i+1)}:=\pi_{\text{update}}^{(i)}$. We then attempt to
learn $p=\text{dim}\left(\Lambda\right)$ number of QoI and denote the learned
map at this iteration as $Q^{(i+1)}$ (not to be confused by the notation used
to denote any potential part of the QoI map associated with a cluster), and
compute the weights $r^{(i+1)}(Q^{(i+1)}(\lambda))$ at this iteration using
the ratio of observed to predicted densities constructed and evaluated from
this learned QoI map as described in Section 2. This results in the $i+1$
update defined as
$\pi_{\text{update}}^{(i+1)}(\lambda):=r^{(i+1)}(Q^{(i+1)}(\lambda))\pi_{\text{init}}^{(i+1)}(\lambda)$.
Denoting $\pi_{\text{init}}$ as the original initial density for the first
iteration, it follows that

$\pi_{\text{update}}^{(i+1)}(\lambda)=\left(\prod_{j\leq
i+1}r^{(j)}(Q^{(j)}(\lambda))\right)\pi_{\text{init}}(\lambda).$ (9)

This is summarized in Algorithm 1 below where we emphasize how the DCI
diagnostic is utilized at each iteration to systematically determine the
number of components of each $Q^{(i+1)}$ map to utilize. Note the possibility
that no components are utilized. This is a feature of the algorithm that
avoids two undesirable situations. One situation involves a violation of the
predictability assumption, which can occur, for example, due to an inability
to sufficiently filter noise at a particular time step. Another situation
involves data that are insensitive to parameter value perturbations, which can
result, for example, from dynamics converging to steady-state or equilibrium
values resulting in data for which it is fundamentally not possible to
quantitatively learn a useful QoI map for DCI.

Algorithm 1 $i+1$ Iteration of Iterative DCI

1:apply kPCA to spatial data at current time step to learn new $d$-dimensional
QoI map $Q^{(i+1)}(\lambda)$

2:evaluate $r=\prod_{j\leq i}r^{(j)}(Q^{(j)}(\lambda))$ at initial parameter
samples

3:estimate $\pi_{\text{obs}}$ from learned QoI samples and $\pi_{\text{pred}}$
as a weighted KDE of learned QoI samples weighted by $r$

4:evaluate
$r^{(i+1)}(Q^{(i+1)}(\lambda))=\pi_{\text{obs}}(Q^{(i+1)}(\lambda))/\pi_{\text{pred}}(Q^{(i+1)}(\lambda))$
at learned predicted QoI samples

5:if $|1-\text{mean}(r^{(i+1)}r)|>tol$ then

6: restart using $(d\leftarrow d-1)$ or discard time step if $d=1$

While the DCI diagnostic in this iterative approach informs us when it is
possible to utilize data from a particular time step to update the
distribution for the next iteration, it does not state whether an iterative
update is computationally worthwhile. The impact of a potential $i+1$
iteration to improve the overall solution requires the contour structure of
$Q^{(i+1)}$ to differ from that of $Q^{(i)}$. This naturally leads to
questions involving optimal experimental design (OED) for learning QoI maps.
While defining any formal OED metrics based on quantification of geometric
structures is outside the scope of this work, it is the subject of ongoing
research in DCI. Here, we note that the iterative approach does permit, a
simple statistic that quantifies, in a sense, the degree of change in a
proposed updated density in terms of the standard deviation of the $r$-values.
It is straightforward to verify that the variance of the $r$-values for a
given QoI map $Q$ with respect to $\mathbb{P}_{\text{init}}$ is given by

$\displaystyle\mathbb{V}_{\text{init}}\left(r(Q(\lambda))\right):=\mathbb{E}_{\text{init}}\left(r^{2}(Q(\lambda))\right)-\mathbb{E}_{\text{init}}^{2}\left(r(Q(\lambda))\right)=\mathbb{E}_{\text{update}}\left(r(Q(\lambda))\right)-1$

where $\mathbb{E}_{\text{init}}$ and $\mathbb{E}_{\text{update}}$ are the
expected value operators with respect to $\mathbb{P}_{\text{init}}$ and
$\mathbb{P}_{\text{update}}$, respectively. We note that
$\mathbb{V}_{\text{init}}\left(r(Q(\lambda))\right)=\mathbb{E}_{\text{update}}\left(r(Q(\lambda))\right)-1$
implies that $\mathbb{E}_{\text{update}}\left(r(Q(\lambda))\right)\geq 1$.
Moreover, $\mathbb{E}_{\text{update}}\left(r(Q(\lambda))\right)$ quantifies
the expected weight given to samples of $\mathbb{P}_{\text{pred}}$ that are
aligned with $\mathbb{P}_{\text{obs}}$.

At a computational level, if the $i+1$ iteration results in
$\mathbb{V}_{\text{init}}\left(\prod_{j\leq
i+1}r^{(j)}(Q^{(j)}(\lambda))\right)>\mathbb{V}_{\text{init}}\left(\prod_{j\leq
i}r^{(j)}(Q^{(j)}(\lambda))\right)$, then the $i+1$ iteration is providing a
statistically substantial change to the iterative DCI solution. Additionally,
since
$\mathbb{V}_{\text{init}}\left(r(Q(\lambda))\right)=\mathbb{E}_{\text{init}}\left(r^{2}(Q(\lambda))\right)-1$,
we can simply use the $r^{2}$-values of iid samples drawn from the initial
distribution and compute the sample average analogous to the DCI diagnostic
that computes the sample average of the $r$-values. That is, an optional
“significance” check can be added to Algorithm 1, which we summarize in
Algorithm 2 below.

Algorithm 2 $i+1$ Iteration of Iterative DCI with Significance Check

1:apply kPCA to spatial data at current time step to learn new $d$-dimensional
QoI map $Q^{(i+1)}(\lambda)$

2:evaluate $r=\prod_{j\leq i}r^{(j)}(Q^{(j)}(\lambda))$ at initial parameter
samples

3:estimate $\pi_{\text{obs}}$ from learned QoI samples and $\pi_{\text{pred}}$
as a weighted KDE of learned QoI samples weighted by $r$

4:evaluate
$r^{(i+1)}(Q^{(i+1)}(\lambda))=\pi_{\text{obs}}(Q^{(i+1)}(\lambda))/\pi_{\text{pred}}(Q^{(i+1)}(\lambda))$
at learned predicted QoI samples

5:if $|1-\text{mean}(r^{(i+1)}r)|>tol$ then

6: restart using $(d\leftarrow d-1)$ or discard time step if $d=1$

7:else if $\text{mean}((r^{(i+1)}r)^{2})\leq\text{mean}(r^{2})$ then

8: discard time step

An important feature of the iterative DCI framework is that it allows
immediate utilization of new observed data to update the DCI solution without
requiring any re-computation of prior information. In other words, once an
iterative solution is computed, the only values required to store in memory
are the original samples of $\pi_{\text{init}}$, the current $r$-values at
these samples, and any model information associated with the initial samples
required to generate predicted data at the next potential iteration. It is
also possible to setup either algorithm to utilize LUQ with different
filtering options (in terms of centering data, utilizing different depths of
NNs, and/or including polynomial trend fitting) at each iteration as the
complexity of the underlying spatial signal at each time evolves.

###  5.2 Illustrative Example: Randomly Generated Waves (Part IV)

Revisiting the wave equation example, we show how spatio-temporal data can be
fully utilized. Suppose that there are sensors at each $(0.5i,0.5j)$ for
$i,j=1,...,9$ that record the wave height every 0.5 time units from $t=2.5$ to
$t=7.0$. As in Part III of this illustrative example, “white” noise is
simulated by a normal distribution cetered at 0.0 with standard deviation of
2.5E-3. The spatial data at each time step are filtered by iteratively fitting
1-7 RBFs to approximate the wave’s deviation from its approximate mean value
by starting with 1 RBF and increasing to 7 RBFs stopping if a relative error
of less than 0.01 is reached. We use the same filtered data coordinates as in
Part III of this illustrative example. At each time step, kPCA is performed
and the corresponding predicted and observed densities are estimated from the
transformed samples corresponding to the first two principal components. We
consider both Algorithms 1 and 2 with $tol$ to $0.095$. The resulting updated
joint densities using both algorithm are compared in Figure 9, and the
marginals are compared in Figure 10.

Figure 9: Color plots of the DCI solutions using $10E3$ uniform samples
colored by the updated densities for part IV of the illustrative example with
data-generating samples overlaid in black. The left plot is from using
Algorithm 1 that does not utilize the $r^{2}$ values, and the right plot is
from using Algorithm 2 that does utilize the $r^{2}$ values. The resulting TV
metrics are 0.161 for Algorithm 1 and 0.160 for Algorithm 2. For reference,
the TV between the initial uniform distribution and the exact data-generating
distribution is about 0.808 while the TV between the KDE of the data-
generating distribution and the exact data-generating distribution is about
0.119.

(a) Algorithm 1

(b) Algorithm 2

Figure 10: Updated density marginals for wave equation example part IV along
with the original data-generating marginals. Top is from using Algorithm 1;
bottom is from using Algorithm 2. The resulting TV metrics for parameter $a$
are 0.034 and 0.037 for Algorithms 1 and 2, respectively; for parameter $b$,
the metrics are 0.071 and 0.037 for Algorithms 1 and 2, respectively. These
values are within the expected range of metric values of about 0.01 to 0.08
for parameter $a$ and 0.01 to 0.10 for parameter $b$ as seen in Figure 3.

The TV metrics between the final updated density and the exact DG density is
about 0.161 using Algorithm 1 and about 0.160 using Algorithm 2. For the
marginal densities, the TV metrics between the marginals for $a$ are about
0.034 and 0.037 and for $b$ are about 0.071 and 0.037 when using Algorithms 1
and 2, respectively. Plots of the TV metrics at each iteration for both
algorithms are shown in Figure 11 (for the joint densities) and Figure 12 (for
the marginals). Note that the final values are within the range of TV values
expected from a KDE estimate of the DG density and its marginals as shown in
Figure 3.

Figure 11: TV metrics between updated iterative densities and the exact data-
generating distribution for the illustrative wave equation example part IV.
Left is from using Algorithm 1; right is from using Algorithm 2. The time of
each iteration used is annotated for each data point. The dotted lines show
the TV between the KDE of the data-generating distribution and the exact data-
generating distribution.

Figure 12: TV metrics between updated iterative marginal densities and the
exact data-generating distribution for the illustrative wave equation example
part IV. Left is from using Algorithm 1; right is from using Algorithm 2. The
time of each iteration used is annotated for each data point. The dotted lines
show the TV metrics between the KDEs of the data-generating marginal
distributions and the exact data-generating marginal distributions.

##  6 Learned QoI as Functions of Kernels and Numerical Sufficiency of
Filtered Data

In this section, we derive a closed-form representation of a learned QoI map
as a function of the kernel utilized within the kPCA approach applied in the
final step of the LUQ framework. We then apply the mathematical insight
provided by this representation to develop a quantitative procedure for
identifying when a sufficient number of filtered data are computed for
learning the QoI.

###  6.1 QoI maps as functions of kernels

At a high-level, a kernel function denoted by $\mathcal{K}(\cdot,\cdot)$ is a
measure of similarity, i.e., $\mathcal{K}(a,b)>\mathcal{K}(a,c)$ if objects
$a$ and $b$ are considered “more similar” than objects $a$ and $c$ [62]. In
the context of this work, the objects we consider are the filtered datasets.
Let $\mathcal{D}$ denote the space of unfiltered data (either observed or
predicted) and $\mathcal{F}(\mathcal{D})$ the space of filtered data obtained
by the filtering step of LUQ. With this notation, we assume that
$\mathcal{K}:\mathcal{F}(\mathcal{D})\times\mathcal{F}(\mathcal{D})\to[0,\infty)$
is symmetric and positive definite. With these assumptions, the kernel is both
a measure of similarity between filtered data vectors as well as a
representation of an inner product in the so-called reproducing kernel Hilbert
space (RKHS) with reproducing kernel $\mathcal{K}$. This Hilbert space,
denoted by $\mathcal{H}$, is the span of the collection
$\left\\{\mathcal{K}\left(\cdot,\bm{d}\right)\right\\}_{\bm{d}\in\mathcal{F}(\mathcal{D})}$.
In this context, for every $\bm{d}\in\mathcal{F}(\mathcal{D})$,
$\mathcal{K}\left(\cdot,\bm{d}\right)$ represents the unique element in
$\mathcal{H}$ such that for every $h\in\mathcal{H}$,

$h(\bm{d})=\langle\mathcal{K}\left(\cdot,\bm{d}\right),h(\cdot)\rangle_{\mathcal{H}}$
(10)

where $\langle\cdot,\cdot\rangle_{\mathcal{H}}$ is the inner product of
$\mathcal{H}$.

The above equation is referred to as the reproducing property of
$\mathcal{K}$. The Moore-Aronszajn theorem [51] implies that for each
symmetric positive definite kernel $\mathcal{K}$ defined on
$\mathcal{F}(\mathcal{D})$, there exists a unique RKHS $\mathcal{H}$ with
reproducing kernel $\mathcal{K}$. Working directly with the kernel in a RKHS
has advantages in many contexts. For instance, [63] presents analysis to
clarify the “dilemma” that directly working with the kernel “seemingly
def[ies] the curse of dimensionality” in the context of SVMs, ridge
regression, and regularization networks. The application of the kernel to
significantly reduce computational complexity for kernel ridge regression is
analyzed in [64]. Below, we work in the RKHS as a means of deriving the
closed-form representation of learned QoI maps as functions of the kernel
utilized in a kPCA.

Given a set of $N$ iid parameter samples drawn from the initial distribution,
let $\bm{d}^{(1)},\dots,\bm{d}^{(N)}\in\mathcal{F}(\mathcal{D})$ denote the
associated filtered data vectors. Given an arbitrary error function
$E:\left(\mathcal{F}(\mathcal{D})\times\mathbb{R}\right)^{N}\to\mathbb{R}\cup\left\\{\infty\right\\}$,
consider the optimization problem with a regularization constraint,

minimize: $\displaystyle
E\left(\left(\bm{d}^{(1)},h\left(\bm{d}^{(1)}\right)\right),\dots,\left(\bm{d}^{(N)},h\left(\bm{d}^{(N)}\right)\right)\right)$
subject to: $\displaystyle\lVert h\rVert_{\mathcal{H}}\leq c$

where $c$ is some positive constant and $\lVert
h\rVert_{\mathcal{H}}=\sqrt{\langle h,h\rangle_{\mathcal{H}}}$. Consider the
decomposition of $h\in\mathcal{H}$ as $h=h_{\parallel}+h_{\perp}$ where
$h_{\parallel}$ is in the span of
$\mathcal{K}\left(\cdot,\bm{d}^{(1)}\right),\dots,\mathcal{K}\left(\cdot,\bm{d}^{(N)}\right)$
and $h_{\perp}$ is in the orthogonal complement of this span. By the
reproducing property,
$h\left(\bm{d}^{(i)}\right)=\left\langle\mathcal{K}\left(\cdot,\bm{d}^{(i)}\right),h(\cdot)\right\rangle_{\mathcal{H}}$
for each $1\leq i\leq N$. It follows that
$h\left(\bm{d}^{(i)}\right)=h_{\parallel}\left(\bm{d}^{(i)}\right)$ for each
$1\leq i\leq N$. Subsequently, $E$ is unaffected by $h_{\perp}$. Furthermore,
since $\lVert h_{\parallel}\rVert_{\mathcal{H}}\leq\lVert
h\rVert_{\mathcal{H}}$, the optimal feasible solution, denoted by
$h^{*}\in\mathcal{H}$, is in the span of
$\left\\{\mathcal{K}\left(\cdot,\bm{d}^{(1)}\right),\dots,\mathcal{K}\left(\cdot,\bm{d}^{(N)}\right)\right\\}$.
In other words, the optimal feasible solution is given by choosing $h_{\perp}$
to be the zero function and determining the unique constants
$\alpha_{1},\dots,\alpha_{N}\in\mathbb{R}$ such that

$h^{*}=\sum_{i=1}^{N}\alpha_{i}\mathcal{K}\left(\cdot,\bm{d}^{(i)}\right)\in\mathcal{H}\
\Longrightarrow\
h^{*}(\bm{d})=\left\langle\mathcal{K}(\cdot,\bm{d}),\sum_{i=1}^{N}\alpha_{i}\mathcal{K}(\cdot,\bm{d}^{(i)})\right\rangle_{\mathcal{H}},\
\forall\bm{d}\in\mathcal{F}(\mathcal{D}).$

The reproducing property of $\mathcal{K}$ implies that
$\mathcal{K}\left(\bm{d},\bm{d}^{(i)}\right):=\left\langle\mathcal{K}\left(\cdot,\bm{d}\right),\mathcal{K}\left(\cdot,\bm{d}^{(i)}\right)\right\rangle$
for all $\bm{d}\in\mathcal{F}(\mathcal{D})$ so that $h^{*}$ can be constructed
by computing the $\alpha_{i}$’s without working in the RKHS directly.
Furthermore, the argument above implies that while $\mathcal{H}$ may be
infinite dimensional, we are only concerned with the finite dimensional
subspace spanned by
$\left\\{\mathcal{K}\left(\cdot,\bm{d}^{(1)}\right),\dots,\mathcal{K}\left(\cdot,\bm{d}^{(N)}\right)\right\\}$
allowing us to perfom linear PCA in $\mathcal{H}$. We can therefore interpret
a kPCA as applying a standard PCA on the data vectors after they are mapped to
$\mathcal{H}$ via the so called feature map $\phi$ defined as
$\phi\left(\bm{d}\right):=\mathcal{K}\left(\cdot,\bm{d}\right)$. With this
interpretation, the kPCA seeks the unit vector pointing in the direction that
maximizes the variance of the transformed _centered_ data. Denote this vector
by $h^{*}$. This vector subsequently defines the first component of the
learned QoI map. We make this mathematically precise below.

Define $D$ as the $N\times\text{dim}\left(\mathcal{F}(\mathcal{D})\right)$
matrix $D:=\left[\bm{d}^{(1)},\dots,\bm{d}^{(N)}\right]^{\top}$ and
$\phi\left(D\right):=\left[\phi\left(\bm{d}^{(1)}\right),\dots,\phi\left(\bm{d}^{(N)}\right)\right]^{\top}$.
The transformed centered data is denoted as
$\tilde{\phi}\left(D\right):=\phi\left(D\right)-\frac{1}{N}\boldsymbol{1}\sum_{i=1}^{N}\phi\left(\bm{d}^{(i)}\right)$
where $\boldsymbol{1}$ denotes the
$N\times\text{dim}\left(\mathcal{F}(\mathcal{D})\right)$ matrix of all ones.
With this notation, the kPCA seeks to solve the following optimization
problem:

maximize: $\displaystyle\mathbb{V}\left(\tilde{\phi}\left(D\right)h\right)$
subject to: $\displaystyle\lVert h\rVert_{\mathcal{H}}\leq 1$

Note that the variance operator
$\mathbb{V}\left(\tilde{\phi}\left(D\right)h\right)=\langle
h,\mathcal{C}h\rangle_{\mathcal{H}}$ where
$\mathcal{C}:=\frac{1}{N-1}\tilde{\phi}\left(D\right)\tilde{\phi}\left(D\right)^{\top}$
denotes the sample covariance matrix in the finite dimensional subspace of
$\mathcal{H}$ spanned by
$\left\\{\mathcal{K}\left(\cdot,\bm{d}^{(1)}\right),\dots,\mathcal{K}\left(\cdot,\bm{d}^{(N)}\right)\right\\}$.
Since the variance will increase with scaling of $h$, we take $\lVert
h\rVert_{\mathcal{H}}=1$ and by formulating the Lagrangian of the above
optimization problem in order to set the partial derivatives of the Lagrangian
with respect to $h$ and Lagrange multiplier $\ell$ equal to zero, one obtains
the eigenvalue equation

$\mathcal{C}h=\ell h$

and $\mathbb{V}\left(\tilde{\phi}\left(D\right)h\right)=\ell$. The above
optimization problem is therefore equivalent to solving the following least
squares projection problem:

minimize:
$\displaystyle\left\lVert\tilde{\phi}\left(D\right)-\tilde{\phi}\left(D\right)hh^{\top}\right\rVert_{\mathcal{H}}^{2}$
subject to: $\displaystyle\lVert h\rVert_{\mathcal{H}}\leq 1$

From the discussion above, we can rewrite the solution $h^{*}$ as
$h^{*}=\sum_{i=1}^{N}\alpha_{i}\tilde{\mathcal{K}}\left(\cdot,\bm{d}^{(i)}\right)$
where $\tilde{\mathcal{K}}$ is the centered kernel defined as

$\tilde{\mathcal{K}}\left(\bm{d}^{(i)},\bm{d}^{(j)}\right):=\left\langle\phi\left(\bm{d}^{(i)}\right)-\frac{1}{N}\sum_{n=1}^{N}\phi\left(\bm{d}^{(n)}\right),\phi\left(\bm{d}^{(j)}\right)-\frac{1}{N}\sum_{n=1}^{N}\phi\left(\bm{d}^{(n)}\right)\right\rangle_{\mathcal{H}}.$

Substituting this expression of $h^{*}$ into the eigenvalue equation
$\mathcal{C}h^{*}=\ell h^{*}$, we obtain the matrix equation
$\tilde{K}^{2}\bm{\alpha}=\left(N-1\right)\ell\tilde{K}\bm{\alpha}$ where
$\tilde{K}$ is the centered Gram matrix, i.e. the $N\times N$ matrix with
$\tilde{K}_{i,j}:=\tilde{\mathcal{K}}\left(\bm{d}^{(i)},\bm{d}^{(j)}\right)$,
and $\bm{\alpha}=\left[\alpha_{1},\dots,\alpha_{N}\right]^{\top}$. Since
$\tilde{K}$ is symmetric positive definite, then the matrix equation reduces
to

$\tilde{K}\bm{\alpha}=\left(N-1\right)\ell\bm{\alpha}.$

Note that $\lVert h\rVert_{\mathcal{H}}=1$ is equivalent to
$\lVert\bm{\alpha}\rVert_{2}=\frac{1}{\sqrt{(N-1)\ell}}$. Subsequently,
solving the eigenvalue equation
$\tilde{K}\bm{\alpha}=\left(N-1\right)\ell\bm{\alpha}$ for $\bm{\alpha}$ and
rescaling according to

$\bm{\alpha}^{(1)}\mapsfrom\dfrac{1}{(N-1)\ell}\dfrac{\bm{\alpha}}{\lVert\bm{\alpha}\rVert_{2}}$

produces the vector of constants defining
$h^{*}=\left[\mathcal{K}\left(\cdot,\bm{d}^{(1)}\right),\dots,\mathcal{K}\left(\cdot,\bm{d}^{(N)}\right)\right]\bm{\alpha}^{(1)}$.
We have thus derived the first component of the learned QoI map associated
with this kernel:

$Q_{1}\left(\cdot\right):=\sum_{j=1}^{N}\bm{\alpha}^{(1)}_{j}\tilde{\mathcal{K}}\left(\cdot,\bm{d}^{(j)}\right).$

For the other components of $Q$, we organize the similarly rescaled
eigenvectors $\bm{\alpha}^{(i)}$ of the eigenvalue equations
$\tilde{K}\bm{\alpha}=(N-1)\ell\bm{\alpha}$ in descending order of
corresponding eigenvalues $\ell_{i}$ (which is equivalent to the descending
order of variance). Utilizing the first $q$ of these eigenvalues and
eigenvectors, we construct the learned QoI map
$Q=\left[Q_{1},\dots,Q_{q}\right]^{\top}$ where

$Q_{i}\left(\cdot\right)=\sum_{j=1}^{N}\bm{\alpha}_{j}^{(i)}\tilde{\mathcal{K}}\left(\cdot,\bm{d}^{(j)}\right),\
1\leq i\leq q.$ (11)

###  6.2 Numerical Sufficiency Test for Filtered Data

In the above analysis, the data vectors
$\bm{d}^{(1)},\ldots,\bm{d}^{(N)}\in\mathcal{F}(\mathcal{D})$ are associated
with the $N$ iid parameter samples drawn from the initial distribution. From
(11), we see that the $\bm{\alpha}^{(i)}\in\mathbb{R}^{N}$ for each $1\leq
i\leq q$. In other words, the vectors learned from the kPCA that define each
component of the QoI map have dimension that is equal to the number of samples
drawn from the initial distribution. Consequently, the dimension of these
vectors is invariant to the dimension of $\mathcal{F}(\mathcal{D})$. This
observation is key to developing a quantitative diagnostic for analyzing the
sufficiency of filtered data.

Let $A$ and $B$ denote two potential sets of coordinates for which filtered
data are constructed by evaluation of the learned filtered functions in Step 1
of LUQ. In the context of how the following ideas are applied, it is helpful
to think of $B$ as a potential refinement of $A$. For instance, if $A$ is
defined by evaluating the filtered functions on an $n\times n$ grid of spatial
coordinates, suppose $B$ is defined by evaluating the filtered functions on an
$m\times m$ grid with $m>n$. Denote by $\mathcal{F}(\mathcal{D})^{(A)}$ and
$\mathcal{F}(\mathcal{D})^{(B)}$, respectively, the associated potential
filtered data spaces for these two sets of coordinates. Finally, for a given
component of the QoI vector, denote by $\bm{\alpha}^{(A)}$ and
$\bm{\alpha}^{(B)}$, respectively, the associated learned $\alpha$ vectors
from the kPCA analysis as given by (11).

Viewing the components of $\bm{\alpha}^{(A)}$ as predictors for the components
of $\bm{\alpha}^{(B)}$, we can use simple linear regression on these
components to quantify what is learned from refining the filtered coordinate
system. However, since the eigenvalues $\ell$ associated with the Gram matrix
that are used to rescale the $\bm{\alpha}$ vectors may vary significantly from
the one filtered coordinate system to another, we first normalize the
$\bm{\alpha}$ vectors with respect to the standard $2$-norm in
$\mathbb{R}^{N}$ before applying regression. For a given QoI component learned
from the two filtered coordinate systems $A$ and $B$, the filtered coordinates
of $A$ are deemed sufficient if the associated slope and coefficient of
determination are both near unity. In other words, the additional filtered
data obtained by the $B$ coordinate system do not significantly alter the QoI
component to justify the additional computational burden of increasing the
dimension. Below, we illustrate the application of the concept in the context
of evaluating sufficiency of filtered data for the example of randomly
generated waves. A future work will incorporate this perspective into a more
comprehensive OED framework that evaluates optimal unfiltered data collection
in conjunction with optimally learning the QoI from such data.

###  6.3 Illustrative Example: Randomly Generated Waves (Part V)

Revisiting the wave equation example from Part III in Section 4.4, recall that
the original mesh over the spatial domain is a uniformly spaced $99\times 99$
grid at spacial locations $\left(0.05i,0.05j\right)$, $i,j=1,\dots,99$. We
investigate potential filtered coordinates defined by sub-meshes of the
computational mesh. Specifically, we consider grid of sizes $4\times 4$,
$9\times 9$, $19\times 19$, $49\times 49$, and finally the original $99\times
99$ mesh. While we focus our attention on uniformly spaced meshes that
coincide with the original mesh for computational simplicity in filtering the
predicted data, this is not necessary due to the filtering method allowing the
data to be evaluated at any set of data coordinates.

Starting with the two coarsest potential filtered coordinate sets, we compute
the desired QoI maps, normalize the $\bm{\alpha}$ vectors for each component,
and perform the regression analysis as described above in Section 6.2 to
quantitatively compare the coefficient of determination denoted by $R^{2}$ and
the absolute value of the slope against values of unity. In the interest of
space, we only show results for comparing the $4\times 4$ grid with the
$9\times 9$ grid and comparing the $49\times 49$ grid with the $99\times 99$,
which are shown in Figure 13. The interested reader can follow the
instructions in A to obtain all the data and scripts required to generate the
intermediate comparisons. Note that going from the $4\times 4$ grid (i.e.,
using $16$ data points) to the $9\times 9$ grid (i.e., using $81$ data points)
produces significantly different QoI components with coefficients of
determination 0.496 and 0.243 and slopes -0.705 and -0.493 for each component,
respectively. Note that the sign of the slopes is arbitrary since multiplying
a given $\alpha$ obtained in the kPCA analysis by $-1$ does not alter the QoI
learned from this kPCA. By comparison, once we use the $49\times 49$ grid
(i.e.,using $2401$ data points) of potential filtered coordinates, we see that
going to a $99\times 99$ grid (i.e., using $9801$ data points) produces
minimal change in QoI components with coefficients of determination 0.813 and
0.783 and slopes -0.902 and -0.885 for each component, respectively. In other
words, there is strong statistical evidence for using the $49\times 49$
filtered coordinate system to learn the QoI map since any further refinement
of this mesh does not significantly alter the learned $\bm{\alpha}$ vectors
defining the QoI components. This does not imply that using one of the coarser
filtered data coordinate sets will provide a lower quality QoI map for DCI
than using the data points from the $49\times 49$ grid, but it is clear that
there is a limit to how much filtered data obtained at more data coordinates,
either in space or time, can effect the kPCA components. In fact, many of the
updated densities shown in prior parts of this illustrative example are
similar both qualitatively and quantitatively if we utilize the $9\times 9$
grid for the filtered coordinates, which the interested reader can easily
verify by adjusting the filtered data coordinates variables within the
provided data and scripts. At a practical level, the analysis discussed in
this section is useful for establishing a numerical method for identifying an
upper bound on the amount of potentially useful filtered data for learning the
QoI.

Figure 13: Plots of normalized $\boldsymbol{\alpha}$ vectors from performing
kernel PCA on the 4x4 and 9x9 grids (top) and from 49x49 and 99x99 grids
(bottom). For each row, the left plots are of $\boldsymbol{\alpha}^{(1)}$ from
each grid and the right plots are of $\boldsymbol{\alpha}^{(2)}$ from each
grid.

##  7 Numerical Experiments

Our numerical examples concern shells of revolution. Specifically, we seek to
use displacement data from a population of shells to infer distributions of
parameters defining material properties of the shells. All structures are
considered to be thin with dimensionless thickness of $1/100$ and a Poisson
number $\nu=0.3$ approximating steel. All realizations are solved using
$hp$-FEM. For details on the underlying elasticity equations see B.

###  7.1 Shell of Revolution: Random Profile and Material

Mechanical characteristics of shells of revolution are strongly dependent on
their shapes, or more precisely, types of midsurface geometry. Given an
interval $I$, we define a positive profile function $f(x)>0$, $x\in I$, which
revolves around the $x$-axis. In this example, we take $I=[-1,1]$. The
different kinds of shells of revolution formed in this way can be classified
in terms of Gaussian curvature (see for instance [65]) which can be described
in terms of the derivatives of the profile function $f(x)$. Assuming that the
second derivative has a constant sign over the whole interval $I$, one has

  * 1.

_Parabolic [Zero Gaussian curvature shells]_ if $f^{\prime\prime}(x)=0$ (see
left plot in Figure 14),

  * 2.

_Elliptic [Positive Gaussian curvature shells]_ if $f^{\prime\prime}(x)<0$
(see mid plot in Figure 14),

  * 3.

_Hyperbolic [Negative Gaussian curvature shells]_ if $f^{\prime\prime}(x)>0$
(see right plot in Figure 14).

(a) Parabolic: $f(x)=1$. (b) Elliptic: $f(x)=1-x^{2}/4$. (c) Hyperbolic:
$f(x)=1+x^{2}/4$.

Figure 14: Three shells of revolution: $x\in[-1,1]$.

In this example, both the Young’s modulus $E=a_{0}$ and profile function
$f(x)$ are considered uncertain. We assume based on the types of materials
used in the manufacturing as well as the machining of the shells that it is
known that $a_{0}$ should be within $[0.8,1.2]$ (dimensionless) and
$f(x)=1+a_{1}(x-1)(x+1)$ with $a_{1}$ varying within $[-0.2,0.2]$ to yield
shells that vary from elliptic/positive Gaussian curvature to
hyperbolic/negative Gaussian curvature for $x\in[-1,1]$. It follows that
$\Lambda=[0.8,1.2]\times[-0.2,0.2]$ in this example.

In the context of this example, it is useful to consider three classes of
shapes defined as “clearly elliptic”, “almost parabolic with mild
deformation”, or “clearly hyperbolic” that are associated, respectively, with
$a_{1}$ belonging to the interval $[-0.2,-0.075]$, $(-0.075,0.075)$, or
$[0.075,0.2]$.

Loading is a periodic unit “wind load” acting on the tranverse direction. In
the notation of B, this implies that the $w$-component is formally given by
$f_{w}(x,y)=\cos y$. This is intended to simulate the strong directional winds
that wind turbine towers are subjected to in a wind farm. As further described
in B, in the reduced 1D formulation, the load becomes $f_{w}(x)=1$, with the
angular wave number $k=1$.

####  7.1.1 Example highlights

This experiment demonstrates several features of the LUQ and DCI frameworks
not previously highlighted with the illustrative examples. First, we
demonstrate that updated distributions can reasonably approximate DG
distributions with relatively few samples from both the DG and initial
distributions. In doing this, we instantiate/re-use the same LUQ object on
several distinct observation data sets. This exploits a mathematical feature
of the updated density obtained within the DCI framework where any observed
density satisfying the predictability assumption immediately defines an
updated density via direct substitution into the numerator of either (3) or
(6) once the predicted density is estimated for a given initial density. In
other words, it is trivial to solve multiple stochastic inverse problems for
the same initial density once a predicted density is estimated. To properly
exploit this fact, we designed the LUQ software to allow for the instantiation
and saving of a LUQ object from initial/predicted sample data in cases where
either observable data are not yet available or multiple observed data sets
are to be analyzed. Additionally, to demonstrate the flexibility of the
filtering step in the LUQ code, we assume that the observable and predicted
data are both obtained via experiments with differing levels of signal-to-
noise (SNR) ratios defined by the ratio of the variance of “exact” data (the
signal) to variance in the noise and emphasize that it is unnecessary to
specify the SNR value to perform the filtering as this step is agnostic to
such details.

####  7.1.2 DCI and LUQ Setup

Input | Details
---|---
DG distributions |  |  $a_{0}\sim\text{Beta}(2,6)$ on $[0.8,1.2]$ |  $a_{1}\sim\text{Beta}(3,4)$ on $[-0.2,0.2]$
---|---
Init. distributions |  | $a_{0}\sim\mathcal{U}([0.8,1.2])$ | $a_{1}\sim\mathcal{U}([-0.2,0.2])$
---|---
DG sample sets | 3 sets of 75 samples from DG distributions
Init. sample sets | 1 set of 200 samples from init. distributions
Predicted data labels |  $0$ if $a_{1}\in[-0.2,-0.075]$, $1$ if $a_{1}\in(-0.075,0.075)$, and $2$ if $a_{1}\in[0.075,0.2]$
Observable data | 901 noisy data with SNR of 5 for each DG sample
Predicted data | 901 noisy data with SNR of 10 for each init. sample
# Filtered data | 60 filtered data for each observation of prediction sample
Table 1: Summary of setup details for DCI and LUQ frameworks in shell of
revolution example.

Table 1 summarizes the relevant setup details for both the DCI problem and LUQ
frameworks in this example. We assume (noisy) displacement data are obtained
via an electronic “strip” equipped with 901 equally spaced sensors that is
placed on the opposite side of each shell from the direction of forcing for
$x\in[-0.9,0.9]$. The noise is assumed to be additive Gaussian and iid at each
measurement location. The goal is to learn the distributions on $a_{0}$ and
$a_{1}$ from this noisy observed displacement data. For the purposes of this
example, the DG distributions are taken to be $a_{0}\sim\text{Beta}(2,6)$
scaled and shifted to $[0.8,1.2]$ and $a_{1}\sim\text{Beta}(3,4)$ scaled and
shifted to $[-0.2,0.2]$ while the initial distributions on $a_{0}$ and $a_{1}$
are assumed uniform across $\Lambda$. To simulate the difference between
controlled experiments versus field studies, we assume samples taken from the
initial distribution have measurement data producing a SNR of ten while
samples taken from the observed distribution have measurement data producing a
SNR of five. In other words, the “field” data is twice as noisy as the
“predicted” data. Figure 15 illustrates representative results for the
different levels of noise in these types of data and the estimated filtered
responses from LUQ utilizing splines. See A for details on obtaining the full
scripts and datasets utilized for creating this and the remaining figures in
this example.

Figure 15: Noisy data for shell of revolution example with different SNR
levels for predicted and observed data. The filtering step in LUQ can flexibly
approximate the underlying response even if each sample within a data set has
its own unique SNR level.

We use 200 initial samples and three different sets of 75 samples from the DG
distribution to demonstrate the simplicity of solving multiple stochastic
inverse problems and the LUQ code as mentioned in the example highlights. The
initial samples are drawn from a set of ${1}\mathrm{E}{4}$ simulations
associated with the initial distribution while the three sets of 75 samples
are drawn from a set of ${3}\mathrm{E}{3}$ simulations associated with the DG
distribution, and the interested reader can change both the number of samples
and how they are drawn from these sets in the provided code. The provided code
also considers two different LUQ setups that utilize either two or three
classes based on the shape characteristics (i.e., the $a_{1}$ values). In the
interest of space, we only present the results for the three classes case
since the diagnostics obtained from a pilot study involving the noise-free
data indicated better performance than in the two classes case (these pilot
studies are included in the code). As mentioned above, we define the three
classes as being for $a_{1}$ in the intervals $[-0.2,-0.075]$,
$(-0.075,0.075)$, and $[0.075,0.2]$, which roughly correspond to the shapes
being either “clearly elliptic”, “almost parabolic with mild deformation”, or
“clearly hyperbolic”. The interested reader can alter the ranges of these
intervals or the number of classes defined and use the diagnostics discussed
within this example and provided within the code as guides for assessing the
quality of their choices.

####  7.1.3 Results

Output | Details
---|---
Diagnostic: $\mathbb{E}_{\text{init}}(r(Q(\lambda)))\approx 1$ |  $0.934$ (DG set 1) $0.937$ (DG set 2) $0.928$ (DG set 3)
TV between $a_{0}$ update/DG marginals | 0.086 (DG set 1) 0.096 (DG set 2) 0.099 (DG set 3)
TV between $a_{1}$ update/DG marginals | 0.083 (DG set 1) 0.028 (DG set 2) 0.060 (DG set 3)
Table 2: Summary of output details for DCI in shell of revolution example.

Table 2 summarizes the relevant outputs details for DCI in this example. To
obtain these results requires the filtering of data, training of the SVM
classifier, and learning of QoI in the LUQ framework. To accomplish those
steps, the predicted dataset generated from the initial samples is used to
instantiate a LUQ object where the predicted data are subsequently filtered,
utilized to learn/train an SVM classifier for the three-clusters associated
with $a_{1}$ values, and learn the 2-dimensional QoI maps. While 901 noisy
data are utilized in the filtering, a total of 60 filtered data are produced
in the filtering step to learn the QoI maps. The sufficiency of 60 filtered
data is verified in the code by comparing to 30 filtered data in the manner
described in Section 6.2 and demonstrated on the illustrative example in
Section 6.3. In this case, we observe both slopes and $R^{2}$ values from the
regression analysis are all above 0.95 for each component of the QoI map
across all three clusters indicating sufficiency in the number of filtered
data utilized to learn the QoI.

(a) Data-generating set 1.

(b) Data-generating set 2.

(c) Data-generating set 3.

Figure 16: Marginal densities for $a_{0}$ for shell of revolution example from
three different data-generating sets of 75 samples. The horizontal line
indicates the uniform initial density, the dashed line indicates the KDE of
the DG distribution computed from the samples associated with the given
observation data set, and the dashed-dotted line indicates the KDE of the
updated density.

(a) Data-generating set 1.

(b) Data-generating set 2.

(c) Data-generating set 3.

Figure 17: Marginal densities for $a_{1}$ for shell of revolution example from
three different data-generating sets of 75 samples. The horizontal line
indicates the uniform initial density, the dashed line indicates the KDE of
the DG distribution computed from the samples associated with the given
observation data set, and the dashed-dotted line indicates the KDE of the
updated density. The vertical dashed-dotted line at $a_{1}=0$ indicates where
the shells change from elliptic to hyperbolic in shape.

The LUQ object is then separately applied to each of the three distinct
observation sets. Given the clustering and relatively small number of samples
utilized in the DCI solution overall and the clusters specifically, we compute
a “weighted mean diagnostic” defined by the overall sample average of the
updated density given in (6) to assess the predictability assumption and
quality of construction of the updated density as a whole. For observation
sets 1, 2, and 3, the associated diagnostics are $0.934,0.937$ and $0.928$,
which are all within 10% of unity and indicate that both the predictability
assumption is satisfied and the updated density is well approximated.

Having numerically verified the predictability assumption, we now evaluate how
well the densities for the individual parameters are estimated in each of the
observation sets. Figures 16 and 17 show the associated uniform initial
densities (horizontal line), KDE approximation of the DG distribution
constructed from the samples associated with the observation set (dashed
line), and KDE of the updated density (dashed-dotted line) for $a_{0}$ and
$a_{1}$, respectively. For the $a_{1}$ density plots, we also provide a
vertical dashed-dotted line at $a_{1}=0$ as a visual indicator of where the
shapes change from elliptic (on the left) to hyperbolic (on the right). It is
worth emphasizing that there is clear variation between the KDEs of the DG
distributions within the $a_{0}$ and $a_{1}$ plots due to the relatively small
sample size (75 samples) utilized for each of the observation sets. The
marginal density plots qualitatively indicate that each of the updated
densities, which are all computed by re-weighting the same 200 initial samples
in every plot, are reasonable approximations to the associated KDE of the DG
distribution that generated the observable data. Moreover, the TV metrics are
all under 0.1 for each marginal of each parameter across all sets of
observable data, which as the illustrative examples have demonstrated, is a
value we expect for a “good” KDE approximation of a given distribution.
Specifically, the TV metrics (rounded to three decimal places) between the
KDEs of the updated and DG densities for sets 1, 2, and 3 for $a_{0}$ are
$0.086$, $0.096$, and $0.099$ and for $a_{1}$ are $0.083$, $0.028$, and
$0.060$, respectfully.

###  7.2 Trommel Screen: MultiPanel Construction with Material Variation

The second example is a trommel screen with three sections each with different
regular penetration patterns but constant 35% hole coverage as shown in the
introduction in Figure 1a. The computational domain is a square
$\Omega=[-\pi,\pi]\times[0,2\pi]$. The boundaries at $x=\pm\pi$ are clamped.
The bottom ($y=0$) and the top ($y=2\pi$) are periodic. All holes are free,
that is, there are no kinematical constraints. As with the previous example,
we refer to B for more details related to the underlying elasticity equations
and solution of the model.

In this example, due to different penetration patterns in each of the three
sections and variations in the manufacturing, the Young’s modulus $E$ is
considered uncertain in each section. We let $a_{0}$, $a_{1}$, and $a_{2}$
denote, respectively, the values of Young’s modulus in the left-, middle-, and
right-sections of each trommel screen. We assume that due to the larger holes
in the the left-section that the material is weakest with $a_{0}$ taking
values in $[0.576,0.704]$. The values of $a_{1}$ and $a_{2}$ are assumed,
respectively, to take values between $[0.72,0.88]$ and $[0.9,1.1]$. These
intervals represent 10% variation around nominal design values of $0.6,0.8,$
and $1.0$, respectively.

Loading is assumed to be a linear combination of concentrated loads acting on
specific points $(x_{k},y_{k})$. The model for the concentrated load is simply

$f(x,y)=\exp(-100((x-x_{k})^{2}+(y-y_{k})^{2})).$ (12)

We set $y_{k}=\pi$ (constant), and let $x_{k}=-2\pi/3+(k-1)2\pi/3$, $k=1,2,3$.
Large debris that is separated by the left-section of the trommel will have
impact points associated with $k=1$ values. Moderate sized debris separated by
the middle-section of the trommel will have impacts points associated with
$k=2$ values. Finally, small debris separated by the right-section of the
trommel will have impact points associated with $k=3$.

####  7.2.1 Example highlights

Here, we assume that three separate experiments are carried out across a
random sample of trommel screens in controlled experiments (with known
parameter values) and a different number of trommel screens in field studies
(with unknown parameter values). In each experiment, either “large”,
“moderate”, or “small” debris are separated by the trommel with sensors
recording the displacement from debris impact inside the trommel to finite
precision, which represents noise in terms of “loss of information/precision”
in the data. In each case, these datasets are available synchronously meaning
that the relative order of obtaining the spatial displacement data is
considered unknown. We subsequently utilize this data within the LUQ framework
by bypassing the filtering step for each dataset (since noise is due to
truncation in precision), learn a distinct number of QoI for each experiment,
and then iterate through the DCI procedure to arrive at reasonable updated
estimates of the DG distributions for each parameter. Diagnostics and
statistical analysis of data guided the iteration order and number of QoI
learned for each experiment, but the interested reader can alter all of these
choices and observe the subsequent results and diagnostics.

####  7.2.2 DCI and LUQ Setup

Input | Details
---|---
DG distribution $(a_{0})$ |  $a_{0}\sim\text{Beta}(2,6)$ on $[0.576,0.705]$
DG distribution $(a_{1})$ |  $a_{1}\sim\text{Beta}(3,4)$ on $[0.72,0.88]$
DG Distribution $(a_{2})$ |  $a_{2}\sim\text{Beta}(2,4)$ on $[0.9,1.1]$
Init. distrs. | Uniform on respective domains
DG sample set | 50 samples from DG distributions
Init. sample set | 500 samples from init. distributions
Observable & Predicted data | Sensor data (see Figure 1b) truncated to 2 digits each sample
Table 3: Summary of setup details for DCI and LUQ frameworks in trommel
example.

Table 3 summarizes the relevant setup details for both the DCI and LUQ
frameworks in this example. We assume that the measurements for each
experiment are truncated at two decimal places for both the observed and
predicted data. The goal is to learn the distributions of $a_{0},a_{1}$, and
$a_{2}$ from this finite precision data obtained from the separate
experiments. For the purposes of this example, the DG distributions are
$a_{0}\sim\text{Beta}(2,6)$, $a_{1}\sim\text{Beta}(3,4)$, and
$a_{2}\sim\text{Beta}(2,4)$ appropriately scaled and shifted to their
respective domains. The initial distributions are assumed to be uniform for
each parameter.

We use $500$ initial samples and $50$ DG samples to generate the predicted and
observed datasets. To three decimal places, the variances in the “full
precision” predicted data associated with the large, moderate, and small
debris are $0.101$, $0.181$, and $0.015$, which represents significant loss of
information from the finite precision measurements associated with truncating
all data to two decimal places. Since no filtering nor classification is
required, a LUQ object is instantiated for the finite-precision predicted data
associated with each experiment to learn a 2-dimensional, 3-dimensional, and
1-dimensional QoI map for the large, moderate, and small debris experiments,
respectively. The observed data from each experiment are subsequently
transformed into samples associated with their learned QoI maps.

####  7.2.3 Results

Output | Details
---|---
Diagnostic: $\mathbb{E}_{\text{init}}(r(Q(\lambda)))\approx 1$ |  $0.974$ (iter 1, 2 QoI) $0.938$ (iter 2, 3 QoI) $0.937$ (iter 3, 1 QoI)
TV between $a_{0}$ update/DG marginals | 0.090 (iter 1, 2 QoI) 0.064 (iter 2, 3 QoI) 0.064 (iter 3, 1 QoI)
TV between $a_{1}$ update/DG marginals | 0.062 (iter 1, 2 QoI) 0.018 (iter 2, 3 QoI) 0.020 (iter 3, 1 QoI)
TV between $a_{2}$ update/DG marginals | 0.313 (iter 1, 2 QoI) 0.119 (iter 2, 3 QoI) 0.090 (iter 3, 1 QoI)
Table 4: Summary of output details for DCI in shell of revolution example.

We iterate through the large, moderate, and small debris experiments to
sequentially update the initial distributions. This is similar to the temporal
iteration described in Section 5 where an updated distribution obtained from
one experiment serves as the initial distribution for the next experiment. The
results are summarized in Figure 18 and Table 4. We observe that the first
iteration, for which we learn a 2-dimensional QoI map and is associated with
the large debris filtered through the left-section of the trommel appears to
mostly learn the $a_{0}$ DG distribution while also learning some aspects of
the $a_{1}$ DG distribution and providing no significant update to the $a_{2}$
initial distribution. A general improvement across all three parameters is
observed as we iterate across the experiments. By the final experiment, the TV
metric between the KDE estimate of the updated marginal and DG marginal for
each parameter is under $0.1$ indicating a reasonable approximation across all
parameters. However, it is clearly the case that the updated marginal for the
$a_{2}$ parameter is the worst estimate of its associated DG distribution.
This is easily explained as the $a_{2}$ parameter is associated with the
“stiffest” part of the trommel and the final experiment that involves impact
loads on this part of the trommel involves the small debris that produces the
smallest variation in the data and is most affected by the finite precision of
the measurements. The interested reader can simply improve the precision in
the predicted and observed data in the provided scripts to see how this
improves these estimates (and allows for more QoI to be learned and utilized
for futher improvement as well). This is related to a future work exploring
optimal experimental design questions within the LUQ and DCI frameworks.

(a) $a_{0}$ densities

(b) $a_{1}$

(c) $a_{2}$

Figure 18: Marginal densities for the trommel screen parameters associated
with left- ($a_{0}$), middle- ($a_{1}$), and right-sections ($a_{2})$ of the
screen shown in Figure 1a. The updated marginals are iteratively obtained by a
sequence of experiments involving “large”, “moderate”, and “small” debris that
cause impacts at various points on the trommel screen shown in Figure 1b and
are filtered by the left-, middle-, and right-sections, respectively, of the
trommel.

##  8 Conclusions

Uncertainties in the computational models of engineered or physical systems
governed by principles of mechanics served as the motivation for this work. We
consider both aleatoric (i.e., irreducible) and epistemic (i.e., reducible)
sources of uncertainty. To quantify these uncertainties, we combine two
frameworks. The Learning Uncertain Quantities (LUQ) framework defines a formal
three-step machine-learning enabled process for transforming datasets into
samples of a learned Quantities of Interest (QoI) map to enable solution of a
stochastic inverse problem posed in a data-consistent (DC) framework.

In this work, we develop a robust filtering step in LUQ that can learn the
most useful quantitative information present in spatio-temporal datasets.
Reproducing Kernel Hilbert Space (RKHS) theory is leveraged to mathematically
analyze the learned QoI map and develop a quantitative sufficiency test for
evaluating the filtered data. An iterative DC-based inversion scheme is also
developed and considered including a new diagnostic to identify both the
quality and impact of inversion at each iteration.

The practical utility of these contributions is demonstrated through the
application of three problems: an illustrative example involving random
excitation of waves, the simulated manufacturing of thin shells of revolution
with random shape profiles and uncertain stiffness, and the simulated
manufacturing of trommel screens that are utilized in industrial applications
involving the mechanical separation of materials of different sizes.

##  9 Acknowledgments

T. Butler and T. Roper’s work is supported by the National Science Foundation
under Grant No. DMS-2208460. T. Butler’s work is also supported by NSF IR/D
program, while working at National Science Foundation. However, any opinion,
finding, and conclusions or recommendations expressed in this material are
those of the author and do not necessarily reflect the views of the National
Science Foundation.

## References

  * Chapelle and Bathe [2003] D. Chapelle, K. J. Bathe, The Finite Element Analysis of Shells, Springer–Verlag, 2003.
  * Bai and Jin [2016] Y. Bai, W.-L. Jin,  Chapter 33 - random variables and uncertainty analysis,  in: Y. Bai, W.-L. Jin (Eds.), Marine Structural Design (Second Edition), second edition ed., Butterworth-Heinemann, Oxford, 2016, pp. 615–625. URL: https://www.sciencedirect.com/science/article/pii/B9780080999975000332. doi:https://doi.org/10.1016/B978-0-08-099997-5.00033-2.
  * Roy and Oberkampf [2011] C. J. Roy, W. L. Oberkampf,  A comprehensive framework for verification, validation, and uncertainty quantification in scientific computing,  Computer Methods in Applied Mechanics and Engineering 200 (2011) 2131 – 2144. URL: http://www.sciencedirect.com/science/article/pii/S0045782511001290. doi:https://doi.org/10.1016/j.cma.2011.03.016.
  * Tran and Wildey [2021] A. Tran, T. Wildey,  Solving stochastic inverse problems for property–structure linkages using data-consistent inversion and machine learning,  JOM 73 (2021) 72–89.
  * Butler and Hakula [2020] T. Butler, H. Hakula,  What do we hear from a drum? A data-consistent approach to quantifying irreducible uncertainty on model inputs by extracting information from correlated model output data,  Computer Methods in Applied Mechanics and Engineering 370 (2020) 113228\. doi:10.1016/j.cma.2020.113228.
  * Mateker et al. [2013] W. R. Mateker, J. D. Douglas, C. Cabanetos, I. T. Sachs-Quintana, J. A. Bartelt, E. T. Hoke, A. El Labban, P. M. Beaujuge, J. M. J. Fréchet, M. D. McGehee,  Improving the long-term stability of pbdttpd polymer solar cells through material purification aimed at removing organic impurities,  Energy Environ. Sci. 6 (2013) 2529–2537. URL: http://dx.doi.org/10.1039/C3EE41328D. doi:10.1039/C3EE41328D.
  * Liu and Song [2013] M. Liu, G.-L. Song,  Impurity control and corrosion resistance of magnesium–aluminum alloy,  Corrosion Science 77 (2013) 143 – 150. URL: http://www.sciencedirect.com/science/article/pii/S0010938X13003478. doi:https://doi.org/10.1016/j.corsci.2013.07.037.
  * Lu and Hellawell [1987] S.-Z. Lu, A. Hellawell,  The mechanism of silicon modification in aluminum-silicon alloys: Impurity induced twinning,  Metallurgical Transactions A 18 (1987) 1721–1733. URL: https://doi.org/10.1007/BF02646204. doi:10.1007/BF02646204.
  * Ramesh et al. [2000a] R. Ramesh, M. Mannan, A. Poo,  Error compensation in machine tools — a review: Part i: geometric, cutting-force induced and fixture-dependent errors,  International Journal of Machine Tools and Manufacture 40 (2000a) 1235 – 1256. URL: http://www.sciencedirect.com/science/article/pii/S0890695500000092. doi:https://doi.org/10.1016/S0890-6955(00)00009-2.
  * Ramesh et al. [2000b] R. Ramesh, M. Mannan, A. Poo,  Error compensation in machine tools — a review: Part ii: thermal errors,  International Journal of Machine Tools and Manufacture 40 (2000b) 1257 – 1284. URL: http://www.sciencedirect.com/science/article/pii/S0890695500000109. doi:https://doi.org/10.1016/S0890-6955(00)00010-9.
  * Abellan-Nebot et al. [2012] J. V. Abellan-Nebot, J. Liu, F. R. Subirón, J. Shi,  State Space Modeling of Variation Propagation in Multistation Machining Processes Considering Machining-Induced Variations,  Journal of Manufacturing Science and Engineering 134 (2012). URL: https://doi.org/10.1115/1.4005790. doi:10.1115/1.4005790. arXiv:https://asmedigitalcollection.asme.org/manufacturingscience/article-pdf/134/2/021002/5637745/021002_1.pdf, 021002\.
  * Kostic et al. [2022] S. Kostic, J. Miljojkovic, G. Simunovic, D. Vukelic, B. Tadic,  Uncertainty in the determination of elastic modulus by tensile testing,  Engineering Science and Technology, an International Journal 25 (2022) 100998\. URL: https://www.sciencedirect.com/science/article/pii/S2215098621001105. doi:https://doi.org/10.1016/j.jestch.2021.05.002.
  * Glaub et al. [1982] J. C. Glaub, D. B. Jones, G. M. Savage,  Design and use of trommel screens for processing municipal solid waste,  1982\. URL: https://api.semanticscholar.org/CorpusID:4181072.
  * Kapteyn et al. [2021] M. G. Kapteyn, J. V. R. Pretorius, K. E. Willcox,  A probabilistic graphical model foundation for enabling predictive digital twins at scale,  Nature Computational Science 1 (2021) 337–347. URL: https://doi.org/10.1038/s43588-021-00069-0. doi:10.1038/s43588-021-00069-0.
  * Calvetti et al. [2014] D. Calvetti, J. Kaipio, E. Somersalo,  Inverse problems in the Bayesian framework,  Inverse Problems 30 (2014) 110301. URL: http://stacks.iop.org/0266-5611/30/i=11/a=110301.
  * Gelman et al. [2013] A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, D. B. Rubin, Bayesian Data Analysis, Third Edition, Chapman and Hall/CRC, 2013.
  * Cotter et al. [2010] S. Cotter, M. Dashti, A. Stuart,  Approximation of bayesian inverse problems,  SIAM Journal of Numerical Analysis 48 (2010) 322–345.
  * Kennedy and O’Hagan [2001] M. C. Kennedy, A. O’Hagan,  Bayesian calibration of computer models,  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 63 (2001) 425–464. URL: http://dx.doi.org/10.1111/1467-9868.00294. doi:10.1111/1467-9868.00294.
  * Berger et al. [1994] J. O. Berger, E. Moreno, L. R. Pericchi, M. J. Bayarri, J. M. Bernardo, J. A. Cano, J. De la Horra, J. Martín, D. Ríos-Insúa, B. Betrò, A. Dasgupta, P. Gustafson, L. Wasserman, J. B. Kadane, C. Srinivasan, M. Lavine, A. O’Hagan, W. Polasek, C. P. Robert, C. Goutis, F. Ruggeri, G. Salinetti, S. Sivaganesan,  An overview of robust bayesian analysis,  Test 3 (1994) 5–124. URL: http://dx.doi.org/10.1007/BF02562676. doi:10.1007/BF02562676.
  * Fitzpatrick [1991] B. Fitzpatrick,  Bayesian analysis in inverse problems,  Inverse Problems 7 (1991) 675. URL: {http://stacks.iop.org/0266-5611/7/i=5/a=003}.
  * Burger and Lucka [2014] M. Burger, F. Lucka,  Maximum a posteriori estimates in linear inverse problems with log-concave priors are proper bayes estimators,  Inverse Problems 30 (2014) 114004. URL: https://doi.org/10.1088/0266-5611/30/11/114004. doi:10.1088/0266-5611/30/11/114004.
  * Petra et al. [2014] N. Petra, J. Martin, G. Stadler, O. Ghattas,  A computational framework for infinite-dimensional bayesian inverse problems, part ii: Stochastic newton mcmc with application to ice sheet flow inverse problems,  SIAM Journal on Scientific Computing 36 (2014) A1525–A1555. URL: https://doi.org/10.1137/130934805. doi:10.1137/130934805. arXiv:https://doi.org/10.1137/130934805.
  * Alexanderian et al. [2016] A. Alexanderian, N. Petra, G. Stadler, O. Ghattas,  A fast and scalable method for a-optimal design of experiments for infinite-dimensional bayesian nonlinear inverse problems,  SIAM Journal on Scientific Computing 38 (2016) A243–A272. URL: https://doi.org/10.1137/140992564. doi:10.1137/140992564. arXiv:https://doi.org/10.1137/140992564.
  * Vaart [1998] A. W. v. d. Vaart, Asymptotic statistics, Cambridge University Press, Cambridge, UK;New York, NY, USA;, 1998.
  * Wikle et al. [1998] C. K. Wikle, L. M. Berliner, N. Cressie,  Hierarchical bayesian space-time models,  Environmental and Ecological Statistics 5 (1998) 117–154. URL: https://doi.org/10.1023/A:1009662704779. doi:10.1023/A:1009662704779.
  * Gelfand and Kottas [2002] A. E. Gelfand, A. Kottas,  A computational approach for full nonparametric bayesian inference under dirichlet process mixture models,  2002\.
  * Teh [2010] Y. W. Teh, Dirichlet Process, Springer US, Boston, MA, 2010, pp. 280–287. URL: https://doi.org/10.1007/978-0-387-30164-8_219. doi:10.1007/978-0-387-30164-8\\_219.
  * Butler et al. [2014] T. Butler, D. Estep, S. Tavener, C. Dawson, J. Westerink,  A Measure-Theoretic Computational Method for Inverse Sensitivity Problems III: Multiple Quantities of Interest,  SIAM/ASA Journal on Uncertainty Quantification 2 (2014) 174–202. URL: http://dx.doi.org/10.1137/130930406. doi:10.1137/130930406. arXiv:http://dx.doi.org/10.1137/130930406, doi:10.1137/130930406.
  * Butler et al. [2018] T. Butler, J. Jakeman, T. Wildey,  Combining Push-Forward Measures and Bayes’ Rule to Construct Consistent Solutions to Stochastic Inverse Problems,  SIAM Journal on Scientific Computing 40 (2018) A984–A1011. URL: https://doi.org/10.1137/16M1087229. doi:10.1137/16M1087229. arXiv:https://doi.org/10.1137/16M1087229.
  * Bingham et al. [2024] D. Bingham, T. Butler, D. Estep,  Inverse problems for physics-based process models,  Annual Review of Statistics and Its Application 11 (2024) null. URL: https://doi.org/10.1146/annurev-statistics-031017-100108. doi:10.1146/annurev-statistics-031017-100108. arXiv:https://doi.org/10.1146/annurev-statistics-031017-100108.
  * Change and Pollard [1997] J. Change, D. Pollard,  Conditioning as disintegration,  Statistica Neerlandica 51 (1997) 287–317. doi:10.1111/1467-9574.00056.
  * Saghafi et al. [2023] S. Saghafi, T. Rumbell, V. Gurev, J. Kozloski, F. Tamagnini, K. C. A. Wedgwood, C. O. Diekman,  Inferring parameters of pyramidal neuron excitability in mouse models of alzheimer’s disease using biophysical modeling and deep learning,  bioRxiv (2023). URL: https://www.biorxiv.org/content/early/2023/04/19/2023.04.18.537149. doi:10.1101/2023.04.18.537149. arXiv:https://www.biorxiv.org/content/early/2023/04/19/2023.04.18.537149.full.pdf.
  * Zhang and Mikelsons [2023] Y. Zhang, L. Mikelsons,  Solving stochastic inverse problems with stochastic bayesflow,  in: 2023 IEEE/ASME International Conference on Advanced Intelligent Mechatronics (AIM), 2023, pp. 966–972. doi:10.1109/AIM46323.2023.10196190.
  * Rumbell et al. [2023a] T. Rumbell, J. Parikh, J. Kozloski, V. Gurev, Novel and flexible parameter estimation methods for data-consistent inversion in mechanistic modeling, 2023a. arXiv:2009.08267.
  * Rumbell et al. [2023b] T. Rumbell, C. Wanjiru, I. O. Mulang’, S. Obonyo, J. Kozloski, V. Gurev,  Sequential data-consistent model inversion,  in: NeurIPS 2023 Workshop on Deep Learning and Inverse Problems, 2023b. URL: https://openreview.net/forum?id=5NoVk8nBWc.
  * Butler et al. [2022] T. Butler, T. Wildey, W. Zhang,  $l^{p}$ convergence of approximate maps and probability densities for forward and inverse problems in uncertainty quantification,  International Journal for Uncertainty Quantification 12 (2022). URL: https://par.nsf.gov/biblio/10341097. doi:10.1615/Int.J.UncertaintyQuantification.2022038086.
  * Bruder et al. [????] L. Bruder, M. Gee, T. Wildey,  Data-consistent solutions to stochastic inverse problems using a probabilistic multi-fidelity method based on conditional densities,  International Journal for Uncertainty Quantificaiton 10 (????) 399–424.
  * Poole and Raftery [2000] D. Poole, A. E. Raftery,  Inference for deterministic simulation models: The bayesian melding approach,  Journal of the American Statistical Association 95 (2000) 1244–1255. doi:10.1080/01621459.2000.10474324.
  * Butler et al. [2018] T. Butler, J. Jakeman, T. Wildey,  Convergence of probability densities using approximate models for forward and inverse problems in uncertainty quantification,  SIAM Journal on Scientific Computing 40 (2018) A3523–A3548. URL: https://doi.org/10.1137/18M1181675. doi:10.1137/18M1181675. arXiv:https://doi.org/10.1137/18M1181675.
  * Pilosov et al. [2023] M. Pilosov, C. del Castillo-Negrete, T. Y. Yen, T. Butler, C. Dawson,  Parameter estimation with maximal updated densities,  Computer Methods in Applied Mechanics and Engineering 407 (2023) 115906\. URL: https://www.sciencedirect.com/science/article/pii/S0045782523000294. doi:https://doi.org/10.1016/j.cma.2023.115906.
  * Butler et al. [2020] T. Butler, T. Wildey, T. Y. Yen,  Data-consistent inversion for stochastic input-to-output maps,  Inverse Problems 36 (2020) 085015. doi:10.1088/1361-6420/ab8f83.
  * Mattis et al. [2022] S. Mattis, K. Steffen, T. Butler, C. Dawson, D. Estep,  Learning quantities of interest from dynamical systems for observation-consistent inversion,  Computer Methods in Applied Mechanics and Engineering 388 (2022) 114230\. URL: https://www.sciencedirect.com/science/article/pii/S0045782521005582. doi:https://doi.org/10.1016/j.cma.2021.114230.
  * Pearson [1901] K. Pearson,  LIII. On lines and planes of closest fit to systems of points in space,  The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 2 (1901) 559–572. doi:10.1080/14786440109462720.
  * Schölkopf et al. [1997] B. Schölkopf, A. Smola, K.-R. Müller,  Kernel principal component analysis,  in: W. Gerstner, A. Germond, M. Hasler, J.-D. Nicoud (Eds.), Artificial Neural Networks — ICANN’97, Springer Berlin Heidelberg, Berlin, Heidelberg, 1997, pp. 583–588. doi:10.1007/BFb0020217.
  * Mika et al. [1999] S. Mika, B. Schölkopf, A. J. Smola, K.-R. Müller, M. Scholz, G. Rätsch,  Kernel PCA and de-noising in feature spaces,  in: M. J. Kearns, S. A. Solla, D. A. Cohn (Eds.), Advances in Neural Information Processing Systems 11, MIT Press, 1999, pp. 536–542. URL: https://papers.nips.cc/paper/1491-kernel-pca-and-de-noising-in-feature-spaces.
  * Abdi and Williams [2010] H. Abdi, L. J. Williams,  Principal component analysis,  Wiley Interdisciplinary Reviews: Computational Statistics 2 (2010) 433–459. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/wics.101. doi:10.1002/wics.101. arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/wics.101.
  * Mattis et al. [2024] S. Mattis, T. Roper, T. Butler, Cu-denver-uq/luq: v2.0.0, 2024\. URL: https://doi.org/10.5281/zenodo.10680749. doi:10.5281/zenodo.10680749.
  * Unser et al. [2017] M. Unser, J. Fageot, J. P. Ward,  Splines are universal solutions of linear inverse problems with generalized TV regularization,  SIAM Review 59 (2017) 769–793. URL: https://doi.org/10.1137/16M1061199. doi:10.1137/16M1061199. arXiv:https://doi.org/10.1137/16M1061199.
  * Hornik et al. [1989] K. Hornik, M. Stinchcombe, H. White,  Multilayer feedforward networks are universal approximators,  Neural Networks 2 (1989) 359–366. URL: https://www.sciencedirect.com/science/article/pii/0893608089900208. doi:https://doi.org/10.1016/0893-6080(89)90020-8.
  * Jiang et al. [2022] Q. Jiang, L. Zhu, C. Shu, V. Sekar,  An efficient multilayer rbf neural network and its application to regression problems,  Neural Computing and Applications 34 (2022) 4133–4150. URL: https://doi.org/10.1007/s00521-021-06373-0. doi:10.1007/s00521-021-06373-0.
  * Aronszajn [1950] N. Aronszajn,  Theory of reproducing kernels,  Transactions of the American Mathematical Society 68 (1950) 337–404. URL: http://www.jstor.org/stable/1990404. doi:10.2307/1990404, full publication date: May, 1950\.