# Levitation, oscillations, and wave propagation in a stratified fluid

Marina Carpineti1, Fabrizio Croccolo2, and Alberto Vailati1 1Dipartimento di
Fisica Aldo Pontremoli, Università degli Studi di Milano, I-20133 Milano,
Italy  2Universite de Pau et des Pays de l’Adour, E2S UPPA, CNRS, TOTAL, LFCR
UMR5150, Anglet, France<EMAIL_ADDRESS>
###### Abstract

We present an engaging levitation experiment that students can perform at home
or in a simple laboratory using everyday objects. A cork, modified to be
slightly denser than water, is placed in a jug containing tap water and coarse
kitchen salt delivered at the bottom without stirring. The salt gradually
diffuses and determines a stable density stratification of water, the bottom
layers being denser than the top ones. During the dissolution of salt, the
cork slowly rises at an increasing height, where at any instant its density is
balanced by that of the surrounding water. If the cork is gently pushed off
its temporary equilibrium position, it experiences a restoring force and
starts to oscillate. Students can perform many different measurements of the
phenomena involved and tackle non-trivial physical issues related to the
behaviour of a macroscopic body immersed in a stratified fluid. Despite its
simplicity, this experiment allows to introduce various theoretical concepts
of relevance for the physics of the atmosphere and stars and offers students
the opportunity of getting acquainted with a simple system that can serve as a
model to understand complex phenomena such as oscillations at the Brunt-
Väisälä frequency and the propagation of internal gravity waves in a
stratified medium.

##  1 Introduction

The levitation of macroscopic objects has fascinated humankind for a very long
time. Its history is intertwined with the research performed by several Nobel
laureates, because beyond its mystical aspects levitation can be profitably
used both as a tool in fundamental research and in technological applications
[1, 2]. Beyond its historical role in science, the investigation of levitation
is still full of puzzling results, such as the apparent reversal of the
gravity acceleration recently reported for a solid body at the surface of a
levitating liquid [3, 4].

The fascination determined by levitation comes from the fact that we rarely
experience it, as objects close to the surface of the Earth do not float
freely into the air, due to the strong gravitational attraction and the small
buoyancy force exerted by air onto solid bodies. When the gravitational force
is not balanced, the object is in free fall, and in the frame of reference of
the falling object, the acceleration of gravity is not felt. This determines
weightlessness levitation conditions that can be employed profitably to
perform experiments in the absence of gravity on platforms such as drop
towers, parabolic flights, and artificial satellites like the International
Space Station [5]. On Earth, the gravitational force acting on a body can be
also balanced by the buoyancy force acting onto it when the body is immersed
in a fluid medium. When the density of the body equals that of the surrounding
fluid, the weight of the body is exactly balanced by the buoyancy force and
the body is in equilibrium inside the fluid neither rising nor sinking, in a
condition similar to levitation. However, the phenomenon of levitation does
not involve simply the balance of forces, as occurs in buoyancy phenomena, but
the mechanical stability of the levitating body [6]. The particularity of a
levitation phenomenon is that restoring forces drive back the body when it is
moved away from its equilibrium position.

A levitation condition can occur in a stable stratified medium, where the
density decreases as a function of height. A body can therefore be in
equilibrium only at a certain height where its density is exactly matched and
a displacement in a direction parallel to the density gradient determines a
restoring force that brings back the object to the layer of fluid with the
same density. This force gives rise to damped harmonic oscillations that occur
at a frequency, known as Brunt-Väisälä frequency [7], which is determined by
the local density gradient present in the stratified fluid. The oscillations
gradually damp both through viscous dissipation and through the emission of
transverse internal gravity waves propagating in the plane perpendicular to
the direction of oscillation. Examples of natural systems where a density
stratification occurs include giant planets and stars [8], where the mass is
large enough to determine the confinement of a fluid phase. A typical example
is represented by the atmosphere of a planet, where the stratification of a
gaseous phase is due to the variation of pressure with altitude, which
determines a gradual decrease of density.

In this work, we discuss a simple experiment that can be realized by students
using kitchenware to investigate levitation, oscillations at the Brunt-Väisälä
frequency, and internal gravity waves in a stratified medium. The experiment
makes use of a modified cork that sinks into water but suddenly starts
levitating when a large enough saline gradient is generated by pouring coarse
salt inside water. Its simplicity makes this demonstration suitable for
students starting from elementary school up to undergraduate school, due to
the opportunity of achieving different levels of understanding and modelling
of the physical processes involved.

For undergraduate students, in particular, the experiment can be used to
convey the hydrostatic principles behind the balance of the forces involved.
Moreover, an external perturbation of the cork allows the investigation of the
physics of oscillations in a stratified medium and the radiative transfer of
the energy of oscillations determined by internal gravity waves.

The paper is organized as follows. The first paragraph reviews the
experimental contexts where levitation is observed in nature, in order to give
a general framework where the experiment may be inserted. In the following, we
present the experiment, discuss its relevance in an educational context, and
describe the experimental materials and methods. We then provide the basic
theoretical concepts needed to describe the hydrostatic equilibrium of a fluid
and the conditions that can give rise to either buoyancy or levitation of a
body embedded in it. We discuss the case of the oscillations of a body in
equilibrium in a stratified fluid and the internal gravity waves originated by
the dissipation of energy during the oscillations. Finally, we describe,
analyze and discuss the experimental results.

##  2 Levitation in science

Although it hardly occurs spontaneously in nature, levitation has been
profitably used in science to perform several challenging experiments. A
notable example of its usage in fundamental research is represented by the
celebrated experiment by Robert Millikan for the measurement of the charge of
the electron, where a negatively charged droplet of oil was levitated by using
an electrostatic field [9]. Millikan was awarded the Nobel prize in Physics in
1923 for this experimental demonstration of the discrete nature of the
electric charge [10]. Levitation can also be achieved by using a magnetic
field, as it happens for example by placing a superconducting disk above a
strong magnet [11, 12, 1]. Magnetic levitation has important technological
applications for the development of frictionless railways [13, 14]. A strong
magnetic field can also be used for challenging applications like the
levitation of a living being, as it has been done in a popular experiment by
Andre Geim and Michael Berry. By using a strong magnetic field, they were able
to levitate a living frog, and this result allowed them to win the Ig-Nobel
prize in Physics in the year 2000 [15]. A few years later, in 2010, Geim
shared the Nobel prize in Physics with Konstantin Novoselov for the discovery
of graphene [16]. Beyond the use of external fields, levitation can be also
achieved by using radiation [1]. The phenomenon of acoustic levitation was
studied by Rayleigh. Currently, acoustic levitation can be even achieved by
amateurs at home or in the classroom, by building affordable instrumentation
based on off-the-shelf components such as Tiny-Lev [17], a device that allows
levitating small objects with a maximum size of 4 mm and density below 2
g/cm3. The levitation of particles using light can be achieved by using
optical tweezers, developed by Arthur Ashkin [18], an invention that allowed
him to win the Nobel Prize in Physics in 2018 [19]. Optical tweezers allow the
non-invasive manipulation of nanoparticles by taking advantage of a focused
Laser beam that exerts a non-isotropic radiation pressure, resulting in the
trapping of a particle in a potential well. The introduction of optical
tweezers determined a revolution in the fields of Biophysics and Soft Matter,
where they allow the investigation of the behaviour of colloidal particles,
molecular motors, and cells, to mention just a few representative applications
[20, 21].

Neutral buoyancy represents a condition similar to levitation, very common in
the case of animals populating the waters of the Earth. Meaningful strategies
adopted by marine animals to achieve neutral buoyancy are represented by the
swim bladder of teleost fish, by the squalene oil contained in large
quantities into the liver of sharks, and by the large hollow chambers built by
the nautilus into its shell [22]. In the case of marine animals, neutral
buoyancy is achieved in a medium of nearly constant density, and this feature
determines an indifferent equilibrium condition, where the weight of the
animal can be balanced by buoyancy irrespectively of the position of the
animal. The lack of a unique equilibrium position and of a restoring force
indicates that neutral buoyancy cannot properly be considered a levitation
phenomenon.

At variance, the equilibrium condition realized in the presence of a stable
density stratification can be well described as levitation. In particular,
when the density gradient evolves in time, it is possible to observe that an
object, initially at rest, suddenly starts to rise under the action of
invisible forces. At any instant, its equilibrium is at a well-defined height
where it is subjected to restoring forces that make it oscillate in case it is
slightly displaced. Several macroscopic systems exhibit these conditions and
share the common feature that the mechanical energy associated with the local
oscillation can be transferred at a distance by internal gravity waves
generated by the oscillation itself. Although internal gravity waves in a
fluid cannot be visualized, the effect of the energy transferred by them at a
distance can be appreciated from the effect that it has on macroscopic bodies
embedded in the fluid. An important example is represented by the Earth’s
atmosphere, where mountain reliefs and thunderstorms can interact with gravity
waves to give rise to complex turbulent phenomena [7]. Another notable example
is represented by oceans, where the density stratification is both determined
by changes in salinity and temperature. Internal gravity waves have been
reported also at the mesoscale in laboratory experiments on non-equilibrium
fluids, where the displacement of parcels of fluids is determined by the
thermal energy $k_{B}T$ [23].

##  3 The experiment

The idea of this experiment came during the lockdown period due to the
COVID-19 pandemic when forced to stay far from our laboratories we were
stimulated to invent simple and catching experiments, easy to be performed at
home but able at the same time to teach physics. The experiment is ideally
presented to students through an inquiry based path, where the phenomenology
is first demonstrated to the students by a facilitator. Students are then
instructed about the details of the experiment so that they can perform it
autonomously. Eventually, the students are lectured about the physical
principles behind the experiment and about the methods for the analysis of
results.

###  3.1 Inquiry-driven presentation of the experiment

Levitation is a surprising, intriguing, and counter-intuitive theme that
offers the chance to discuss many physical phenomena. It is extremely
difficult to observe it in practice and all the experiments involving it
require a great technological effort. In this work, we describe a very simple
levitation experiment, where a weighed-down cork, initially located at the
bottom of a glass jug containing water, suddenly starts levitating inside the
fluid and gradually rises (Fig. 1).

Figure 1: Shots of the cork during its rise. The marked point and the
reference axis used for tracking are shown. Coordinates are fixed so that the
origin coincides with the position of the tracked point at the starting time.
For each figure, the time of acquisition from the beginning of the measurement
is shown.

Apparently, no source of energy is involved in the process, and this leaves
the students with a first scientific puzzle: i) where does the mechanical work
needed to raise the cork come from? Another striking aspect of this experiment
is that when the levitating cork is slightly displaced vertically from its
position, it undergoes damped harmonic oscillations (Fig. 2). This aspect
leaves the students with a second puzzling question: ii) what is the physical
origin of the restoring force that drives the cork towards its equilibrium
position during oscillations?

Figure 2: Shots of the cork during oscillation induced by a slight initial
push. A video is included as supplemental material (video1).

When two levitating corks are hosted inside the same container, pushing one of
the corks vertically determines, after a latency time, the onset of
oscillations in the second cork as well (Fig. 3). This feature raises a third
puzzling question for the students: iii) what is the origin of the interaction
between the two corks?

Figure 3: Two corks in a tank. After some oscillations of the right cork, also
the left cork starts to move although with an oscillation amplitude much lower
than the first cork. A video is included as supplemental material (video2).

###  3.2 Materials and methods

For the experiment we had in mind, we needed an easy-to-find object, with a
density slightly larger than water. The choice was to weigh-down the cork of a
wine bottle. The cork was cut in two halves digging a hole in each of them to
host a weight, for which we chose another item frequently available at home: a
steel marble with a diameter of $10\ mm$ borrowed from a toy. We soaked the
sphere in vinyl glue, pasted back the two halves of the cork (see Fig. 4), and
used ethyl cyanoacrylate glue along the junction. Finally, we cut thin slices
of the cork, until it had a density slightly larger than water and was able to
sink into it.

Figure 4: The picture on the left shows the two halves of the cork, with the
hemispherical dimples carved into them, both filled with vinyl glue and one of
them hosting the metal sphere. The picture on the right shows the assembled
cork. In this configuration, the cork density was still smaller than that of
water and it was necessary to cut thin slices of the cork to make it slightly
denser than water.

