# New links between invariant dynamical structures and uncertainty
quantification.

Guillermo García-Sánchez  Instituto de Ciencias Matematicas,CSIC, C/Nicolas
Cabrera 15, Campus Cantoblanco, Madrid, 28049, Spain  Escuela Técnica Superior
de Ingenieros de Telecomunicación, Universidad Politécnica de Madrid, Av.
Complutense, 30, Madrid,28040,Spain  Ana Maria Mancho  Instituto de Ciencias
Matematicas,CSIC, C/Nicolas Cabrera 15, Campus Cantoblanco, Madrid, 28049,
Spain  Makrina Agaoglou  Departamento de Matemática Aplicada a la Ingeniería
Industrial, Escuela Técnica Superior de Ingenieros Industriales, Universidad
Politécnica de Madrid, c/ José Gutiérrez Abascal, 2, Madrid, 28006,Spain
Stephen Wiggins  School of Mathematics, University of Bristol, Fry Building,
Woodland Road, Bristol, BS8 1UG, United Kingdom.  Department of Mathematics,
United States Naval Academy, Chauvenet Hall, 572C Holloway Road, Annapolis, MD
21402-5002

###### Abstract

This paper proposes a new uncertainty measure, appropriate for quantifying the
performance of transport models in assessing the origin or source of a given
observation. It is found that in a neighbourhood of the observation the
proposed uncertainty measure is related to the invariant dynamical structures
of the model. The paper illustrates the implementation of the proposed
definition to quantify the performance of ocean data sets in the context of a
real oil spill event in the Eastern Mediterranean in 2021.

##  I Introduction

Models in real-world applications are subjected to uncertainty in their
results, and this is particularly notorious in models that describe the ocean
or the atmosphere state. In these geophysical contexts, the availability of
tools that quantify the uncertainty of model outputs is very much desired as
this allows the evaluation of their reliability and the assessment of their
use in real applications. Uncertainty quantification is a very broad question
and it has been addressed from different perspectives. One approach to this
problem involves performing a sensitivity analysis, which determines what
inputs in the model, which typically is a partial differential equation,
affect the outputs the most. A related issue, referred to as forward
propagation of uncertainty, consists of examining how uncertainty in the model
parameter inputs affects uncertainty in the model outputs. In this framework
it may happen that experimental data are available. In this case, it is
possible to infer the set of parameters that better fit the outputs for the
proposed model. This setting can be precisely formulated as a Bayesian inverse
problem, which for a given observation, and an assumed noise model, determines
the parameters that are most likely to have produced the data. Typically,
Bayesian analyses provide the framework for many inverse uncertainty
quantification applications [1, 2]. These methods employ Bayesian inference
theory along with exploring the posterior probability density function (PDF)
by Markov Chain Monte Carlo (MCMC) sampling. Other methods include statistical
methods based on maximum likelihood estimation (MLE) [3] and data assimilation
methods. Data assimilation methods combine experimental observations with code
predictions and their respective errors to provide an improved estimate of the
system state and of the associated uncertainty [4]. All of these approaches
assume that the model is known, with undetermined parameters, and perform
adaptations on the model inputs. However, there is a different type of problem
for which only an approximate model is available and its exact expression is
unknown. We delve into an explanation of such problems.

Recently, García-Sánchez et al. [5] have proposed to quantify uncertainty
associated with transport by currents in settings in which the exact model
connecting two successive observations (or the currents producing the observed
transport) was not known. They have quantified uncertainty motivated by the
need to judge the performance of ocean model outputs to describe oil spill
events such as that described in [6]. In this reported event an identified
spill occurs at a certain time and it is required to describe its evolution to
determine if it will affect critical areas, and there is a need to assess
ocean models to this end. García-Sánchez et al. [5] approach evaluates the
reliability of the ocean currents with respect to the transport that they
produce. Indeed, transport in the ocean surface is produced by fluid parcels
that follow trajectories ${\bf x}(t)$ that evolve according to the dynamical
system:

$\frac{d{\bf x}}{dt}={\bf v}({\bf x},t),$ (1)

This system is a nonlinear non-autonomous dynamical system in which there are
uncertainties in the velocities ${\bf v}({\bf x},t)$, because they are not
exactly known given that they come from solving highly sensitive partial
differential equations. This is what we mean by not knowing the exact problem:
the right-hand side of Equation (1) is only an approximation. Typically, in
other problems where uncertainty has been studied, the system obeys an
explicitly known partial differential equation where the unknowns are not
within the model itself but rather in adjustable parameters of the equation or
in the observations. However, our case is of a different nature because the
right-hand side of Eq. (1) is a key part of the equation and not exactly
known. The uncertainties in the velocities ${\bf v}({\bf x},t)$ produce
uncertainties in the solutions ${\bf x}(t)$. García-Sánchez et al. [5]
quantified this type of uncertainty by means of error measures in settings
where the initial observation was known, ${\bf x_{0}}$, and a target state at
a later final time, ${\bf x^{*}}$. This approach was adequate to quantify the
suitability of the model given in Eq. (1) to represent this sequence of
observations in forward time. In doing so they found links between the
uncertainty and the stable invariant manifolds of hyperbolic trajectories
present in the system (1).

This paper aims to complete the described picture as follows. While results in
[5] consider the uncertainty of a model based on ocean currents in order to
describe the final fate of an observation as it evolves forward in time, the
new setting considers the uncertainty of a similar model in order to identify
the original location of an observation in a past time, once the position,
${\bf x_{1}}$, of this observation is known at the present time. In a
nutshell, previous results were about: where do trajectories end up? How could
it be measured if they are consistently described by a model constructed from
approximate velocities? Now the results we desire are about: where do
trajectories come from? how could it be measured if their origin is
consistently described by a model that is constructed from approximate
velocities? In this context, this paper provides a new approach to measuring
uncertainty using errors adequate for the backward time setting. We find that
similar to what was found in García-Sánchez et al. [5], this uncertainty has a
structure in a neighborhood of the observation ${\bf x_{1}}$, and this
structure is related to the unstable invariant manifolds of the hyperbolic
trajectories present in the vector field of equation (1).

Once this new type of uncertainty is introduced, we proceed to assess ocean
data capabilities in real events that demand this type of perspective. In
particular, a recent contamination event that occurred in 2021 in the Eastern
Mediterranean is considered. Specifically in this event, thoroughly described
in [7], several beaches were affected by a spill of unknown origin. In
contrast to the results by [5], now observations are done at the arrival point
and the identification of their source is the goal. Equipped with the new tool
we describe in detail its use to quantify the performance of different ocean
data sets for describing observations adequately. There are other attempts to
identify the best dataset for predicting the drift of oil spills [8, 9].
However, their approach consists of using uncertainty estimators in forward
time.

In summary, the goal of this article is to propose a new uncertainty measure
appropriate for quantifying the suitability of a model to describe the origin
of trajectories. It is found that the considered uncertainty has a structure
correlated with the unstable manifolds of the hyperbolic trajectories of the
model. The capacity of this new measure is exploited for characterizing the
performance of ocean data sets in real events. Accordingly, the structure of
this article is as follows. Section 2 provides a new definition for
uncertainty quantification suitable to describe the new setting and, in a
simple example, discusses links with unstable manifolds of hyperbolic
trajectories. Section 3 presents the application of this tool to the mentioned
oil spill event that affected the Eastern Mediterranean, and shows how it
helps to discriminate between two data sets, and supports the identification
of the most suitable one for searching the spill source. In Section 4, the
paper concludes with a discussion of the implications and potential
applications of this approach.

Figure 1: Panels a), b) and c) display a graphical representation of two
sequential observations and their forward time evolution. a) The initial
observation at time $t_{0}$ is expressed by the red circle with initial
condition ${\bf x_{0}}$. The final observed state ${\bf x^{*}}$ at time
$t^{*}$ is referred to as the "target". In general, the evolution law linking
these observations is unknown, but in our setting, it is approximated by Eq.
(2), a system with a hyperbolic fixed point at the origin with stable and
unstable manifolds in blue and red, respectively. The forward evolution of
${\bf X_{0}}$ (the pink blob) to ${\bf X^{*}}$ (the light blue blob) according
to this model is illustrated; b) a representation of $L_{UQ}$, in a domain
beyond the neighbourhood $X_{0}$ for $\tau=3$, where singular features are
identified aligned with the stable manifold (in blue); c) a representation of
the stable manifold. Black arrows point out the forward time evolution for
initial condition on the stable manifold towards the fixed point at the
origin; Panels d), e) and f) display a graphical representation of two
sequential observations and their backward time evolution, d) the final
observation at time $t_{1}$ is expressed with the red circle at the position
${\bf x_{1}}$. Although both the origin and the model of this observation are
unknown, its evolution is approached by Eq. (2), a system with a hyperbolic
fixed point at the origin with stable and unstable manifolds in blue and red
respectively. The backward evolution of ${\bf X_{1}}$ (the pink blob) to ${\bf
X^{*}}$ (the light blue blob) according to this model is illustrated; e) a
representation of $L_{BUQ}$, in a domain beyond the neighborhood $X_{1}$ for
$\tau=3$, where singular features are identified aligned with the unstable
manifold (in red); f) a representation of the unstable manifold. Red arrows
point out the evolution in backward time towards the fixed point for initial
conditions on the unstable manifold. Black arrows point out the evolution
forward time moving away from the fixed point. The color bar at the top
indicates a color scale for $L_{UQ}$ and $L_{BUQ}$ in panels b) and e).

