Modeling a Hidden Dynamical System Using Energy Minimization and Kernel Density Estimates

TL;DR

This paper introduces a kernel density estimation method for modeling and forecasting trajectories of hidden dynamical systems on manifolds, effectively handling sparse and noisy data through energy minimization.

Contribution

It presents a novel KDE-based approach for trajectory modeling that minimizes an energy function, addressing sparse sampling and noise in hidden dynamical systems.

Findings

Estimator minimizes a specific energy function as data samples increase.

Method effectively reconstructs trajectories from sparse, noisy data.

Applicable to modeling recurrent trajectories on manifolds.

Abstract

In this paper we develop a kernel density estimation (KDE) approach to modeling and forecasting recurrent trajectories on a compact manifold. For the purposes of this paper, a trajectory is a sequence of coordinates in a phase space defined by an underlying hidden dynamical system. Our work is inspired by earlier work on the use of KDE to detect shipping anomalies using high-density, high-quality automated information system (AIS) data as well as our own earlier work in trajectory modeling. We focus specifically on the sparse, noisy trajectory reconstruction problem in which the data are (i) sparsely sampled and (ii) subject to an imperfect observer that introduces noise. Under certain regularity assumptions, we show that the constructed estimator minimizes a specific energy function defined over the trajectory as the number of samples obtained grows.

Tables

Table 1. Table 1: Table of results for cruise ship data

Ship	Search Radius	Average Error (NM)	Standard Deviation of Error (NM)	% in HDR
Dream	10	35.9	43.9	90.9
Dream	20	65.6	125.5	89.7
Dream	40	74.8	128.5	88.7
Freedom	10	93.9	66.7	71.7
Freedom	20	90.5	68.3	76.3
Freedom	40	85.9	69.9	78.9

Table 2. Table 2: A ​ P ​ E 𝐴 𝑃 𝐸 APE results for forecasts of first and second trajectories

Time	Prediction Error 1st Test	In 70% HDR?	Prediction Error 2nd Test	In 95% HDR?
0.5	0.277079	Yes	0.0748012	Yes
1.0	0.192594	Yes	0.056595	Yes
1.5	0.271879	Yes	0.333663	Yes
2.0	0.39765	Yes	0.415056	Yes
2.5	0.305548	Yes	0.50537	No
3.0	0.194307	Yes	0.272652	Yes
3.5	0.202291	Yes	0.786428	No
4.0	0.207912	Yes	1.33602	No
4.5	0.122324	Yes	1.73251	No
5.0	0.0112425	Yes	2.33992	No
5.5	0.0681667	Yes	3.02051	No
6.0	0.0993839	Yes	3.21997	No
6.5	0.128976	Yes	3.64753	No
7.0	0.139858	Yes	3.65211	No
7.5	0.333929	Yes	3.16792	No
8.0	0.107408	Yes	3.27594	No
8.5	0.274433	Yes	2.86188	Yes
9.0	2.10894	Yes	2.9664	Yes
9.5	0.0949259	Yes	3.07647	Yes
10.0	2.52883	Yes	3.28763	No
10.5	2.32441	Yes	4.15615	No

Table 3. Table 3: Error between predicted point at time t 𝑡 t and the nearest true point ignoring time for the trial in Fig.   4b

Time	Distance	Less than $ϵ=.25$
0.5	0.0748012	Yes
1.0	0.056595	Yes
1.5	0.25832	No
2.0	0.211505	Yes
2.5	0.159474	Yes
3.0	0.0399122	Yes
3.5	0.0406217	Yes
4.0	0.0733248	Yes
4.5	0.0401314	Yes
5.0	0.0283953	Yes
5.5	0.0472536	Yes
6.0	0.119545	Yes
6.5	0.250742	No
7.0	0.493783	No
7.5	0.388359	No
8.0	0.544281	No
8.5	0.397532	No
9.0	0.458194	No
9.5	0.328007	No
10.0	0.211558	Yes
10.5	0.193224	Yes

Table 4. Table 4: Selected APE APE \mathrm{APE} for Lorenz system forecast

Time	Prediction Error	In 70% HDR?
5	5.419	Yes
10	2.376	Yes
15	15.731	Yes
20	4.621	Yes
25	7.712	Yes
30	20.881	Yes
35	7.834	Yes
40	10.752	Yes
45	12.306	Yes
50	4.908	Yes
55	2.394	Yes
60	1.709	Yes
65	6.154	Yes
70	8.984	Yes
75	6.244	Yes
80	2.603	Yes
85	18.719	Yes
90	11.649	Yes
95	18.328	Yes
100	13.352	Yes

Equations