The experiment is prepared by pouring $500\ ml$ of tap water into a
cylindrical glass jug, and adding a defined amount of course kitchen salt
without applying any stirring, so that the salt grains rapidly sediment to the
bottom of the jug. We decided to use coarse rather than fine salt to avoid its
partial dissolution while settling at the bottom of the jug, and therefore
guarantee to start from better-defined initial conditions. The experiment is
started by rapidly putting the cork at the bottom before the salt has much
time to dissolve. The cork, initially sitting at the bottom of the jug, is
expected to start to levitate when the surrounding density becomes higher than
its own and to move towards a height where an iso-density condition is
reached. The process occurs because the salt undergoes a diffusion process
that determines a density stratification in the water: the salt-rich layers of
fluid at the bottom of the jug are initially denser than the salt-poor layers
at the top. Diffusion determines the spread of the salt ions across water,
until eventually they become uniformly distributed inside the available
volume.

Once the density of the cork is chosen, measurements can be performed with
different amounts of salt that determine the behaviour of the cork. We span
densities varying from 10 to 46 g/l.

In our experiment, for the smaller amounts used (typically 10-15 g/l), the
cork rises after some time but eventually goes back to the bottom of the jug,
because the equilibrium density (i.e. for infinite time) of the aqueous
solution is smaller than that of the cork. For larger amounts of salt, the
cork rises at the beginning and eventually reaches the top of the liquid,
finally floating at its surface. In the latter case, the equilibrium density
of the liquid is larger than that of the cork. Refining the quantity of salt
dissolved in the tap water around the values obtained for the two different
conditions, one can get a somewhat precise measurement of the density of the
cork. We could estimate its density to be approximately $1,030\ g/cm^{3}$,
starting from an estimation of the critical amount of salt of 40 g dissolved
in 1 kg of tap water, although the cork density can slightly change from one
experiment to another due to its moisture content.

The cork ascent lasts usually many hours and we recorded it using a Fujifilm
XT20 camera programmed to shoot pictures at regular time intervals, typically
every five minutes. The framing is chosen so as to capture the entire height
of the liquid. In order to study the oscillation process, we apply a punctual
acceleration to the cork and record its oscillations at 24 frames/second by
using the same camera. Should this experiment be proposed to students, even a
readily available smartphone could be used to make many of the measurements
and this was indeed the procedure followed at the beginning. However, when
available, an automated system is definitely preferable.

Figure 5: Picture of the starting time of a typical measurement. The cork lies
at the bottom of the jug where the coarse salt has settled.

Both cork’s rise and oscillations were characterized using the Open Source
software ”Tracker” [24], a Java-based application that allows to track the
different positions of an object in a video. To analyse the cork’s ascent, the
sequence of collected images has been first converted in an AVI video, using
the Fiji Open Source software [25], and then imported into Tracker, setting
the correct frame rate. Spatial calibration is done for each video by
assigning to the water height the real value previously measured with a ruler
($6.8\ cm$). To track the different positions of the cork, a mark is then
placed at a recognizable point of it, such as a corner or a small detail, and
the origin of the reference system is fixed on it. Due to a large number of
frames to analyse, the choice of a well-identifiable and non-ambiguous point
is extremely important to allow the use of the auto-tracking procedure, much
faster than the manual one. At the end of the process, the different positions
of the point as a function of time are obtained.

##  4 Theory review

In this section, we will provide the basic physical concepts needed for the
understanding of the hydrostatic equilibrium of a fluid, which can give rise
to the levitation of a body embedded into it and, consequently, to its
oscillations when it is gently pushed off its equilibrium position. The theory
of internal gravity waves will be also briefly outlined.

###  4.1 Hydrostatic equilibrium, buoyancy, and levitation

Let us consider a very general configuration suitable to describe fluids such
as the atmosphere and oceans: a layer of fluid of density $\rho$, distributed
in a spherical shell of thickness $L$, located ar a distance $R$ from the
centre of mass of a distribution of total mass M.

A parcel of fluid of volume V is under the action of the gravity force $F=\rho
Vg$ determined by the Law of Universal Gravitation $F=G\frac{mM}{R^{2}}$,
where the acceleration of gravity is given by $g(R)=G\frac{M}{R^{2}}$. At a
height $z$, the weight of the column of fluid above it determines a
hydrostatic pressure $p$ inside the fluid, governed by the equation:

$\nabla p=\rho\,\mathbf{g}$ (1)

In the case of an incompressible fluid like a liquid, the density is constant
and the hydrostatic pressure is governed by Stevino’s equation

$p(z)=p_{0}-\rho gz$ (2)

which is a particular solution of Eqn. 1, where $p_{0}$ is the pressure at the
surface of the fluid $z=0$.

A consequence of the hydrostatic pressure of Eqn. 1 is that the surface $S$ of
a body of volume $V_{0}$ completely immersed in the fluid is subjected to
forces acting perpendicularly to its surface. As the horizontal components are
balanced, the body is under the action of a buoyancy force $\mathbf{F}_{A}$
directed vertically, which can be obtained by integrating the hydrostatic
pressure provided by Eqn. 2 across the surface $S$ of the body:

${F}_{A}=-\left|\oint_{S}p(z)\mathop{}\\!\mathrm{d}\mathbf{S}\right|=\int_{V_{0}}\frac{dp}{dz}\mathop{}\\!\mathrm{d}V=\rho\,{g}\,V_{0}$
(3)

where $\mathop{}\\!\mathrm{d}\mathbf{S}$ is an element of the surface with
outgoing normal, and we have applied the divergence theorem to convert a
surface integral to a volume integral.

The total force acting on the body can be obtained by adding to the buoyancy
force the weight of the body $\rho_{0}\,V_{0}\,{g}$ with the proper sign:

${F}_{T}=\left(\rho-\rho_{0}\right)\,{g}\,V_{0}$ (4)

where $\rho_{0}$ is the density of the body. An immediate consequence of Eqn.
4 is that when the density of the body is matched with that of the fluid,
$\mathbf{F}_{T}=0$. Under these conditions, the body is in an indifferent
equilibrium condition called neutral buoyancy, i.e the equilibrium condition
is achieved at each point of the fluid. Such condition is achieved by several
aquatic organisms, which have devised a series of strategies to compensate for
the variation of pressure as a function of depth [22]. A major limitation of
neutral buoyancy is that it requires a very close matching between the density
of the fluid and that of the body. Even a small difference between these
quantities would give rise to a condition where the body is under the action
of a constant acceleration, which brings it progressively farther from its
initial position.

A stable equilibrium condition can be achieved by immersing the body in a
stratified compressible fluid, like for example a gas stratified by the
gravity force. Under this condition, the body, after a transient, migrates to
the layer of fluid matching its density. This condition is difficult to obtain
in practice because the average density of most solid bodies is much larger
than that of a gas. Matching the density of the solid would require the use of
a liquid, which is intrinsically not compressible and, therefore, cannot be
stratified by the gravity force.

Notwithstanding that, a gravitationally stable density stratification can be
obtained in a liquid under non equilibrium conditions. In fact, if the liquid
is under the action of a vertical temperature gradient $\nabla T$ its thermal
expansion determines a density gradient
$\nabla\rho=\frac{\partial\rho}{\partial T}\nabla T$. Moreover, if the liquid
is a mixture of two components, a variation of the concentration $c$ with
height determines a density gradient $\nabla\rho=\frac{\partial\rho}{\partial
c}\nabla c$. In general, combining these two terms, the overall vertical
density gradient is given by:

$\nabla\rho=\rho\left(\alpha\nabla T+\beta\nabla c\right)$ (5)

where $\alpha$ is the thermal expansion coefficient, and $\beta$ is the
solutal expansion coefficient.

Under the condition that the overall density gradient points downwards, the
density profile is gravitationally stable, and a body of density $\rho_{0}$
such as $\rho_{min}<\rho_{0}<\rho_{max}$ can levitate in a layer of a
stratified fluid.

###  4.2 Oscillations in a stratified fluid

Let us now consider the case where a body of density $\rho_{0}$ and volume
$V_{0}$ is immersed in a layer, with the same density, of a gravitationally
stable stratified fluid. If the volume $V_{0}$ undergoes a small displacement
$z$ in the vertical direction, the density gradient can be assumed to be
locally constant so that the density changes linearly with height:

$\rho(z)=\rho_{0}+\frac{\partial\rho}{\partial z}\,z$ (6)

Combining Eqns. 1, 3, and 6, one can calculate the buoyancy force acting on
volume $V_{0}$:

${F}_{B}=-\oint_{S}p(z)\mathop{}\\!\mathrm{d}\mathbf{S}=\int_{V_{0}}\frac{dp}{dz}\mathop{}\\!\mathrm{d}V=\rho_{0}\,{g}\,V_{0}+\frac{\partial\rho}{\partial
z}\,{g}\,V_{0}\,z$ (7)

and recalling Eqn. 4 the total force acting on $V_{0}$ is:

${F}_{T}=\frac{\partial\rho}{\partial z}\,{g}\,V_{0}\,z$ (8)

By applying the second principle of dynamics
$\rho_{0}\,V_{0}\,\ddot{z}=F_{T}$, one obtains the equation of motion of a
harmonic oscillator (notice that $\frac{\partial\rho}{\partial z}$ is
negative) :

$\ddot{z}=\frac{g}{\rho_{0}}\frac{\partial\rho}{\partial z}\,z$ (9)

where oscillations occur at the Brunt-Väisälä frequency [7]:

$N=\sqrt{-\frac{g}{\rho_{0}}\frac{\partial\rho}{\partial z}}$ (10)

Oscillations at the Brunt-Väisälä frequency also occur when a parcel of fluid
is displaced vertically with respect to the surrounding fluid. This process
occurs in natural systems such as the atmosphere and the oceans and in
astrophysical systems such as stars. Interestingly, the same physical
mechanism drives oscillations at mesoscopic length scales during a thermal
diffusion process occurring in a mixture of fluids [23]. In that case, the
density stratification stems from temperature and concentration gradient, thus
resulting in an additional contribution related to the temperature gradient.
Anyway, the equation can be written in the same form of Eq.10 of the present
paper if one considers the dependence from the density gradient.

###  4.3 Internal gravity waves

When the volume $V_{0}$ oscillates, its mechanical energy is gradually
transferred to the surrounding fluid. Part of the energy undergoes viscous
dissipation, while another part gets transferred at a distance by the layered
fluid in the form of an internal gravity wave. Internal gravity waves are
qualitatively similar to the surface waves in a liquid. However, their density
gradient is, of course, much smaller than that of surface waves. In the
following we will describe the emission of transverse waves by a body that
undergoes vertical harmonic oscillations at the Brunt-Väisälä frequency in a
stratified fluid . For a more rigorous and general treatment of internal
gravity waves see Tritton [26] and Nappo [7].

We consider an incompressible inviscid fluid in the presence of a stable
density stratification described by Eqn. 6. The Eulerian hydrodynamic
equations for a fluid of velocity $\mathbf{u}$, pressure $p^{\prime}$ and
density $\rho^{\prime}$ are:

$\displaystyle\nabla\cdot\mathbf{u}=0$ (11)
$\displaystyle\rho^{\prime}\,\frac{\partial\mathbf{u}}{\partial
t}+\rho^{\prime}\,\mathbf{u}\cdot\nabla\mathbf{u}=-\nabla
p^{\prime}+\rho^{\prime}\,\mathbf{g}$ (12)
$\displaystyle\frac{\partial\rho^{\prime}}{\partial
t}+\mathbf{u}\cdot\nabla\rho^{\prime}=0$ (13)

Equation 11 is the continuity equation and expresses the conservation of mass,
while Eqn. 12 is the second law of dynamics and Eqn. 13 expresses the fact
that the fluid is not compressible.

We assume that the variables for the fluid in motion with velocity
$\mathbf{u}$ can be written as a superposition of an hydrostatic contribution
and a small perturbation:

$\displaystyle p^{\prime}(x,y,z,t)=p(z)+\delta p(x,y,z,t)$ (14)
$\displaystyle\rho^{\prime}(x,y,z,t)=\rho(z)+\delta\rho(x,y,z,t)$ (15)

where $p$ and $\rho$ obey the hydrostatic Eqns. $\nabla p=\rho g$ and
$\partial\rho/\partial t=0$, respectively. By inserting Eqns. 14 and 15 into
Eqns 11-13, we obtain:

$\displaystyle\nabla\cdot\mathbf{u}=0$ (16)
$\displaystyle\rho\,\frac{\partial\mathbf{u}}{\partial
t}+\delta\rho\,\frac{\partial\mathbf{u}}{\partial
t}+\rho\,\mathbf{u}\cdot\nabla\mathbf{u}+\delta\rho\,\mathbf{u}\cdot\nabla\mathbf{u}=-\nabla\delta
p+\delta\rho\,\mathbf{g}$ (17)
$\displaystyle\frac{\partial\delta\rho}{\partial
t}+\mathbf{u}\cdot\nabla\rho+\mathbf{u}\cdot\nabla\delta\rho=0$ (18)

We now assume that the perturbations of density and pressure are small with
respect to the hydrostatic quantities, $|\delta p/p|\ll 1$ and
$|\delta\rho/\rho|\ll 1$, and linearize the equations in $\mathbf{u}$, $\delta
p$, and $\delta\rho$:

$\displaystyle\nabla\cdot\mathbf{u}=0$ (19)
$\displaystyle\rho\,\frac{\partial\mathbf{u}}{\partial t}=-\nabla\delta
p+\delta\rho\,\mathbf{g}$ (20)
$\displaystyle\frac{\partial\delta\rho}{\partial
t}+\mathbf{u}\cdot\nabla\rho=0$ (21)

We look for solutions of the equations in the form of propagating harmonic
waves of wave vector $\mathbf{k}$ and frequency $\omega$:

$\displaystyle\mathbf{u}=\mathbf{U}\,exp\left[i\left(\mathbf{k}\cdot\mathbf{x}-\omega
t\right)\right]$ (22)
$\displaystyle\delta\rho=R\,exp\left[i\left(\mathbf{k}\cdot\mathbf{x}-\omega
t\right)\right]$ (23) $\displaystyle\delta
p=P\,exp\left[i\left(\mathbf{k}\cdot\mathbf{x}-\omega t\right)\right]$ (24)

For simplicity, we restrict ourselves to a two-dimensional system, of
coordinates $x$ (horizontal), and $z$ (vertical). We assume that the parcel of
fluid is located in the origin of the coordinate system at time $t=0$, and
undergoes harmonic motion in the vertical direction. This implies that the
amplitude $\mathbf{U}$ of the oscillation in Eqn. 22 must be directed in the
vertical direction, that is $\mathbf{U}=(0\ ;U)$.

By imposing that Eqn. 22 must be a solution of the continuity Eqn. 19 we
obtain that the velocity must be perpendicular to the wave vector:

$\mathbf{U}\cdot\mathbf{k}=0\;\Rightarrow k_{z}=0$ (25)

Therefore, the velocity wave defined by Eqn. 22 is transverse, and its wave
vector is directed horizontally, $\mathbf{k}=(k;0)$.

By imposing that Eqns. 23 and 24 are solutions of Eqns. 20 and 21 we get:

$\displaystyle k\,P=0\Rightarrow P=0$ (26) $\displaystyle
i\,\omega\,\rho\,U=g\,R$ (27) $\displaystyle
i\,\omega\,R=\frac{\mathop{}\\!\mathrm{d}\rho}{\mathop{}\\!\mathrm{d}z}\,U$
(28)

Equation 26 implies that the pressure does not propagate and the dissipation
of energy does not occur through the emission of a sound wave. Multiplication
of Eqn. 27 and 28 allows to determine the frequency of the internal gravity
waves, which coincides with the Brunt-Väisälä frequency of the oscillating
parcel of fluid:

$\omega=\sqrt{-\frac{g}{\rho_{0}}\frac{\partial\rho}{\partial z}}=N$ (29)

Notice that, although this theoretical model confirms that the harmonic motion
of a parcel of fluid in the vertical direction gives rise to transverse
internal gravity waves propagating horizontally, the wave number and the
velocity of the wave are not predicted by it, and in principle, there are no
constraints on the values that $k$ and the velocity can assume.

##  5 Results and discussion

We performed various experiments using different amounts of salt, ranging
between 10 and 46 g/l. In each measurement, we tracked the cork’s rise until
its position was approximately at half of the water height. Then, we took
movies of the oscillations in order to measure the Brunt-Väisälä frequency.

In Figure 1 we show a sequence of images taken when the cork had a low
density. Thanks to the relatively large speed of the process – it lasted
approximately three hours – it was possible to follow the rise until the cork
reached the water surface. Both the marked point and the reference system used
with Tracker are shown. In Figure 6 the corresponding plot of the position of
the cork versus time is shown. In the beginning, shots were taken every five
minutes but, due to the high rise speed, the time intervals were soon reduced
to one minute. The trend shown in Fig. 6 exhibits the same features observed
in all cases, also when the cork density was increased. There is an initial
latency time during which the cork starts to randomly slip at the bottom of
the jug, then the rise starts and the height increases almost linearly in
time. What is extremely different from one measure to another, and not easily
predictable, are the typical latency times, which do not seem to depend
trivially on the salt amount poured in water. There are probably multiple
reasons for this behaviour. One is that the cork’s density can slightly vary
from one measure to the other.

Figure 6: Cork’s rise obtained by tracking the position of the marked point
shown in Fig. 1.

It is also very difficult to guarantee that the salt fills uniformly the
bottom of the jug. Moreover, the finite size of the cork is not negligible
with respect to the size of the jug. Therefore, the presence of the cork
itself perturbs salt dissolution and therefore the typical times of the
process. Finally, the temperature of the water is not controlled at all,
therefore a change in the transport properties of the mixture is to be
expected from one measure to the other.

It is interesting to focus on the energy balance of the system in order to
understand the relationship between the different mechanisms acting on the
fluid and the cork. If one writes down the total mechanical energy of the cork
as the sum of its gravitational potential energy plus its kinetic energy, when
the cork is lifted from the bottom of the jug, the conservation of mechanical
energy is obviously not satisfied, since the cork potential energy increases
without any reason. In order to solve this apparent energy mystery, one should
seek the origin of the force lifting the cork. The cork is subjected to the
weight force, which is constant in time, and to the buoyancy determined by the
surrounding liquid, which changes in time. The work needed to raise the cork
is provided by the gradual increase of the density of the liquid surrounding
the cork, determined by the diffusive dissolution of the salt. Therefore, the
cork rises because the center of mass of the fluid surrounding it gradually
shifts upwards. This leaves us with another open question: what is the origin
of the mechanical energy needed to lift the center of mass of the fluid? The
dissolution of salt provides an energy contribution to the system in terms of
enthalpy and is driven by the thermal agitation of the component molecules
$k_{B}T$. The energy is, in the end, taken from the thermal energy of the
system, so that an accurate measurement of the fluid temperature should reveal
a decrease of the temperature during the rise of the salt and the cork. This
would, of course, be possible only in an ideal adiabatic system not exchanging
energy with other systems. In the experiment described here, the jug and the
water are actually in thermal equilibrium with the surrounding environment,
and the temperature increase is then not measurable since the environment acts
as the ultimate thermal reservoir. Nevertheless, it is interesting to observe
that in the initial condition, with the salt and the cork at the bottom of the
jug, the potential energy of the system is at a minimum, while after some
time, both the salt and the cork have been lifted increasing their potential
energy at the expense of the enthalpy of the salt solution and ultimately to
thermal energy. This is one of the cases of conversion of thermal energy to
mechanical energy.

At any time during its rise, the cork floats in a stratified medium and, if
set in motion with a slight tap, it starts to oscillate at the Brunt-Väisälä
frequency that is directly related to the local density gradient and gravity
as shown in Eqn. 10. In Fig. 2, we show a sequence of shots of the cork that
oscillates after a slight tap, given when it had reached approximately half
height within the liquid. At a visual inspection, it was evident that the
regular oscillations of the cork were periodically slowed down until the cork
almost stopped, to start again to oscillate after a while. A video is included
as supplemental material (video1). In Fig. 7 we show a sequence of tracking of
the oscillations taken at different times during the cork rise for a slow
process where we used the cork with its final density. The periodical regular
slowing down of the cork that can be observed in almost all the figures
follows a complex dynamics, which cannot be trivially related to that of a
damped harmonic oscillator. However, by fitting the first oscillations with
the formula for a damped harmonic oscillator

$y(t)=A*cos(\omega*t)*exp(-(t/\tau))+y_{0}$ (30)

one gets a typical oscillation frequency $\omega$ in the range from 2.67 to
3.45 $rad/s$ and a damping constant $\tau$ from 2.38 to 4.87. The values of
the frequency resulting from this simple fit are reported in Fig. 12.

Figure 7: Tracking of the cork oscillation taken at different times of the
same measure with 30 g/l NaCl.

A more accurate scrutiny of Fig. 7 reveals that the dynamics of the
oscillations is similar to that observed for a beating phenomenon that occurs
when two waves of slightly different frequencies interfere. As discussed at
the end of Section 4.3, the oscillations due to the gradient of density are
expected to damp emitting transverse internal gravity waves that propagate
perpendicularly to the direction of oscillation. Due to the finite size of the
container, we expect that these waves are reflected back by the jug walls and
finally interfere with the cork oscillation.

To test this hypothesis we performed a further experiment using a rectangular
plexiglass tank as a container. We prepared another cork similar to the first
one and immersed both of them in water and salt. As soon as the corks had
reached a proper height, we slightly pushed one of them and verified that
after a while also the other one started to oscillate. This is a piece of
evidence that a signal propagated horizontally from the first cork to the
second, as expected for internal gravity waves. In Figure 3 an image of the
two corks is shown, while a video is included as supplemental material
(video2). Gravity waves in a transparent fluid are almost invisible, and this
feature makes them hard to investigate, but some of their features can be
determined from their effect on bodies immersed in the fluid, such as the two
corks that we used. A remarkable example is represented by the velocity of the
wave, which can be obtained from the ratio $\frac{d_{c}}{t_{m}}=v=6,6\pm 1,4\
cm/s$, where $d_{c}=8,3\pm 0,5\ cm$ is the distance between the two corks,
measured from the images, and $t_{m}\simeq 1,25\ \pm 0,2s$ is the time
interval between the push on the first cork and the time at which the second
one starts to oscillate itself. A more thorough investigation of the process
requires the possibility of visualizing the wave. This is made possible by the
fact that the waves occur through density variations, which are associated to
changes of the index of refraction that could be visualized by observing the
deformations of a reference object placed behind the liquid. In this
experiment, we used a chequered sheet placed in contact with the back window
of the tank. To analyse the gravity waves, we processed the images of the two
corks with Fiji by expanding the images in the vertical direction by a factor
of 16 to achieve better visualization of the horizontal lines of the chequered
reference pattern. Images were then converted to grayscale and a gamma
correction of exponent 3 was applied to increase the contrast of the
horizontal lines. A video is available as supplemental material (video3). In
Fig. 8 we show a sequence of processed images collected at different times
after the start of oscillations of the first cork. The deformation of the
horizontal line can be clearly observed and used to calculate the speed and
frequency of the internal gravity wave. The positions of the wave peaks were
followed in time with the Tracker program, and in Fig. 9, a plot of the
displacement along the x-direction vs. time for one of them is shown. The
continuous straight line is the best linear fit of the experimental points and
its slope gives an estimate of the wave speed $v_{w}\approx 6,0\pm 0,1cm/s$.

Figure 8: Sequence of images taken with a frame rate of 8 img/s at different
times of a portion of the sheet of paper placed behind the tank between the
two corks. Time increases from left to right and from top to bottom. The
horizontal line is deformed due to the variation of the refractive index
produced by the transit of the internal gravity wave. A video is available as
supplemental material (video3)

We have tracked three different peaks spanning the space between the two corks
obtaining an average value for $v_{w}$, $v_{av}=6,7\pm 0,6\ cm/s$.

We also measured the frequency of the internal gravity wave by using the line
deformations along the vertical direction, induced by the index of refraction
inhomogeneities shown in Fig. 8. Using the Tracker program, fixed points of
the line were chosen and followed in time during the passage of the wave to
measure their displacement along the y-axis. A typical result is shown in Fig.
10.

Figure 9: Peak position of the wave detectable as a deformation of the line
shown in Fig. 8. The solid line is the best linear fit of the experimental
points and its slope gives an estimate of the wave speed that, for this curve
is $6,0\pm 0,1\ cm/s$. Figure 10: Tracking of the position of a point on the
line shown in Fig. 8 during the passage of the gravitational wave. The
continous line is the fit obtained using Eqn. 30.