##  II Uncertainty Quantification in Backward Time

This section is focused on proposing a new definition of uncertainty
quantification to describe the uncertainty in the origin or source of certain
observations based on predictions by a model constructed from velocities. This
section is limited to present the definition and to discuss its implementation
in one simple example. This simplified scenario enables us to establish
connections between the new definition of uncertainty and the unstable
manifolds of hyperbolic trajectories. The objective, however, is to apply this
definition to distinguish the most appropriate velocity field among various
options, the one that best matches the observations.

We start by recalling results reported in [5] in the context of the simple
example illustrated in Figure 1. Panel a) within this figure shows an initial
observation, ${\bf x}_{0}$, at time $t_{0}$ marked with a red circle. At a
later time, $t^{*}$, this observation is at a target position, ${\bf x}^{*}$,
marked with a red asterisk. This evolution is assumed to be following the
model provided by Eq. (2), with a vector field in the background, which is
represented at time $t_{0}$, in the neighborhood of ${\bf x}_{0}$ and at time
$t^{*}$ in the neighborhood of ${\bf x}^{*}$. This vector field is chosen to
be stationary, therefore its representation at later times, $t^{*}$, is the
same as at $t_{0}$. However, in general, this vector field will be time-
dependent and therefore its representation at $t^{*}$ is different from that
at $t_{0}$. The explicit expression for this vector field is:

$\begin{cases}\frac{dx}{dt}=x,\\\ \frac{dy}{dt}=-y,\end{cases}$ (2)

This system has a stable manifold aligned with the vertical axis at $x=0$
(blue line in Figure 1) and an unstable manifold aligned with the horizontal
axis at $y=0$ (red line in figure 1).

García-Sánchez et al. [5] proposed several expressions to measure the error in
a neighborhood ${\bf X}_{0}$ of ${\bf x}_{0}$. Among those we consider this
one:

$L_{UQ}({\bf
X}_{0},t_{0},\tau,p)=\sum_{i=1}^{n}\left|x_{i}(t_{o}+\tau)-x_{i}^{*}\right|^{p},\
p\leq 1,\ {\bf X}_{0}\in\mathbb{R}^{n},\tau>0.$ (3)

Here, ${\bf X}_{0}$ represents a set of initial conditions whose coordinates
$(x_{1},x_{2},..,x_{n})$, belong to a $n$ dimensional system (1) and are in
the neighbourhood of the observation ${\bf x}_{0}$. The coordinates of the
target, towards which the system is assumed to evolve in a time interval
$\tau$, are $(x_{1}^{*},x_{2}^{*},..,x_{n}^{*})$. Panel a) in Figure 1
illustrates the evolution from the set ${\bf X}_{0}$ (the pink blob) to ${\bf
X}^{*}$ (the light blue blob). Panel b) displays $L_{UQ}$ evaluated in the
whole domain. This representation expresses how much uncertainty would be
associated with the system (2) beyond the neighborhood ${\bf X}_{0}$, assuming
that initial observations start at different points of the representation
domain and that they evolve towards the assumed observed target position at
the red asterisk in $\tau=3$. It is clear from the representation, that
$L_{UQ}$ presents a structure that consists of singular features aligned with
the stable manifold of the hyperbolic fixed point present in the vector field,
as discussed and proven in [5].

Figure 2: A representation in the plane of $L_{BUQ}$ for the dynamical system
(2) with target $x^{*}=(0.5,1)$. a) $\tau=1$; b) $\tau=2$; c) $\tau=15$. The
red line represents the position of the unstable manifold.

In this paper, our interest is in quantifying the uncertainty of a model for
identifying a target source ${\bf x}^{*}$, which is consistent with a later
observation ${\bf x}_{1}$ at time $t_{1}$. In this way, the target source,
${\bf x}^{*}$, is located at an earlier time $t^{*}=t_{1}-\tau$. Similarly to
the previous definition, for this setting, we propose the uncertainty to be:

$L_{BUQ}({\bf
X}_{1},t_{1},\tau,p)=\sum_{i=1}^{n}\left|x_{i}(t_{1}-\tau)-x_{i}^{*}\right|^{p},\
p\leq 1,\ {\bf X}_{0}\in\mathbb{R}^{n},\tau>0.$ (4)

where now ${\bf X}_{1}$ represents a neighbourhood around the observation
${\bf x}_{1}=(x_{1},x_{2},..,x_{n})$ at time $t_{1}$ being $n$ the dimension
of the system (1), which in our specific examples is $n=2$. The coordinates of
the target are $(x_{1}^{*},x_{2}^{*})$, and uncertainty is provided in terms
of a distance metric between the target and the backward evolution of points
near to ${\bf x}_{1}$ for a period $\tau$.

Panel d) in Figure 1 illustrates the backward in time evolution of ${\bf
X}_{1}$ (the pink blob) to ${\bf X}^{*}$ (the light blue blob) In the
background, the panel also illustrates the vector field in a neighborhood of
${\bf x}_{1}$ at time $t_{1}$, and in the neighborhood of ${\bf x}^{*}$ at
time $t^{*}$. Given that this vector field is stationary, its representation
at an earlier time, $t^{*}=t_{1}-\tau$, will be the same as at $t_{1}$,
although in general this vector field is time-dependent and therefore its
representation at $t^{*}$ is different from that at $t_{1}$.

Panel e) displays in the whole domain, that is beyond the neighborhood ${\bf
X}_{1}$, $L_{BUQ}$ evaluated for $p=0.5$, assuming $\tau=3$. This
representation expresses how much uncertainty is associated to the system (2),
assuming that observations ${\bf x}_{1}$ could be anywhere in the domain and
that they have evolved there from the target position ${\bf x}^{*}$ at
$t^{*}=t_{1}-\tau$. It is clear from the representation that $L_{BUQ}$
presents a structure that consists of singular features aligned with the
unstable manifold of the hyperbolic fixed point present in the vector field. A
formal proof of this observation is given in the appendix. In the proof it is
shown how a sufficiently large $\tau$ is required for the singular structures
obtained from Eq. (4) to become aligned with the unstable manifolds. This is
more clearly visualized in Figure 2. Panel a) shows the evaluation of
$L_{BUQ}$ for $\tau=1$, where clearly the alignment with the unstable manifold
is not yet achieved. Panels b) and c) confirm the convergence of the alignment
towards the unstable manifold for, respectively, $\tau=2$ and $\tau=15$. This
result completes the link between invariant dynamical structures and
uncertainty quantification.

Panels c) and f) in Figure 1 graphically illustrate why the forward time
uncertainty for system (2) is related to the stable manifold of the hyperbolic
fixed point, while the backward time uncertainty is linked to the unstable
manifold of that fixed point. Indeed, the forward time integration for initial
conditions that stay on the stable manifold, for a sufficiently large $\tau$,
approach the fixed point at the origin, maintaining a fixed distance to the
target, while all other initial conditions in the plane will increase their
distance to the target for increasing $\tau$. This is the origin of the
singular features aligned with the stable manifold in forward time.
Analogously, the backward time integration of initial conditions that stay on
the unstable manifold, for a sufficiently large $\tau$, approaches the fixed
point at the origin, maintaining a fixed distance to the target. All other
initial conditions in the plane will increase their distance to the target for
increasing backward integration. This is the origin of the singular features
aligned with the unstable manifold in backward time. Analogously this
explanation also illustrates why both $L_{UQ}$ and $L_{BUQ}$ take minimum
values along the unstable manifold. Indeed, they express the finite distance
between the fixed point and the target, while the value of $L_{UQ}$ and
$L_{BUQ}$ diverge for growing $\tau$ for all other points in the plane.

Figure 3 illustrates further these ideas. Panel a) represents the evolution
versus $\tau$ of $L_{UQ}$ at points: (0,2) in red – which is on the stable
manifold– and (0,0) in black –the hyperbolic fixed point– as a function of
$\tau$, both with a target ${\bf x}^{*}=(1.5,0.2)$ and $p=0.5$. There is a
transient behavior in $\tau$ for the point on the stable manifold, while in
the limit of large $\tau$ both values converge to $|1.5|^{0.5}+|0.2|^{0.5}\sim
1.672$, which is the distance, in the defined metric, between the fixed point
and the target. Similarly, panel b) represents the evolution of $L_{BUQ}$ at
points: (2, 0) in red –which is on the unstable manifold– and (0,0) in black
–the hyperbolic fixed point– as a function of $\tau$, both with a target ${\bf
x}^{*}=(0.5,1)$ and $p=0.5$. In the limit of large $\tau$, both converge to
the constant $\sim 1.71$.

Figure 3: For the system given by Eq. (2), a) evolution of $L_{UQ}$ versus
$\tau$ on the fixed point (black line) and on a point on the stable manifold
(red line) with target ${\bf x}^{*}=(1.5,0.2)$; b) evolution of $L_{BUQ}$
versus $\tau$ on the fixed point (black line) and on a point on the unstable
manifold (red line) with target ${\bf x}^{*}=(0.5,1)$.  Figure 4: A graphical
representation of the spills observed along the coastline of the Eastern
Mediterranean and satellite observations matching the sources.

##  III Results

The backward-time uncertainty definition provided in Eq. (4), or its analog
for the Euclidean distance given by:

$L_{BUQ}({\bf
X}_{1},t_{1},\tau)=\left[\sum_{i=1}^{n}\left|x_{i}(t_{1}-\tau)-x_{i}^{*}\right|^{2}\right]^{\frac{1}{2}},\
{\bf X}_{1}\in\mathbb{R}^{n},\tau>0$ (5)

