Machine Learning of Time Series Using Time-delay Embedding and Precision Annealing

TL;DR

This paper introduces a method combining time-delay embedding and precision annealing to improve machine learning predictions of time series, ensuring better generalization and identifying optimal training data requirements.

Contribution

It presents a novel approach that integrates time-delay embedding with a precision annealing training method for enhanced time series prediction.

Findings

Effective identification of training data needed for good predictions

Use of nonlinear time series analysis to create embedding spaces

Application of multi-layer perceptrons in high-dimensional embedding space

Abstract

Tasking machine learning to predict segments of a time series requires estimating the parameters of a ML model with input/output pairs from the time series. Using the equivalence between statistical data assimilation and supervised machine learning, we revisit this task. The training method for the machine utilizes a precision annealing approach to identifying the global minimum of the action (-log[P]). In this way we are able to identify the number of training pairs required to produce good generalizations (predictions) for the time series. We proceed from a scalar time series $s(t_n); t_n = t_0 + n \Delta t$ and using methods of nonlinear time series analysis show how to produce a $D_E > 1$ dimensional time delay embedding space in which the time series has no false neighbors as does the observed $s(t_n)$ time series. In that $D_E$-dimensional space we explore the use of feed forward multi-layer perceptrons as network models operating on $D_E$-dimensional input and producing $D_E$-dimensional outputs.

Equations

$$\mbox{$\bf{s}$}(n)=[s(n),s(n+\tau),s(n+2\tau),...,s(n+(D_{E}-1)\tau)],$$

$$S_q(n)=s(n+(q-1)\tau );q=1,2,...,D_E;$$

$$AMI(\tau)=\sum_{\{s(n),s(n+\tau)\}}P(s(n),s(n+\tau))\log\biggl{[}\frac{P(s(n),s(n+\tau))}{P(s(n))\,P(s(n+\tau))}\biggr{]}.$$

$$\displaystyle\mbox{$\bf{Y}$}^{(k)}(l_{0})=\{y(k),y(k+\tau),y(k+2\tau),y(k+3\tau),y(k+4\tau)\},$$

$$\displaystyle\mbox{$\bf{Y}$}^{(k)}(l_{F})=\{y(k+1),y(k+1+\tau),y(k+1+2\tau),y(k+1+3\tau),y(k+1+4\tau)\},$$

$$C(k)$$

$$x_{q}(l+1)=f_{q}\biggl{(}\sum_{v=1}^{D_{hl}}\mbox{\bf W}(l)_{qv}x_{v}(l)\biggr{)}\;\;q=1,...D_{h(l+1)},$$

$$A(x(l);k)=C(k)+\frac{R_{f}}{2}\frac{1}{\sum_{l=l_{1}}^{l_{F}}D_{hl}}\sum_{l=l_{0}}^{l=l_{F}}\sum_{q=1}^{D_{h(l+1)}}\biggl{(}x^{(k)}_{q}(l+1)-f_{q}(\sum_{v=1}^{D_{hl}}\mbox{\bf W}(l)_{qv}x^{(k)}_{v}(l))\biggr{)}^{2},$$

$$A(x_q^{(k)}(l),\text{mbox}\mathbf{W}(l))$$

$$E[G(\mbox{{\bf X}})|\mbox{$\bf{Y}$}]=\langle G(\mbox{{\bf X}})\rangle=\frac{\int d\mbox{{\bf X}}\,G(\mbox{{\bf X}})\exp[-A(\mbox{{\bf X}})]}{\int d\mbox{{\bf X}}\,\exp[-A(\mbox{{\bf X}})]}$$

$$y(n)\to \frac{2y(n)-(y_{max}+y_{min})}{y_{max}-y_{min}},$$

$$MS{E}_{\text{mbox}TrainingError}(M)=\frac{1}{M}\sum _{k=1}^M\frac{1}{D_E}\sum _{q=1}^{D_E}(x_q^{(k)}(l_F)-y(k+1+(q-1)\tau ){)}^2$$

$$MS{E}_{\text{mbox}ValidationError}(M)=\frac{1}{M_{\text{mbox}total}-M}\sum _{k=M}^{M_{\text{mbox}total}}\frac{1}{D_E}\sum _{q=1}^{D_E}(x_q^{(k)}(l_F)-y(k+1+(q-1)\tau ){)}^2.$$

Full text

**Machine Learning of Time Series

Using Time-delay Embedding and Precision Annealing**

Alexander J. A. Ty, Zheng Fang, Rivver A. Gonzalez,