$$\{\begin{array}{rl}{\phi }_{{\mathbf{x}}_0}(t-t_0)\stackrel{\Delta }{=}\arg \underset{\phi }{\min }& {\int }_{t_0}^{t_f}\mathcal{L}(\phi ,\dot{\phi },t)dt\\ s.t.& {\mathbf{g}}_t(\phi ,\dot{\phi })\le 0\\ & {\phi }_{{\mathbf{x}}_0}(t_0)={\mathbf{x}}_0\\ & {\phi }_{{\mathbf{x}}_0}(t_f-t_0)={\mathbf{x}}_f\end{array}$$

$${\mathbf{x}}_i=\eta (\phi (t_i))+{\mathbf{ϵ}}_{t_i}$$

$$\hat{f}(x)=\frac{1}{Nh}\sum _{i=1}^{N}K\left(\frac{x-x_i}{h}\right),$$

$$\hat{f}(\mathbf{x})=\frac{1}{N{h}_1{h}_2\cdots h_m}\sum _{i=1}^N\left(\prod _{j=1}^{m}K\left(\frac{x_j-x_{i_j}}{h_j}\right)\right),$$

$$K(x)=\{\begin{array}{ll}\frac{3}{4}(1-x^2)& x\in [-1,1]\\ 0& \text{otherwise}\end{array},$$

$${\iint }_{{\mathbb{R}}^2}\mathbb{E}[\hat{f}(\mathbf{x})-f(\mathbf{x}){]}^{2}d\mathbf{x},$$

$$\begin{gathered} \hat{f}(x,y)= \\ \frac{9}{16N{h}_1{h}_2}\sum _{i=1}^N\left(1-\frac{(x-x_i{)}^2}{h_1^2}\right)\left(1-\frac{(y-y_i{)}^2}{h_2^2}\right). \end{gathered}$$

$$f(\lambda ,\varphi )=c(\kappa ,\beta {)}^{-1}\exp \left(\kappa \cos \lambda +\beta {\sin }^2\lambda \cos 2\varphi \right).$$

$$c(\kappa ,\beta )=\underset{0}{\overset{\pi }{\int }}\underset{0}{\overset{2\pi }{\int }}\exp \left(\kappa \cos x+\beta {\sin }^{2}x\cos 2y\right)\sin xdydx.$$

$$\tilde{\mathbf{x}}=\arg \underset{\mathbf{x}\in \Omega }{\max }f(\mathbf{x}),$$

$$\tilde{\mathbf{y}}=\arg \underset{\mathbf{y}\in M}{\max }f(\mathbf{y})$$

$$R(c)=\{\mathbf{x}\in X:f(\mathbf{x})\ge c\}.$$

$$c_{\alpha }=\arg \underset{c}{\max }{\text{idotsint}}_{R(c)}f(\mathbf{x})d\mathbf{x}\ge (1-\alpha ).$$

$$H=\{{\mathbf{x}}_{i_j}\in P:j\in J\}$$

$$\overline{\mathcal{P}}=\left\{{\left\{{\overline{\mathbf{x}}}_i^{(j)}\right\}}_{i=i_j}^{i_j+\overline{q}}:j\in \{1,\dots ,|H|\}\right\}$$

$$X_i=\left\{{\overline{\mathbf{x}}}_{i_j+(i-N)}^{(j)}:j\in \{1,\dots ,|H|\},i_j+i\le N\right\}$$

$$\delta ({\mathbf{v}}_i,{\mathbf{v}}_j)\stackrel{\Delta }{=}1-\frac{⟨{\mathbf{v}}_i,{\mathbf{v}}_j⟩}{\Vert {\mathbf{v}}_i\Vert \Vert {\mathbf{v}}_j\Vert }$$

$${\mathbf{v}}_N\approx \frac{{\mathbf{x}}_N-{\mathbf{x}}_{N-1}}{t_N-t_{N-1}}.$$

$$P({\mathbf{x}}_i,T)=\left\{{\mathbf{x}}_k\in P:k\ge i\wedge t_k-t_i\le T\right\}$$

$$\{\begin{array}{rl}\underset{\phi }{\min }& {\int }_{t_j}^{t_{j+1}}\hat{\mathcal{L}}(\phi ,\dot{\phi },t)dt\\ s.t.& \hat{\mathbf{g}}(\phi ,\dot{\phi })\le 0\\ & \phi (t_j)={\mathbf{x}}_j,\phi (t_{j+1})={\mathbf{x}}_{j+1}\end{array}$$

$${\hat{\mathcal{L}}}_G(\mathbf{x})\stackrel{\Delta }{=}\sum _{j=1}^N\left(1-\frac{1}{\sqrt{2\pi {\sigma }_i^2}}\exp \left(-\frac{d(\mathbf{x},{\mathbf{x}}_i{)}^2}{2{\sigma }_i^2}\right)\right)$$

$${\hat{\mathcal{L}}}_{LS}(\mathbf{x})\stackrel{\Delta }{=}\sum _{j=1}^N\frac{1}{{\sigma }_i}d(\mathbf{x},{\mathbf{x}}_i{)}^2$$

$$\{\begin{array}{rl}\underset{\phi }{\min }& {\int }_{\phi }{\hat{\mathcal{L}}}_{\ast }(\mathbf{x})d\mathbf{x}\\ s.t.& \hat{\mathbf{g}}(\phi ,\dot{\phi })\le 0\\ & \phi (t_j)={\mathbf{x}}_j,\phi (t_{j+1})={\mathbf{x}}_{j+1}\end{array}$$

$${\mathbf{x}}^{(i)}(t)={\phi }_{{\mathbf{x}}_0}(t-t_0)+{\mathbf{ϵ}}_t^{(i)}=\phi (t_i)+{\mathbf{ϵ}}_t^{(i)}$$

$$\underset{\gamma }{\min }{\int }_{t_0}^{t_f}\sum _{i=1}^{n}d{\left(\gamma (t),\mathbf{x}(t)\right)}^2\Vert {\gamma }^{\prime }(t)\Vert dt,$$

$${\gamma }^{\ast }(t)=\frac{1}{n}\sum _{i=1}^n{\mathbf{x}}^{(i)}(t)=\phi (t)+\frac{1}{n}\sum _{i=1}^n{ϵ}_t^{(i)}$$

$$\frac{1}{n}\sum _{i=1}^n{\mathbf{ϵ}}_t^{(i)}\to 0$$

$$\rho (S_1,S_2)=\inf \{d(\mathbf{x},\mathbf{y}):\mathbf{x}\in S_1,\mathbf{y}\in S_2\}$$

$${\mathbf{p}}_i\in \Omega \cup \{\mathbf{m}\in M:\rho (\partial \overline{Y},\mathbf{m})\le \Vert \mathbf{h}\Vert +\underset{1\le i\le N}{\max }\{\Vert {\mathbf{ϵ}}_{t_i}\Vert \}\}$$

$$K_{\mathbf{h}}(⟨x_1,x_2,...,x_d⟩)=K\left(\frac{x_1-x_1^j}{h_1},\frac{x_2-x_2^j}{h_2},...,\frac{x_d-x_d^j}{h_d}\right)$$

Full text

I Introduction

In this paper we propose an algorithm for modeling and forecasting a sparse, noisy, recurrent trajectory that lies entirely on a smooth Riemannian manifold embedded in an arbitrary dimensional Euclidean space. By sparse, we mean the signal may be subject to long gaps in observation; by noisy we mean the signal is sampled by an (unknown) imperfect observer. We will define recurrent precisely in the context of the underlying mathematical model, however in general we mean the trajectory visits a neighborhood (or collection of neighborhoods) infinitely often. Examples of these trajectories include vehicle (ship, plane and car) tracks, migration data (e.g., in birds, whales and sharks), and some economics data subject to seasonality (e.g., detrended annual sales).

Our approach uses a combination of kernel density estimation (KDE) and energy minimizing inference for sparse trajectory reconstruction prior to model learning. Our goal in using a KDE is to construct distribution estimators, rather than pointwise estimators with confidence intervals. That is, rather than using a traditional pointwise time series forecasting method, our objective is to generate a sequence of probability distributions that can be used to generate an optimized pointwise estimator on demand. The methods proposed in this paper will generalize to smooth Riemannian manifolds in arbitrary dimensions, however we will focus specifically on examples from compact subsets of $\mathbb{R}^{2}$ and the $2$-sphere $\mathbb{S}^{2}$ as a representation of the Earth. By way of comparison with similar work Berry et al. (2015), we will show that the proposed approach works well on the attracting set of Lorenz63 Lorenz (1963), thus illustrating the proposed approach can generalize outside the assumed smooth Riemannian manifold assumption.

I.1 Related Work

Our work extends and is related to the basic statistical problem of time series modeling. Linear and non-linear time series modeling is a well established field of statistics Box et al. (2008) and statistical process control del Castillo (2002). Basic linear regression Hogg and Tanis (2006) (Chapter 8) and non-linear regression Serber and Wild (2003) (Chapter 2) attempt to model observations $\{\mathbf{x}_{t_{i}}\}_{i=1}^{N}$ as functions of a variable $t\in\mathbb{R}$. In one dimension, Autoregressive Integrated Moving Average Models (ARIMA) extend these notions by allowing the model $\mathbf{x}_{t}$ to vary as a function of past values and past shocks Box et al. (2008) (Chapters 3 and 7). Seasonal ARIMA models extend this notion by adding seasonal periodicity Ghysels and Osborn (2001) (Chapter 5). Fractional ARIMA Hosking (1981) models add short and long range dependence, not expressible with classical ARIMA techniques. In particular, these non-linear models are better able to express persistence and anti-persistence. Finally (Generalized) Autoregressive Heteroskedastic Models (ARCH/GARCH) add heteroskedastic behavior to the error components of the time series, allowing globally stationary and locally non-stationary error terms to be analyzed Engle (2001). Many of these techniques can be extended to vector valued functions (of the type we consider). In particular, Vector Autoregressive models (VAR and VARIMA) Hacker and Hatemi-J (2008) can be used to model time series of vector valued functions. The most general models are the stochastic differential and difference equations that use Weiner and Lévy processes to model stochasticity Oksendal (2007) (Chapter 1). Kernel based approaches for forecasting stochastic dynamical systems modeled by (hidden) stochastic differential equations are considered extensively by Giannakis et al. Berry et al. (2015); Giannakis (2019). In particular, in Berry et al. (2015); Giannakis (2019) the authors use a diffusion forecasting approach. The shift map of the stochastic process is expressed in a smooth basis of eigenfunctions. This is used to estimate the semigroup solution of the unknown stochastic differential equation without specific parameter estimation.

Grid-based methods that approximate the trajectory can be employed when the space of the time series is continuous but can be partitioned into a collection of discrete grid points and the trajectory modeled as a time series of these grid points. The work in Martinerie and Forster (1992); Rabiner (1989); Streitt and Barrett (1990) describes methods of using hidden Markov models (and in the case of Martinerie and Forster (1992), dynamic programming) to identify optimal estimators for the behavior of trajectories passing through the discretized state space. Griffin et al. (2011) uses a multi-resolution grid model and a continuous time model to construct a hybrid track estimator that attempts to retain the simplicity of a grid-based model without sacrificing the accuracy of a continuous model. We note that many (but not all) of the approaches discussed are designed to generate pointwise forecasts with confidence regions, while our objective in using a KDE-based approach is to generate a sequence of probability distributions, which can be used to generate a pointwise forecast.

Forecasting dynamical systems, especially non-linear dynamical systems, is a well known problem in physics with applications to noise reduction and experimental data smoothing. Molecular trajectory modeling is considered in N uske et al. (2016, 2014) using a variational approach with user supplied basis functions. This approach is in contrast to the the standard Markov process approaches, which are more reminiscent of grid-based methods, already discussed. Kostelich and Yorke (1988); Kostelich and Schreiber (1993); Schreiber (1993) consider noise reduction in dynamical systems with Schreiber (1993) providing a fitting approach that is qualitatively similar to the work presented in this paper. Anomaly detection is considered in Cao et al. (2004) with stated goals similar to those in Pallotta et al. (2013), but applied to one-dimensional chaotic signals. Forecasting and non-linear modeling is considered in McSharry and Smith (1999); Voss (2001); Danforth and Yorke (2006); de S. Cavalcante et al. (2013); Berry et al. (2015); Giannakis (2019). In addition to this work, Li and Cheng (2010) applies stochastic hidden Markov models to fuzzy time series forecasting. Fuzzy time series forecasting is also considered in Huarng et al. (2007). Han et al. (2018) considers the problem of non-uniform state space reconstruction of chaotic time series. Chaotic time series forecasting is also considered in Shen et al. (2013), which uses an ant colony optimization algorithm to optimally embed a time series in an appropriate space. Joint continuous and discrete forecasting is considered in Ramasso et al. (2013), while outlier detection of time series is considered in Rasheed and Alhajj (2014). More recent work has applied multi-layer perceptron neural networks to time series forecasting Rahman et al. (2016).

Using a KDE for the purpose of modeling and forecasting recurrent trajectories has been studied by other authors in more restricted contexts. Pallotta, Vespe, and Bryan Pallotta et al. (2013) use a kernel density estimation technique to model shipping routes using Automatic Identification System (AIS) data. They use the resulting distributions to identify anomalous behavior in ship routes. As they note, AIS data is exceptionally dense, and can be used in real-time to track ships.

Additionally, it is well known that a kernel density estimate can be used as a convolutional filter on noisy data. This was done in Lampe and Hauser (2011) in order to visualize streaming data from an aircraft. This is a trajectory in $\mathbb{R}^{3}$, although Lampe and Hauser (2011) only considers the projection onto $\mathbb{R}^{2}$. In both Pallotta et al. (2013) and Lampe and Hauser (2011), the data are highly dense with minimal noise. This is not realistic in antagonistic situations or in cases where the trajectory cannot be observed with high fidelity. This occurs naturally when biologists observe animals in their natural habitat (e.g., see OC1 (2018)). This paper considers situations in which the sampled trajectory is neither dense nor exhibits high signal-to-noise ratio (SNR). We contrast this to the work in Giannakis (2019); Berry et al. (2015), where the data are assumed to be more dense.

I.2 Paper Organization

The remainder of this paper is organized as follows: In Section II we introduce notation and the underlying mathematical model to be used throughout the rest of the paper. In Section III we discuss our proposed modeling and forecasting algorithms. Theoretical results on the algorithms are provided in Section IV. We present empirical results using synthetic and real-world data sets in Section V. Finally conclusions and future directions are presented in Section VI.

II Notation and Preliminaries

II.1 Notation and Assumptions

Let $\mathbb{R}$ denote the real numbers, and ${\mathbb{R}^{+}=\{x\in\mathbb{R}:x\geq 0\}}$. Let $\langle{M,\mathbb{R}^{+},\varphi}\rangle$ be a (hidden) dynamical system describing the motion over time ($\mathbb{R}^{+}$) of an (autonomous) particle on a $d$-dimensional, smooth Riemannian manifold $M$. The manifold $M$ may be embedded in Euclidean space of dimension $m\leq 2d$, and such an embedding is guaranteed to exist by Whitney’s Strong Embedding Theorem. Throughout this paper, boldface symbols will indicate positions on the manifold in an appropriate coordinate system; e.g., $\mathbf{x}\in\mathbb{R}^{d}$ and $\mathbf{x}=\langle{x_{1},\dots,x_{d}}\rangle$. We will often (without explicitly stating) identify the set of points in $M$ with their image in $\mathbb{R}^{d}$ under an appropriate chart.

If $M$ is known, then our approach may be taken using $M$ itself, e.g. using the KDE theory developed in le Brigant and Puechmorel (2019). In the case $M$ is unknown, but the image of the embedding $M\hookrightarrow\mathbb{R}^{m}$ is known, our approach may be used taking $\mathbb{R}^{m}$ to be the manifold of interest (even though the data may be drawn from a different underlying manifold). Determining $M$ given data in $\mathbb{R}^{m}$ is the fundamental problem in topological data analysis Ghrist (2007); Carlsson (2009) and will not be explored further here.

Since our manifold is Riemannian, it may be endowed with an appropriate metric. For example, when $M\equiv\mathbb{R}^{d}$, the standard Euclidean metric is used; when $M\equiv\mathbb{S}^{2}$ (the 2-sphere), the Haversine metric is applicable. We denote distance between two points $\mathbf{x},\mathbf{y}$ in $M$ as $d(\mathbf{x},\mathbf{y})$. Again by our assumption of a Riemannian manifold, we have the existence of an inner product (positive-definite metric tensor), denoted $\langle\mathbf{x},\mathbf{y}\rangle$. This should not be confused with a two-dimensional vector as the entries are vectors in bold typeface rather than coordinates in standard typeface. We will use the inner product to quantify the degree to which (e.g., velocity) vectors located at $\mathbf{x},\mathbf{y}$ have similar heading. In Euclidean space, we would choose the usual dot product. Throughout the paper, we use $\overset{\Delta}{=}$ to denote equality by definition rather than derivation.

The dynamical system we study is hybrid in the following sense: Fix a finite set $\mathscr{O}\subseteq M$. At any time $t$, either:

There are positions $\mathbf{x}_{0},\mathbf{x}_{f}\in\mathscr{O}$ and the function $\varphi$ defines a sub-trajectory $\mathbf{x}_{t}=\varphi_{\mathbf{x}_{0}}(t-t_{0})$ in $M$ so that

[TABLE]

where $\mathscr{L}:(\varphi,\dot{\varphi},t)\mapsto r\in\mathbb{R}$ is a hidden energy function (Lagrangian) and $\mathbf{g}_{t}:(\varphi,\dot{\varphi})\mapsto\mathbf{b}\in\mathbb{R}^{m}$ are hidden (possibly time parameterized) constraints. 2. 2.

We have $\varphi(t)=\mathbf{x}_{0}\in\mathscr{O}$. At some time $t+\tau$, a new $\mathbf{x}_{f}\in\mathscr{O}$ is chosen (possibly at random).

We assume the dynamical system is recurrent in the sense that the choice of $\mathscr{O}$ is governed by an ergodic or periodic (hidden) Markov chain with no transient states. Therefore, if $\varphi(t)=\mathbf{x}_{0}$, there is some $T<\infty$ so that $\varphi(t+T)=\mathbf{x}_{0}$.

In the problem we discuss, all relevant information about the dynamics, including $\mathscr{O}$, the Lagrangian $\mathscr{L}$ and some (perhaps all) of the constraints $\mathbf{g}_{t}$ are hidden. The assumption that $\varphi(t)$ is constructed from piecewise optimal paths is used to justify our method of inferring missing information in sparsely sampled data. For simplicity, in the remainder of this paper, we will assume that $\mathbf{g}_{t}(\varphi,\dot{\varphi})$ are time invariant and denote the constraint functions by $\mathbf{g}$. In the sequel, we assume data are sampled discretely via a sampling function $\eta:M\rightarrow\mathbb{R}^{d}$ (or $\eta:M\rightarrow\mathbb{R}^{m}$ as appropriate) and with unbiased noise to produce a sparse noisy signal:

[TABLE]

here the $\bm{\epsilon}_{t_{i}}$ are unbiased noise vectors of appropriate dimension. In the sequel, we will elide the observation function $\eta$ for the sake of clarity and identify $\varphi(t_{i})$ with $\eta\circ\varphi(t_{i})$. Whether we are using $\varphi$ to mean a trajectory on $M$ or its image in $\mathbb{R}^{d}$ or $\mathbb{R}^{m}$ will be clear from context.

We note, our approach is a parameter-free approximation method and our focus is not on estimating the distribution that describes $\bm{\epsilon}_{t_{i}}$, unlike (e.g.) in the traditional Kalman filter estimation (see Kalman (1960)).

In addition to the recurrence assumption, we assume that $\varphi(t)$ is piecewise smooth, and in particular at $t_{0}$ an instantaneous velocity can be constructed using initial conditions. In practice velocity is numerically approximated by a difference quotient. Finally, since we assume that $\varphi_{\mathbf{x}_{0}}(t-t_{0})$ obeys a set of externally imposed constraints defined by $\mathbf{g}(\varphi,\dot{\varphi})$ in Eq. 1, we assume there is a feasible region $\Omega\subseteq M$ defined by $\mathbf{g}(\varphi,\dot{\varphi})$ and for all $t$, $\varphi(t)\in\Omega$. In particular, $M$ may be convex as a set in $\mathbb{R}^{m}$, but $\Omega$ may not be, making the problem more challenging.

II.2 Techniques

We provide a brief overview of Kernel Density Estimates (KDE) in Euclidean space, which are a foundational element of our proposed algorithm. The interested reader may consult Scott (2015) for a more detailed overview of Euclidean KDE methods, or le Brigant and Puechmorel (2019) for KDE methods on a Riemannian manifold.

The KDE is a non-parametric estimate of the probability distribution of a set of data $\{x_{i}\}_{i=1}^{N}$. Formally, the univariate KDE $\hat{f}$ is given by

[TABLE]

with the requirement that the kernel $K$ be a valid probability density function (PDF). Intuitively, one may think of a KDE as a mean of many probability distribution functions with a normalizing constant, $\frac{1}{h}$. The parameter $h>0$ is referred to as the bandwidth and there are many well established rules-of-thumb for choosing $h$ given arbitrary data (see, e.g., Jones et al. (1996)). The bandwidth controls the variance of the kernel and acts as a smoothing control on the resulting KDE.

The KDE can be generalized to arbitrary-dimensional Euclidean space. Let $\{\mathbf{x}_{i}\}_{i=1}^{N}$ be $N$ data points in $\mathbb{R}^{m}$. Define $\mathbf{x}_{i}=\langle{x_{i_{1}},\dots,x_{i_{m}}}\rangle$. The KDE generalizes to:

[TABLE]

where $K$ is a PDF defined on $\mathbb{R}^{d}$ and $h_{j}$ is the bandwidth parameter for the coordinate $x_{j}$. We will write the bandwidth vector $\mathbf{h}={\langle h_{1},h_{2},...,h_{d}\rangle}\in(\mathbb{R}^{+})^{d}$. The product defined inside the sum is called the product kernel

Choice of kernel: One consideration which must be taken into account when using KDE methods is the choice of kernel. While the only restriction on $K$ is that it is a valid PDF, there are a few canonical choices of kernel. The first, often used in image processing applications as the kernel of a convolutional filter on images, is the Gaussian/normal PDF. Moreover, the Gaussian kernel is used in both Pallotta et al. (2013) and Lampe and Hauser (2011).

Another canonical choice is the Epanechnikov kernel. Defined as

[TABLE]

the Epanechnikov kernel is a compactly supported parabola. In Epanechnikov (1969), Epanechnikov shows that the kernel minimizes relative global error in the following sense: Let $f$ be the true distribution and assume that $f$ is analytic. Suppose $f$ is approximated by $\hat{f}$ using Epanechnikov kernel. Then $\hat{f}$ minimizes the functional:

[TABLE]

where $\mathbb{E}$ is the classical expectation operator. The preceding result is derived using the calculus of variations with constraints and is independent both of distribution (subject to the requirement of analyticity) and constraints on the kernel function.

Eq. 3 defines the Mean Integrated Square Error (MISE). Therefore, the Epanechnikov kernel minimizes MISE. Scott Scott (2015) notes “…the MISE error criterion has two different, though equivalent, interpretations: it is a measure of both the average global error and the accumulated point-wise error.” The fact that the Epanechnikov kernel is a MISE minimizer makes it the natural choice of non-Gaussian kernel.

Moreover, the fact that the Epanechnikov kernel has compact support also implies that the inferred distribution has compact support (as the union of finitely many compact sets), which is sometimes desirable. In the context of our work we may implicitly constrain the velocity of the forecast by requiring the KDE have a compact support. The Gaussian kernel does not respect such a constraint, as it gives non-zero probability to the whole ambient space, rather than restricting the non-zero probability to a region that satisfies velocity constraints. As we also show in Section IV, a kernel with bounded support can also provide meaningful results regarding feasible regions.

KDE in the Plane: In $\mathbb{R}^{2}$, which is a sufficient for approximation of the earth’s surface over small regions, we may choose the Epanechnikov product kernel, given by the expression:

[TABLE]

This is our choice of kernel for the experimental results of Section V.2.

KDE on the Sphere: In Section V, we apply the proposed algorithm to cruise ships traveling on the surface of the earth. While for this paper we consider a small enough region to approximate by the Earth as a plane, on larger scales this may lead to significant distortion. In such a case one might wish to approximate the Earth as a 2-sphere $\mathbb{S}^{2}$. In this case, one could use use the Kent distribution, the analogue to the Gaussian distribution on a sphere, as the choice of kernel. Let $\lambda$ be the longitude in degrees, and $\phi$ the the latitude in degrees. The general formulation of the Kent distribution in spherical coordinates is

[TABLE]

Here $\kappa$ is a parameter representing the concentration of the distribution, $\beta$ is an analogue of covariance, which Kent describes as the “ovalness,” and $c(\kappa,\beta)$ is a normalizing constant given by

[TABLE]

We make the simplifying assumption that the covariance of $x,y$ is [math], which is equivalent to setting $\beta=0$. Then we have $c(\kappa,0)=4\pi\kappa^{-1}\sinh\kappa$, simplifying the double integral above. A full description of the Kent distribution can be found in Kent (1982). It is also possible to find a KDE on other compact Riemannian manifolds without boundary (see Pelletier (2005); le Brigant and Puechmorel (2019)). We do not further discuss this case, as Euclidean space is sufficient for our (and most other) applications.

II.3 Deriving a Point Estimator and Uncertainty Regions with a KDE

Given a PDF (potentially a KDE) $f$ on $M$, there are several ways to find a point estimator. The most obvious method is to compute the argument of expectation, $\mathbf{e}\in M$ such that $f(\mathbf{e})=\mathbb{E}[f]$. However, if $f$ is multimodal on $M$ and there are constraints on the dynamics (see Eq. 1), then it is possible $\mathbf{e}\not\in\Omega$, where $\Omega\subseteq M$ is the feasible subset of $M$. In this case it is more useful to to compute:

[TABLE]

as the point estimator. Depending on the numerical complexity, one can also define the conditional distribution $f_{\mathbf{x}|\Omega}$ and compute $\mathbf{e},\tilde{\mathbf{x}}$ accordingly. We note that Eq. 4 may be a non-convex optimization problem. In this case, it might be simpler to compute the unconstrained optimization:

[TABLE]

and either (i) accept that $\tilde{\mathbf{y}}\not\in\Omega$ or (ii) project $\tilde{\mathbf{y}}$ onto the boundary of $\Omega$. We discuss the limits of this approach in Section IV.

To find uncertainty regions, we compute the highest density region, following the Monte Carlo technique established in Hyndman (1996). Intuitively, the highest density region (HDR) of a distribution function $f:M\rightarrow\mathbb{R}^{+}$ is the subset of $M$ corresponding to the preimage of a horizontal “slice” of ${\mathbb{R}^{+}}$ where the slice includes the maximum of the PDF, and continues extending down until the probability measure of the preimage of the slice meets a predetermined threshold.

More formally, given a PDF $f$ with support $X\subseteq M$ define:

[TABLE]

The $100\%\times(1-\alpha)$ highest density region is the set $R(c_{\alpha})$ where:

[TABLE]

It is clear from this definition that $\tilde{\mathbf{y}}\in R(c_{\alpha})$ for any $0\leq\alpha<1$. If we compute the HDR of the conditional distribution $f_{\mathbf{x}|\Omega}$, then the probability measure of $M\backslash\Omega$ is zero, and so $\tilde{\mathbf{x}}\in R(c_{\alpha})$. This is also the case if the unconditional probability given to $M\backslash\Omega$ is sufficiently small.

Hyndman Hyndman (1996) proposes a Monte Carlo algorithm that samples the computed PDF (in our case the KDE) and uses the $\alpha$ quantile as an estimator for $c_{\alpha}$. We use this algorithm in the sequel.

Lastly, since we consider a sequence of KDE’s that advance in time, there is an implicit inclusion of the velocity of the object. Thus, we do not need to actually approximate or compute $\dot{\varphi}_{\mathbf{x}_{0}}(t-t_{0})$ explicitly while generating a forecast. This stands in direct contrast to alternate approaches (e.g., Giannakis (2019); Berry et al. (2015); N uske et al. (2014, 2016)), which view forecasting as finding and solving a system of stochastic differential equations describing the motion of a particle.

III Algorithm Description

In this section we motivate Algorithm 1, which forecasts a finite sequence of triples $\mathscr{F}=\{(\mathbf{p}_{i},\hat{f}_{i},R_{i})\}_{i=N+1}^{N+\overline{q}}$, with $\overline{q}\in\mathbb{N}$, where $\mathbf{p}_{i}$ estimates the value $\varphi_{\mathbf{x}_{0}}(t_{i}-t_{0})$ and $\hat{f}_{i}:M\to\mathbb{R}$ is a distribution used to construct the HDR $R_{i}$ that acts as an uncertainty region for $\mathbf{p}_{i}$. The algorithm takes as input the observation set $P=\{\mathbf{x}_{i}\}_{i=1}^{N}$, a time-indexed sequence of observed positions. Recall from Eq. 2, $\mathbf{x}_{i}$ is observed with noise. Note that we do not require $t_{i}-t_{i-1}=t_{i+1}-t_{i}$ for the points in $P$.

The algorithm is broken into four stages:

Identify a collection of points

[TABLE]

similar (defined by the metric and/or inner product on $M$) to the last known state of the particle. The set $J$ is an index set of consecutive integers which respect the time series $P$, that is $i_{j}<i_{j+1}$ for all $j$. 2. 2.

Extract sub-trajectories of $P$ corresponding to each identified point in $H$ that (i) begin at the identified point and (ii) span an input forecast time. This set of sub-trajectories is denoted $\mathscr{P}$. 3. 3.

Densify the observed sub-trajectories in $\mathscr{P}$ to obtain $\overline{\mathscr{P}}$ using a line integral minimization on an estimator $\hat{\mathscr{L}}$ of $\mathscr{L}$. Each densified trajectory in $\overline{\mathscr{P}}$ is composed of points on $M$ that are equally spaced in time. 4. 4.

Let:

[TABLE]

be the densified trajectories. For each $i\geq N$ (time index) use the set:

[TABLE]

to construct a KDE $\hat{f}_{i}$. Use the KDE to construct $\mathbf{p}_{i}$ and an associated HDR representing an uncertainty region.

Metrics and Tolerances: In Section II, we have already defined the distance $d(\mathbf{x},\mathbf{y})$ and inner product $\langle\mathbf{x},\mathbf{y}\rangle$ on the manifold $M$. Given two velocity vectors $\mathbf{v}_{i}$ and $\mathbf{v}_{j}$, let the angle metric be:

[TABLE]

This is just the standard cosine distance when $M\equiv\mathbb{R}^{d}$. In the absence of velocity data, we can use a finite difference method so that (e.g.) in $\mathbb{R}^{d}$, the velocity vector of the last observed point in $P$ is:

[TABLE]

More generally, if we are working in $M$, then we may use the language of velocity vectors on smooth manifolds (see e.g. Boumal (2010) or Chapter 3 of Lee (2013)) to approximate $\mathbf{v}_{N}$. It is in these details that we wish for $M$ to be smooth. The details of this computation obfuscate the presentation of the proposed algorithm, thus we omit them and refer the readers to the provided references.

As input to Algorithm 1, we take two parameters $\epsilon>0$, which is a a tolerance on $d(\cdot,\cdot)$ and $\vartheta$, which is a tolerance on $\delta(\cdot,\cdot)$. These function as hyper-parameters in our proposed algorithm.

Forecast Duration: Define $P(\mathbf{x}_{i},T)$ to be the forward time restriction of $P$ beginning with $\mathbf{x}_{i}\in P$ and including all points $\mathbf{x}_{k}$ so that $t_{k}-t_{i}\leq T$. That is:

[TABLE]

In Algorithm 1, $T$ is the duration of the forecast and is an input parameter.

Sampling Period: For sparse track reconstruction and forecast generation, we require a sampling frequency. The sampling frequency is a value $\Delta t$ so that if $\overline{Q}=\left\{\overline{\mathbf{x}}_{i}^{(j)}\right\}_{i=i_{j}}^{i_{j}+\overline{q}}\in\overline{\mathscr{P}}$ is a densified trajectory corresponding to some subtrajectory $Q\in\mathscr{P}$, then for all ${i\in\{i_{j},\dots,i_{j}+\overline{q}\}}$: $t_{i+1}-t_{i}=\Delta t$, where $\overline{\mathbf{x}}_{i}^{(j)}\in\overline{Q}$. This sampling period gives resolution to intermediate points of prediction but does not affect predictions made at any given point. It is now straightforward to see that $\overline{q}=\lceil T/\Delta t\rceil$.

Track Densification: Suppose $\mathbf{x}_{i}\in H$ and $P(\mathbf{x}_{i},T)$ must be densified; i.e., there is some pair $\mathbf{x}_{j},\mathbf{x}_{j+1}\in P(\mathbf{x}_{i},T)$ so that so that $t_{j+1}-t_{j}>\Delta t$. (Note: if $P(\mathbf{x}_{i},T)$ is too dense, it is trivial to downsample it to make it sparser.) If approximations $\hat{\mathscr{L}}$ and $\hat{\mathbf{g}}$ are available, then it is trivial to solve (numerically):

[TABLE]

The resulting solution can be used to provide an estimated track of arbitrary density. If $\hat{\mathscr{L}}$ is not already available, then it is straight-forward to define a Gaussian well function:

[TABLE]

or a least squares function:

[TABLE]

and solve the constrained line integral minimization problem:

[TABLE]

Here, $*$ indicates the choice of Lagrangian. Using Eq. 9 in Eq. 8 causes inferred trajectories to follow historical trends, since the center of the Gaussian well yields the minimal energy, while using Eq. 10 causes the path to minimize the square error. One benefit to Eq. 10 is it has simpler theoretical properties as we show in the sequel. On the other hand, when using Eq. 9, if $\sigma_{i}$ is an increasing function of $t_{i}$, then this is a continuous variant of pheromone routing Sim and Sun (2003); Hu et al. (2010).

Constraint inference is more difficult. In practical situations (e.g., ship tracks) there are obvious constraints in play, like land avoidance (see Section V for examples). For the remainder of this paper, we assume that the constraint function $\mathbf{g}$ (or at least $\hat{\mathbf{g}}$, and hence the feasible region $\Omega$) is known and consider constraint estimation as future work. Algorithm 1 shows the pseudo-code for the proposed algorithm.

IV Theoretical Results

If the Lagrangian $\mathscr{L}$ is stationary, we can show that the optimal solution to the problem given in Eq. 11 is asymptotically $\varphi(t)$ when $t\in[t_{0},t_{f}]$ and $\varphi(t)\equiv\varphi_{\mathbf{x}_{0}}(t-t_{0})$ as the sampling rate increases. Assume $\varphi_{\mathbf{x}_{0}}(t_{f}-t_{0})=\mathbf{x}_{f}\in\mathscr{O}$. Further assume we have $n\in\mathbb{N}$ observations of the continuous path connecting $\mathbf{x}_{0}$ with $\mathbf{x}_{f}$ denoted as $\{\mathbf{x}^{(i)}(t)\}_{i=1}^{n}$ with $t_{i}$ representing the time at which position is observed. We are considering the asymptotic case when the sampling rate is infinite (i.e., $\mathbf{x}^{(i)}$ can be thought of as a function from $[t_{0},t_{f}]\rightarrow M$), so:

[TABLE]

Assuming we use $\hat{\mathscr{L}}_{LS}$ as our estimation for $\mathscr{L}$ then we solve:

[TABLE]

subject to $\mathbf{g}$. Here $\left\|\gamma^{\prime}(t)\right\rVert$ accounts for the length of $\gamma$ on $M$ so that geodesic trajectories are preferred. At any time instant $t$, the value $\gamma(t)$ minimizing $d\left(\gamma(t),\mathbf{x}^{(i)}(t)\right)^{2}$ is:

[TABLE]

from Eq. 12. (To see this, note the integrand is simply the energy function for a mechanical equilibrium point.) We assumed $\bm{\epsilon}_{t}^{(i)}$ was unbiased. Therefore, as $n\to\infty$,

[TABLE]

and $\gamma^{*}(t)\to\varphi(t)$. We ignored constraints $\mathbf{g}(\gamma,\dot{\gamma})\leq\mathbf{0}$ only because we can see that $\varphi$ must satisfy these constraints and therefore, asymptotically so will $\gamma$. A similar argument can be made for $\mathscr{L}_{G}$, but it is not as clean, owing to the additional exponential function in the Gaussian.

Using the above results, we see that the proposed technique for filling in missing information in our discretely sampled noisy signal is (in some sense) an optimal one, assuming a stationary Lagrangian. In the case of non-stationarity, the problem is more difficult, and hence the use of heteroskedastic variances $\sigma_{i}$ (see Eq. 9) related to the time of the observation.

The inferred point predictor given in Eq. 5 is simple to implement but does not take constraints into consideration. Supposing we know the true feasible region $\Omega$, we quantify how far outside $\Omega$ a point predictor $\mathbf{p}_{i}$ could be. This can be used to determine whether it is appropriate to go through the effort of computing Eq. 4 or to simply use Eq. 5.

As before, let $\Omega\subseteq M$ be the feasible region for the trajectory $\varphi_{\mathbf{x}_{0}}(t-t_{0})$. Without loss of generality, assume $\Omega$ is a proper subset of $M$, so that the feasible region is non-trivial. Let $Y=M\setminus\Omega\neq\emptyset$ be the infeasible or forbidden region. Denote by $\overline{\,\,\cdot\,\,}$ the topological closure of a set and denote by $\partial$ the topological boundary of a set.

We show that feasible regions (and hence forbidden regions) are (partially) inferred as a part of Algorithm 1. To do this, we will use the Hausdorff distance, defined on the power set $2^{M}$ of $M$ by ${\rho:(2^{M})^{2}\to\mathbb{R}^{+}}$ by

[TABLE]

That is, $\rho$ is the smallest distance between points in $S_{1}$ and $S_{2}$. When we write $\rho$ with a set and a single point $\mathbf{x}\in M$, we will abuse notation and understand ${\rho(S_{1},\mathbf{x})\overset{\Delta}{=}\rho(S_{1},\{\mathbf{x}\})}$ so that the singleton $\{\mathbf{x}\}\in 2^{M}$.

Let $Y$ be a fixed forbidden region with closure and boundary denoted as above. Assume the prediction point $\mathbf{p}_{i}$ is computed with the unconstrained Eq. 5 and the Epanechnikov kernel $K(\mathbf{x})$ with bandwidth vector $\mathbf{h}$. Recall $d$ is the dimension of $M$, and let $\left\|\cdot\right\rVert$ be the Euclidean metric on $\mathbb{R}^{d}$. Then:

[TABLE]

In other words, the distance from any prediction point to the boundary (of the closure) of the forbidden region is at most the magnitude of the worst-case noise plus the magnitude of the bandwidth $\left\|\mathbf{h}\right\rVert$. If $\mathbf{p}_{i}\in\Omega$ this is trivial, so we consider the case when $\mathbf{p}_{i}\in Y\subseteq M$.

Let $K_{\mathbf{h}}(\mathbf{x})$ be the shifted Epanechnikov product kernel

[TABLE]

centered at $\mathbf{x}^{(j)}=\left\langle{x_{1}^{j},x_{2}^{j},...,x_{d}^{j}}\right\rangle$. Then the support of $K_{\mathbf{h}}(\mathbf{x})$ is the parallelepiped:

[TABLE]

The support of $\hat{f}_{i}$ (the $i^{\text{th}}$ estimated distribution) is the union of the supports of the individual kernels centered at $\mathbf{x}_{i}^{(j)}$ for $1\leq j\leq|H|$, hence there is some point $\mathbf{x}\in\{\mathbf{x}_{i}^{(j)}\}_{j=1}^{|H|}$ such that ${d(\mathbf{x},\mathbf{p}_{i})\leq||\mathbf{h}||}$. For the right discrete time point $t$, $\mathbf{x}=\varphi(t)+\bm{\epsilon}_{t}$ is perturbed by at most $\max_{1\leq i\leq N}\{{\left\|\bm{\epsilon}_{t_{i}}\right\rVert}\}$, then

[TABLE]

since $\varphi(t)\in\Omega$.

To maximize

[TABLE]

that is, to have conditions which allow $\mathbf{p}_{i}$ to be as far, away from the boundary of (the closure) of $Y$, while not being in $\Omega$, we need $\mathbf{x}$ to minimize $\rho(Y\backslash\Omega_{\epsilon},\mathbf{x})$. The closest that ${\mathbf{x}}$ could be to $Y\backslash\Omega_{\epsilon}$ without being in it is if ${\mathbf{x}}\in(\partial\overline{Y\backslash\Omega_{\epsilon}})\cap\Omega_{\epsilon}$, the boundary of $\overline{Y\backslash\Omega_{\epsilon}}$. This, of course, assumes that $Y\backslash\Omega_{\epsilon}$ is open - if it were closed, then $\rho(Y\backslash\Omega_{\epsilon},{\mathbf{x}})$ is strictly positive and our argument still holds. We now see that $\rho\left((\partial\overline{Y\backslash\Omega_{\epsilon}}),\mathbf{p}_{i}\right)\leq\left\|\mathbf{h}\right\rVert$, which implies Eq. 13 via the triangle inequality. Assuming that the noise is sufficiently small to allow all observations to be in the feasible region $\Omega$, then $\left\|\mathbf{h}\right\rVert$ alone serves as an upper bound on the distance inside $Y$ at which $\mathbf{p}_{i}$ may appear.

It should be noted that the essence of this argument extends to any Kernel whose support is bounded. We also note that if the topological diameter of a forbidden region $Y^{\prime}$ is smaller than $||\mathbf{h}||$, then this property does not prohibit $\mathbf{p}_{i}$ from being at any point of $Y^{\prime}$.

V Experimental Results

We discuss two sets of experiments to test Algorithm 1. In the first experiment, we forecast two cruise ships over several days (Carnival line’s Freedom and Dream ) to evaluate the performance of Algorithm 1 in a real world context. In the second experiment, we evaluate the efficacy of Algorithm 1 with three synthetic data sets. Use of a synthetic data set allows us to more closely control the underlying dynamical system and provides a method for exploring potential limitations of the proposed technique.

We use two error metrics to evaluate the algorithm: absolute pointwise error ($\mathrm{APE}$), and percent in highest density region ($\mathrm{\%HDR}$). Let $\mathbf{x}_{i}$ be the true position of particle $s$ at time $t_{i}$. We create a forecast $\mathscr{F}$ using Algorithm 1 (including information prior to $t_{i}$), and obtain prediction point $\mathbf{p}_{i}$ for time $t_{i}$. The $\mathrm{APE}$ function is defined as the distance $\mathrm{APE}(t_{i})\overset{\Delta}{=}d(\mathbf{x}_{i},\mathbf{p}_{i})$. As noted, we can construct HDR $R_{i}$ at $t_{i}$. Let:

[TABLE]

be an indicator function, then define

[TABLE]

$\mathrm{APE}$ tells us how far off the pointwise forecast is while $\mathrm{\%HDR}$ tells if the true position is in the derived uncertainty region. We compute mean and standard deviation of $\mathrm{APE}$ for an entire forecast, to give a global error metric for the forecast as a whole.

V.1 Ship Track Forecasts

Cruise ships exhibit highly recurrent behavior as they travel from port to port, according to a list of destinations which appeal to tourists. Cruise ships also use AIS to give their positions at a high sampling rate with low noise. This makes them excellent subjects on which to test our algorithm, as we can downsample a portion of a known track and generate noise to create training data. After creating a forecast from this sparse, noisy training data, we can then use the remainder of the track as high-resolution ground truth for an error analysis of the forecast.

In order to test our algorithm with this data, we used (approximately) two-years of positional data on the Carnival line cruise ships Dream (from December 2011 - July 2012) and Freedom (from March 2012 - June 2013). The data was taken from sai , under fair use. This data was densely sampled, giving a location for each ship on average about once every fifteen minutes. The first 80% of the historical trajectory of each ship was used as the historical data for “training” the KDE model and the last 20% was used as an unseen track on which to test.

We downsampled the training data (the first 80% for each ship) to give one position every day with exactly 24 hour intervals, while retaining the resolution of the unseen track. Since there were usually not AIS positions at exactly 24 hours of time-difference, we linearly interpolated between the nearest known position before the 24-hour mark and nearest known position after the 24-hour mark. The choice of linear interpolation is “wrong” in the sense that we are working on an oblate spheroid as the manifold, but served the purpose of introducing noise into the training data. The result was a sparse noisy track; this was the desired condition for our historical data.

We generated several forecasts with different parameters. In particular, we considered forecast windows of one week with 15 minute resolution, and with input search radii of 10NM, 20NM, and 40 NM. In each case we chose a bandwidth of 1.5 degrees of latitude/longitude, or 90 nautical miles. The bandwidth was chosen by trial-and-error to yield a smooth forecast. We consider the problem of automated bandwidth selection as a problem for future work.

After performing the energy minimization step of Algorithm 1, we have densely sampled paths. Error metrics for forecasting are are plotted in Fig. 1. The trajectories of cruise ships change frequently (e.g., as a result of stop over at ports of call). By way of comparison, we note that course change dynamics would have to be known a priori when using (e.g.) a Kalman Filter.

Table 1 shows summary statistics for the cruise ship forecasting experiments. For Dream, as we increase the search radius, we see an increase in the average error and standard deviation of error, and a decrease in percent in HDR. For Freedom, the average error went down from 10NM to 20NM to 40NM. The standard deviation also only marginally increased and %HDR increased significantly.

It is curious that the examples exhibit different behavior as the search radius grows. One possible explanation for this phenomenon is that for Dream with radius 10NM, the prediction in the first several hours is less than 10NM, while for Freedom with radius 10NM, there are predictions made in the first hour with error greater than 10NM. By declaring a search radius we are saying two points are in the same location if they are sufficiently close; i.e., we are creating an equivalence class on the observed data. Thus, we cannot expect our error to be smaller than our search radius. If the error is initially smaller than the search radius, then expanding the search radius would add extra data, which would contribute to a less accurate prediction. On the other hand, if the error is greater than the search radius, increasing the radius does not add data that is further away than the error, and so it might improve the forecast.

At a higher level, this example provides relevant information about the proposed method. First, there is not necessarily a single best initial search radius - it is context dependent; i.e., a parameter of the model that must be fit. Secondly, it validates our choice of HDR as uncertainty region, because our worst prediction (Freedom, 10NM), was still within the HDR $71.7\%$ of the time. Moreover, Table 1 shows that greater error and standard deviation of error corresponds to a lower $\%HDR$. It is true that the HDR depends on the bandwidth parameter, but with a properly chosen bandwidth, we have a reasonable uncertainty region.

Another positive aspect of our forecast is how it treats land. For the most part, the forecast respects the fact that it must remain in the water, and in the few cases were the forecast does go over land, it is only over small islands or tips of peninsulas, which may not be considered by the energy minimization step 111This step was carried out using a gridded representation of the Earth using a numerical approximation of the line integral.. This is important to note because we did not give as an input the location of landmasses in the statistical model. Not only does it respect navigational constraints in this manner, but we also see two interesting patterns in the Freedom forecast. When the ship goes around the Bahamas, the true track goes west of the islands, while the forecast goes mostly east of the islands. In this sense, we have a valid navigational pattern given by the forecast. A similar occurrence happens as the ship proceeds southeasterly, north of Puerto Rico. The forecast goes south in between Hispaniola and Puerto Rico avoiding both landmasses, before proceeding northwest at which point the error goes down to around 50NM. While the forecast was wrong, it gave valid outputs with all knowledge of land contained in the estimate $\hat{\mathscr{L}}$ and upsampling step, and no knowledge of land in computing $\hat{f}_{i}$ or $\mathbf{p}_{i}$. This indicates if the data was not sparse (e.g. streaming AIS positions), then the forecast would have avoided land without explicitly computing an approximation $\hat{\mathscr{L}}$.

Examining the plots in Fig. 1 more closely, we do not see a general trend in error with respect to time. For Dream the worst error occurs within the first two days of prediction and is almost entirely temporal. For Freedom the worst error is in the last two days and is almost entirely spatial. A traditional forecaster such as a Kalman or particle filter would be expected to have increasing error over time. Finally, we note that the error seems to be periodic. More specifically, it seems to cycle roughly in relation to days. It is possible that this is an artifact of our daily downsampling while preprocessing the historical data. This claim is supported by computing the auto-correlation function for the error time-series, which can be seen in Fig. 2. The plot shows statistically significant auto-correlation around the 48 hour period for both forecast, as well as for Dream at around the 80 hour and 125 hour marks.

To summarize: Algorithm 1 gives a very reasonable forecast, with only minimal information about the manifold of interest. It respects navigational constraints, and does not appear to loose accuracy over time in any general way.

V.2 Synthetic Data Forecasts

In order to have more control over both the training and testing data set (thereby guaranteeing the data meets our assumptions), we created three synthetic histories of trajectories. The first two were trajectories on $\mathbb{R}^{2}$, and included loiter points and followed the same underlying model, with a different geometry. The third was the Lorenz system in $\mathbb{R}^{3}$. This enabled us to compare our results with those in Berry et al. (2015).

V.2.1 Plane trajectories

The first two synthetic histories of trajectories on $\mathbb{R}^{2}$ consisted of 10,000 data points with Gaussian noise, shown in Fig. 3a and Fig. 3b. The generated tracks included six and five loiter points respectively, including a bifurcating trajectory where, after reaching the top leftmost loiter point the particle uniformly chose to go towards the center of the system or proceed due south. The purpose of the loiter points is to understand how the algorithm treats speed implicitly, while the purpose of the bifurcation is to see how the forecasting algorithm deals with such phenomena.

In both cases, we forecasted a length $\overline{q}=20$ trajectory $\{\mathbf{x}_{i}\}_{i=N+1}^{N+\overline{q}}$ for the dynamical system that generated the first $N=10,000$ data points, with $\mathbf{x}_{N+1}=(3.5,7.5)$. In reporting (see Table 2) we consider each index to be a half-unit of time, so that $\mathbf{x}_{N+1}$ corresponds to reporting time $t=0.5$ and $\mathbf{x}_{N+\overline{q}}$ corresponds to time $t=10.5$. We note that $\{\mathbf{x}_{i}\}_{i=N+1}^{N+\overline{q}}$ was not used to train the model (in the sense of contributing to the data used to build a KDE). These tracks are shown in blue in Fig. 4a and Fig. 4b, where the points represent the actual $\{\mathbf{x}_{i}\}$ and the points are connected linearly in both space and time to give a position for any $t\in[t_{N+1},t_{N+\overline{q}}]\subseteq\mathbb{R}^{+}$.

After generating training data and a ground truth trajectory for comparison to predictions, we compute a prediction $\{\mathbf{p}_{i}\}_{i=N+1}^{N+\overline{q}}$. We plot the predictions in red in Fig. 4a and Fig. 4b.

In Fig. 4a, we see that the overall shape of the prediction and ground truth are roughly similar, and they are extremely similar in Fig. 4b. In particular, the predicted path and ground truth are very close up through the predictions made after the ground truth leaves the loiter regions at $(1,5)$, and $(8,6)$ respectively.

For the first forecast, the average $\mathrm{APE}$ is $0.495$ with a standard deviation of $0.773$ units. It is clear that this relatively large deviation comes from only $t=9,10,10.5$. Looking only from $t=0.5$ to $t=8.5$, we see the $\mathrm{APE}$ is less than $0.4$ units, and is as low as $0.0112$ units at $t=5$. Moreover, every observation is within a 70% HDR.

For the second forecast, the average $\mathrm{APE}$ is $2.104$ with a standard deviation of $1.414$ units. Interestingly while the error is not monotonically increasing, the error is below average for every point in the first half of the trial, and above average for the second half of the trial. We also see drastically different performance with regards to the HDR. Only $8$ of the $20$ points lie within a 95% HDR. We chose to report the 95% HDR to test the importance of the choice of $\alpha$. The choice of 95% is also in keeping with the standard 95% confidence interval for a test-statistic.

The first notable failure in the proposed method can be observed in the first trajectory between times $t=8.5$ and $t=9$. The predictions linger at the loiter point for a slightly longer time than ground truth. After there is a prediction made outside of the loiter region ($t=9$), the distance to $\mathbf{p}_{N+17}$ from $\mathbf{p}_{N+16}$ is much larger than the distance between any two points on the test path. We see this phenomenon to a smaller degree between $\mathbf{p}_{N+13}$ and $\mathbf{p}_{N+14}$ in the second prediction.

To some extent this phenomenon is negative, because in some cases (as in the first trajectory) this may result in a prediction that is traveling too fast; i.e., the resulting predicted trajectory may not be physically realizable. On the other hand, the observed behavior is positive because the algorithm self-corrects, in the sense that after a loiter point, it catches up to where it should be, rather than just proceeding on from the loiter point at a constant velocity. If it did proceed from the loiter point at a constant velocity, it would then be temporally behind the correct position regardless of how accurate the prediction was spatially.

We do note that in the second trajectory, error is almost entirely temporal. For both cases the spatial nearness threshold set for the algorithm was $\epsilon=0.25$. In Table 3 we have computed the distance between every predicted point $\mathbf{p}_{i}$ and the closest point on the true trajectory. This distance ignores the time at which a forecast expects the particle to be in a certain location.

In Table 3, we see that mean distance between a forecast point and any true point is less than $\epsilon$, and 61.9% of the individual points were less than $\epsilon$ away from the nearest point on the forecast. We draw attention to this fact because by choosing $\epsilon$, as we noted, we are essentially constructing an equivalence relation on the data. Thus, the majority of the forecast points were in the “same place” as at least one of the true points.

The second notable failure occurs $t=9.5$ in the first trajectory and $t=10.5$ in the second trajectory, where the prediction back-tracks toward the loiter point, rather than remaining on its present trajectory. This is not consistent with any behavior exhibited by the ground-truth particle at any point on the trajectory. As with the first failure, there are both positive and negative aspects to this phenomenon. On the one hand, for the first trajectory the $\mathrm{APE}=0.095$ at $t=9.5$ is quite low. On the other hand, in both trajectories, the back-tracking movement is not physically representative of any pattern occurring in the data.

We note that the first predicted path correctly used the right fork of the bifurcating path. This is coincidental as the underlying model chose uniformly to move to the right or down. The prediction points happened to follow the correct path because there was a slight skew in the data toward the correct fork and then the test path also went right. We note that the underlying data used to train the KDE contained information on the bifurcation and the uncertainty included both possible paths. In the proposed algorithm, the uncertainty regions are far more important than specific pointwise predictors. This can be seen concretely in the multi-modal distribution of Fig. 5a for $t=10$.

The third notable failure also occurs at time $t=10$ for the first trajectory. We see the prediction is made at the loiter point $(1,1)$ corresponding to the left fork. This makes sense because the distance between the loiter point at $(1,5)$ and $(1,1)$ is less than the distance between $(1,5)$ and $(3,1)$. In other words, at the time when a particle following the left fork would be at the loiter point, a particle following the right fork would along the line between $(1,5)$ and $(3,1)$ (approximately). The concentration of mass of historical points (corresponding to the greatest peak in the KDE) at this time should then be expected to occur, as is the case, at $(1,1)$, because the historical points on the right fork are more spread out. We can see this shown in the contour plot of the KDE at time $10$ in Fig. 5a. This phenomenon also occurs (switching the role of the two forks) at time $t=10.5$, where we see the prediction is made near $(3,1)$, another loiter point. This phenomenon does not occur in the second trajectory because the region in which the trajectory was forecat does not include any bifurcating portion of the track.

We note, it is relatively straightforward to construct a multi-hypothesis pointwise forecast Griffin et al. (2016) that would provide branching paths following multiple local maxima of the KDE. This is considered in future work.

With the exception of the three errors noted, the predictions made were quantitatively and qualitatively good. The prediction was close in space (if not always in time) to the ground truth prior to reaching loiter points and along bifurcations in the path. Considering that we never explicitly compute speed, this is good performance. We also see that near the loiter points, the error becomes quite low. This makes sense, because we are most sure of the position of the particle when it is at a loiter point. Using the membership in the HDR as a measure of performance, we see that it is a good way to measure the temporal accuracy of a forecast, but not necessarily the spatial accuracy. This fact is exemplified by trajectory two. We also note that the HDR does not represent a “confidence region” in the frequentist sense of the term, but instead has a more Bayesian flavor.

V.2.2 Lorenz System

We use the classic Lorenz system defined by the system of differential equations in $\mathbb{R}^{3}$:

[TABLE]

The Lorenz system is the standard example of a chaotic dynamical system Hirsch et al. (2013) with a fractal strange attractor. Work on forecasting a particle whose underlying dynamics are governed by the Lorenz system has been considered in Berry et al. (2015). Our goal is this section is to create a forecast using Algorithm 1 that can be compared to the work in Berry et al. (2015). We used a Mathematica code snippet Wikipedia (2019) to generate 30000 points from a Lorenz system starting at the point $\langle 1,1,1\rangle\in\mathbb{R}^{3}$. The parameters were fixed to those in Lorenz’s original paper, with $\sigma=10,\,\rho=28,\,\beta=\frac{8}{3}$.

To the sample of 30000 data points, we added unbiased Gaussian noise. We set the spatial tolerance to $\epsilon=3$ and for the directional tolerance required the dot product to be positive. This is equivalent to partitioning the space with the plane whose normal vector is the velocity vector at $t_{0}$. Given these conditions we attempted to forecast the movement of the dynamics for the middle $\frac{1}{3}$ of the original data (prior to adding noise) given the noisy data. Since the Lorenz system is chaotic but deterministic, this is equivalent to to generating a new forecast. We used the naïve interpretation of manifold and simply used the embedding in $\mathbb{R}^{3}$. We note that there may be a better coordinate system to choose. However, finding an ideal manifold in which to embed the strange attractor may be difficult due to its fractal nature.

A visualization of the true path (blue) and our forecast (red) can be seen in Fig. 6. In Table 4 we report the error of every fifth time step for our forecast of the Lorenz system. In Fig. 7 we show a plot of the full results for all 101 points in time. The mean $\mathrm{APE}$ of the forecast was 9.01 units, with a standard deviation of 7.27. Moreover, 92 out of 101 data points were in the 70% HDR, even though Table 4 has no points outside of the 70% HDR.

It is difficult to see in the static image, but by animating the forecast reveals that this forecast has errors which are consistent with the planar forecast in Section V.2.1, namely that (i) error is often temporal in nature, (ii) the forecast regresses towards the system’s center of mass, and (iii) the forecast spends some amount of time on the wrong “lobe” of the Lorenz system. We also do not see any strong trend in the error as a function of time. At first glance, the error appears somewhat periodic. This is a reasonable suggestion, due to the pseudo-periodic nature of the movement of a single lobe of the orbit of the Lorenz system. However a discrete Fourier analysis does not support this observation; the error is mostly likely chaotic consistent with the motion of the underlying dynamical system. Moreover, $12.8\%$ of the forecast points have $\mathrm{APE}\leq 3$. This is notable, because we used a spatial tolerance of $3$. As we have previously observed, we cannot reasonably ask that the forecast has error less than that. We also note that the performance of the algorithm with respect to the highest density region computation is quite good.

Finally, we make one surprising observation. Consider the error curve shown in Fig. 7 as a function of time, which we will write $\mathrm{APE}(t)$. We integrate the curve with respect to time; i.e., we numerically compute:

[TABLE]

The plot of $g(t)$ is shown in Fig. 8, along with the line of best fit when the $y$-intercept is forced to zero. We see that a line $g_{\ell}(t)$ is a plausible model for $g$, and moreover that the approximate slope of $g_{\ell}$ is 8.28. Using all default settings from Mathematica to compute the fit, we have a 95% confidence interval on the slope of $g_{\ell}$ of $(8.131,8.425)$. This appears similar in value to the invariant measure computed for the Lorenz system in Berry et al. (2015). Since the invariant measure is the asymptotic error, were it true that a theoretical value of $\dot{g}_{\ell}$ were the same as the invariant measure $\mathrm{IM}$, then it would be that

[TABLE]

However, we only conjecture this relationship through experimental evidence rather than any theoretical result.

VI Conclusions and Future Directions

The previous literature on modeling and forecasting using Kernel Density Estimates has focused on well sampled data. In this paper, we develop Algorithm 1 to generate forecasts for the cases when we have a sparse, noisy data drawn from a recurrent trajectory. This approach offers an alternative to existing forecasting methods in cases when data are sparse and noisy, since it requires significantly less data. The energy minimization technique we use allows us to “connect the dots” between existing data points in an intelligent way, so that we can sample data as finely as we wish.

Algorithm 1 has several useful features. Traditional forecasters use some approximation or knowledge of speed, which is then projected forward in time. Our forecasting algorithm treats speed implicitly as we develop time-indexed probability distributions over a smooth Riemannian manifold. This simplifies the computation, and makes the algorithm path-independent. In practice, this allows us to pick specific times $t_{i}$ at which to provide a forecast position $\mathbf{p}_{i}$ and uncertainty function $f_{i}$, rather than having the requirement of computing intermediate steps to predict the whole path.

The algorithm also has nice theoretical properties. In particular, the optimal solution of a minimal energy trajectory is the asymptotic limit of the prediction as the sampling rate goes to infinity. Moreover, with a reasonable choice of parameters, it respects the existence of forbidden regions, albeit with a “fuzzy” boundary. With reasonable assumptions on (or knowledge of) the noise, we can easily quantify the fuzziness of the boundary.

When we assume but cannot guarantee that the governing energy function is the same, (as in the cruise ship forecasts) we still get reasonably good results. When we can guarantee that the underlying energy function is the same (as in our synthetic forecast) the algorithm performs quite well. It gives forecasts, with error that is mostly temporal, rather than spatial. When the error is spatial, it is as a result of a valid path which was not taken. In both the spatial and temporal error cases we still have an uncertainty region for each point, given by the highest density region of the KDE $\hat{f}$, and experimental evidence supports our claim that this makes sense as an uncertainty region.

While the algorithm performs well, there are certainly improvements which can be made. In particular, we note that often times the error is temporal rather than spatial. That is, the predicted point is not off because the predicted trajectory has high error in space, but rather because the time is off. A next step would be to consider how to improve this weakness in the algorithm. Additionally, for our ship forecast, we heuristically tuned the bandwidth to give a smooth path. It would be desirable to find a non-heuristic way to find the appropriate bandwidth, as the existing bandwidth selection rules (particularly Scott’s and Silverman’s Rules of Thumb) give a track with high error. Finally, improving our handling of bifurcating tracks when constructing $\mathbf{p}_{i}$ may also improve the resulting forecasts.

VII Acknowledgements

This work was supported in part by the Office of Naval Research through Naval Sea Systems Command DO 0451 (Task 23875). We would like to thank Brady Bickel, Eric Rothoff, and Douglas Mercer for helpful conversations during the development of this paper.