The continuous line is the best fit of the experimental points with Eq. 30. By
averaging the oscillation frequencies obtained by tracking different points of
the line, we obtained an average value for $\omega=4.4\pm 0.4\ rad/s$ that
corresponds to a frequency $f=0,64\pm 0.06Hz$ for the internal gravity wave.
This value is close to the oscillation frequency of the cork pushed in motion,
$\omega_{c}=4.26\ rad/s$, as measured by tracking the cork position. This
result confirms that the internal gravity wave has the same frequency as the
Brunt-Väisälä wave that generates it, as shown in Eqn. 29, although it must be
stressed that the image quality and the strong damping of the waves make it
difficult to obtain a precise estimation of the frequencies.

To have better quality measurements of Brunt-Väisälä frequencies in different
conditions, we finally fitted the individual cork’s oscillations shown in
Figs. 2 and 7. In order to fit the entire curves and not only the first
oscillations, we used a phenomenological fitting function more refined with
respect to Eq. 30 that takes into account both the oscillation damping and the
beatings between the cork oscillation and the back-reflected internal gravity
wave:

$y(t)=A*cos(\omega_{1}*t+\phi_{1})*cos(\Delta\omega*t+\phi_{2})*exp(-(t/\tau)^{\gamma})+y_{0}$
(31)

where $\omega_{1}$, $\Delta\omega$, $\phi_{1}$, $\phi_{2}$, $\tau$, $\gamma$
and $y_{0}$ are the free parameters. Eq. 31 is the product of two cosine
functions which typically describes beatings that occur when two waves of
similar frequencies interfere. In beating phenomena, $\omega_{1}$ is the
average of the frequencies of the interfering waves, while $\Delta\omega$ is
their difference. The stretched exponential term $exp(-(t/\tau)^{\gamma})$
with $\gamma=2$ turned out to be the best function to account for the severe
damping of the two oscillating terms that have different time constants. In
Fig. 11 we show the remarkable result of the fit. It can be noticed that while
Eqn. 30 is able to fit only the first oscillations, Eqn 31 is able to account
for all the oscillations features.

Figure 11: Fit of the oscillation tracking for a measurement performed with 30
g/l NaCl. Data are taken 4 hours after the beginning of data acquisition, when
the cork was almost 1,5 cm from the bottom of the jug.

As shown in Fig. 12, for all the data shown in Fig. 7, $\Delta\omega$ is much
smaller than $\omega_{1}$ and this result confirms that the two beating waves
have very close frequencies. As a consequence, $\omega_{1}$ gives a good
estimation of the Brunt-Väisälä frequency and turns out to be in very good
agreement with the frequency $\omega$ that is obtained by fitting the data
with the much simpler Eqn. 30. In Fig. 12, it can be appreciated that the
frequency of oscillation decreases in time as expected because the density
gradient diminishes while the concentration approaches the final equilibrium
value.

Figure 12:  The plot shows the values of $\omega_{1}$ (blue squares) and
$\Delta\omega$ (red squares), as obtained by fitting the data in Fig. 7 with
Eqn. 31. Open squares are the values of $\omega$ obtained by fitting the first
oscillations of the same data with Eqn. 30. Both $\omega$ and $\omega_{1}$
give a good estimation of the Brunt-Väisälä frequency.

From the Brunt-Väisälä frequency, it is possible to calculate the value of the
density gradient $\frac{\partial\rho}{\partial z}$ from Eqn. 10. Having fitted
all the oscillations measured with all the different salt concentrations, we
found for the frequency $N=\omega$ values between $1,5$ and $4,5$ rad/s
depending on the initial salt concentration and the hight of the cork in the
fluid and on the time of the measurement. Using Eqn. 10 the corresponding
values for $\frac{\partial\rho}{\partial z}$ ranges from $236\ kg/m^{4}$ to
$2140\ kg/m^{4}$ corresponding to concentration variations $\frac{\partial
c}{\partial z}=\frac{\partial\rho}{\partial z}/\frac{\partial\rho}{\partial
c}$ ranging from $0.32\ m^{-1}$ to $3.0\ m^{-1}$, where for a solution of NaCl
in water at small concentration $\frac{\partial\rho}{\partial c}=720kg/m^{4}$
[27].

##  6 Conclusions

We have discussed an engaging, inexpensive, and simple experiment that allows
introducing undergraduate students to many arguments on hydrostatic, buoyancy,
and levitation in a stratified fluid. The general scientific framework of
levitation, and the theoretical description of the phenomena that are
experimentally accessible, are also presented in this work. Brunt-Väisälä
oscillations are observed and the measurement of their frequency may be used
to estimate the concentration gradient inside the fluid. Moreover, the effect
of gravity waves can be observed and measured. We propose a phenomenological
function to fit the data, able to keep into account both the damping of the
oscillations and the beatings with the internal gravity wave they generate.
The motion of the cork determines a partial remixing of the stratified fluid,
thus reducing the concentration gradient that drives the oscillations. In this
respect, the experiment represents a remarkable example of a classical system
where performing a measurement significantly affects the state of the system.
This feature makes a long series of repeated measurements not useful, in
contrast with the usual physical systems usually encountered by students
during laboratory activities. The experiment is exploitable by students of any
age as the level of deepening can be tuned as a function of the students’
knowledge. However, it is particularly suitable for undergraduate students
that can appreciate all its implications.

## Acknowledgements

Work partially supported by the European Space Agency, CORA-MAP TechNES
Contract No. 4000128933/19/NL/PG. This research was carried under the
framework of the E2S UPPA Hub Newpores and Industrial Chair CO2ES supported by
the ‘Investissements d’Avenir’ French programme managed by ANR
(ANR–16–IDEX–0002).

## References

## References

  * [1] Brandt E H 1989 Science 243 349
  * [2] N Callens et al 2011 Microgravity Sci. and Tech. 23 181–189
  * [3] Sorokin V and Blekhman I I 2020 Nature 585 31–32
  * [4] Apffel B, Novkoski F, Eddi A and Fort E 2020 Nature 585 48–52
  * [5] M Braibanti et al 2019 Eur. Phys. J. E 42 86
  * [6] Jones T B, Washizu M and Gans R 1997 J. Appl. Phys. 82 883
  * [7] Nappo C J 2012 An introduction to Atmospheric Gravity Waves (Amsterdam: Elsevier)
  * [8] M Le Bars D L 2020 Fluid Mechanics of Planets and Stars (Switzerland: Springer Nature)
  * [9] Millikan R A 1913 Phys. Rev. 2 109
  * [10] Millikan R A Nobel lecture https://www.nobelprize.org/prizes/physics/1923/millikan/lecture/ accessed: 26 Jan 2021
  * [11] Landau L D and Lifshitz E M 1960 Electrodynamics of Continuous Media (Oxford: Pergamon)
  * [12] Saslow W M 1991 Am. J. Phys. 59 16
  * [13] Powell J and Danby G 1971 Cryogenics 11 192
  * [14] Coffey H T 1994 In: Status of Maglev: Opportunities in Cryogenics and Superconductivity (Boston, MA: Kittel P. (eds))
  * [15] Berry M and Geim A 1997 Eur. J. Phys. 18 307
  * [16] Geim A Nobel lecture https://www.nobelprize.org/prizes/physics/2010/geim/lecture/ accessed: 25 Jan 2021
  * [17] Marzo A, Barnes A and Drinkwater B W 2017 Rev. Sci. Instrum. 88 085105
  * [18] Ashkin A 1970 Phys. Rev. Lett. 24 156
  * [19] Ashkin A Nobel lecture https://www.nobelprize.org/prizes/physics/2018/ashkin/lecture/ accessed: 25 Jan 2021
  * [20] Moffitt J R, Chemla Y R, Smith S B and Bustamante C 2008 Annu. Rev. Biochem. 77 205
  * [21] McGloin D and Reid J P 2010 Opt. Photonics News 21 20
  * [22] Vogel S 2013 Comparative Biomechanics ISBN 9780691155661
  * [23] Croccolo F, García-Fernández L, Bataller H, Vailati A and Ortiz de Zárate J M 2019 Phys. Rev. E 99 012602
  * [24] Tracker app: https://physlets.org/tracker/ accessed: 2020-12-29
  * [25] Fiji app: https://imagej.net/Fiji
  * [26] Tritton D J 1988 Physical Fluid Dynamics (Oxford: Clarendon Press)
  * [27] Lide D R 2005 CRC Handbook of Chemistry and Physics, 86th Edition


# Bayesian meta learning for trustworthy uncertainty quantification

Zhenyuan Yuan1 Thinh T. Doan1 1Zhenyuan Yuan and Thinh T. Doan are with
Bradley Department of Electrical and Computer Science, Virginia Tech,
Blacksburg, VA 24061, USA (email: {zyuan18,thinhdoan}@vt.edu). This work was
supported by the National Science Foundation under ECCS-CAREER Grant No.
2339509.

###### Abstract

We consider the problem of Bayesian regression with trustworthy uncertainty
quantification. We define that the uncertainty quantification is trustworthy
if the ground truth can be captured by intervals dependent on the predictive
distributions with a pre-specified probability. Furthermore, we propose,
Trust-Bayes, a novel optimization framework for Bayesian meta learning which
is cognizant of trustworthy uncertainty quantification without explicit
assumptions on the prior model/distribution of the functions. We characterize
the lower bounds of the probabilities of the ground truth being captured by
the specified intervals and analyze the sample complexity with respect to the
feasible probability for trustworthy uncertainty quantification. Monte Carlo
simulation of a case study using Gaussian process regression is conducted for
verification and comparison with the Meta-prior algorithm.

##  I Introduction

Engineering systems operating in the real world are usually subject to unknown
uncertainties. Examples include autonomous cars driving in urban scenarios,
unmanned aerial robots for outdoor package delivery and mobile robots for
search-and-rescue. In order to ensure mission success and the safety of these
systems while maintaining high autonomy, it is necessary to learn and quantify
these uncertainties in a trustworthy manner such that these systems can
operate with minimal human supervision/intervention.

Bayesian learning [1] is a class of statistical learning frameworks, including
but not limited to Gaussian process regression (GPR) [2], Kalman filtering
[3], Bayesian neural network [4] and particle filtering [5]. In general,
Bayesian learning first models a target (e.g., system state, parameter or
function) as a sample from a distribution a priori, then given a set of data
it utilizes the Bayesian inference framework to compute a posterior
distribution of the target for prediction. With proper choice of the prior
distribution and mild assumptions on the target, Bayesian learning is able to
consistently approximate the target [6][7]. Established analysis, such as the
well-known PAC-Bayesian theorems [8], has shown that the generalization error
for the performances of Bayesian learning methods decreases at the rate of
$\mathcal{O}(\frac{1}{\sqrt{n}})$, where $n$ is the number of data samples.
Furthermore, the predictive distributions, i.e., the prior and the posterior
distributions, inherently allow Bayesian models to predict with uncertainty
quantification for each input. These aforementioned advantages make Bayesian
learning a powerful tool in a variety of applications, e.g., optimization [9],
learning-based control [10, 11, 12], motion planning[13][14], state estimation
[3] and system identification [15].

Based on the Bayes rule [1], the posterior distribution reflects predictive
uncertainties accurately only if the likelihood function and the prior
distribution are correctly specified, which are also a common assumption in
the analysis for uncertainty quantification [16][17][18]. However, in some
cases, the aforementioned information may not be obtained a priori accurately.
To relax this assumption, existing works instead assume the likelihood and/or
the prior distribution are specified structurally such that the
hyperparameters can be learned through data [2][19][20]. While there are
powerful function approximation models, such as deep neural networks and
finite-dimensional basis functions, which are able to consistently approximate
a wide range of possible models, in some cases it can be hard to ensure the
selected class of approximation models fully capture the target when there is
not much prior information about the structure of the target. Furthermore,
when the data is scarce, the hyperparameters may not be well learned [21]. As
a result, the uncertainty quantification obtained from the predictive
distributions may not be trustworthy, and the subsequent operations, e.g.,
synthesis of safety controller [22][23] and selection of safe decisions
[24][25], leveraging the uncertainty quantification can be unreliable.