Department of Physics

University of California, San Diego

9500 Gilman Drive

La Jolla, CA 92093-0357

Paul. J. Rozdeba,

Institut für Mathematik

Universität Potsdam,

Karl-Liebknecht-Str. 24-25, 14476 Potsdam, Germany

and

Henry D. I. Abarbanel

Department of Physics

and

Marine Physical Laboratory (Scripps Institution of Oceanography)

University of California, San Diego

9500 Gilman Drive

La Jolla, CA 92093-0357

Abstract

Tasking machine learning to predict segments of a time series requires estimating the parameters of a ML model with input/output pairs from the time series. We borrow two techniques used in statistical data assimilation in order to accomplish this task: (1) time-delay embedding to prepare our input data, and (2) precision annealing as a training method. The precision annealing approach identifies the global minimum of the action ($-\log[P]$). In this way we are able to identify the number of training pairs required to produce good generalizations (predictions) for the time series. We proceed from a scalar time series $s(t_{n});t_{n}=t_{0}+n\Delta t$ and using methods of nonlinear time series analysis show how to produce a $D_{E}>1$ dimensional time delay embedding space in which the time series has no false neighbors as does the observed $s(t_{n})$ time series. In that $D_{E}$-dimensional space we explore the use of feed forward multi-layer perceptrons as network models operating on $D_{E}$-dimensional input and producing $D_{E}$-dimensional outputs.

1 Background

Machine learning methods for capturing the structure of a time series with the goal of predicting future segments of that time series have been analyzed for many years  Frank \BOthers. (\APACyear2001); Kajitani \BOthers. (\APACyear2005); Goodfellow \BOthers. (\APACyear2016) . We revisit this problem using analysis tools allowing one to explore questions such as: if we are given a time series data set and a network architecture with which to predict a future segment of the time series, how many distinct samples of input/output pairs used in training the network are required to achieve very good prediction (generalization)? Ascertaining the number of training examples in order to attain a given performance metric, classification error for example, have been limited to the study of learning curves in the current ML literature.

In Abarbanel \BOthers. (\APACyear2018) two of the present authors recognized for the first time the equivalence between supervised machine learning (ML) and statistical data assimilation (SDA) as widely utilized in large Physics, Geophysical, and Biophysical modeling. This recognition opens up a variety of opportunities to use methods from SDA in tasks asked of ML with the possibility of both improving the performance of ML solutions as well as gaining insight as to how these solutions work. This paper builds on Abarbanel \BOthers. (\APACyear2018) using our knowledge of how variational principles in ML may be implemented using methods not often practiced in that literature. Another of the insights in Abarbanel \BOthers. (\APACyear2018) that we called ‘deepest learning’ when the number of layers in ML networks becomes continuous will be further visited in future publications Abarbanel \BBA Durstewitz (\APACyearFall, 2019).

We work within a setting where we are presented with scalar time series data $s(t_{n})=s(t_{0}+n\Delta t)=s(n)$, sampled every $\Delta t$. A sample of these data is shown in Fig. (1). We wish to present segments of these data to a multi-layer perceptron network and train the network to learn subsequent segments of the time series. The task asked of the network in this paper is to predict one step forward in time, namely $s(n+1)$ given $s(n)$. There is no barrier to training this class of network to answer other questions about the data series. To train the selected network using the given data, we use a precision annealing (PA) method Ye \BOthers. (\APACyear2014, \APACyear2015).

In this paper we explore the ability of a feed forward multi-layer perceptron (MLP) to accomplish learning this task. We show how PA allows us to answer questions about how many input/output pairs are required to achieve good generalization, namely, allowing the trained network to reliably predict from inputs not seen in the network training phase. In cases where there are practical limitations to the number of training data are available, (e.g. cost, ethical considerations, rarity, etc.) it is of interest to determine this. Our networks have only a few hidden layers, though there seems to be no barrier to making the network much deeper. The method we present can be used with other network architectures, for example, recurrent networks, with no fundamental change in approach. We address this configuration in the later parts of this paper.

2 Preparing the Data

We are presented with a time series, part of which is shown in Fig. (1). The data set is comprised of a large number of data values uniformly sampled in time at times $t_{n}=t_{0}+n\Delta t$. We do not know $\Delta t$. We are not given any further information about the sequence $\{s(n);\;n=0,1,...N\}$. We were given $\approx 10^{5}$ data points, and we discarded about 104 of them to eliminate potential ‘transients’. Only 2048 of the data points are shown in Fig. (1).