are very suitable for assessing the performance of ocean models for describing
events where the source of a certain observable, for instance, a contaminant,
is unknown. This section describes how to implement the use of the definitions
given in (4) or (5) to assess and determine the model that better achieves the
observations in a recent oil spill of this type. This event is of a different
nature to those described in [6], where after a reported oil spill accident
the interest was in describing and controlling its evolution in order to
minimize its impact. In contrast, now we are interested in assessing models
that provide information about when and where did the spill originate. Good
models in the setting described in Figure 1d) are expected to have small
distances between ${\bf X^{*}}$ (the blue domain) and ${\bf x^{*}}$ (the
target) and we will exploit this idea to measure the quality of available
ocean data sets.

We describe next the specific event to be considered under this perspective.
At the beginning of 2021, the coastline of several Middle Eastern countries in
the Eastern Mediterranean was affected by the presence of oil from an unknown
source(s) [10]. No accidents were reported for these spills previous to their
arrival at the coast, and therefore there was no hint of their origin. Israeli
authorities estimated that more than 1000 tons of tar [11] landed along 180 km
of the Israeli and Lebanese shoreline [12, 13, 14, 15, 16, 17] in mid-
February. They combined the use of satellite images with simulations by
Garcia-Sánchez et al. [7], which provided possible sources for this spill.
Their work studied this event both from the remote sensing and modeling
perspectives. Figure 4 illustrates some affected geographical areas. In this
figure, orange and green colors along the coastline mark the oil arrivals. The
systematic search of satellite images identified potential sources of these
observations and a careful examination concluded that those sources surrounded
by green and orange polygons, which correspond to satellite observations, were
located in position and time in a manner consistent with the position and
dates of the oil that subsequently arrived at the coast. According to [7],
orange polygons labeled as A and B, are successive earlier locations to the
orange arrival labeled as C. Analogously green polygons labeled as D and E,
correspond to the green arrival labeled as F. The successive dates for each
observation are given in the figure.

Simulations in [7] demonstrated that the spill was closely following fluid
parcel trajectories ${\bf x}(t)$ that obeyed the equation:

$\begin{cases}\dfrac{d\lambda}{dt}=\dfrac{u(\lambda,\phi,t)}{R\cos\phi}\\\\[8.5359pt]
\dfrac{d\phi}{dt}=\dfrac{v(\lambda,\phi,t)}{R}\end{cases}$ (6)

where the position of the fluid parcel at the ocean surface is given in
longitude ($\lambda$) and latitude ($\phi$), and $R$ is the Earth’s radius.
The system (6) is the particular form for our study of the general expression
given by Eq. (1), where ${\bf x}=(\lambda,\phi)$ and its right-hand side
represent the two components of the vector field ${\bf v}$, which are
determined by the zonal ($u$) and meridional ($v$) velocities. These velocity
components are obtained as data sets from the Copernicus Marine Monitoring
Environmental Service (CMEMS).

Equation (6) is a nonlinear non-autonomous dynamical system in 2D and it is
required for implementing the backward evolution of fluid parcels in
expression (5). $L_{BUQ}$ in Eq. (5) will have different evaluations for
different CMEMS data sets. The goal of this section is to use these outputs to
discriminate models. CMEMS provides two data sets describing currents in the
Mediterranean Sea. One is the operational Mercator global ocean analysis and
forecast system that provides 10 days of 3D global ocean forecasts updated
daily. This product includes daily and monthly mean fields of variables such
as temperature, salinity, currents, sea level, mixed layer depth, and ice
parameters from the top to the bottom of the global ocean. It also includes
hourly mean surface fields for sea level height, temperature, and currents.
The global ocean output files are displayed with a 1/12 degree horizontal
resolution with an equirectangular longitude/latitude projection and 50
vertical levels ranging from 0 to 5500 meters.

The second data set is the physical component of the Mediterranean Forecasting
System (Med-Currents), which is a coupled hydrodynamic-wave model, including
tides, implemented over the whole Mediterranean Basin, with a higher
resolution than the global model. Indeed, its horizontal grid resolution is
1/24 degree with 141 unevenly spaced vertical levels. The model solutions are
corrected with a variational data assimilation scheme (3DVAR) of vertical
profiles of temperature and salinity along the satellite track that provides
observations of sea level anomalies (SLA).

Given that both models provide information for the same area, a natural
question is to determine whether one is more suitable than the other for
describing transport in the type of event described above. Models such as that
given in Eq. (6), possess a transport signature based on invariant manifolds
associated with hyperbolic trajectories, which following Poincaré’s idea,
constitute geometrical structures on the ocean surface that organize particles
schematically by regions corresponding to qualitatively different types of
trajectories. These geometrical features allow a more robust analysis of the
transport capacity of ocean currents than that based on individual
trajectories. Indeed, there exist mathematical results that discuss the
persistence of these geometrical structures in the presence of small
perturbations, while individual trajectories may be affected much more
drastically by perturbations [18, 19]. In the past these features have been
used to assess the transport performance of data sets [20, 21, 22, 23, 24,
25]. In the context of geophysical flows, these geometrical structures are
referred to as Lagrangian Coherent Structures (LCSs). The use of LCS allows a
qualitative assessment of the performance of data sets, however, definitions
such as that given in Eq. (5) allow a quantitative analysis that we will
implement next in the context of the oil spill event described in [7].

In this work, we compute LCS using Lagrangian Descriptors [26, 27], a method
based on computing the function defined as:

$M(\mathbf{x}_{0},t_{0},\tau)=\int^{t_{0}+\tau}_{t_{0}-\tau}||\mathbf{v}(\mathbf{x}(\mathbf{x}_{0},t),t)||\;dt,$
(7)

where the Euclidean norm, $||.||$, of the vector field $\mathbf{v}$ is
evaluated along fluid parcel trajectories $\mathbf{x}(\mathbf{x}_{0},t)$ of
Eq. (6), and therefore the integral in Eq. (7) evaluates the arclength of this
trajectory on the latitude-longitude plane at the ocean surface. In expression
(7), the forward integration from $t_{0}$ reveals by means of singular
features the stable manifolds associated to hyperbolic trajectories, whereas
the backward integration highlights their unstable manifolds. This is
explicitly expressed by splitting the integral (7) into two terms as follows:

$M(\mathbf{x}_{0},t_{0},\tau)=\int_{t_{0}-\tau}^{t_{0}}\mathbf{v}(\mathbf{x}(\mathbf{x}_{0},t),t)\;dt\;+\int_{t_{0}}^{t_{0}+\tau}\mathbf{v}(\mathbf{x}(\mathbf{x}_{0},t),t)\;dt=M^{(b)}(\mathbf{x}_{0},t_{0},\tau)+M^{(f)}(\mathbf{x}_{0},t_{0},\tau)$
(8)

a)  b)

c)  d)

Figure 5: Evaluation on the 16th of February 2021 of $L_{BUQ}$ using the
target, ${\bf x^{*}}=(34.36^{\circ}\textrm{E},31.78^{\circ}\textrm{N})$ and of
$M^{(b)}$ using $\tau=16$ days; a) $L_{BUQ}$ with the CMEMS global product; b)
$M^{(b)}$ with the CMEMS global product; c) $L_{BUQ}$ with the CMEMS
Mediterranean product; d) $M^{(b)}$ with the CMEMS Mediterranean product.

We explore next the implementation of the ideas discussed in the previous
section, to evaluate the performance of the two CMEMS products in the context
of this environmental damage and their links to the unstable manifolds of
hyperbolic trajectories present in the system Eq. (6). Following oil
detections displayed in figure 4, we start the discussion considering the
first arrival to the coast on the 16th of February 2021, which is the green
mark at the Israeli coast (F) and is considered to be the observation ${\bf
x_{1}}$. For it, we consider that the position (D) on the 31st of January 2021
is its backward time target ${\bf
x^{*}}=(34.355734^{\circ}\textrm{E},31.777414^{\circ}\textrm{N})$. Figure 5a)
shows the evaluation of $L_{BUQ}$ according to Eq. (5) in a large area,
considering that the eastward ($u$) and northward ($v$) velocity components
are provided by the CMEMS global model. It is remarkable that its structure
presents important analogies with the $M^{(b)}$, visible in panel b), which
has been computed for $\tau=16$ days. This image highlights, by means of
singular features, the unstable manifolds of hyperbolic trajectories present
in the vector field supplied by Eq. (6). This result confirms that the
uncertainty associated with the model has a structure closely related to
invariant dynamical structures. The red square in panel a) highlights a
neighborhood of the observation whose origin is conjectured to be at target
${\bf x^{*}}$. This domain ${\bf X_{1}}$ is used to characterize the
performance of the model with respect to the hypothesis, assuming that the
model is good if the target is consistent with arrivals in all this domain
${\bf X_{1}}$. The model is characterized through mean values of $L_{BUQ}$ in
this box. Panels c) and d) reproduce results for the CMEMS Mediterranean
product. It is again clear the close relationship between the structure of
$L_{BUQ}$ and the unstable manifolds of hyperbolic trajectories present in Eq.
(6).

a)  b)

c)  d)