Meta learning is a machine learning framework which aims to utilize the data
of a collection of tasks to identify a good initialization/prior for the
learning algorithm in new tasks such that fast learning can be achieved using
a small amount of data [26]. The process of identifying a good
initialization/prior is known as the meta training procedure, and the learning
in a new task is known as the adaptation procedure. The problem is usually
formulated as an optimization problem, where the objective function is the
expected performance of the adapted model in a new task. Most of the works do
not consider uncertainty quantification [26, 27, 28, 29]. Generalization
errors, which can be used to derive for uniform uncertainty quantification,
are considered in [30][31] for adapted models. Bayesian meta learning is
considered in [32, 33, 34, 35, 36, 37, 38], where the posterior distributions
provide input-dependent uncertainty quantification. In particular, papers
[32][33] consider meta learning for Bayesian neural networks, where the prior
distributions of the parameters in the neural networks are meta trained to
optimize the performances of the neural networks with parameters sampled from
the posterior distributions. Papers [34, 35, 36, 37, 38] consider Gaussian
process or Bayesian linear regression, where the hyperparameters in the prior
covariance and/or mean (functions) are meta trained provided pre-specified
structures. Prediction accuracy is considered in these methods, however,
whether the learned prior and posterior provide trustworthy uncertainty
quantification remains an open question.

Contribution statement. In this paper, we consider the problem of Bayesian
regression with trustworthy uncertainty quantification. We propose a Bayesian
meta learning framework which is cognizant of trustworthy uncertainty
quantification without explicit prior assumptions on the model/distribution of
the functions. Specifically, we define that the uncertainty quantification is
trustworthy if the ground truth can be captured by intervals dependent on the
predictive distributions with a pre-specified probability. We then propose
Trust-Bayes, a novel optimization framework for Bayesian meta learning with
constraints on trustworthy uncertainty quantification using the meta-trained
prior distribution and the posterior distribution. We characterize the lower
bounds of the probabilities of the ground truth being captured by the
specified intervals in terms of the empirical estimates from meta training. We
further analyze the sample complexity with respect to the feasible pre-
specified probability for trustworthy uncertainty quantification. In summary,
our major contributions are threefold:

  * •

We mathematically formulate trustworthy uncertainty quantification for
Bayesian regression.

  * •

We propose Trust-Bayes, a novel optimization framework for trustworthy
uncertainty quantification.

  * •

We characterize the lower bounds of the probabilities of the ground truth
being captured by the specified intervals and analyze the sample complexity
with respect to the feasible probability for trustworthy uncertainty
quantification.

We conduct Monte Carlo simulation and consider a case study using GPR for
verification of trustworthy uncertainty quantification by Trust-Bayes and for
comparison against Meta-prior [34, 35, 36] for its necessary.

Notation. Let $P_{x}\Big{(}\mathcal{E}\Big{)}$ return the probability of event
$\mathcal{E}$ with respect to the probability measure of $x$ and
$\mathbb{E}_{x}[\cdot]$ return the expected value with respect to the
probability measure of $x$. Define indicator function $\bm{1}[\mathcal{E}]=1$
if event $\mathcal{E}$ is true and $\bm{1}[\mathcal{E}]=0$ otherwise.

##  II Problem formulation

Observation model. Consider a distribution of functions $\mathcal{P}_{f}$,
where each unknown function $f^{i}\sim\mathcal{P}_{f}$,
$f^{i}:\mathcal{R}^{n_{x}}\to\mathbb{R}$, can be observed as

$\displaystyle y^{i}_{t}=f^{i}(x^{i}_{t}),$ (1)

where $x^{i}_{t}\in\mathcal{X}\subset\mathcal{R}^{n_{x}}$ is the $t$th input
to function $f^{i}$, $n_{x}$ is the dimensionality of $x_{t}$, and
$y^{i}_{t}\in\mathbb{R}$ is the corresponding output. In this paper, we
consider noiseless observation or the noise is inherent in $f^{i}$, which is a
typical model for system identification [17][39][40]. Here, $f^{i}$ can be
different dynamic models of a system when operating in different environments,
such as the dynamics of an autonomous car when operating in different weather
conditions.

Regression with uncertainty quantification. Denote the dataset of observations
for function $f^{i}$ as
$\mathcal{D}^{i}_{tr}\triangleq\\{(y^{i}_{t},x^{i}_{t})\\}_{t=1}^{t^{i}_{tr}}$,
where $t^{i}_{tr}$ depends on $f^{i}$. Denote prior mean function
$m:\mathcal{R}^{n_{x}}\to\mathbb{R}$, prior covariance function
$k:\mathcal{R}^{n_{x}}\times\mathcal{R}^{n_{x}}\to\mathbb{R}_{>0}$ which is
positive semidefinite. Given dataset $\mathcal{D}^{i}_{tr}$, we denote the
corresponding posterior mean function
$\mu_{m,k}^{i}:\mathcal{R}^{n_{x}}\to\mathbb{R}$ and posterior standard
deviation function $\sigma^{i}_{k}:\mathcal{R}^{n_{x}}\to\mathbb{R}$, which
are output from some algorithm $\mathcal{ALG}(m,k,\mathcal{D}^{i}_{tr})$. The
goal of this paper is to identify prior functions $m$ and $k$ which minimize
some objective function $J(m,k)$ and meanwhile capture the ground truth
$f^{i}(x)$ by intervals

$\displaystyle\mathcal{I}_{m,k}(x)$
$\displaystyle\triangleq[m(x)-q\sqrt{k(x,x)},m(x)+q\sqrt{k(x,x)}]$
$\displaystyle\mathcal{I}^{i}_{m,k}(x)$
$\displaystyle\triangleq[\mu^{i}_{m,k}(x)-q\sigma^{i}_{k}(x),\mu^{i}_{m,k}(x)+q\sigma^{i}_{k}(x)]$

with a pre-specified probability at least $1-\delta$, $\delta\in[0,1]$ for any
$x\in\mathcal{X}$. Notice that constants $q$ and $\delta$ are fixed a priori,
and they do not need to satisfy the relation between reliability factor and
level-of-confidence in the settings of confidence interval. Formally, the
above specifications for regression with uncertainty quantification is
formulated as the optimization problem below

$\displaystyle\min_{m,k}$ $\displaystyle\leavevmode\nobreak\ J(m,k)$ (2a) s.t.
$\displaystyle
P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}_{m,k}(x)\Big{)}\geqslant 1-\delta$
(2b) $\displaystyle
P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}^{i}_{m,k}(x)\Big{)}\geqslant
1-\delta$ (2c)
$\displaystyle\mu^{i}_{m,k},\sigma^{i}_{k}=\mathcal{ALG}(m,k,\mathcal{D}^{i}_{tr}).$
(2d)

Constraints (2b) and (2c) aim to ensure trustworthy uncertainty
quantification, and one example of the objective function $J(m,k)$ can be
negative marginal log likelihood (MLL) [2] evaluated over a dataset as in
[34][35]. Notice that in most cases given an arbitrary function $m(\cdot)$,
intervals $\mathcal{I}_{m,k}(\cdot)$ and $\mathcal{I}^{i}_{m,k}(\cdot)$ can be
arbitrarily large, and hence constraints (2b) and (2c) can always be
satisfied, by choosing function $k$ such that $k(x,x)$ is sufficiently large
if $\mathcal{ALG}$ follows a Bayesian inference framework, e.g., the prior
distribution and the likelihood function are Gaussian [5].

We do not assume any additional structure on distribution $\mathcal{P}_{f}$ or
on function $f^{i}$, e.g., $f^{i}\sim\mathcal{GP}(m,k)$, as in [16][17][18].
Instead, we assume we have access to a meta dataset
$\mathcal{D}^{meta}\triangleq\\{\mathcal{D}^{i}\\}_{i=1}^{n}$, where
$\mathcal{D}^{i}$ contains the observations of $f^{i}\sim\mathcal{P}_{f}$
i.i.d. drawn offline, which can be potentially used to estimate and optimize
$J(m,k)$ as well as the left hand sides of (2b) and (2c). For testing, we
consider functions $f^{j}\sim\mathcal{P}_{f}$, which are not necessarily
observed in $\mathcal{D}^{meta}$.

######  Remark II.1.

(Motivation of formulation (2)). Problem (2) is motivated by the problems of
safe learning/exploration [17] [14][34][24][41], where the system is required
to safely explore and online learn about an unknown environment using the data
collected along the system’s operation. During the early stage of exploration
where there are only a few (or even no) collected data, the control of the
system is mainly based on the prior knowledge, i.e., $(m,k)$, of the
environment. Therefore, it is crucial that the selected prior is trustworthy
such that the ground truth of function $f^{i}$ can be captured a priori by
$\mathcal{I}_{m,k}$, to ensure system safety when the system just starts to
explore the environment and no data is collected, and a posteriori by
$\mathcal{I}^{i}_{m,k}$ after collecting a (small) amount of data of $f^{i}$.
$\hfill\blacksquare$

Parameterization. Note that Problem (2) is a functional optimization problem
and can be hard to solve in general. To make the problem tractable, we
approximate the spaces of $m$ and $k$ using spaces of parameterized functions,
e.g., neural networks and finite-dimensional basis functions. Specifically, we
consider parameterized prior mean function $m_{\theta}$, with parameters
$\theta\in\mathbb{R}^{n_{\theta}}$, e.g., the weights of a deep neural
network, and parameterized prior mean function $k_{\phi}$, with parameters
$\phi\in\mathbb{R}^{n_{\phi}}$. Then Problem (2) can be rewritten as

$\displaystyle\min_{\theta,\phi}$ $\displaystyle\leavevmode\nobreak\
J(\theta,\phi)$ (3a) s.t. $\displaystyle
P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}_{\theta,\phi}(x)\Big{)}\geqslant
1-\delta$ (3b) $\displaystyle
P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}^{i}_{\theta,\phi}(x)\Big{)}\geqslant
1-\delta$ (3c)
$\displaystyle\mu^{i}_{\theta,\phi},\sigma^{i}_{\phi}=\mathcal{ALG}(m_{\theta},k_{\phi},\mathcal{D}^{i}_{tr}),$
(3d)

with the terms rewritten accordingly. In particular, we require $k_{\phi}$ to
be scalable, i.e., $k_{\phi}\triangleq\phi_{1}\kappa_{\phi_{2}}$, where
$\phi_{1}\in\mathbb{R}_{>0}$, $\phi_{2}\in\mathbb{R}^{n_{\phi}-1}$,
$\phi=[\phi_{1},\phi_{2}]$, and
$\kappa_{\phi_{2}}:\mathcal{X}\times\mathcal{X}\to[0,1]$ is a covariance
function, such that the sizes of the intervals $\mathcal{I}_{\theta,\phi}(x)$
and $\mathcal{I}^{i}_{m,k}(x)$ can be arbitrarily large by increasing
$\phi_{1}$, if $\mathcal{ALG}$ follows some Bayesian inference frameworks such
as GPR [2]. Similar to (2), this implies that there always exists $\phi$ which
satisfies (3b) and (3c), and hence Problem (3) is always feasible. One example
of covariance function which satisfies the above requirement is the widely-
used square-exponential kernel
$k(x,x^{\prime})=\sigma_{f}^{2}\exp(-\frac{\|x-x^{\prime}\|^{2}}{2\ell})+\sigma^{2}_{e}\delta_{xx^{\prime}}$.

##  III The Trust-Bayes framework

In this section, since distributions $\mathcal{P}_{f}$ and possibly
$\mathcal{P}_{x}$ are unknown, we introduce the framework Trust-Bayes, to
solve problem (3) leveraging the empirical estimates from meta dataset
$\mathcal{D}^{meta}$.

Define 0-1 loss functions

$\displaystyle c^{i}_{1}(x)$
$\displaystyle\triangleq\bm{1}[f^{i}(x)\in\mathcal{I}_{\theta,\phi}(x)]$
$\displaystyle c^{i}_{2}(x)$
$\displaystyle\triangleq\bm{1}[f^{i}(x)\in\mathcal{I}^{i}_{\theta,\phi}(x)].$

For each $\mathcal{D}^{i}\in\mathcal{D}^{meta}$, it is split into
$\mathcal{D}^{i}_{tr}\subset\mathcal{D}^{i}$ for obtaining the posterior
functions $(\mu^{i}_{\theta,\phi},\sigma^{i}_{k})$ for predictions and an
evaluation dataset
$\mathcal{D}^{i}_{eval}\triangleq\mathcal{D}^{i}\setminus\mathcal{D}^{i}_{tr}$
for evaluating the performances of the posterior functions and estimating (3b)
and (3c). We write
$\mathcal{D}^{i}_{eval}=\\{(y^{i}_{t},x^{i}_{t})\\}_{t=1}^{t^{i}_{eval}}$ and
$\mathcal{D}^{i}_{tr}=\\{(y^{i}_{t},x^{i}_{t})\\}_{t=1}^{t^{i}_{tr}}$.For any
$\gamma\in(0,0.5]$, define