Our goal here is to train a feedforward MLP network architecture to give as output $s(n+1)$ when presented with input $s(n)$. We could have used the method described here to train the network to predict $s(n+K)$ for any integer $K\geq 1$; we restrict our discussion here to $K=1$. For larger values of $K$, extra caution would be required to assure that $K$ is not so large that the input and output are not correlated.

Without further knowledge of the signal $s(n)$, we assume that although it is a sequence of scalars, it might have come from projection onto the $s$-axis from the operation of a higher dimensional dynamical system. To examine this we seek a ‘proxy space’ which carries the essential properties of the original higher dimensional source of the observed signal $s(n)$. For this purpose we turn to techniques of nonlinear time series analysis Abarbanel (\APACyear1996); Kantz \BBA Schreiber (\APACyear2003).

3 Time-Delay Vectors

If the observed time series $s(n)$ comes from projecting onto the $s$-axis, then points which appear to be nearby in time $t_{n}=t_{0}+n\Delta t$ may be neighbors due to the projection rather than due to the dynamics that moves the actual system of interest forward in time in a higher dimensional space. Nonlinear time series methods for unfolding the scalar time series Aeyels (\APACyear1981\APACexlab\BCnt1, \APACyear1981\APACexlab\BCnt2); Takens (\APACyear1981) use the data $s(t_{n})=s(n)$ along with the time delays of the data at time points $t_{n}+q\tau\Delta t=t_{0}+(n+q\tau)\Delta t:\;s(n+q\tau)$. $\tau$ and $q$ are integers.

The idea here is that $s(n+\tau)$ contains information on how the dynamics of the source of the time series $s(n)\to s(n+\tau)$ moves the system of interest forward in time. This information is not available in $s(n)$ alone.

This leads us to form the $D_{E}$-dimensional vector extending from each time $t_{n}$:

[TABLE]

or, in components,

[TABLE]

$D_{E}$ is also an integer.

To use this idea in a practical sense, we must estimate the time delay $\tau$ and the dimension $D_{E}$ of the vectors \mbox{\bf{S}}(n) containing the properties of the original state space from which $s(n)$ is projected. The value of $\tau$ should not be too small or the system will not have revealed the new information coming from the operation of the underlying dynamics, and $\tau$ should also not be too large or noise and intrinsic instabilities of the (nonlinear) dynamics will erase the utility of information at time $n+(D_{E}-1)\tau$.

3.1 Selecting $\tau$

To estimate $\tau$ we use an information theoretic ‘correlation function’, the average mutual information (AMI) Fano (\APACyear1961); Fraser \BBA Swinney (\APACyear1986). This function is nonnegative Fano (\APACyear1961), and the Fraser Fraser \BBA Swinney (\APACyear1986) criterion is to select the first minimum of the AMI as a balance between $\tau$ being too large or too small. The minimum means that the coordinates $s(n)$ and $s(n+\tau)$ are correlated, but not so strongly correlated that no new information on the origin of the time series results from knowing both $s(n)$ and $s(n+\tau)$.

$AMI(\tau)$ requires the joint distribution of $\{s(n),s(n+\tau)\},\;P(s(n),s(n+\tau))$, as well as $P(s(n))$ and $P(s(n+\tau))$. The latter come from the marginal distributions of $P(s(n),s(n+\tau))$.

[TABLE]

This quantity, $AMI(\tau)$, answers the question: how much information (in bits if the logarithm is to base 2) do we learn from $s(n)$ about $s(n+\tau)$ on the average over all joint values of $\{s(n),s(n+\tau)\}$.

In the Left Panel of Fig. (2) we display $AMI(\tau)$ evaluated from 215 samples of the time series of $s(n)$. The first minimum of this is near $\tau\approx\,$ 7 or 8.

3.2 Selecting $D_{E}$

Once $\tau$ has been selected, the estimation of $D_{E}$ is made by systematically asking when neighbors in dimension $D$ for \mbox{\bf{S}}(n) remain neighbors when \mbox{\bf{S}}(n) is expressed in dimension $D+1$. This method of false nearest neighbors reveals a global property of the source of the time series; namely, the minimum dimension $D_{E}$ within which the vectors \mbox{\bf{S}}(n) can represent trajectories that contain no neighbors arriving through projection from higher dimensions. The analysis of these global false nearest neighbors $FNN(D)$ is shown in the Right Panel of Fig. (2) where $\tau=7$ is used. It is good to check that our results are robust against selecting $\tau=7$; using $\tau=6$ or $\tau=8$ each yield $D_{E}=5$. This is discussed in Abarbanel (\APACyear1996) in more detail.