Figure 6: Evaluation on the 19th of February 2021 of $L_{BUQ}$ using the
target, ${\bf x^{*}}=(34.32^{\circ}\textrm{E},31.35^{\circ}\textrm{N})$ and of
$M^{(b)}$ using $\tau=25$ days; a) $L_{BUQ}$ with the CMEMS global product; b)
$M^{(b)}$ with the CMEMS global product; c) $L_{BUQ}$ with the CMEMS
Mediterranean product; d) $M^{b}$ with the CMEMS Mediterranean product.

We continue the discussion considering that the observation ${\bf x_{1}}$ is
the orange mark at the Gaza coast (C) on the 19th of February 2021. For it, we
consider that the position (A) on the 25th of January 2021 is its backward
time target: ${\bf
x^{*}}=(32.319813^{\circ}\textrm{E},31.350392^{\circ}\textrm{N})$. Figure 6
displays the results both for CMEMS global and Mediterranean products. Panels
a) and c) display the outputs of $L_{BUQ}$, and the red box highlights the
neighborhood ${\bf X_{1}}$ used to characterize the performance of the model.
Panels b) and d) again confirm the similarities between the unstable manifolds
of hyperbolic trajectories in Eq. (6) and the structure of $L_{BUQ}$.

It is remarkable that for both cases examined in figures 5 and 6, $L_{BUQ}$
takes, in the neighborhood ${\bf X_{1}}$, considerably larger values for the
CMEMS Mediterranean product than for the global one. Table 1 confirms this
point by providing the averages in both domains in the column labeled as
“Remote Targets", suggesting the better performance of the global model for
these cases.

a)  b)

c)  d)

Figure 7: Evaluation $L_{BUQ}$ on the Israel and Gaza coast for the close
targets using the global and the Mediterranean data. $L_{BUQ}$ on 16th of
February 2021 using the target, ${\bf
x^{*}}=(34.72^{\circ}\textrm{E},32.59^{\circ}\textrm{N})$; a) $L_{BUQ}$ on the
Israel coast with the CMEMS global product; b) $M^{(b)}$ on the Israel coast
with the CMEMS Mediterranean product. $L_{BUQ}$ on 19th of February 2021 using
the target, ${\bf x^{*}}=(34.28^{\circ}\textrm{E},32.59^{\circ}\textrm{N})$;
c) $L_{BUQ}$ on the Gaza coast with the CMEMS Global product; d) $M^{(b)}$ on
the Gaza coast with the CMEMS Mediterranean product.

Figure 4 shows that related to the oil arrivals marked in green (F) and orange
(C), there exist additional targets ${\bf x^{*}}$ nearer and closer to the
observation ${\bf x_{1}}$, than those previously discussed. These backward-in-
time targets, labeled (B) and (E), were identified from satellite imagery
between the dates of the original sources and the dates of arrival on shore.
These spills, after several days have elapsed since their leak, spread out on
the ocean surface due to the ocean movement, which stretches and folds them in
a chaotic dynamic. Indeed, this description matches well the appearance of
targets (B) and (E), as discussed in [7]. Now we consider the evaluation of
$L_{BUQ}$ using these new ${\bf x^{*}}$ to assess both models. We start
considering the green observation (F) (${\bf x_{1}}$) on the Israeli coast.
For it, central coordinates of the spill contour on the 12th of February 2021
are ${\bf x^{*}}=(34.72008^{\circ}\textrm{E},32.59311^{\circ}\textrm{N})$.
Other target positions within the extended spill could have been considered,
but they do not lead to substantial changes either in the results or in the
discussions. Panels a) and b) in figure 7 display the evaluation $L_{BUQ}$ for
this target, considering the CMEMS global and Mediterranean products
respectively. Regarding the orange observation (C) coordinates that correspond
to its backward-in-time target on the 12th of February 2021 are ${\bf
x^{*}}=(34.2891^{\circ}\textrm{E},32.59311^{\circ}\textrm{N})$. As before,
given that this is a rather spread contour, other target positions could have
been considered, but we do not discuss them as they provide no changes in the
conclusions. For these assumptions panels c) and d) of figure 7 display
$L_{BUQ}$ considering the CMEMS global and Mediterranean products
respectively. In these pictures, the red box delimits the boundaries of the
neighborhood ${\bf X_{1}}$ used to evaluate the model performance. In this
case the column “Close Targets" in Table 1 shows the results for both models.

A discussion of these results requires highlighting the sizes involved in the
different settings. The distance between the remote targets and the Gaza
(orange) and Israeli (green) shore impact points are, respectively, $\sim$ 215
km and 127 km. These size ranges involve ocean mesoscale structures, which
seem to be better represented in the global model. The close targets are at
distances $\sim$ 28 km and 33 km respectively of the Gaza and Israel arrival
points (C) and (F). These sizes involve ocean submesoscale structures, which
should be better represented by the higher resolution model, which is the
CMEMS Mediterranean product. Indeed, this is the one that on average performs
better for the closer targets, although the global model shows a good
performance as well for the Israel case, despite its lower resolution.

| Remote Targets | Close Targets
---|---|---
Case | Global | Mediterranean | Global | Mediterranean
Israel | 0.2989 | 0.9434 | 0.1642 | 0.3375
Gaza | 0.34074 | 1.5426 | 0.7810 | 0.24666
Table 1: Averages of $L_{BUQ}$ within the ${\bf X_{1}}$ domains highlighted
with red boxes.

##  IV Conclusions

This article proposes a new definition for uncertainty quantification that
extends those recently proposed in oceanic contexts in [6, 28]. The new
definition is referred to as Lagrangian Uncertainty Quantification in Backward
time ($L_{BUQ}$), and it is suited to measure if a given transport model
consistently describes the source of an observation whose origin is unknown.

We have found that the defined quantity $L_{BUQ}$, when evaluated in the
neighborhood of an observation with respect to its backward-in-time target,
has a structure which has been linked to the unstable invariant manifolds of
the hyperbolic trajectories present in the model vector field. This link has
been rigorously proven in a simple example, and numerically verified in a real
inspired case: the oil spill accident that affected the Eastern Mediterranean
in 2021, which was studied in [7].

The new definition, $L_{BUQ}$, has been exploited to quantify the performance
of different CMEMS products to describe the sequence of events reported in [7]
regarding the oil spill accident in the Eastern Mediterranean in 2021. In this
event, some of the backward-in-time targets are very far from the impact
point, and their evolution involves mesoscale structures, which our analysis
shows are better represented in the CMEMS global model. On the contrary for
targets that are close, involving submesoscale structures, on average the
CMEMS Mediterranean product performs better, although in some cases CMEMS
global product also performs well.

Our results are set within a different framework compared to other inverse
uncertainty studies. In those studies, the model typically consists of a
dissipative partial differential equation with adjustable parameters, and the
uncertainty analysis aims to determine the best-fitting parameters based on
observations. However, in our framework, the model comprises a two-dimensional
set of volume-preserving (non-dissipative) non-autonomous differential
equations, and their exact form is unknown as they rely on approximated
velocity fields. Many classical inverse problems, such as the inverse heat
equation, are ill-posed, meaning that inferring a previous temperature
distribution from final data is highly sensitive to changes in the final data.
Although our framework differs from such cases, we still observe a high
sensitivity to the uncertainty measured around a final observation, as
demonstrated in Figures 5-7. In the neighborhood surrounding the final
observation, the uncertainty exhibits a complex structure with high
oscillations, and we have established a connection with the unstable manifolds
of hyperbolic trajectories. Thus, we illustrate how concepts from dynamical
systems are linked to uncertainty quantification, whereas a statistical
perspective has traditionally been the main focus in this context. Finally, we
have shown how to implement these concepts for the practical purpose of
discriminating the best data set performer on specific environmental
applications.

## Acknowledgements

GGS and AMM acknowledge the support of a CSIC PIE project Ref. 202250E001.
AMM, GGS and MA acknowledge the support from grant PID2021-123348OB-I00 funded
by MCIN/ AEI /10.13039/501100011033/ and by FEDER A way to making Europe. MA
acknowledges the support from the grant CEX2019-000904-S and IJC2019-040168-I
funded by: MCIN/AEI/10.13039/501100011033. AMM, GGS, and MA are active members
of the CSIC Interdisciplinary Thematic Platforms TELEDETECT. SW gratefully
acknowledges the support of the William R. Davis ’68 Chair at the United
States Naval Academy and the support of EPSRC Grant No. EP/P021123/1.

## Appendix

Our proof follows the spirit of the work by [28, 29, 30, 31]. We assume the
definition of singular features given there, by considering that these are
features of $L_{BUQ}$ (we consider for simplicity the expression for $L_{BUQ}$
given in Eq. 4) on which the transversal derivative is not defined. We prove
that, for the simple case of the linear saddle, the unstable manifold of the
hyperbolic fixed point is aligned with singular features of $L_{BUQ}$.

We study the vector field of the autonomous saddle case, where the equations
of motion are the following:

$\begin{cases}\frac{dx}{dt}=\lambda x,\\\ \frac{dy}{dt}=-\lambda
y,\end{cases}\quad\lambda>0$ (9)

The system has a unique solution for a given condition $(x_{1},y_{1})$ at time
$t=t_{1}$:

$\begin{cases}x(t)=x_{1}e^{\lambda(t-t_{1})}\\\
y(t)=y_{1}e^{-\lambda(t-t_{1})},\end{cases}\quad\lambda>0$ (10)

Notice that the origin $(0,0)$ is a hyperbolic fixed point with stable and
unstable manifolds which are the following:

$W^{s}(0,0)=\\{(x,y)\in\mathbb{R}^{2}\colon x=0,y\neq 0\\},$ (11)
$W^{u}(0,0)=\\{(x,y)\in\mathbb{R}^{2}\colon x\neq 0,y=0\\},$ (12)

Since in this paper we study the backwards uncertainty quantification method
we will consider according to (4):

$\displaystyle L_{BUQ}({\bf X}_{1},t_{1},\tau,p)$
$\displaystyle=\sum_{i=1}^{n}\left|x_{i}(t_{1}-\tau)-x_{i}^{*}\right|^{p}$
$\displaystyle=|x_{1}e^{\lambda\tilde{t}}-x^{*}|^{p}+|y_{1}e^{-\lambda\tilde{t}}-y^{*}|^{p}$
$\displaystyle=|x_{1}|^{p}|\omega-a|^{p}+|y_{1}|^{p}\omega^{-p}|1-b\omega|^{p}$
(13)

where $\tilde{t}=t_{1}-\tau-t_{1}=-\tau<0$,
$a=\frac{x^{*}}{x_{1}},b=\frac{y^{*}}{y_{1}}$ and
$\omega=e^{\lambda\tilde{t}}=e^{-\lambda\tau}$, $p\leq 1$. In order to compute
equation (Appendix) we will calculate, for simplicity, each of the terms
separately. In the first term we notice that the sign of the expression
$(\omega-a)$, for $\tau\gg 1$, will be positive when $a<0$ and could be
negative when $a>0$.

  * •

For $a<0$:

$\displaystyle|\omega-a|^{p}$ $\displaystyle=(\omega-a)^{p}$ (14)
$\displaystyle=(-a)^{p}+p\omega(-a)^{p-1}+\frac{1}{2}(p-1)p\omega^{2}(-a)^{p-2}+\frac{1}{6}(p-2)(p-1)p\omega^{3}(-a)^{p-3}+O\left(\omega^{4}\right)$
(15)
$\displaystyle=(-a)^{p}+O\left(e^{-\lambda\tau}\right),e^{-\lambda\tau}\ll 1$
(16)

  * •

For $a>0$ let be $t_{L}$ such that for $\tau>t_{L}$, we get that
$(\omega-a)<0$:

$\displaystyle|\omega-a|^{p}$ $\displaystyle=(-\omega+a)^{p}=a^{p}+p\omega
a^{p-1}+\frac{1}{2}(p-1)p\omega^{2}a^{p-2}+\frac{1}{6}(p-2)(p-1)p\omega^{3}a^{p-3}+O\left(\omega^{4}\right)$
(17) $\displaystyle=a^{p}+O\left(\omega\right)$ (18)

Therefore

$|\omega-a|^{p}=|a|^{p}+O\left(e^{-\lambda\tau}\right),\,\omega\ll 1.$ (19)

Note that $a=x^{*}/x_{1}$. Thus

$|x_{1}|^{p}|\omega-a|^{p}=|x_{1}|^{p}|e^{-\lambda\tau}-a|^{p}=|x^{*}|^{p}+O\left(|x_{1}|^{p}e^{-\lambda\tau}\right),\,e^{-\lambda\tau}\ll
1.$

For the factor $|1-b\omega|^{p}$ of the second term there exists a $t_{L}$
such that if $\tau>t_{L}$ , then $b\omega<<1$ and $(1-b\omega)>0$. This is
always the case for $b<0$ and is a plausible assumption for $b>0$, if
$b<<\omega^{-1}$. We recall that $b=y^{*}/y_{1}$ and that therefore such
$t_{L}$ exists if $y_{1}\neq 0$. In this case positiveness is guaranteed for
sufficiently large $\tau$, i.e, a Taylor series around $\omega=0$, attained if
$\tau\gg 1$ and $\tau>t_{L}$, is performed for the binomial:

$(1-b\omega)^{p}=1-bp\omega+\frac{1}{2}b^{2}(p-1)p\omega^{2}-\frac{1}{6}\omega^{3}\left(b^{3}(p-2)(p-1)p\right)+O\left(\omega^{4}\right).$

Therefore,

$\frac{1}{\omega^{p}}(1-b\omega)^{p}=\omega^{-p}-bp\omega^{(1-p)}+\frac{1}{2}b^{2}(p-1)p\omega^{(2-p)}-O\left(\omega^{(3-p)}\right).$

Finally we get:

$\displaystyle|y_{1}|^{p}\omega^{-p}|1-b\omega|^{p}$
$\displaystyle=|y_{1}|^{p}|\omega^{-p}-b|^{p}$ (20)
$\displaystyle=|y_{1}|^{p}|e^{\lambda\tau p}-b|^{p}$ (21)
$\displaystyle=|y_{1}|^{p}e^{\lambda\tau
p}+O\left(|y_{1}|^{p}be^{-\lambda\tau(1-p)}\right)$ (22)

Thus, we can approximate $L_{BUQ}$ as

$L_{BUQ}(\textbf{X}_{1},t_{1},\tau,p)\approx|x^{*}|^{p}+|y_{1}e^{\lambda\tau}|^{p}=|x^{*}|^{p}+|y_{1}|^{p}e^{\lambda\tau
p}$ (23)

the dominant term is the term $|y_{1}|^{p}e^{\lambda\tau p}$. Therefore the
unstable manifold at $y=0$, for ‘sufficiently large’ $\tau$, is aligned with a
singular feature of $L_{BUQ}$, as long as $y_{1}\neq 0$, and for any
$|y_{1}|>0$ for a sufficient large $\tau$ satisfying, $\tau>t_{L}$.

## Declaration of Generative AI and AI-assisted technologies in the writing
process

During the preparation of this work, AMM, MA used Google Translator,
Grammarly, and Chat GPT to improve the use of English in writing. After using
these tools/services, AMM, MA reviewed and edited the content as necessary and
take full responsibility for the content of the publication.