$\displaystyle
p_{1}(\gamma)\triangleq\Big{(}1-2\gamma\Big{)}\cdot\Big{(}\frac{1}{n}\sum_{i=1}^{n}\frac{1}{t^{i}_{eval}}\sum_{t=1}^{t^{i}_{eval}}c^{i}_{1}(x^{i}_{t})$
$\displaystyle\qquad\qquad-\sqrt{\frac{\log(2/\gamma)\sum_{i=1}^{n}\frac{1}{t^{i}_{eval}}}{2n^{2}}}-\sqrt{\frac{\log(2/\gamma)}{2n}}\Big{)},$
$\displaystyle
p_{2}(\gamma)\triangleq\Big{(}1-2\gamma\Big{)}\cdot\Big{(}\frac{1}{n}\sum_{i=1}^{n}\frac{1}{t^{i}_{eval}}\sum_{t=1}^{t^{i}_{eval}}c^{i}_{2}(x^{i}_{t})$
$\displaystyle\qquad\qquad-\sqrt{\frac{\log(2/\gamma)\sum_{i=1}^{n}\frac{1}{t^{i}_{eval}}}{2n^{2}}}-\sqrt{\frac{\log(2/\gamma)}{2n}}\Big{)}.$

Then the following theorem characterizes the lower bounds of
$P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}_{\theta,\phi}(x)\Big{)}$ and
$P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}^{i}_{\theta,\phi}(x)\Big{)}$.

######  Theorem III.1.

Suppose $\\{f^{i}\\}_{i=1}^{n}$ are sampled i.i.d. from a latent distribution
$\mathcal{P}_{f}$. Suppose the training dataset $\mathcal{D}^{i}_{tr}$ for
each function $f^{i}$ is generated through a latent conditional distribution
$\mathcal{P}_{\mathcal{D}^{i}_{tr}\mid f^{i}}$, i.e.,
$\mathcal{D}^{i}_{tr}\sim\mathcal{P}_{\mathcal{D}^{i}_{tr}\mid f^{i}}$.
Suppose for each $f^{i}$, $\\{x^{i}_{t}\\}_{t=1}^{t^{i}_{eval}}$ are sampled
i.i.d. from a latent function $\mathcal{P}_{x}$. Then the following
inequalities hold:

$\displaystyle
P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}_{\theta,\phi}(x)\Big{)}\geqslant\max_{\gamma\in(0,0.5]}p_{1}(\gamma),$
$\displaystyle
P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}^{i}_{\theta,\phi}(x)\Big{)}\geqslant\max_{\gamma\in(0,0.5]}p_{2}(\gamma).$

$\hfill\blacksquare$

The proof of the theorem is deferred to Section IV. The lower bounds in
Theorem III.1 indicates that
$P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}_{\theta,\phi}(x)\Big{)}$ and
$P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}^{i}_{\theta,\phi}(x)\Big{)}$ can be
lower bounded by their corresponding empirical estimates of (i.e.,
$\frac{1}{n}\sum_{i=1}^{n}\frac{1}{t^{i}_{eval}}\sum_{t=1}^{t^{i}_{eval}}c^{i}_{1}(x^{i}_{t})$
and
$\frac{1}{n}\sum_{i=1}^{n}\frac{1}{t^{i}_{eval}}\sum_{t=1}^{t^{i}_{eval}}c^{i}_{2}(x^{i}_{t})$)
using the evaluation datasets $\mathcal{D}^{i}_{eval}$, $i=1,\cdots,n$, in the
meta training dataset $\mathcal{D}^{meta}$, and the error terms diminish by
$t^{i}_{eval}$ and $n$, the sizes of $\mathcal{D}^{i}_{eval}$ and the number
of functions in $\mathcal{D}^{meta}$. This provides empirical underestimates
of $P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}_{\theta,\phi}(x)\Big{)}$ and
$P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}^{i}_{\theta,\phi}(x)\Big{)}$ but
sufficient verification for constraints (3b) and (3c).

By the lower bounds in the theorem above, we can approximate Problem (3) with

$\displaystyle\min_{\theta,\phi}$ $\displaystyle\leavevmode\nobreak\
J(\theta,\phi)$ (4a) s.t.
$\displaystyle\max_{\gamma\in(0,0.5]}p_{1}(\gamma)\geqslant 1-\delta,$ (4b)
$\displaystyle\max_{\gamma\in(0,0.5]}p_{2}(\gamma)\geqslant 1-\delta,$ (4c)
$\displaystyle\mu^{i}_{\theta,\phi},\sigma^{i}_{\phi}=\mathcal{ALG}(m_{\theta},k_{\phi},\mathcal{D}^{i}_{tr}),$
(4d)

which is our proposed formulation of Trust-Bayes.

######  Remark III.1.

(Feasibility and sample complexity). As discussed below Problem (3), the sizes
of intervals $\mathcal{I}_{\theta,\phi}(x)$ and
$\mathcal{I}_{\theta,\phi}^{i}(x)$ can be arbitrarily large by choosing $\phi$
accordingly. This implies that
$\frac{1}{n}\sum_{i=1}^{n}\frac{1}{t^{i}_{eval}}\sum_{t=1}^{t^{i}_{eval}}c^{i}_{j}(x_{t})$,
$j\in\\{1,2\\}$, can be as large as one. Therefore, by the forms of (4b) and
(4c), the feasibility of (4) can be checked by verifying whether

$\displaystyle\max_{\gamma\in(0,0.5]}$
$\displaystyle(1-2\gamma)\Big{(}1-\sqrt{\frac{\log(2/\gamma)\sum_{i=1}^{n}\frac{1}{t^{i}_{eval}}}{2n^{2}}}-\sqrt{\frac{\log(2/\gamma)}{2n}}\Big{)}$
$\displaystyle\geqslant 1-\delta$ (5)

holds for a given $\delta$, which can be done through numerically solving the
left hand side using, e.g., gradient ascent. Note that for any
$\gamma\in(0,0.5]$,
$\sqrt{\frac{\log(2/\gamma)\sum_{i=1}^{n}\frac{1}{t^{i}_{eval}}}{2n^{2}}}$ and
$\sqrt{\frac{\log(2/\gamma)}{2n}}$ can be arbitrarily small by increasing
sample sizes $t^{i}_{eval}$ and $n$. Therefore, inequality (III.1) provides
the sample complexity for the lower bound of feasible $\delta$ and implies
that the freedom of selecting $\delta$ can be increased by increasing
$t^{i}_{eval}$ and $n$. Note that (4b) and (4c) approximate the probability
constraints (3b) and (3c) by underestimating the probabilities with the
corresponding empirical estimates. This is reminiscent of the no-free-lunch
property between estimation error and the number of data in statistical
learning theory [21]. $\hfill\blacksquare$

##  IV Proof of Theorem III.1

The proof of the theorem leverages several properties of sub-Gaussian random
variables defined as follows.

######  Definition IV.1.

A random variable $X$ with mean $\mu=\mathbb{E}[X]$ is sub-Gaussian if there
is a positive number $\sigma$ such that
$\mathbb{E}[e^{\lambda(X-\mu)}]\leqslant e^{\sigma^{2}\lambda^{2}/2}$ for all
$\lambda\in\mathbb{R}$. $\hfill\blacksquare$

Based on the above definition, we can denote a sub-Gaussian random variable
$X$ as $X\sim\textsf{sub-Gauss}(\mu,\sigma)$ The properties leveraged in the
proof are included below for completeness, and they can be found in [42].

######  Lemma IV.2.

(Bounded random variables are sub-Gaussian). A random variable $X$ with mean
$\mu=\mathbb{E}[X]$ bounded as $a\leqslant X\leqslant b$ is sub-Gaussian with
$\sigma^{2}=(\frac{b-a}{2})^{2}$. $\hfill\blacksquare$

######  Lemma IV.3.

(Preservation under superposition). If $X_{1}\sim\textsf{sub-
Gauss}(\mu_{1},\sigma_{1}^{2})$ and $X_{2}\sim\textsf{sub-
Gauss}(\mu_{2},\sigma_{2}^{2})$, then $X_{1}+X_{2}\sim\textsf{sub-
Gauss}(\mu_{1}+\mu_{2},\sigma_{1}^{2}+\sigma_{2}^{2})$ and
$aX_{1}\sim\textsf{sub-Gauss}(a\mu_{1},a^{2}\sigma_{1}^{2})$.
$\hfill\blacksquare$

######  Lemma IV.4.

(Hoeffding’s inequality). Suppose that the variables $X_{i}\sim\textsf{sub-
Gauss}(\mu,\sigma^{2})$, $i=1,\cdots,n$, are independent. Then for all
$\epsilon\geqslant 0$, we have

$P\Big{(}|X-\mu|\geqslant\epsilon\Big{)}\leqslant
2\exp\Big{(}-\frac{\epsilon^{2}}{2\sigma^{2}}\Big{)}.$

$\hfill\blacksquare$

Before proving the theorem, we first introduce the following expectations and
empirical estimates:

$\displaystyle C_{1}$
$\displaystyle\triangleq\mathbb{E}_{f^{i},x}[c^{i}_{1}(x)]=P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}_{\theta,\phi}(x)\Big{)}$
$\displaystyle C_{2}$
$\displaystyle\triangleq\mathbb{E}_{f^{i},x}[c^{i}_{2}(x)]=P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}^{i}_{\theta,\phi}(x)\Big{)}$
$\displaystyle C^{i}_{1}$
$\displaystyle\triangleq\mathbb{E}_{x}[c^{i}_{1}(x_{t})\mid f^{i}],\quad
C^{i}_{2}\triangleq\mathbb{E}_{x}[c^{i}_{2}(x_{t})\mid f^{i}],$
$\displaystyle\hat{C}^{i}_{1}$
$\displaystyle\triangleq\frac{1}{t^{i}_{eval}}\sum_{t=1}^{t^{i}_{eval}}c^{i}_{1}(x^{i}_{t}),\quad\hat{C}^{i}_{2}\triangleq\frac{1}{t^{i}_{eval}}\sum_{t=1}^{t^{i}_{eval}}c^{i}_{2}(x^{i}_{t}).$

Then the following lemma characterizes the lower bounds for expectations
$C_{1}$ and $C_{2}$ in terms of empirical estimates $\hat{C}^{i}_{1}$ and
$\hat{C}^{i}_{2}$.

######  Lemma IV.5.

Each of the followings hold with probability at least $1-2\gamma$:

$\displaystyle C_{1}$
$\displaystyle\geqslant\frac{1}{n}\sum_{i=1}^{n}\hat{C}^{i}_{1}-\sqrt{\frac{\log(2/\gamma)\sum_{i=1}^{n}\frac{1}{t^{i}_{eval}}}{2n^{2}}}-\sqrt{\frac{\log(2/\gamma)}{2n}},$
$\displaystyle C_{2}$
$\displaystyle\geqslant\frac{1}{n}\sum_{i=1}^{n}\hat{C}^{i}_{2}-\sqrt{\frac{\log(2/\gamma)\sum_{i=1}^{n}\frac{1}{t^{i}_{eval}}}{2n^{2}}}-\sqrt{\frac{\log(2/\gamma)}{2n}}.$

Proof: By the definitions of $c^{i}_{1}$ and $C^{i}_{1}$, we have
$C^{i}_{1}\in[0,1]$ and $\mathbb{E}_{f^{i}}[C^{i}_{1}]=C_{1}$. Then Lemma IV.2
renders that $C^{i}_{1}\sim\textsf{sub-Gauss}(C_{1},\frac{1}{4})$, and hence
Lemma IV.3 renders that $\frac{1}{n}\sum_{i=1}^{n}C^{i}_{1}\sim\textsf{sub-
Gauss}(C_{1},\frac{1}{4n})$. Therefore, by Lemma IV.4, we have
$P\Big{(}|\frac{1}{n}\sum_{i=1}^{n}C^{i}_{1}-C_{1}|\geqslant\epsilon\Big{)}\leqslant
2\exp\Big{(}-2n\epsilon^{2}\Big{)}.$ After some simple algebraic
transformations, we have

$\displaystyle
C_{1}\geqslant\frac{1}{n}\sum_{i=1}^{n}C^{i}_{1}-\sqrt{\frac{\log(2/\gamma)}{2n}}.$
(6)

with probability at least $1-\gamma$.