Using $D_{E}=5$, we can evaluate the Lyapunov exponents of the dynamical system at the source of the data $s(n)$. The methods for doing this are described in Abarbanel (\APACyear1996); Kantz \BBA Schreiber (\APACyear2003). Briefly summarized: one uses the development of trajectories in the $D_{E}$-dimensional space, and following a few trajectories nearby each other in this space from one location in space to another construct a local map from one set of points to the location where they go in one step. This permits one to read off the local Jacobian matrix step by step through out the trajectories in $D_{E}$-dimensional space. According to the Oseledec theorem Abarbanel (\APACyear1996) the sequential products of these $D_{E}\times D_{E}$ matrices when diagonalized yields the Lyapunov exponents. To accurately estimate all $D_{E}$ Lyapunov exponents, one uses a recursive Q****R  decomposition.

Following this protocol, we find two positive exponents. There is one zero exponent telling us that the source of the data is some (unknown) differential equation for five state variables. There are two negative exponents. The sum of the Lyapunov exponents is negative, and the associated information dimension of the attractor is about 4.4 Abarbanel (\APACyear1996); Kantz \BBA Schreiber (\APACyear2003).

Until this point, we have withheld information about the method in which the data was generated for this experiment. We now reveal that it is a $D=5$ Lorenz 96 model with forcing term $F=8.15$. Information was kept from the reader in order to illustrate that in order to perform such an analysis, information about the original system is not required.

Now that we have established estimates for $\tau$ and $D_{E}$, we want to use the observed time evolution \mbox{\bf{S}}(n)\to\mbox{\bf{S}}(n+1) to train a network to implement this discrete time map in $D_{E}$-dimensions.

4 Training a Selected Neural Network

The idea is now to work in $D_{E}=5$ dimensional space on vectors \mbox{\bf{S}}(n)=[s(n),s(n+\tau),s(n+2\tau),...,s(n+(D_{E}-1)\tau)], and build a machine that learns the discrete time mapping \mbox{\bf{S}}(n)\to\mbox{\bf{S}}(n+1).

To the scalar data \mbox{\bf{s}}(n) we add noise of mean zero and rms error $\sigma$ to form noisy scalar data $y(n)=s(n)+\mbox{Noise}(0,\sigma)$. In our calculations we selected $\sigma$ to be 2% of the dynamical range of the observed data. With higher noise levels, the information content of each individual training pair is reduced, leading to a need of additional training pairs to achieve similar prediction capabilities.

Using these noisy scalar data we form a data library of many input/output $D_{E}$-dimensional vector pairs to be used at the input port at layer $l_{0}$ and the output port at layer $l_{F}$, to train a neural network. We use $k=1,2,...,M$ members of this library as our training set.

[TABLE]

$k=1,2,...,M$ here and $D_{E}=5$.

The network we choose is a Multi-layer Perceptron, and we wish to train it to take at the input vectors \mbox{\bf{Y}}^{(k)}(l_{0};n) and produce at the output vectors \mbox{\bf{Y}}^{(k)}(l_{F};n). At the input layer $l_{0}$ we have one input port with $D_{E}$ slots. At the output layer $l_{F}$ we have one output port with $D_{E}$ slots. The network has $l_{F}-2$ hidden layers $l=\{l_{1},l_{2},...,l_{F}-1\}.$ At the hidden layers we have $D_{hl}$ active units (’neurons’) at layer $l$.

As a function of the three quantities $\{l_{F},D_{hl},M\}$: $l_{F}$, the number of layers or the ‘depth’ of the network; $D_{hl}$, the number of active units in layer $l$ or the breadth of the network; and $M$, the number of distinct input/output pairs containing the information presented to the network for training, we wish to analyze, using statistical Physics methods, the quality of the training, the accuracy of the operation of the trained network on input/output pairs not used in training, and the ability of the trained network to represent the information in the M data pairs. In the networks we develop here, we take $D_{hl}$ to be independent of $l$.

4.1 The Action

In much of machine learning one seeks to minimize a cost function evaluated at the input and the output layers of a selected network. We call the activity variables (‘neurons’) at layer $l$ $x_{q}(l)$ for active unit $q=1,2,...,D_{hl}$ in layer $l$. The cost function for each input/output pair is at time k

[TABLE]