## References

  * Wu et al. [2018] X. Wu, T. Kozlowski, H. Meidani, K. Shirvan,  Inverse uncertainty quantification using the modular bayesian approach based on gaussian process, part 2: Application to trace,  Nuclear Engineering and Design 335 (2018) 417–431.
  * Domitr et al. [2022] P. Domitr, M. Włostowski, R. Laskowski, R. Jurkowski,  Comparison of inverse uncertainty quantification methods for critical flow test,  Energy (2022) 125640\.
  * de Crècy [1996] A. de Crècy, Determination of the uncertainties of the constitutive relationships in the CALTHARE 2 code, Technical Report, American Society of Mechanical Engineers, New York, NY (United States), 1996.
  * Petruzzi [2019] A. Petruzzi,  The casualidad method for uncertainty evaluation of best-estimate system thermal-hydraulic calculations,  Nuclear Technology 205 (2019) 1554–1566.
  * García-Sánchez et al. [2022] G. García-Sánchez, A. Mancho, S. Wiggins,  A bridge between invariant dynamical structures and uncertainty quantification,  Commun Nonlinear Sci Numer Simulat 104 (2022) 106016. doi:10.1016/j.cnsns.2021.106016.
  * García-Sánchez et al. [2021] G. García-Sánchez, A. M. Mancho, A. G. Ramos, J. Coca, B. Pérez-Gómez, E. Álvarez-Fanjul, M. G. Sotillo, M. García-León, V. J. García-Garrido, S. Wiggins,  Very high resolution tools for the monitoring and assessment of environmental hazards in coastal areas,  Frontiers in Marine Science 7 (2021).
  * García-Sánchez et al. [2022] G. García-Sánchez, A. M. Mancho, R. A. G., C. J., S. Wiggins,  Structured pathways in the turbulence organizing recent oil spill events in the eastern mediterranean,  Scientific Reports 12 (2022) 3662.
  * Brushett et al. [2011] B. A. Brushett, B. A. King, C. Lemckert,  Evaluation of met-ocean forecast data effectiveness for tracking drifters deployed during operational oil spill response in australian waters,  Journal of Coastal Research (2011) 991–994.
  * Zhang et al. [2020] X. Zhang, L. Cheng, F. Zhang, J. Wu, S. Li, J. Liu, S. Chu, N. Xia, K. Min, X. Zuo, et al.,  Evaluation of multi-source forcing datasets for drift trajectory prediction using lagrangian models in the south china sea,  Applied Ocean Research 104 (2020) 102395.
  * Wikipedia [2021] Wikipedia, 2021 Mediterranean Oil Spill., Wikipedia, https://en.wikipedia.org/wiki/2021_Mediterranean_oil_spill, 2021\.
  * Z. Rinat and A. Ben Zikri [2021] Z. Rinat and A. Ben Zikri, Oil Spill Off Israel’s Coast Is Its Worst Maritime Pollution in Decades, and Cleanup ’Could Take Years’, Haaretz, https://www.haaretz.com/israel-news/.premium-oil-spill-off-israel-s-coast-is-its-worst-maritime-pollution-in-decades-1.9553528, 2021\.
  * Tercatin [2021] R. Tercatin, Damage to Israeli marine environment from tar spill extreme, experts say, Jerusalem Post, https://www.jpost.com/israel-news/damage-to-israeli-marine-environment-from-tar-spill-extreme-experts-say-659716, 2021\.
  * Kaplan-Zantopp [2021] M. Kaplan-Zantopp, Israel oil spill: How did it happen and what will we do going forward?, Jerusalem Post, https://www.jpost.com/israel-news/israel-oil-spill-how-did-it-happen-and-what-will-we-do-going-forward-660242, 2021\.
  * T. Joffre and R. Tercatin [2021] T. Joffre and R. Tercatin , Israel investigates tar spill calamity, places inquiry under gag order, Jerusalem Post, https://www.jpost.com/israel-news/israel-oil-spill-disaster-investigation-details-placed-under-censor-659765, 2021\.
  * AL JAZEERA AND NEWS AGENCIES [2021] AL JAZEERA AND NEWS AGENCIES, Lebanon begins cleaning beaches after oil spill, Aljazeera, https://www.aljazeera.com/news/2021/2/27/lebanon-begins-cleaning-beaches-after-oil-spill, 2021\.
  * Rutgers Staff [2021] Rutgers Staff, Oil spill off Israel reaches south Lebanese beaches, Rutgers, https://www.reuters.com/article/us-israel-environment-oil-spill-lebanon-idUSKBN2AM19V, 2021\.
  * Polidura [2021] A. Polidura, Israel’s oil spill now affects entire Lebanese coastline, Atalayar, https://atalayar.com/en/content/israels-oil-spill-now-affects-entire-lebanese-coastline, 2021\.
  * García-Garrido et al. [2018] V. J. García-Garrido, J. Curbelo, A. M. Mancho, S. Wiggins, C. R. Mechoso,  The application of lagrangian descriptors to 3d vector fields,  Regul Chaotic Dyn 23 (2018) 551–568. doi:10.1134/S1560354718050052.
  * Haller [2002] G. Haller,  Lagrangian coherent structures from approximate velocity data.,  Phys. Fluids 14 (2002) 1851–1861.
  * Kuznetsov et al. [2002] L. Kuznetsov, M. Toner, A. D. Kirwan, C. K. R. T. Jones, L. H. Kantha, J. Choi,  The Loop Current and adjacent rings delineated by Lagrangian analysis of the near-surface flow,  J. Mar. Res. 60 (2002) 405–429.
  * Shadden et al. [2009] S. C. Shadden, F. Lekien, J. D. Paduan, F. P. Chavez, J. E. Marsden,  The correlation between surface drifters and coherent structures based on high-frequency radar data in Monterey Bay,  Deep Sea Res. II 56 (2009) 161–172.
  * Haza et al. [2007] A. C. Haza, A. Griffa, P. Martin, A. Molcard, T. M. Özgökmen, A. Poje, R. Barbanti, J. Book, P. Poulain, M. Rixen, P. Zanasca,  Model-based directed drifter launches in the Adriatic Sea: Results from the DART experiment,  Geophys. Res. Lett 34 (2007) L10605.
  * Haza et al. [2010] A. C. Haza, T. Özgökmen, A. Griffa, A. Molcard, P. M. Poulain, G. Peggion,  Transport properties in small-scale coastal flows: relative dispersion from VHF radar measurements in the Gulf of La Spezia,  Ocean Dynam. 60 (2010) 861–882.
  * Mendoza et al. [2014] C. Mendoza, A. Mancho, S. Wiggins,  Lagrangian descriptors and the assessment of the predictive capacity of oceanic data sets,  Nonlinear Processes in Geophysics 21 (2014) 677–689.
  * Garcia-Garrido et al. [2015] V. J. Garcia-Garrido, A. M. Mancho, S. Wiggins, C. Mendoza,  A dynamical systems approach to the surface search for debris associated with the disappearance of flight MH370.,  Nonlin. Processes Geophys. 22 (2015) 701–712. doi:10.5194/npg-22-701-2015.
  * Mendoza and Mancho [2010] C. Mendoza, A. M. Mancho,  Hidden geometry of ocean flows,  Phys Rev Lett 105 (2010) 038501. doi:10.1103/PhysRevLett.105.038501.
  * Mancho et al. [2013] A. M. Mancho, S. Wiggins, J. Curbelo, C. Mendoza,  Lagrangian descriptors: A method for revealing phase space structures of general time dependent dynamical systems,  Commun. Nonlinear Sci. Numer. Simul. 18 (2013) 3530–3557. doi:https://doi.org/10.1016/j.cnsns.2013.05.002.
  * García-Sánchez et al. [2022] G. García-Sánchez, A. M. Mancho, S. Wiggins,  A bridge between invariant dynamical structures and uncertainty quantification,  Commun. Nonlinear Sci. Numer. Simul. 104 (2022).
  * Lopesino et al. [2017] C. Lopesino, F. Balibrea-Iniesta, V. J. García-Garrido, S. Wiggins, A. M. Mancho,  A theoretical framework for lagrangian descriptors,  International Journal of Bifurcation and Chaos 27 (2017) 1730001.
  * Lopesino et al. [2015] C. Lopesino, F. Balibrea, S. Wiggins, A. M. Mancho,  Lagrangian descriptors for two dimensional, area preserving, autonomous and nonautonomous maps,  Communications in Nonlinear Science and Numerical Simulation 27 (2015) 40–51.
  * García-Garrido et al. [2018] V. J. García-Garrido, F. Balibrea-Iniesta, S. Wiggins, A. M. Mancho, C. Lopesino,  Detection of phase space structures of the cat map with lagrangian descriptors,  Regular and Chaotic Dynamics 23 (2018) 751–766.


* [68] P. Appell and M. Kampé de Fériet, “Fonctions hypergéométriques et hypersphériques,” Gauthier-Villars (1926) .
  * [69] H. Bateman, A. Erdélyi, and U. of Naval Research, “Bateman Manuscript Project, Higher transcendental functions,” No. v. 1 in Higher Transcendental Functions. McGraw-Hill (1953) .
  * [70] A. Bzowski, P. McFadden, and K. Skenderis, “Implications of conformal invariance in momentum space,” JHEP 03 (2014) 111, arXiv:1304.7760 [hep-th].
  * [71] H. Exton, “On the system of partial differential equations associated with appell’s function f4,” Journal of Physics A: Mathematical and General 28 no. 3, (Feb, 1995) 631. https://dx.doi.org/10.1088/0305-4470/28/3/017.
  * [72] A. Prudnikov, Y. Brychkov, and O. Marichev, “Integrals and series, vol 2,” Gordon and Breach (1992) .
  * [73] J. M. Maldacena and G. L. Pimentel, “On graviton non-Gaussianities during inflation,” JHEP 09 (2011) 045, arXiv:1104.2846 [hep-th].
  * [74] P. Creminelli, J. Noreña, and M. Simonović, “Conformal consistency relations for single-field inflation,” JCAP 07 (2012) 052, arXiv:1203.4595 [hep-th].
  * [75] D.-P. Du, B. Wang, and R.-K. Su, “Quasinormal modes in pure de Sitter space-times,” Phys. Rev. D 70 (2004) 064024, arXiv:hep-th/0404047.
  * [76] D. L. Jafferis, A. Lupsasca, V. Lysov, G. S. Ng, and A. Strominger, “Quasinormal quantization in de Sitter spacetime,” JHEP 01 (2015) 004, arXiv:1305.5523 [hep-th].


# Automated Audio Captioning using Transfer Learning and Reconstruction Latent
Space Similarity Regularization

###### Abstract

In this paper, we examine the use of Transfer Learning using Pretrained Audio
Neural Networks (PANNs) [1], and propose an architecture that is able to
better leverage the acoustic features provided by PANNs for the Automated
Audio Captioning Task [2]. We also introduce a novel self-supervised
objective, Reconstruction Latent Space Similarity Regularization (RLSSR). The
RLSSR module supplements the training of the model by minimizing the
similarity between the encoder and decoder embedding. The combination of both
methods allows us to surpass state of the art results by a significant margin
on the Clotho dataset [3] across several metrics and benchmarks.

##  1 Introduction

Automated Audio Captioning is an emerging research area [2] which aims to
describe in textual form, the events taking place in the foreground and
background of an audio clip. Automated Audio Captioning is the amalgamation of
both Audio Processing and Natural Language Processing research, much like
Automated Speech Recognition. However, unlike Automated Speech Recognition,
Automated Audio Captioning does not require an alignment between the source
audio and target text sequence. Furthermore, Automated Audio Captioning
requires the model to detect both high level events in the background and the
more intricate details that occur in shorter time-frames in the foreground.

Convolutional Neural Networks [1] are a popular choice of encoder for encoding
information from the audio input. Convolutional Neural Networks has done well
for tasks such as sound event detection, audio tagging, by being able to
extract local features using convolutional filters. This has also led to
previous work in Automated Audio Captioning also similarly using a single
convolutional neural network as the encoder. However, we argue that this is
insufficient for Automated Audio Captioning, as this task requires the model
to capture not just local events (e.g a loud honk) that occur in narrower
timeframes, but also overarching events (e.g. car engine running in the
background) that can span the whole duration of the audio. To combat this, we
propose to append the transformer encoder [4] to the convolutional encoder.
The transformer encoder consists of several layers of multi-head self
attention and feedforward layers, and many variants of the transformer have
broken the state of the art on Natural Language Processing benchmarks [5]. We
hypothesise that the combination of both the convolutional encoder and
transformer encoder, along with transfer learning, will allow the model to
perform much better.

Additionally, we propose a simple self-supervised learning objective,
Reconstruction Latent Space Similarity Regularization (RLSSR), to regularize
the training of model. As the Automated Audio Captioning is arguably less
fine-grained and more abstract in its objective compared to Automated Speech
Recognition, we argue that it is insufficient to simply use a cross attention
mechanism to bridge information from the audio encoder to the text decoder.
RLSSR provides a straightforward way to regularize the training of the decoder
by reducing the euclidean distance between the reconstruction of the decoder
layer in the latent space and the latent output of the encoder.