Recall the definitions of $c^{i}_{1}(x_{t})$ and $C^{i}_{1}$. Notice that
$c^{i}_{1}(x_{t})\in[0,1]$ for all $i=1,\cdots,n$ and $t=1,2,\cdots$. Then
Lemma IV.2 renders that $c^{i}_{1}(x_{t})\sim\textsf{sub-
Gauss}(C^{i}_{1},\frac{1}{4})$, and Lemma IV.3 renders that
$\hat{C}^{i}_{1}=\frac{1}{t^{i}_{eval}}\sum_{t=1}^{t^{i}_{eval}}c^{i}_{1}(x_{t})\sim\textsf{sub-
Gauss}(C^{i}_{1},\frac{1}{4t^{i}_{eval}})$. Furthermore, by Lemma IV.3, we
also have $\frac{1}{n}\sum_{i=1}^{n}\hat{C}^{i}_{1}\sim\textsf{sub-
Gauss}(\frac{1}{n}\sum_{i=1}^{n}C^{i}_{1},\frac{1}{n^{2}}\sum_{i=1}^{n}\frac{1}{4t^{i}_{eval}})$.
Then following similar logic to (6), we have

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}C^{i}_{1}\geqslant\frac{1}{n}\sum_{i=1}^{n}\hat{C}^{i}_{1}-\sqrt{\frac{\log(2/\gamma)\sum_{i=1}^{n}\frac{1}{t^{i}_{eval}}}{2n^{2}}}.$
(7)

with probability at least $1-\gamma$. Combining (7) with (6) renders

$\displaystyle
C_{1}\geqslant\frac{1}{n}\sum_{i=1}^{n}\hat{C}^{i}_{1}-\sqrt{\frac{\log(2/\gamma)\sum_{i=1}^{n}\frac{1}{t^{i}_{eval}}}{2n^{2}}}-\sqrt{\frac{\log(2/\gamma)}{2n}}$

with probability at least $(1-\gamma)^{2}\geqslant 1-2\gamma$. Following the
same logic, we also have, with probability at least $1-2\gamma$,

$\displaystyle
C_{2}\geqslant\frac{1}{n}\sum_{i=1}^{n}\hat{C}^{i}_{2}-\sqrt{\frac{\log(2/\gamma)\sum_{i=1}^{n}\frac{1}{t^{i}_{eval}}}{2n^{2}}}-\sqrt{\frac{\log(2/\gamma)}{2n}}.$

$\hfill\blacksquare$

The the proof of Theorem III.1 is given as follows.

Let
$G_{1,\gamma}\triangleq\frac{1}{n}\sum_{i=1}^{n}\hat{C}^{i}_{1}-\sqrt{\frac{\log(2/\gamma)\sum_{i=1}^{n}\frac{1}{t^{i}_{eval}}}{2n^{2}}}-\sqrt{\frac{\log(2/\gamma)}{2n}}$.
Then by Lemma IV.5, we have

$\displaystyle
P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}_{\theta,\phi}(x)\Big{)}$
$\displaystyle=P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}_{\theta,\phi}(x)|C_{1}\geqslant
G_{\hat{C}^{i}_{1},\gamma}\Big{)}P\Big{(}C_{1}\geqslant
G_{\hat{C}^{i}_{1},\gamma}\Big{)}$ $\displaystyle\leavevmode\nobreak\
+P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}_{\theta,\phi}(x)|C_{1}<G_{\hat{C}^{i}_{1},\gamma}\Big{)}P\Big{(}C_{1}<G_{\hat{C}^{i}_{1},\gamma}\Big{)}$
$\displaystyle\geqslant
P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}_{\theta,\phi}(x)|C_{1}\geqslant
G_{\hat{C}^{i}_{1},\gamma}\Big{)}P\Big{(}C_{1}\geqslant
G_{\hat{C}^{i}_{1},\gamma}\Big{)}$
$\displaystyle\geqslant\Big{(}\frac{1}{n}\sum_{i=1}^{n}\hat{C}^{i}_{1}-\sqrt{\frac{\log(2/\gamma)\sum_{i=1}^{n}\frac{1}{t^{i}_{eval}}}{2n^{2}}}-\sqrt{\frac{\log(2/\gamma)}{2n}}\Big{)}$
$\displaystyle\qquad\cdot\Big{(}1-2\gamma\Big{)}.$ (8)

for any $\gamma\in(0,0.5]$. Following the same logic, we have

$\displaystyle
P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}^{i}_{\theta,\phi}(x)\Big{)}$
$\displaystyle\geqslant\Big{(}\frac{1}{n}\sum_{i=1}^{n}\hat{C}^{i}_{2}-\sqrt{\frac{\log(2/\gamma)\sum_{i=1}^{n}\frac{1}{t^{i}_{eval}}}{2n^{2}}}-\sqrt{\frac{\log(2/\gamma)}{2n}}\Big{)})$
$\displaystyle\qquad\cdot\Big{(}1-2\gamma\Big{)}.$ (9)

for any $\gamma\in(0,0.5]$. To ensure the tightness of the bounds, we maximize
the right hand sides of (IV) and (IV) with respect to $\gamma$. This completes
the proof of the theorem. $\hfill\blacksquare$

##  V Case study

In this section, we aim to verify whether the Trust-Bayes formulation in (4)
can provide trustworthy uncertainty quantification, i.e., satisfying the
constraints in (3), and whether it is necessary. In this case study, we
consider GPR [2] as the algorithm $\mathcal{ALG}$. Then the prior predictive
distribution for $f^{i}(x)$ is given by
$\mathcal{N}(m_{\theta}(x),k_{\phi}(x,x))$ and the posterior predictive
distribution is given by
$\mathcal{N}(\mu^{i}_{\theta,\phi}(x),(\sigma^{i}_{\phi})^{2}(x))$, where

$\displaystyle\mu^{i}_{\theta,\phi}(x)=m_{\theta}(x)+k_{\phi}(x,X^{i})k_{\phi}^{-1}(X^{i},X^{i})(Y^{i}-m_{\theta}(X^{i}))$
$\displaystyle(\sigma^{i}_{\phi})^{2}(x)=k_{\phi}(x,x)+k_{\phi}(x,X^{i})k_{\phi}^{-1}(X^{i},X^{i})k_{\phi}(X^{i},x).$

where $X^{i}$ aggregates all the inputs $x^{i}_{t}$ in $\mathcal{D}^{i}$,
$Y^{i}$ aggregates all the outputs.

Experiment setup. In this experiment, we let $x\in[0,1]$ and for each
$f^{i}\in\mathcal{P}_{f}$,
$f^{i}(x)=d^{i}x^{2}+\sum_{m=1}^{10}\alpha^{i}a^{i}_{m}\sin(w^{i}_{m}x+\beta^{i}_{m})+(1-\alpha^{i})b^{i}_{m}\sin(u^{i}_{m}x+\beta^{i}_{m})$,
where

$\displaystyle a^{i}_{m}$ $\displaystyle\sim
0.5\mathcal{N}(-20,5)+0.5\mathcal{N}(10,2)$ (10a) $\displaystyle b^{i}_{m}$
$\displaystyle\sim 0.5\mathcal{N}(-1,0.1)+0.5\mathcal{N}(1,0.1)$ (10b)
$\displaystyle w^{i}_{m}$ $\displaystyle\sim
0.5\mathcal{N}(-10,10)+0.5\mathcal{N}(10,10)$ (10c) $\displaystyle u^{i}_{m}$
$\displaystyle\sim 0.5\mathcal{N}(-100,10)+0.5\mathcal{N}(100,10)$ (10d)
$\displaystyle\beta^{i}_{m}$ $\displaystyle\sim\mathcal{N}(0,1)$ (10e)
$\displaystyle d^{i}$ $\displaystyle\sim
0.5\mathcal{N}(-10,1)+0.5\mathcal{N}(10,1)$ (10f) $\displaystyle\alpha^{i}$
$\displaystyle\sim\mathcal{B}(0.5).$ (10g)

$\mathcal{N}$ denotes normal distribution and $\mathcal{B}$ denotes Bernoulli
distribution. For learning using GPR, we consider constant
$m_{\theta}(x)=\theta$ for the prior mean function and kernel
$k_{\phi}(x,x^{\prime})=\phi_{1}^{2}\exp(-\|x-x^{\prime}\|^{2}\phi_{2})$. For
trustworthy uncertainty quantification, we select $\delta=0.1$ and
$q_{\delta}=1.64$ for intervals $\mathcal{I}_{\theta,\phi}(x)$ and
$\mathcal{I}^{i}_{\theta,\phi}(x)$. That is, we require the predicted
intervals with 90% nominal confidence level, according to the Z-score table,
to indeed include at least 90% of the true values. Note that according to
(10), there are parameters following Gaussian mixtures distributions and
Bernoulli distribution, and therefore $f^{i}$ does not follow a Gaussian
process. Furthermore, each $f^{i}$ has a trend, and therefore
$k_{\phi}(x,x^{\prime})$ is not a suitable prior covariance function. The
purpose of this setup is to demonstrate that Trust-Bayes is able to provide
trustworthy uncertainty quantification even under such mis-specification of
the prior distribution, which can happen when there is no much prior knowledge
on the distribution and structure on the target function $f^{i}$.

Training. For the training dataset $\mathcal{D}^{tr}$, we sample 2000
functions $f^{i}$, i.e., $n=2000$, following (10), and for each corresponding
training dataset $\mathcal{D}^{i}$ we uniformly sample over interval $[0,1]$
for inputs $t^{i}_{tr}=20$ for obtaining the posterior and $t^{i}_{eval}=100$
for evaluating the posterior for (4a) to (4c). Note that with the above
specification of $\delta$, $n$, and $t^{i}_{eval}$, (III.1) holds and can be
verified by plugging in $\gamma=0.001$. For the objective function
$J(\theta,\phi)$, we consider negative MLL over the whole meta dataset
$\mathcal{D}^{tr}$ as in [34, 35, 36].

Testing. Notice that the probabilities in (3b) and (3c) can be challenging to
exactly obtain. As an estimation, the trained hyperparameters $(\theta,\phi)$
are tested against a dataset composed of $10,000$ functions $f^{i}$ randomly
sampled over (10) with each $f^{i}$ has testing $10,000$ inputs $x$ uniformly
sampled over $[0,1]$. The corresponding posterior distribution is obtained
after observing each function over 20 inputs without further tuning the
hyperparameters.

Comparison. Our results are compared against the Meta-prior method used in
[34, 35, 36], where the hyperparameters in the prior distributions are meta-
trained using negative MLL. Note that the Bayesian regression approach used in
[34][35] is a special case of GPR [2]. The difference between Trust-Bayes and
Meta-prior in this case study is the addition of (4b) and (4c).

Results. Figure 1 shows the convergence of
$\frac{1}{n}\sum_{i=1}^{n}\frac{1}{t^{i}_{eval}}\sum_{t=1}^{t^{i}_{eval}}c^{i}_{1}(x^{i}_{t})$
and
$\frac{1}{n}\sum_{i=1}^{n}\frac{1}{t^{i}_{eval}}\sum_{t=1}^{t^{i}_{eval}}c^{i}_{2}(x^{i}_{t})$,
the empirical estimates of
$P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}_{\theta,\phi}(x)\Big{)}$ and
$P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}^{i}_{\theta,\phi}(x)\Big{)}$,
respectively, by the evaluation dataset
$\\{\mathcal{D}^{i}_{eval}\\}_{i=1}^{n}$ in the meta dataset
$\mathcal{D}^{meta}$. From Figure 1, we can see that the training of Meta-
prior converges when the empirical estimates are still at a level much lower
than the required rate of inclusion $1-\delta=0.9$. In contrast, the training
of Trust-Bayes converges much faster and beyond the required rate of
inclusion.

Table I provides a comparison between the $(\theta,\phi)$ trained till
convergence by Trust-Bayes and that by Meta-prior over the estimated
$P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}_{\theta,\phi}(x)\Big{)}$,
$P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}^{i}_{\theta,\phi}(x)\Big{)}$, their
empirical estimates by $\mathcal{D}^{i}_{eval}$ as well as the mean-squared
error (MSE) for the posterior predictions over the testing dataset. From Table
I, we can see that Trust-Bayes performs better than Meta-prior over all the
metrics. Specifically,
$P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}_{\theta,\phi}(x)\Big{)}$ and
$P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}^{i}_{\theta,\phi}(x)\Big{)}$ are
higher than the required inclusion rate (i.e, 0.9) when $(\theta,\phi)$ are
trained by Trust-Bayes, which verifies that Trust-Bayes can provide
trustworthy uncertainty quantification. In contrast, when $(\theta,\phi)$ are
trained by Meta-prior, the rates of inclusion are much lower than the required
inclusion rate although the intervals are constructed with 90% nominal level
of confidence. This comparison with Meta-prior shows the necessary of using
Trust-Bayes to provide trustworthy uncertainty quantification when it is
uncertain whether the prior is properly specified or not.