where the noise or errors in the input and output data have been taken to be Gaussian with zero mean and diagonal precision matrix $R_{m}$. $D_{h0}=D_{hF}=D_{E}$ for us.

This is to be minimized subject to a layer-to-layer connection rule

[TABLE]

with $\mbox{\bf W}(l)$ a matrix of weights to be determined in the minimization of $C(k)$.

If Gaussian errors with precision matrix $R_{f}$ are accepted in the layer-to-layer rule Eq. (4.3), then the full cost function is

[TABLE]

and we call this the ‘action’, after its usage in statistical Physics, for a single input/output data pair chosen at time $k$.

When we have many input/output pairs, we add a label to the active states in the network $x_{q}(l)\to x^{(k)}_{q}(l)$, and our goal is to minimize the action

[TABLE]

with respect to the connection weight matrices $\mbox{\bf W}(l)$ and the activities $x^{(k)}_{q}(l)$. Minimizing this action recognizes that for each input/output training pair, the activity of the network nodes may differ, but averaging over all M presentations of pairs from the library will train a possible generalizable network characterized by the $\mbox{\bf W}(l)$ and any other fixed parameters in the nonlinear functions $f_{q}$.

5 Use of the Action

The action $A(\mbox{{\bf X}})$, where X is the collection of all $x^{(k)}_{q}(l)$ in the network as well as the $\mbox{\bf W}(l)$ and other fixed parameters, is proportional to the negative of the logarithm of the conditional probability of the full state X conditioned on the M members of the input/output library, collected into a quantity $\bf{Y}$, used in the training set: P(\mbox{{\bf X}}|\mbox{\bf{Y}})\propto\exp[-A(\mbox{{\bf X}})]. An important use of this conditional probability density is the evaluation of expected values of functions $G(\mbox{{\bf X}})$ on the variables X, and this is evaluated by doing the integral

[TABLE]

It is here that the connection of machine learning with statistical Physics becomes apparent.

Estimating this integral can always be done with various Monte Carlo methods, and depending on the action $A(\mbox{{\bf X}})$ surfaces in X may be accomplished by finding the maxima of P(\mbox{{\bf X}}|\mbox{\bf{Y}}), or equivalently the minima of $A(\mbox{{\bf X}})$. The latter method Laplace (\APACyear1774, \APACyear1986) is why we are interested in the paths X which yield minima of $A(\mbox{{\bf X}})$.

6 Precision Annealing

We have developed a precision annealing (PA) approach Ye \BOthers. (\APACyear2014, \APACyear2015) for the minimization of the action Eq. (4.1) directed to finding the path with the smallest value of the action. The problem of finding the global minimum of the action, a nonlinear objective function of X, is NP-complete Murty \BBA Kabadi (\APACyear1987). PA is a continuation method  Allgower \BBA Georg (\APACyear1990) in $R_{f}$ that begins at very small $R_{f}$ where the global minimum is a solution to minimizing a quadratic form; this can be done in a straightforward manner, and moves adiabatically in $R_{f}$ to quite high values. Formally as $R_{f}\to\infty$, the layer-to-layer rule used in constructing the network becomes precise and deterministic.

While we have no mathematical proof that the global minimum is found, our numerical results indicate this may be the case. The PA method produces a set of minima of the action giving a numerical clue as to the roughness of the surface in path X space. It also finds low magnitude action minima with much higher rates of success than starting directly with large $R_{f}$.

The action surface $A(\mbox{{\bf X}})$ depends, among other items, on the number of measurement pairs $M$, on the hyper-parameter $R_{f}$, and on the number of model layers between $l_{0}$ and $l_{F}$. As the number of hidden layers increases, the model architecture deepens.

At the first step of PA we choose a solution to the optimization problem at $R_{f}=0$ and select the states at the hidden layers as drawn from a uniform distribution with ranges known from the dynamical range of the input/output state variables. One can learn that dynamical range well enough by solving the underlying model forward for various initial conditions. We make this draw $N_{I}$ times, and now have $N_{I}$ paths $\mbox{{\bf X}}^{0}$ as candidates for the PA procedure.

Now we select a small value for $R_{f}$, call it $R_{f0}$, and use the previous $N_{I}$ paths $\mbox{{\bf X}}^{0}$ as $N_{I}$ initial choices in our minimization algorithm. After using that minimization procedure we find $N_{I}$ new paths $\mbox{{\bf X}}^{1}$ for the minimization problem with $R_{f}=R_{f0}$. This gives us $N_{I}$ values of the action $A(\mbox{{\bf X}}^{1})$ associated with the new paths $\mbox{{\bf X}}^{1}$.