Our contribution are as follows:

  1. 1.

We examine the usefulness of the transformer encoder in the encoder, and show
that including the transformer encoder is instrumental for the model
performance when using transfer Learning.

  2. 2.

We introduce the novel Reconstuction Latent Space Similarity Regularization
(RLSSR) objective and show that simply using this objective allows a model
trained from scratch to achieve competitive state of the art results.

  3. 3.

We demonstrate that using both Transfer Learning and RLSSR in tandem
consistently beats state of the art results on the Clotho Dataset.

##  2 Related Work

Automated Audio Captioning has been tackled by various authors.

###  2.1 Audio Captioning Datasets

The Clotho Dataset [3] is by far the dataset most researchers have examined
and experimented on. The Clotho Dataset consists of 4981 audio sample of
duration from 15 to 30 seconds. Each audio sample are matched to 5
crowdsourced English audio annotations, or captions, that describe the events
in that audio. The Clotho Dataset was created to address the critiques of the
larger Audiocaps dataset [6], which were found to contain biases in the
annotation process. There are also niche datasets like the Hospital Scene and
Car Scene dataset [7] which are Mandarin annotated.

In this work, we focus only on the Clotho Dataset, as it is the most popular
dataset.

###  2.2 Transfer Learning

Several work has applied transfer learning to Automated Audio Captioning.
Fine-tune PreCNN Transformer[8] used a 3 stage pipeline of first pretraining
an encoder, then applying transfer learning to the pretrained model by
freezing parts of the model for training, and finally finetuning the entire
model.

AT-CNN [9], the current best performing model on Clotho, follows a similar
approach with a 2 pipeline which first pretrains the CNN10 [1] / CRNN5 [10]
encoder on an audio tagging task, and then performs transfer learning on the
pretrained model on the audio tagging task. Similarly, the Audio Captioning
Transformer [11] also pretrains the encoder on audio tagging before training
on audio tagging.

Other approaches uses specialized embeddings for transfer learning. [12]
applied transfer learning to a pretrained GPT2 model, which was pretrained in
the text modality. This approach also requires several other prerequisites,
such as a BERTScore [13] model, and several preprocessing steps for caption
retrieval. In another vein, [14] used a combination of pretrained audio
embeddings and subject-verb embeddings to generate captions.

In this work, we show that applying transfer learning using publicly available
pretrained models in conjunction with the transformer encoder allows for
significant performance improvements.

##  3 Proposed Method

We use a typical encoder-decoder model in all our experiments. In addition to
the encoder-decoder structure, there is an external Reconstruction Latent
Space Similarity Regularization (RLSSR) module. A detailed overview of the
model is shown in Figure 1.

Fig. 1: Detailed Overview of the model. The encoder consists of a CNN10 and a
Transformer Encoder, and the decoder consists of a single Transformer decoder.
The output of the CNN10, $A$, and the output embedding of the transformer
decoder, $T$, is passed to the RLSSR module, where the L1/L2 loss is used as
objective function in the RLSSR Module.

###  3.1 Encoder

The encoder consists of a 10 layer CNN (CNN10) and the encoder of the
transformer, and takes in log-mel spectrograms as input. CNN10 is adapted from
different CNN architectures availables in the Pretrained Audio Neural Networks
collection [1]. The output of CNN10 is passed to a transformer encoder. The
transformer encoder [4] consists of two layers with a hidden size of 192. We
hypothesize that the multihead attention mechanism in the transformer encoder
is instrumental in learning higher-level representations from the features
provided by the CNN. In our experiments, we show that using the transformer
encoder in transfer learning provides a huge boost in performance. Where
transfer learning is applied, only the CNN10 encoder is preloaded with
pretrained weights from, while the transformer encoder is always trained from
scratch.

###  3.2 Decoder

The decoder is a typical Transformer decoder with multihead-attention. It
takes a sequence of word tokens, along with the output from the encoder, as
decoder input and outputs the next most probable word token. In our
experiments, we used a 2 layer Transformer decoder with a hidden size of 192.
While we did experiment on decoders of different layer sizes, we found that
using 2 layers worked the best.

###  3.3 Reconstruction Latent Space Similarity Regularization Module

The aim of the Reconstruction Latent Space Similarity Regularization (RLSSR)
module is to maximize the similarity between the latent space of the decoder
and the latent space of the encoder. Intuitively, this module has to
reconstruct the latent space of the encoder output in the audio domain, $A$,
from the latent space of the decoder output in the text domain, $T$. This is
done by introducing a very small module of parameters that bridges the decoder
output to the encoder.

$\begin{array}[]{l}A_{rec}=\text{FFN}(\text{pool}(T))\\\
A_{pooled}=\text{pool}(A)\\\ L_{RLSR}=\text{sim}(A_{rec},A_{pooled})\\\
\end{array}$ (1)

The RLSSR module takes the output text embedding, $T$ from the last layer of
the decoder and the output audio embedding, $A$, from the CNN from the encoder
as input. A average max pooling is applied to $A$ to extract the most
prominent features to be reconstructed. After which, an average global pooling
is applied to standardize the length of the audio features to obtain Apooled.
The text embeddings $T$ are first passed through a linear layer and then an
average global pooling layer to obtain the reconstructed representation from
the text, Arec. Finally, a similarity metric is used to minimize the distance
between Apooled and Arec. These operations are shown in the Equation 1. While
we did try other similarity metrics such as cosine similarity, Huber loss
[15], we found that the L1 distance, L2 distance worked the best as the
similarity metric.

$\begin{array}[]{l}L=\alpha L_{ce}+\beta L_{RLLSR}\\\ \end{array}$ (2)

The model minimizes the sum of the cross entropy loss from the generated
caption, along with the loss from the RLLSR module. We tried various
combinations of $\alpha$ and $\beta$ and found that weighing both loss values
equally worked the best.

Model | BLEU1 | BLEU2 | BLEU3 | BLEU4 | ROUGEL | METEOR | CIDEr | SPICE | SPIDEr
---|---|---|---|---|---|---|---|---|---
Base (with transformer encoder) | 0.516 | 0.330 | 0.219 | 0.142 | 0.348 | 0.152 | 0.319 | 0.102 | 0.210
Base (with transformer encoder) + PANN | 0.542 | 0.363 | 0.248 | 0.163 | 0.362 | 0.161 | 0.369 | 0.108 | 0.242
Base - transformer enc | 0.523 | 0.337 | 0.227 | 0.151 | 0.353 | 0.153 | 0.332 | 0.102 | 0.211
Base - transformer enc + PANN | 0.538 | 0.350 | 0.235 | 0.153 | 0.362 | 0.158 | 0.348 | 0.106 | 0.227

Table 1: Ablation experiments to determine usefulness of the transformer
encoder. Base refers to our default model, consisting of the convolutional
encoder, transformer encoder and transformer decoder. PANN refers to loading
the pretrained weights for transfer learning.

Model | BLEU1 | BLEU2 | BLEU3 | BLEU4 | ROUGEL | METEOR | CIDEr | SPICE | SPIDEr
---|---|---|---|---|---|---|---|---|---
Base | 0.516 | 0.330 | 0.219 | 0.142 | 0.348 | 0.152 | 0.319 | 0.102 | 0.210
Base + L2 loss | 0.518 | 0.335 | 0.228 | 0.152 | 0.352 | 0.151 | 0.326 | 0.102 | 0.214
Base + L1 loss | 0.515 | 0.338 | 0.234 | 0.159 | 0.351 | 0.151 | 0.325 | 0.098 | 0.212
Base + PANN | 0.542 | 0.363 | 0.248 | 0.163 | 0.362 | 0.161 | 0.369 | 0.108 | 0.242
Base + PANN + L2 loss | 0.552 | 0.370 | 0.251 | 0.166 | 0.369 | 0.163 | 0.375 | 0.112 | 0.240
Base + PANN + L1 loss | 0.551 | 0.369 | 0.252 | 0.168 | 0.373 | 0.165 | 0.38 | 0.111 | 0.246

Table 2: Results of using the RLSSR module. L1 and L2 loss refers to the
distance metric used to optimize the RLSSR module. PANN refers to applying
transfer learning using the pretrained weights from the pretrained audio
neural network.

Model | BLEU1 | BLEU2 | BLEU3 | BLEU4 | ROUGEL | METEOR | CIDEr | SPICE | SPIDEr
---|---|---|---|---|---|---|---|---|---
Baseline [3] | 0.389 | 0.136 | 0.055 | 0.015 | 0.262 | 0.084 | 0.074 | 0.033 | 0.054
Fine-tune PreCNN Transformer [8] | 0.534 | 0.343 | 0.230 | 0.151 | 0.356 | 0.160 | 0.346 | 0.108 | 0.227
AT-CNN10 [9] | 0.556 | 0.363 | 0.242 | 0.159 | 0.368 | 0.169 | 0.377 | 0.115 | -
Base + PANN + L1 loss | 0.551 | 0.369 | 0.252 | 0.168 | 0.373 | 0.165 | 0.380 | 0.111 | 0.246

Table 3: Comparison with other state of the art. Our model consistently beats
the previous state of the art on the BLEUn, ROUGEL, and CIDEr scores.

##  4 Experimental Details

###  4.1 Data

As mentioned in Section 2.1, the experimental dataset in use is the Clotho
dataset. We try to follow as closely as possible the experimental settings of
previous work. These are the specific details. The raw audio files are first
preprocessed into log mel-spectrograms. We use 64 Mel-bands, sampling rate of
44100, FFT window length of 1024, and a hop size of 512.