(a) Trust-Bayes

(b) Meta-prior

Figure 1: Comparison over empirical estimates of $P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}_{\theta,\phi}(x)\Big{)}$ and $P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}^{i}_{\theta,\phi}(x)\Big{)}$ by the evaluation dataset $\\{\mathcal{D}^{i}_{eval}\\}_{i=1}^{n}$ between Trust-Bayes and Meta-prior during meta-training | Trust-Bayes | Meta-prior
---|---|---
$\frac{1}{n}\sum_{i=1}^{n}\frac{1}{t^{i}_{eval}}\sum_{t=1}^{t^{i}_{eval}}c^{i}_{1}(x^{i}_{t})$ | 0.997 | 0.623
$\frac{1}{n}\sum_{i=1}^{n}\frac{1}{t^{i}_{eval}}\sum_{t=1}^{t^{i}_{eval}}c^{i}_{2}(x^{i}_{t})$ | 0.999 | 0.706
$P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}_{\theta,\phi}(x)\Big{)}$ | 0.996 | 0.630
$P_{f^{i},x}\Big{(}f^{i}(x)\in\mathcal{I}^{i}_{\theta,\phi}(x)\Big{)}$ | 0.999 | 0.708
MSE | 91.42 | 125.15
Table I: Comparisons over empirical and expected prior/posterior inclusions
between Trust-Bayes and Meta-prior

Figure 2 provides a visual comparison over the inclusions of the true values
of $f^{i}(x)$ by the prior 90% confidence intervals
$\mathcal{I}_{\theta,\phi}(x)$ trained using Trust-Bayes and Meta-prior,
respectively, and Figure 3 provides a comparison over the inclusions by the
posterior 90% confidence intervals $\mathcal{I}^{i}_{\theta,\phi}(x)$. Ten
functions are sampled from (10) to provide a visualization of the distribution
of the functions. From these two figures, we can see that the 90% intervals
generated by Trust-Bayes are larger and are able to capture most of the values
of the functions including those with larger variation, whereas Meta-prior
generates smaller intervals and fails to capture the values of the functions
with larger variation.

(a) Trust-Bayes

(b) Meta-prior

Figure 2: Ten random functions drawn from (10) and comparison between prior
predictions by Trust-Bayes and Meta-prior

(a) Posterior predictions over ten sample functions by Trust-Bayes

(b) Posterior predictions over ten sample functions by Meta-prior

Figure 3: Comparison between posterior predictions by Trust-Bayes and Meta-
prior over ten random functions from (10)

##  VI Conclusion

We consider trustworthy uncertainty quantification in Bayesian regression
problems. We formulate trustworthy uncertainty quantification as constraints
on capturing the ground truths of the function by intervals depending on the
prior and posterior distributions with a pre-specified probability. We
propose, Trust-Bayes, a Bayesian meta learning framework which is cognizant of
trustworthy uncertainty quantification without explicit assumptions on the
model/distribution of the functions. We characterize the lower bounds of the
probabilities of the ground truth being captured by the specified intervals in
terms of the empirical estimates and analyze the sample complexity with
respect to the feasible pre-specified probability for trustworthy uncertainty
quantification. Monte Carlo simulation is conducted for evaluation and
comparison through a case study using GPR, which verifies the proposed
framework and demonstrates the necessary of Trust-Bayes for trustworthy
uncertainty quantification when the prior is not necessarily correctly
specified.

## References

  * [1] T. M. Mitchell, “Machine learning,” 1997.
  * [2] C. K. Williams and C. E. Rasmussen, _Gaussian Processes for Machine Learning_. MIT Press, 2006.
  * [3] G. Welch, G. Bishop _et al._ , “An introduction to the kalman filter,” 1995\.
  * [4] J. Lampinen and A. Vehtari, “Bayesian approach for neural networks—review and case studies,” _Neural networks_ , vol. 14, no. 3, pp. 257–274, 2001\.
  * [5] C. M. Bishop, _Pattern recognition_. Springer, 2006.
  * [6] T. Choi and M. J. Schervish, “On posterior consistency in nonparametric regression problems,” _Journal of Multivariate Analysis_ , vol. 98, no. 10, pp. 1969–1987, Jan. 2007.
  * [7] G. E. Box and G. C. Tiao, _Bayesian inference in statistical analysis_. John Wiley & Sons, 2011.
  * [8] D. A. McAllester, “Some PAC-Bayesian theorems,” _Machine Learning_ , vol. 37, no. 3, pp. 355–363, 1999.
  * [9] R. Garnett, _Bayesian optimization_. Cambridge University Press, 2023.
  * [10] A. Lederer, Z. Yang, J. Jiao, and S. Hirche, “Cooperative control of uncertain multi-agent systems via distributed gaussian processes,” _IEEE Transactions on Automatic Control_ , 2022.
  * [11] M. Liu, G. Chowdhary, B. C. Da Silva, S.-Y. Liu, and J. P. How, “Gaussian processes for learning and control: A tutorial with examples,” _IEEE Control Systems Magazine_ , vol. 38, no. 5, pp. 53–86, 2018.
  * [12] N. Virani, D. K. Jha, Z. Yuan, I. Shekhawat, and A. Ray, “Imitation of demonstrations using bayesian filtering with nonparametric data-driven models,” _Journal of Dynamic Systems, Measurement, and Control_ , vol. 140, no. 3, 2018.
  * [13] M. Mukadam, X. Yan, and B. Boots, “Gaussian process motion planning,” in _2016 IEEE international conference on robotics and automation (ICRA)_. IEEE, 2016, pp. 9–15.
  * [14] Z. Yuan and M. Zhu, “dSLAP: Distributed safe learning and planning for multi-robot systems,” in _Proc. IEEE Conf. Decision and Control (CDC)_ , 2022, pp. 5864–5869.
  * [15] A. Chiuso and G. Pillonetto, “A bayesian approach to sparse dynamic network identification,” _Automatica_ , vol. 48, no. 8, pp. 1553–1565, 2012.
  * [16] N. Srinivas, A. Krause, S. M. Kakade, and M. W. Seeger, “Information-theoretic regret bounds for Gaussian process optimization in the bandit setting,” _IEEE Trans. Information Theory_ , vol. 58, no. 5, pp. 3250–3265, Jan. 2012\.
  * [17] J. F. Fisac, A. K. Akametalu, M. N. Zeilinger, S. Kaynama, J. Gillula, and C. J. Tomlin, “A general safety framework for learning-based control in uncertain robotic systems,” _IEEE Trans. Automatic Control_ , vol. 64, no. 7, pp. 2737–2752, 2019.
  * [18] Z. Yuan and M. Zhu, “Lightweight distributed gaussian process regression for online machine learning,” _IEEE Transactions on Automatic Control_ , 2024\. To appear.
  * [19] A. G. Wilson, Z. Hu, R. Salakhutdinov, and E. P. Xing, “Deep kernel learning,” in _Artificial intelligence and statistics_. PMLR, 2016, pp. 370–378.
  * [20] Y. Shen, T. Chen, and G. B. Giannakis, “Random feature-based online multi-kernel learning in environments with unknown dynamics,” _The Journal of Machine Learning Research_ , vol. 20, no. 1, pp. 773–808, 2019.
  * [21] V. Vapnik, _The Nature of Statistical Learning Theory_. Springer science & business media, 2013.
  * [22] Z. Zhou, O. S. Oguz, M. Leibold, and M. Buss, “A general framework to increase safety of learning algorithms for dynamical systems based on region of attraction estimation,” _IEEE Trans. Robotics_ , 2020.
  * [23] J. Umlauft, T. Beckers, A. Capone, A. Lederer, and S. Hirche, “Smart forgetting for safe online learning with gaussian processes,” in _Learning for dynamics and control_ , 2020, pp. 160–169.
  * [24] F. Berkenkamp, A. Krause, and A. P. Schoellig, “Bayesian optimization with safety constraints: safe and automatic parameter tuning in robotics,” _Machine Learning_ , vol. 112, no. 10, pp. 3713–3747, 2023.
  * [25] Y. Sui, V. Zhuang, J. Burdick, and Y. Yue, “Stagewise safe bayesian optimization with gaussian processes,” in _Proc. Int. Conf. Machine Learning (ICML)_ , 2018, pp. 4781–4789.
  * [26] C. Finn, P. Abbeel, and S. Levine, “Model-agnostic meta-learning for fast adaptation of deep networks,” in _International conference on machine learning_. PMLR, 2017, pp. 1126–1135.
  * [27] A. Rajeswaran, C. Finn, S. M. Kakade, and S. Levine, “Meta-learning with implicit gradients,” _Proc. Advances in Neural Information Processing Systems (NeurIPS)_ , vol. 32, 2019.
  * [28] M. Khodak, M.-F. F. Balcan, and A. S. Talwalkar, “Adaptive gradient-based meta-learning methods,” _Proc. Advances in Neural Information Processing Systems (NeurIPS)_ , vol. 32, 2019.
  * [29] J. Lee, J. Tack, N. Lee, and J. Shin, “Meta-learning sparse implicit neural representations,” _Advances in Neural Information Processing Systems_ , vol. 34, pp. 11 769–11 780, 2021.
  * [30] S. Xu and M. Zhu, “Online constrained meta-learning: Provable guarantees for generalization,” _Advances in Neural Information Processing Systems_ , vol. 36, 2024.
  * [31] A. Fallah, A. Mokhtari, and A. Ozdaglar, “Generalization of model-agnostic meta-learning algorithms: Recurring and unseen tasks,” _Proc. Advances in Neural Information Processing Systems (NeurIPS)_ , vol. 34, pp. 5469–5480, 2021\.
  * [32] J. Yoon, T. Kim, O. Dia, S. Kim, Y. Bengio, and S. Ahn, “Bayesian model-agnostic meta-learning,” _Proc. Advances in Neural Information Processing Systems (NeurIPS)_ , vol. 31, 2018.
  * [33] X. Zhang, D. Meng, H. Gouk, and T. M. Hospedales, “Shallow bayesian meta learning for real-world few-shot recognition,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2021, pp. 651–660.
  * [34] T. Lew, A. Sharma, J. Harrison, A. Bylard, and M. Pavone, “Safe active dynamics learning and control: A sequential exploration–exploitation framework,” _IEEE Trans. Robotics_ , vol. 38, no. 5, pp. 2888–2907, 2022\.
  * [35] J. Harrison, A. Sharma, and M. Pavone, “Meta-learning priors for efficient online bayesian regression,” in _Algorithmic Foundations of Robotics XIII: Proceedings of the 13th Workshop on the Algorithmic Foundations of Robotics 13_ , 2020, pp. 318–337.
  * [36] J. Rothfuss, V. Fortuin, M. Josifoski, and A. Krause, “Pacoh: Bayes-optimal meta-learning with pac-guarantees,” in _Proc. Int. Conf. Machine Learning (ICML)_ , 2021, pp. 9116–9126.
  * [37] M. Patacchiola, J. Turner, E. J. Crowley, M. O’Boyle, and A. J. Storkey, “Bayesian meta-learning for the few-shot setting via deep kernels,” _Advances in Neural Information Processing Systems_ , vol. 33, pp. 16 108–16 118, 2020.
  * [38] S. Nabi, H. Nassif, J. Hong, H. Mamani, and G. Imbens, “Bayesian meta-prior learning using empirical bayes,” _Management Science_ , vol. 68, no. 3, pp. 1737–1755, 2022.
  * [39] S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” _Proceedings of the national academy of sciences_ , vol. 113, no. 15, pp. 3932–3937, 2016.
  * [40] V. Dhiman, M. J. Khojasteh, M. Franceschetti, and N. Atanasov, “Control barriers in bayesian learning of system dynamics,” _IEEE Trans. Automatic Control_ , 2021.
  * [41] L. Wang, E. A. Theodorou, and M. Egerstedt, “Safe learning of quadrotor dynamics using barrier certificates,” in _Proc. Int. Conf. Robotics and Automation (ICRA)_ , 2018, pp. 2460–2465.
  * [42] R. Vershynin, _High-dimensional probability: An introduction with applications in data science_. Cambridge university press, 2018, vol. 47.