Next we increase the value of $R_{f}$ to $R_{f}=R_{f0}\alpha$ where $\alpha>1$. For this new value of $R_{f}$, we perform the minimization of the action starting with the $N_{I}$ initial paths $\mbox{{\bf X}}^{1}$ from the previous step to arrive at $N_{I}$ new paths $\mbox{{\bf X}}^{2}$. Evaluating the action on these paths $A(\mbox{{\bf X}}^{2})$ now gives us an ordered set of actions that are no longer as degenerate. Many of the paths $\mbox{{\bf X}}^{2}$ may give the same numerical value of the action. However, typically the ‘degeneracy’ lies within the noise level of the data $\approx(1/\sqrt{R_{m}})$.

This procedure is continued until $R_{f}$ is ‘large enough’ which is indicated by at least one of the action levels becoming substantially independent of $R_{f}$ and typically smaller than the others.

Effectively PA starts with a problem ($R_{f}=0$) where the global minimum is apparent and systematically tracks it and many other paths through increases in $R_{f}$. In doing the ‘tracking’ of the global minimum, one must check that the selected value of $\alpha$ is not too large lest one leave the global minimum and land in another minimum. Checking the result using smaller $\alpha$ is always worthwhile.

It is important to note that simply starting with a large value of $R_{f}$, $R_{f}\approx 1$ or larger, places one in the undesirable situation of the action $A(\mbox{{\bf X}})$ having multiple local minima into which any optimization procedure is quite likely to fall.

In the dynamical problems we have examined, one typically finds that as the number of measurement pairs $M$ is increased, more terms are added in the sum in Equation 4.1, thus raising the action levels of minima disproportionately until there is one dominant minimum. This we attribute to the additional information from the augmented set of measurement pairs.

6.1 Smallest Minimum; Not Necessarily a Convex Action

As our goal is to provide accurate estimations of the conditional expected value of functions $G(\mbox{{\bf X}})$ Eq. (5.1) where X, a path in model space, is distributed as $\exp[-A(\mbox{{\bf X}})]$, we actually do not require convexity of $A(\mbox{{\bf X}})$ as a function in path space. From the point of view of accurately estimating expected values, it is sufficient that the lowest action level be much smaller than the second lowest action level. If the action value at the lowest level $A(\mbox{{\bf X}}_{\mbox{lowest}})$ is much smaller than the action value at the next minimum $A(\mbox{{\bf X}}_{\mbox{second lowest}})$, then by a factor $\exp[-\{A(\mbox{{\bf X}}_{\mbox{lowest}})-A(\mbox{{\bf X}}_{\mbox{second lowest}})\}]$, the lowest path $\mbox{{\bf X}}_{\mbox{lowest}}$ dominates the integral to be done and provides a sensible choice for the path at which to evaluate the integral.

We will see in the examples below that when the PA procedure is used we may encounter situations where the action is apparently not convex. However, it may have a distinct smallest action level, much smaller in magnitude than the next lowest action level. That lowest level is expected to give a path which gives an accurate estimation to the expected value of functions $G(\mbox{{\bf X}})$ . This may occur in cases where sufficient information from the data has been transferred to the model, and this can indicate the size model adequate for the problem posed.

7 Action Levels for Our Time Series $\{s(n)\}$

We will now build and train a feedforward MLP Rozdeba (\APACyear2018) to learn the function $s(n)\to s(n+1)$ using $D_{E}$-dimensional data pairs from our library. We examined networks with $D_{E}$ = 5 dimensional input ($l_{0}$) and output ($l_{F}$) layers and 1-5 hidden layers each with the same number $D_{h}$ of active units (‘neurons’). The nonlinear function operating from layer-to-layer was chosen to be $\tanh(\bullet)$. We use the Python based program VarrAnneal Rozdeba (\APACyear2018) to perform the minimization of the action at each value of $R_{f}/R_{m}>0$.

To prepare our data for these network choices, we first scaled all of the noisy inputs $y(n)$ to lie within the range $[-1,1]$ via

[TABLE]

where $y_{max,min}$ are the maximum and minimum values taken by the noisy data. These scaled values were used to construct our data library of input/output pairs.

7.1 Two Hidden Layers; $D_{h}$ = 15; M = 50,…, 1200