###  4.2 Training and Evaluation

These are our training hyperparameters and settings. We use a batch size of 64
with a gradient accumulation steps of 4 for 200 epochs with early stopping.
Learning rate is set to $3\times 10^{-4}$ and SpecAugmentation[16] is applied
to all log mel-spectrogram inputs as an data augmentation tactic. We do not
apply label smoothing. The weights of the cross entropy loss and the
similarity loss from the RLSSR module is weighted and optimized equally. Where
transfer learning is applied, we used the weights of the CNN10 model111This
pretrained model is available publicly for download at
https://github.comqiuqiangkongaudioset_tagging_cnn which scored a mAP of 0.380
on Audioset.

For inference, we perform beam search with beam size of 4 for decoding. As the
Clotho dataset has 5 reference captions for each audio, we follow the COCO
image captioning evaluation process [17] to evaluate the generated caption.
The BLEUn scores [18] calculates the degree of n-gram overlap between the
generated and reference caption. ROUGEL [19] measures the longest common
sequence and gives a score base on the similarity. METEOR [20] is said to
better correlate with human annotators and uses a harmonic mean of unigram
precision and recall to score the generated sequence. CIDEr [21] also measures
similarity using the average cosine similarity between the candidate sentence
and the reference sentences. Like METEOR, SPICE [22] is also said to correlate
with human judgement more by using semantic propositional content and is less
reliant on n-gram overlaps. Finally, SPIDEr [23] is a linear combination of
both SPICE and CIDEr. We use all the aforementioned metrics for evaluation.

Reference |  base-enc |  base+PANN+L1
---|---|---
water flowing at a constant rate then slows down to a drip |  water is running into a sink or tub |  water is running at a constant rate in a sink
a person is arranging objects on a table and pouring water into a glass |  someone is washing dishes and dishes in a sink |  water is being poured into a container while a person speaks
a man slides the file across the wood tapping it a few times at the end |  a door is opened and closed several times |  a hard object is being dragged across a hard surface

Table 4: Table containing the generated examples of our models, along with the
reference captions.

##  5 Experimental Results

###  5.1 Usefulness of the Transformer Encoder

We compare the effectiveness of including the transformer encoder in the
encoder. The results are shown in Table 3. Experiments show that using
pretrained weights along with the transformer encoder outperformed all the
other models by a significant margin. On the other hand, not including the
transformer encoder along with pretrained weights resulted in only marginal
improvements compared to not using pretrained weights. It is also interesting
to note that, without pretrained weights, not using the transformer encoder
boosted the performance. We believe this is because the size of the Clotho
dataset is insufficient to meaningfully train a transformer encoder from
scratch. However, using pretrained weights in the convolutional encoder
provided the initialization required for the transformer encoder to learn much
better.

###  5.2 Effect of the Reconstruction Latent Space Similarity Regularization
Module

We examine the impact of the RLSSR Module. Results are shown in Table 3. From
our experiments, using the L1 loss outperforms the L2 loss consistently.
However, using either the L1 loss or L2 loss in RLSSR module outperforms the
models which did not use the RLSSR module. Models which did not use transfer
learning but used the RLSSR module outperformed the base model by a
significant margin and even performed competitively to state of the art models
using transfer learning. When using transfer learning, the performance of our
models increased even further. This attests to the RLSSR module being an
useful objective to optimize to effectively regularize the model.

###  5.3 Comparison with previous work

The comparison with other work is shown in Table 3. Our best model beats the
best performing model, AT-CNN10 [9], on all metrics except the METEOR and
SPICE metric. In addition to being self-supervised and simple to train, our
approach makes use of publicly available pretrained models, while AT-CNN10
involves a two stage transfer learning process trained from scratch. We have
also included some examples of our generated sequences at Table 4.

##  6 Conclusion

This work introduces the benefit of the transformer encoder when applying
transfer learning. We have also proposed the self-supervised Reconstruction
Latent Space Similarity Regularization objective, which acts as regularization
during training. We have shown that combining both RLSSR and transfer learning
allows us to beat state of the art on the Clotho dataset by a significant
margin across many metrics.

## References

  * [1] Qiuqiang Kong, Yin Cao, Turab Iqbal, Yuxuan Wang, Wenwu Wang, and Mark D. Plumbley,  “Panns: Large-scale pretrained audio neural networks for audio pattern recognition,” 2020.
  * [2] Konstantinos Drossos, Sharath Adavanne, and Tuomas Virtanen,  “Automated Audio Captioning with Recurrent Neural Networks,”  arXiv:1706.10006 [cs], Oct. 2017,  arXiv: 1706.10006.
  * [3] Konstantinos Drossos, Samuel Lipping, and Tuomas Virtanen,  “Clotho: An Audio Captioning Dataset,”  arXiv:1910.09387 [cs, eess], Oct. 2019,  arXiv: 1910.09387.
  * [4] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin,  “Attention is all you need,” 2017.
  * [5] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova,  “Bert: Pre-training of deep bidirectional transformers for language understanding,” 2019.
  * [6] Chris Dongjoo Kim, Byeongchang Kim, Hyunmin Lee, and Gunhee Kim,  “AudioCaps: Generating Captions for Audios in The Wild,”  in Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers), Minneapolis, Minnesota, June 2019, pp. 119–132, Association for Computational Linguistics.
  * [7] Xuenan Xu, Heinrich Dinkel, Mengyue Wu, and Kai Yu,  “Audio Caption in a Car Setting with a Sentence-Level Loss,”  arXiv:1905.13448 [cs, eess], Oct. 2020,  arXiv: 1905.13448.
  * [8] K. Chen, Y. Wu, Z. Wang, X. Zhang, F. Nian, S. Li, and X. Shao,  “Audio captioning based on transformer and pre-trained CNN,”  in Proceedings of the Detection and Classification of Acoustic Scenes and Events 2020 Workshop (DCASE2020), Tokyo, Japan, Nov. 2020, pp. 21–25.
  * [9] Xuenan Xu, Heinrich Dinkel, Mengyue Wu, Zeyu Xie, and Kai Yu,  “Investigating local and global information for automated audio captioning with transfer learning,”  in ICASSP 2021 - 2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2021, pp. 905–909.
  * [10] Heinrich Dinkel and Kai Yu,  “Duration robust weakly supervised sound event detection,”  ICASSP 2020 - 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), May 2020.
  * [11] Xinhao Mei, Xubo Liu, Qiushi Huang, Mark D. Plumbley, and Wenwu Wang,  “Audio captioning transformer,” 2021.
  * [12] Yuma Koizumi, Yasunori Ohishi, Daisuke Niizumi, Daiki Takeuchi, and Masahiro Yasuda,  “Audio Captioning using Pre-Trained Large-Scale Language Model Guided by Audio-based Similar Caption Retrieval,”  arXiv:2012.07331 [cs, eess], Dec. 2020,  arXiv: 2012.07331.
  * [13] Tianyi Zhang, Varsha Kishore, Felix Wu, Kilian Q. Weinberger, and Yoav Artzi,  “Bertscore: Evaluating text generation with bert,” 2020.
  * [14] Ayşegül Özkaya Eren and Mustafa Sert,  “Audio Captioning using Gated Recurrent Units,”  arXiv:2006.03391 [cs, eess], Jan. 2021,  arXiv: 2006.03391.
  * [15] Ross Girshick,  “Fast r-cnn,” 2015.
  * [16] Daniel S. Park, William Chan, Yu Zhang, Chung-Cheng Chiu, Barret Zoph, Ekin D. Cubuk, and Quoc V. Le,  “Specaugment: A simple data augmentation method for automatic speech recognition,”  Interspeech 2019, Sep 2019.
  * [17] Tsung-Yi Lin, Michael Maire, Serge Belongie, Lubomir Bourdev, Ross Girshick, James Hays, Pietro Perona, Deva Ramanan, C. Lawrence Zitnick, and Piotr Dollár,  “Microsoft coco: Common objects in context,” 2015.
  * [18] Kishore Papineni, Salim Roukos, Todd Ward, and Wei-Jing Zhu,  “Bleu: A method for automatic evaluation of machine translation,”  in Proceedings of the 40th Annual Meeting on Association for Computational Linguistics, USA, 2002, ACL ’02, p. 311–318, Association for Computational Linguistics.
  * [19] Chin-Yew Lin,  “ROUGE: A package for automatic evaluation of summaries,”  in Text Summarization Branches Out, Barcelona, Spain, July 2004, pp. 74–81, Association for Computational Linguistics.
  * [20] Michael Denkowski and Alon Lavie,  “Meteor universal: Language specific translation evaluation for any target language,”  in Proceedings of the EACL 2014 Workshop on Statistical Machine Translation, 2014.
  * [21] Ramakrishna Vedantam, C. Lawrence Zitnick, and Devi Parikh,  “Cider: Consensus-based image description evaluation,”  in 2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2015, pp. 4566–4575.
  * [22] Peter Anderson, Basura Fernando, Mark Johnson, and Stephen Gould,  “Spice: Semantic propositional image caption evaluation,” 2016.
  * [23] Siqi Liu, Zhenhai Zhu, Ning Ye, Sergio Guadarrama, and Kevin Murphy,  “Improved image captioning via policy gradient optimization of spider,”  2017 IEEE International Conference on Computer Vision (ICCV), Oct 2017.