Using PA and systematically moving $R_{f}/R_{m}$ from $R_{f0}/R_{m}\approx$ 10*-8* to $R_{f}/R_{m}\approx$ 1011 we evaluated $A(\mbox{{\bf X}})$ for $M=50,100,...,1200$ input/output pairs. $R_{f}$ was slowly increased using $\alpha=1.1$ At each value of $R_{f}/R_{m}$ we used $N_{I}=20$ initializations of the PA procedure.

We first examine the structure of the action levels as a function of $R_{f}/R_{m}$ for $M=50,300,900,1200$ I/O pairs. This is displayed in Fig. (3). Note that the action levels become nearly independent of $R_{f}/R_{m}$ for large values of this hyperparameter. Equally interestingly is the initial rise of $A(\mbox{{\bf X}})$ for large $R_{f}/R_{m}$ as $M$ increases. Then this saturates as the information in the time series $s(n)$ is represented fully in the network. See Fig. (4). Recall we call this the information content as, up to a constant $\langle A(\mbox{{\bf X}})\rangle=\langle-\log[P(\mbox{{\bf X}})]\rangle$.

8 Errors in Training and Validation

Once we have trained the proposed network, we can evaluate its quality when performing the task we have set it. In the example we have discussed here, that task in summarized as: when presented with a $D_{E}$-dimensional vector of inputs, created from time delays of a signal $s(n)$, accurately produce the next element of the time series $s(n+1)$. The quantities $s(n)$ and $s(n+1)$ are the first components of the data vectors.

We have tested (or validated) the operation of the network both on the data used to train the network and on data held aside in our library of I/O pairs. The latter is often called the “test” set or validation set or prediction set portion of the total data available to us Frank \BOthers. (\APACyear2001).

The error on the training set as a function of the number $M$ of I/O pairs used to train the network is given as

[TABLE]

This compares, in a least squares sense, the $D_{E}$-dimensional output $x^{(k)}_{q}(l_{F})$ from the trained network with the data from the training set, $y(k+1+(q-1)\tau)$, that are the output side of the input/output training pairs. The input to the trained network are the values \mbox{\bf{Y}}^{(k)}(l_{0})=\{y(k),y(k+\tau),y(k+2\tau),y(k+3\tau),y(k+4\tau)\}; the trained network operates on this $D_{E}$-dimensional vectors producing the output $x^{(k)}_{q}(l_{F});\;q=1,2,...,D_{E}$. We plot this as a function of $M$, the number of input/output pairs used in the training procedure, and in Fig. (5) we also examine the dependence on the number of active units (“neurons”) in each of the two hidden layers in the network.

We also can determine the accuracy of the trained network when acting on inputs selected from I/O pairs not used in the training of the network. This ‘validation’ error is evaluated as

[TABLE]

This compares the $D_{E}$-dimensional output $x^{(k)}_{q}(l_{F})$ from the trained network with the data from the set of input/output pairs that were not used during the training, $y(k+1+(q-1)\tau)$, that are the output side of the input/output pairs from the data library. We plot this as a function of $M$, the number of input/output pairs used in the training procedure. All of the I/O pairs from the data library not used in training were used in this validation error estimate. $M_{\mbox{Total}}=84971\gg M$.

In the input layer there are five ports into which a vector \mbox{\bf{Y}}(n)=[y(n),y(n+\tau),y(n+2\tau),...,y(n+4\tau)] is presented. In each of the hidden layers there are 15 or 25 or 35 active units (‘neurons’). At the output layer there are five ports within which a vector \mbox{\bf{Y}}(n+1)=[y(n+1),y(n+\tau+1),y(n+2\tau+1),...,y(n+4\tau+1)] is estimated. In Fig. (5) we show the MSEs in the estimation/training window and the prediction/generalization window as a function of the number of active units $D_{h}$ in the hidden layers and as a function of the number of input/output pairs $M$ used in training the network. The results show that for small $M$ the training and validation errors differ substantially, but as $M$ increases, enough information lies in the training set of $M$ training I/O pairs that the overall training error levels out when the network has completed its representation of the information in the data. Similarly, while the validation error is large for small $M$, as the network becomes ‘well trained’ (represents the information in the data series) the prediction MSE is essentially the same as the MSE in training. This result is consistent with the observation that the maximum value of the action levels for large $R_{f}/R_{m}$ becomes independent of $M$; see Fig.(4).

Fig. (6) examines the training and validation MSEs as a function of the number of layers in the network. The number of hidden layers is $l_{F}-2$, and we have evaluated this, using our achitecture, for $l_{F}=\,$ 4,5, and 6.

We display the dependence of the action on the number of input/output samples $M$ and $R_{f}/R_{m}$ relevant in the PA algorithm in Fig. (7) to further illustrate the outcome of our MLP network instantiation.

In Fig. (8) we display the predictions produced by our trained MLP network after the training using $M=400$ input/output pairs of segments of the noisy time series starting at $y(n-1)$ and predicting $y(n)$ in comparison with the known value of $y_{data}(n)$. The training is performed in $D_{E}$-dimensional space and the the output of the network is also in $D_{E}$ dimensional space. We display only the first component of the $D_{E}=5$ dimensional proxy state space vector as that is our (noisy) measured quantity. The predictions are what this network has been trained to do.

If we ask another question of the network: take the trained network as a dynamical system, namely, train the network using $M=400$ input output pairs, then use the trained network to predict $y(n)\to y(n+1)$ forward from the training window, we find the results in Fig. (9). On this task the trained network does not perform as well as on the task it was trained to do

9 Summary and Discussion

Using the interpretation of a familiar machine learning task: using the information flow in a scalar time series to train a rather standard multi-layer perceptron (MLP) to predict one step forward in the time series, as the equivalent of a statistical data assimilation (SDA) task Abarbanel \BOthers. (\APACyear2018), we have shown that using the precision annealing Ye \BOthers. (\APACyear2014, \APACyear2015) training methods of SDA, given the model architecture, leads to a network whose action ($\propto-\log[P(\mbox{{\bf X}})])$ rapidly becomes independent of the precision of the model as well as independent of the number $M$ of input/output model pairs and independent of the number of model layers for two or more hidden layers.

We attribute this independence to the class of models having captured the information content within the time series, and thus the method of training reveals how one may use precision annealing to estimate the number of input/output pairs required for excellent training and accurate prediction/generalization. Efforts to this end have been mainly curve-fitting learning curves, for example: Figueroa \BOthers. (\APACyear2012), Beleites \BOthers. (\APACyear2013). Furthermore, as knowledge of the conditional expected values of model state variables is what we wish to utilize approximations of the conditional probability distribution of model states P(\mbox{{\bf X}}|\mbox{\bf{Y}}) for, X is the collection of all model states at all layers, and $\bf{Y}$ is the collection of all input/output noisy data pairs, we can see how to properly limit the number of data pairs in a training set. This can be important in practical applications.

The data set used in these experiments was generated with the Lorenz96 model equations. Our analysis assumed no knowledge of this to illustrate that the decisions made for preparing the data can be made independent of its source. In curating the data we employed a technique from nonlinear time series analysis Abarbanel (\APACyear1996); Kantz \BBA Schreiber (\APACyear2003) that, while well known in the analysis of time series from nonlinear sources, has been used only once Frank \BOthers. (\APACyear2001), as far as we could tell, in a machine learning context over some decades. Considering its utility, we employed it here in a bit of detail as a friendly suggestion for future time series investigations.

In our earlier paper Abarbanel \BOthers. (\APACyear2018) introducing the analogy between machine learning and SDA, we noted the saturation of actions and prediction quality in a less structured example. We have shown it again here with an attribution to its information theoretic origin. As precision annealing within a Lagrangian training approach from classical methods of variational principles Gelfand \BBA Fomin (\APACyear1963); Marsden \BBA West (\APACyear2001); Kadakia \BOthers. (\APACyear2017) is utilized by us, the success may also be attributed to the capability of precision annealing to follow the global minimum of the action even though it is nonlinear in its variables X Murty \BBA Kabadi (\APACyear1987). This is the value of the action that maximizes the contribution of the conditional expected values of many quantities of interest.

The training method for the MLP network follows that for variational principles in data assimilation Evensen (\APACyear2009); Asch \BOthers. (\APACyear2017); Abarbanel (\APACyear2013); Marsden \BBA West (\APACyear2001) and control theory Kirk (\APACyear1970); Gelfand \BBA Fomin (\APACyear1963) in which ‘backpropagation’ procedures are absent and the methodology is well organized and principled. An additional value of the methods used here and in these references is that the symplectic structure of the variational principles is maintained Gelfand \BBA Fomin (\APACyear1963); Marsden \BBA West (\APACyear2001); Kadakia \BOthers. (\APACyear2017); Abarbanel \BOthers. (\APACyear2018).

Acknowledgment

Participation by PJR in this research has been partially funded by Deutsche Forschungsgemeinschaft (DFG) through grant CRC 1294 ”Data Assimilation” (project A06).