{'fraction_non_alphanumeric': 0.041943244086069406, 'fraction_numerical': 0.01318661736534949, 'mean_word_length': 4.441158647436675, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A non-autonomous equation discovery method for time signal classification

Certain neural network architectures, in the infinite-layer limit, lead to systems of nonlinear differential equations. Motivated by this idea, we develop a framework for analyzing time signals based on non-autonomous dynamical equations. We view the time signal as a forcing function for a dynamical system that governs a time-evolving hidden variable. As in equation discovery, the dynamical system is represented using a dictionary of functions and the coefficients are learned from data. This framework is applied to the time signal classification problem. We show how gradients can be efficiently computed using the adjoint method, and we apply methods from dynamical systems to establish stability of the classifier. Through a variety of experiments, on both synthetic and real datasets, we show that the proposed method uses orders of magnitude fewer parameters than competing methods, while achieving comparable accuracy. We created the synthetic datasets using dynamical systems of increasing complexity; though the ground truth vector fields are often polynomials, we find consistently that a Fourier dictionary yields the best results. We also demonstrate how the proposed method yields graphical interpretability in the form of phase portraits.


Introduction
Time signals, due to their temporal ordering, multiple scales, high dimension, and autocorrelation, require special tools for meaningful analysis. In this work, we propose a non-autonomous dynamical systems framework to address time signal classification, the problem of learning a mapping that assigns a distribution over labels y ∈ R |Y| to a vector-valued, continuous-time signal x : [0, T ] → R n . Here Y is a finite set of labels. We study the proposed method both theoretically and empirically. We find that, with its principled and parsimonious dictionary representation of the dynamical system's vector field, the proposed method approaches and/or exceeds the test set accuracy of competing methods.
We distinguish between time signals and time series; time signals are continuous in time, while time series are discrete in time. Frequently time series are obtained from sampling a time signal at discrete times. Applications of time signal classification include predicting the genre of music based on a sound recording [32], recognizing human activity using mobile sensors [33], diagnosing disease based on electrical biosignals (e.g., EEG, ECG, and EMG) [31,29,30], detecting natural phenomena such as earthquakes or volcanic eruptions using geophysical signals [24], and automatically distinguishing between mosquito species using wing-beat recordings [7].
One of the most promising approaches to time series classification involves deep recurrent neural networks (RNNs); popular methods include gated recurrent units (GRUs) and Long Short-Term Memory (LSTM) networks. The method developed in this paper can be viewed as an infinite-layer or continuum limit of a particular type of RNN; we describe this relationship next.

Deriving non-autonomous Equation Discovery from RNNs
Recurrent Neural Networks (RNN) were first introduced by [27] and have been used for processing sequential data. For input data, x t ∈ R n and a hidden state vector h t ∈ R m (typically initialized with h 0 = 0), a traditional sequence-to-label RNN is given by a discrete-time map together with an output layerŷ = σ(Ah T + b). Here φ is a user-specified activation function, and for classification problems, σ is typically the softmax function. To train an RNN, we learn parameters W , U and b. Typically RNNs are difficult to train due to a loss of long-term memory and suffer from computational issues in backpropagation through time, called exploding or vanishing gradients [23]. Gated RNNs, such as the Long Short Term Memory network (LSTM) were developed in [14] to overcome the challenge of long term dependencies. These more complex RNN architectures can be represented abstractly as a discrete-time map with parameter vector θ-see [11,Eq. (10.5)]: From (2), we derive a continuous-time model as follows. We first insert N − 1 hidden layers between h t−1 and h t and consider the discrete-time map When N = 1, we recover (2). For N ≫ 1, the model has a deep hidden-to-hidden transition [22]; N layers must be traversed to go from h t−1 to h t . Next, we take a near-identity or residual network form of the right-hand side function f . Essentially, f (h, x; θ) = h + δtΦ(h, x; θ) for a non-autonomous vector field Φ and δt = N −1 : Finally, we take the infinite-layer limit N → ∞ and obtain the central equation in the non-autonomous equation discovery (NAED) method : d dt We view the input signal x(t) as a forcing term in a non-autonomous dynamical system governing a hidden variable h : [0, T ] → R m . We introduce a function π : R m → R |Y| to assign a class label to the hidden variable evaluated at the final time,ŷ = π (h(T )). The objective is to learn the right-hand side Φ, parameterized by θ, and the function π, so that given a new time signal x(t), t ∈ [0, T ] we can estimate its class label, y. In the NAED method, we represent the right-hand side function Φ using a predetermined dictionary, a set of candidate functions, which is sufficiently large to capture a wide class of dynamics. There are a variety of choices for dictionaries; here we employ polynomial and Fourier basis functions.
In Section 2, we describe the NAED method in more detail, including efficient gradient computation using the adjoint method. In practice, we are given a time series, which we think of as a discretized time signal and we must also discretize the dynamical system to obtain a discrete-time approximation of the hidden variable. In this paper, we employ the optimize-then-discretize approach, where the gradient is computed analytically (see Theorem 2.2) using the continuous-time hidden variable and input time signal, and then evaluated using the time series and discretized hidden variable. This in contrast to a discretize-then-optimize approach that begins with discrete-time models such as (2) or (3), and then optimizes using gradients computed via backpropagation-through-time.

Motivation and related work
We motivate the NAED method in three ways: 1. As shown above, the NAED method can be viewed as an infinite-layer limit of a particular class of residual RNN architectures. In the context of deep feedforward networks, this residual structure has been shown to improve robustness to noise and generalizability [12]. When we train the NAED method, we avoid common issues with training deep RNNs, such as vanishing/exploding gradients. This allows NAED training to proceed without truncated backpropagation through time, gradient clipping, or other heuristics. 2. While deep RNNs have achieved high accuracy rates on certain time series classification tasks, they are often difficult to interpret. We seek to blend the high accuracy of deep RNN architectures with the interpretability of continuoustime dynamical system methods. In particular, we illustrate in Section 3 that trained NAED models can be interpreted graphically using phase portraits.
Additionally, the differential equation form of the model enables us to prove that the outputs of the classifier are stable with respect to both deterministic and random perturbations. 3. The dictionary representation of the vector field Φ is motivated by the literature on equation discovery. The problem formulation and goal in equation discovery differs from ours; there one assumes that the data consists of observations of the state vector h(t) of a continuous-time dynamical system. Using this data, the goal is to learn the vector field Φ. This is a nonparametric regression problem, equivalent to finding a system of ordinary differential equations that fit the observations h(t). The Sparse Identification of Nonlinear Dynamics (SINDy) method assumes that Φ can be represented as a sparse linear combination of elements from a dictionary Ξ [3]. Representing Φ with a dictionary requires fewer parameters than with a neural network. In SINDy, training proceeds via an iteratively thresholded least squares method whose convergence has been established [34]. In our Algorithm 1, we retain the iterative thresholding step from SINDy. However, we generalize SINDy in the following way: we do not assume access to h(t) at all, but rather the forcing function x(t). Learning Φ is a byproduct of our method, but the goal is to train a model whose predictionŝ y match the true labels y.
Continuous-time RNNs were proposed by Hopfield [15] and studied by many authors-see [9,1] and references therein. Early continuous-time RNNs were proposed as models of associative memory and hence are not directly comparable to the classifiers studied here. Still, early continuous-time RNNs share two features with NAED: the models are expressed as systems of nonlinear differential equations, and inputs are treated as non-autonomous forcing terms. Compared to NAED, early continuous-time RNNs have a rigid right-hand side structure that guarantees Lyapunov stability of the unforced system [15]. In contrast, NAED has a flexible right-hand side Φ that we can often represent as a sparse linear combination of dictionary functions.
More recently, there has been a growing literature that connects deep and recurrent neural networks with ordinary differential equations (ODEs). One branch of this literature seeks to apply ideas from dynamical systems theory to determine stable feedforward architectures [12], RNNs that do not exhibit chaotic dynamics [19], and RNNs that are constrained to be linearly stable [4]. The RNNs considered in these works [19,4] do not involve ODEs.
Another branch stems from Neural ODEs or ODE-Nets [5]. We view both NAED and Neural ODEs as infinite-layer limits of deep networks that are trained via the adjoint method rather than backpropagation. Note that in Neural ODEs, the vector field is typically modeled using a (static) feedforward neural network (rather than with a dictionary), and the input is used as an initial condition (rather than a forcing term) to the ODE system. Recent efforts have sought to make Neural ODE techniques more practical for large-scale problems [8,10,25] and also to better understand the learning of genuinely continuous-time dynamics [21]; we may be able to apply similar ideas to NAED in future work.
Recently, there has been some effort to generalize Neural ODE models to the RNN context. [26] combines ODE-Net and RNN layers, instead of purely relying on ODEs as in NAED. We also find continuous-time versions of GRU and LSTM models [2,16,13]. Compared with NAED, these architectures have more constraints on the right-hand side vector field Φ. Finally, the NAED dynamical system (5) can be viewed as a special case of the recently proposed neural controlled differential equation (NCDE) model [18]. Compared with NAED, the controlled differential equation allows for more general dependency of the hidden state h(t) on the input x(t). While NCDE uses a neural network model of the vector field, NAED uses a dictionary.
Let us briefly outline the present paper. In Section 2, we formally define the NAED method, establish existence/uniqueness of the method's solutions, compute gradients via the adjoint method, and also quantify the method's stability. In Section 3, we report the results of several computational experiments that demonstrate the competitive performance of NAED with respect to widely used and/or related methods for time series classification. We carry out these experiments both for synthetic data and for real data from the UCR Time Series Classification Archive [6]. In these experiments, NAED achieves similar or better accuracy than recurrent neural network methods (including LSTM and CFN architectures) and neural controlled differential equations (NCDE), with orders of magnitude fewer parameters. Additionally, we give examples of how NAED yields graphical interpretability in the form of phase portraits. We conclude in Section 4 with a discussion of the NAED method and ideas for future directions. Appendix A contains proofs of the Theorems given in this paper.
In this section, we describe our proposed non-autonomous equation discovery (NAED) method for time series classification, a gradient-based method for training it, our choice of dictionary in the NAED method, stability of the classifier, and a sparse version of the method.

NAED model for time signal classification
We assume that we are given data of the form {x i , T i , y i } i∈ [N ] , where x i : [0, T i ] → R n is a time signal and y i ∈ R |Y| is a probability mass function over the classes. In practice, y i will be a unit vector and arg max i y i will be the class or label. Note that we allow for the possibility that the time signals have different lengths. We consider the following non-autonomous dynamical system: For each time signal x i (t), we interpret the solution to (6), h i (t) ∈ R m ∀t ∈ [0, T i ] as a time-dependent hidden variable that is being forced by the function x i (t). The solution at time T i is used to make a class predictionŷ i viâ where A ∈ R |Y|×m , b ∈ R |Y| , and σ : We parameterize the vector field Φ : R m × R n → R m using a dictionary D = {ξ j } j∈ [d] , with ξ j : R m → R. We discuss specific choices for the dictionary, D, in Section 2.3, but we have in mind, e.g., multivariate polynomials. Let θ = (β, B). Composing the dictionary elements in a dictionary, where β ∈ R m×d and B ∈ R m×n are unknown coefficients.
To train the classifier, we must learn θ = (β, B), determining Φ via (7), together with the parameters A and b in (6c). We frame this learning problem as one of minimizing the following cross-entropy loss between labels y i and predictionsŷ i : Here, Θ = {β, B, A, b} represents all parameters to be learned. It is understood that h i satisfies (6) for the forcing An important consideration is whether there exists a solution of the dynamical system in (6) with right-hand side given by (7). The following theorem gives a sufficient condition for the existence and uniqueness of a solution.
Then there is an ε > 0 such that the initial value problem in (6) has a unique solution defined on the interval [−ε, ε].
A proof of Theorem 2.1 is provided in Appendix A.

Gradient computation and the adjoint method
subject to (6) where the objective function J(Θ) is defined in (8). To employ a gradient-based optimization method, we need to compute ∇ Θ J. However, directly computing the gradient of J with respect to Θ is complicated and computationally expensive because J involves h i (T ; Θ), the solution to (6) at time t = T . An alternative method to compute ∇ Θ J is to use the adjoint method, as we do in the following theorem.
Theorem 2.2 The gradients of the objective function in (8) with respect to the unknown parameters: β ∈ R m×d , B ∈ R m×n , A ∈ R |Y|×m , and b ∈ R |Y| are given by A proof of Theorem 2.2 is provided in Appendix A. The gradients from Theorem 2.2 are used with an optimization method to minimize the cross-entropy loss function (6c) and thereby train the model.

Dictionary choice
In the NAED model, the right-hand side of the dynamical system is given by Φ(h, x; θ) = βΞ(h)+Bx; see (7). The first term is a linear combination of dictionary functions, There is tremendous freedom in selecting the dictionary functions and this choice is paramount to the model. We tested NAED using two different dictionaries, a polynomial dictionary and a Fourier dictionary, described now in turn.
The polynomial dictionary consists of all possible polynomials of h ∈ R m up to k-th order.
Alternatively, we can consider a Fourier dictionary. Using separation of variables for the function Φ : . Each row of the vector βΞ(h) appearing in the RHS of the dynamical system (6) can be written as m i=1 f i (h i ), where the coefficients a i 0 , a i k , and b i k correspond to entries of β. In other words, our dictionary D consists of functions given by the outer product of harmonic functions, . For instance, for m = 2, and K = 1, the dictionary consists of the following 9 functions: Note that by the Stone-Weierstrass theorem, the polynomial dictionary and Fourier dictionary are dense in the space of continuous functions and L 2 , in the limiting case where k → ∞ and K → ∞, respectively. By choosing these parameters sufficiently large, all smooth dynamical systems can be represented as accurately as is needed.
We would like to apply Theorem 2.1 to guarantee the existence of a unique solution to (6). Assuming that the time signal x is continuous, it is enough to choose a dictionary that satisfies the Lipschitz continuity assumption. If we use the Fourier dictionary, then the Lipschitz constant is approximately L ≈ 2πK L . In this case, our model (6) has a unique solution until time T , provided we initialize β with sufficiently small values. On the other hand, if we use the Polynomial dictionary, there are two cases. If only linear terms are used in the dictionary, then the righthand side is Lipschitz continuous and a unique solution exists on the time interval [0, T ]. However, if we use higher-order polynomials in the dictionary, the righthand side is only locally Lipschitz and Theorem 2.1 can only guarantee a solution on a short time interval; the solution may blow up in finite time. In the numerical results in Section 3, we will observe that models with the Fourier dictionary are generally more accurate and less sensitive to initialization than models with a nonlinear Polynomial dictionary.

Stability of the NAED method
Dynamical systems theory can be used to prove that a given NAED classifier x C − → y is stable to noise; below we do this for both deterministic and stochastic perturbations.
The classifier C is Lipschitz continuous with constant L > 0 defined in the proof. That is, if we have a time signal x(t) and a noise corrupted version, Let W t denote the Wiener process in R d . Consider a time signal x(t) and a version corrupted by Gaussian white noise, Proofs of Theorem 2.3 and 2.4 are provided in Appendix A. Theorems 2.3 and 2.4 can be further interpreted in terms of classification stability as follows. Suppose that for a given time series, x, the NAED classifier gives the estimateŷ = C (x) (a probability vector). Further, suppose that max i C (x) i is uniquely attained so that the distance between C (x) and the decision boundary is positive. Then there exists a positive constant, ε > 0, such that for any corrupted time signalx(t) = x(t)+η(t) with η ≤ ε the two estimates C (x) and C (x) have the same maximum component, and so the assigned class does not change for the corrupted time signal. Here, the corruption can be either deterministic (Theorem 2.3) or stochastic (Theorem 2.4).

Algorithm 1 Sparse NAED method for time signal classification
Input: cut-off value, λ > 0 and initial parameters, Θ = {β, B, A, b}. while not converged: do (descent step) Take a step for the parameters Θ according to a chosen method to minimize J(Θ) in (8). E.g., for the gradient descent method with stepsize t > 0, update the parameters according to Θ ← Θ − t ∇ Θ J.
(threshold step) We threshold the β parameter values by setting

Sparse NAED method
The main task in our proposed learning method is to find the right-hand side (rhs) of the underlying non-autonomous dynamical system in (7), where the rhs is assumed to be a linear combination of dictionary terms. Here we explore the idea of imposing sparsity on the dictionary coefficients, with the goal of finding a simple representation of the underlying dynamics. As in equation discovery methods, we are motivated by the observation that most equations describing physical phenomena involve only a few relevant terms so that the rhs is sparse in the set of all possible functions. Imposing this assumption, we learn a model that balances accuracy and parsimony. Additionally, the sparsity assumption on the dictionary coefficients helps to prevent overfitting on the training dataset, leading to a method that is more robust to noise. Moreover, by assuming sparsity, we also obtain more interpretable dynamical models.
To develop a practical method to promote sparsity in the dictionary coefficients, we adopt the idea of iterative thresholding from [3,34]. The resulting algorithm is given in Algorithm 1. In Algorithm 1, entries of β with magnitude less than λ > 0 are thresholded to zero. This procedure is repeated until β has converged. In general, increasing λ trades accuracy for sparsity. The optimal value of λ will thus depend on the problem and data at hand. In practice, we use cross-validation to tune the value of λ. We note that the convergence of the algorithm depends on the value of λ.

Implementation details
We implemented the NAED method, described in Section 2, using TensorFlow. To solve the optimization problem, we used the ADAM optimizer with gradient computed as in Theorem 2.2. The gradient computation requires us to solve both the forward ODE (6) for h : [0, T ] → R m and the adjoint ODE (11) To approximate the solution of the forward and adjoint ODEs, we used the fourthorder Runge-Kutta (RK4) method, implemented via tfs.integrate.odeint_fixed in the tensorflow_scientific library. To approximate x(t) at times t not in the sampled time series data, we use linear interpolation In all numerical examples we fixed the initial condition for the hidden state in (6b) to be h 0 = 0. For the ADAM method, we used a learning rate ∈ {0.001, 0.005, 0.01, 0.05, 0.1, 0.5}.
For each dataset, we train the NAED model several times for different dimensions, m, of the hidden state, largest degree of polynomials k, or the maximum number of Fourier basis terms K. We report the results for several such models. Other methods. We implemented the RNN and LSTM methods in TensorFlow, using tf.keras.sequential with keras.layers.rnn and keras.layers.LSTM layers. Using the description in [19], we developed our own implementation of the CFN method in TensorFlow. We train the NCDE method using the published code [18]; we implemented this in PyTorch using the torchcde library. We trained each model using the cross-entropy loss function, the ADAM optimization method, and the default initialization. The models were trained until convergence of the loss function. For all considered datasets and all methods, we report the best result among our numerical experiments after varying hyper parameters such as network depth and width.

Forced harmonic oscillator
We consider a forced oscillator with position u(t) satisfying where γ is the damping coefficient, ω 2 is the undamped angular frequency, and x(t) is a specified forcing. To form the ground-truth labels, we record whether the position of u(T ) at the final time t = T is positive or negative, With the above framework, we generate a synthetic dataset as follows. Fix K ∈ N, γ > 0, ω > 0, and T > 0. For a forcing of the form where A k are randomly chosen amplitudes and α k are randomly chosen forcing frequencies, we numerically solve (12) for u(t), t ∈ (0, T ] and compute y via (13). We choose K = 2, γ = 0.2, ω = 1, T = 10, A k iid ∼ N (0, 1), and α k iid ∼ N (0, 1). In this paper, we use N (µ, σ 2 ) to denote the normal distribution with mean µ and variance σ 2 . The process is repeated N = 10000 times to create a dataset with 8000 training examples and 2000 test examples.
In Table 1, we tabulate the accuracy and number of trained parameters for various methods on this dataset. Note that the total number of parameters for the NAED method is given by In the first column of Table 1, additional information about each method is summarized. For the NAED Method with Polynomial dictionary, the parenthetical numbers are (# of units in hidden layer, maximum degree of polynomial in dictionary). The first row block of Table 1 is for the NAED method while varying either the dimension of the hidden units or the maximum degree of the polynomial entries. We observe that the Polynomial dictionary with a two-dimensional hidden state and polynomials up to degree one produces the best accuracy. This model also has the smallest number of parameters of all methods tested. This result might be expected as it agrees with the ground-truth model (harmonic oscillator). For the NAED method with the Fourier dictionary, the parenthetical numbers refer to (# of units in hidden layer, largest multiplier K) where the dictionary consists of units). For the NCDE method, the parenthetical numbers represent (# of units, width-depth of neural network for vector field). Note that repeated runs with different values of these hyperparameters were carried out, but we report only the hyperparameters for the models with the best test accuracy. We observe that all methods performed remarkably well for this simple dataset.
We can visualize our model using phase portraits; examples are given in Figure  1. Here, the black arrows represent the autonomous part of the learned vector field, h → βΞ(h). Also plotted in color are solution trajectories. Note that all samples are initially at the origin and the final positions at T = 10 are indicated by a square. The class associated with each sample is indicated in the legend. The classification decision is made using the final state of a trajectory via the probability vector, we see that the softmax function is being applied to the vector , where a i is the i-th row of A. We can visualize this decision in Figure 1 as follows. At the point h 0 , we draw the two rows of A as green vectors. These vectors partition R 2 into two regions; (each representing a class); we shade the region representing class 0 in red and class 1 in blue. In Figure 1, we observe that the the final states of the chosen trajectories belong to the correctly identified partition component.
We now remark on the identifiability of our model. Recall that a statistical model is said to be identifiable if the parameter values uniquely determine the  probability distribution of the data. For an identifiable model, it is in principle possible to learn the ground-truth parameters used to construct the data. Also recall that the goal of our algorithm is not to learn the mapping h → βΞ(h), but rather the mapping x → y. Since only the solution of the forward ODE at the final time is used to make this prediction, the learned vector field can differ from the ground truth vector field. If we consider Figure 1(b), the learned vector field closely agrees with the ground truth vector field (a), up to conjugation by an orthogonal matrix. The eigenvalues of βΞ(h) for in (a) are λ = −0.1 ± i0.995 which are close to the eigenvalues of βΞ(h) for (b), given by λ = −0.1015 ± i1.001. However, the vector fields in (c) and (d) are seen to differ from (a) considerably.

Forced Van der Pol oscillator
Consider the forced Van der Pol oscillator with position u(t) satisfying where µ = 0.3 controls the strength of nonlinear damping. We choose the forcing x(t) as in Section 3.2.1 and, at time T = 10, we define the label y as in (13). As shown in Table 3.2.2, the best accuracy for the forced Van der Pol dataset is obtained with the Fourier (2,2) dictionary. It is a remarkable result in that NAED only uses 58 parameters. On the other hand, the second best result trains roughly 20 times more numbers of parameters. Since the true system is nonlinear, it is not surprising to see strong performance from the Fourier dictionaries, which contain sums and products of trigonometric functions. Due to the presence of nonlinear polynomials in the Van der Pol system, we might expect that the best dictionary would be the Polynomial (2,3) dictionary. However, as discussed in Section 2.3, the nonlinear entries in the dictionary cause the right-hand side of (6) to be only locally Lipschitz continuous, so that Theorem 2.1 can only guarantee a solution on a short time interval. Since the class prediction is made using (6c), i.e., it depends on the hidden variable h(t) at time t = T , premature blowup of solutions spoils the learning process. The first block in Table 3

Forced Lorenz
Consider the forced nonlinear Lorenz system with positive parameters (σ, ρ, β): The first coordinate is forced by ∼ N (0, 1). Using the position of u 1 (t) at the final time T = 10, we define the label y as in (13). To generate the synthetic data, we choose parameters σ = 5, β = 1.3 and ρ = 10. Note that this dataset is not balanced in each class; the training data consist of 6348 and 1652 instances in classes 0 and 1, respectively, and the test data contains 1566 and 434 instances in classes 0 and 1, respectively.
As shown in Table 3, the highest accuracy for different choices of dictionaries and parameters in the NAED method is obtained by Fourier (3,1). This result demonstrates that complex dictionary entries are required to capture the nonlinearity in the underlying dynamics. It is remarkable that the NAED methods produced comparable results to the other methods using far fewer parameters, although it does not exceeded the classification accuracy of the LSTM method.

Forced Lotka-Volterra equations
Consider the forced Lotka-Volterra system, with initial condition (u 1 (0), u 2 (0)) = (5, 4), and parameters (α, β, δ, γ) = (0.8, 0.1, 0.01, 1.1). We sample A k , a k iid ∼ N (0, 0.5). After numerically solving up to time T = 10, we set the ground truth label y  via the indicator function for argmax(u 1 (T ), u 2 (T )). Note that the x(t) appears as a coefficient in the nonlinear terms. If we introduce the additional variable u 3 (t) = x(t), the forcing occurs linearly, This system suggests that we considerẋ(t) as the time series input data. Hence we train the model using, in turn, either x(t) orẋ(t) as the input. We generatė x(t) using the derivative of x(t) computed by hand. Note that the first block of Table 4 shows results with x(t) as input, while the second block shows results withẋ(t) as input. Comparing these two blocks in Table 4, we see that across all dictionaries and hyperparameters, the NAED method performs better withẋ(t) as input. Note that the NAED method with Fourier dictionary yields similar or better results than other methods regardless of whether x(t) orẋ(t) is used as input.

Stochastic gated partial diffusion equation
Consider the one-dimensional stochastic gated partial diffusion equation [20],  At z = 0, we impose the reflecting (Neumann) boundary condition. At z = 1, we impose the switching (time-dependent Robin) boundary condition, where for all t ∈ [0, 1], we have x(t) ∈ {0, 1}, a switching function. For the initial condition, we use an approximation to the Dirac delta δ(z − 0.5), given by u 0 (z) = , where σ = 0.1. The solution has an interpretation in terms of a particle experiencing Brownian motion on the interval. The probability of finding the particle at time t and position z is given by u(z, t). The initial condition is interpreted as the particles all starting near z = 0.5. The boundary condition at z = 1 has the interpretation that when x(t) = 1 a particle leaves the interval when it reaches the boundary and when x(t) = 0 the particles are reflected. The proportion of particles remaining in the interval at time t, referred to as the survival probability is given by S(t) = 1 0 u(z, t) dz. For a given switching function x(t), we assign a label y based on the survival probability at time t = 1; The classification problem seeks the mapping from the switching function x(t) to the binary class y. We generate a synthetic dataset for this problem with 8000 training examples and 2000 testing examples as follows. To generate each switching function x(t) we choose an integer, q, between zero and ten uniformly. We then randomly select q times in the interval [0, 1] and starting with x(t) = 0, we set x(t) to alternate between 0 and 1 at these times. For each switching function, x(t), we approximately solve the heat equation for u(z, t) as follows. We apply a forward difference in time and a second-order central difference scheme for the space derivative. We use a spatial discretization size of dx = 0.05 and temporal step size of dt = 0.01. To obtain roughly balanced class sizes, we choose the diffusion coefficient to be κ = 0.165. For this choice of parameters, the CFL condition κ (dt) (dx) 2 ≈ 0.66 < 1 is satisfied, so the numerical method is stable. The solution at time t = 1 is used to define the label y as in (17).
A comparison of the accuracy of various methods is given in Table 3.2.5. As shown in the first block of Table 3.2.5, the proposed NAED method works well on this dataset generated using a partial diffusion equation. Among the several choices of entries for the dictionary, we achieve the best accuracy with the Fourier (2-1) dictionary. The NAED method provides comparable accuracy to other methods with substantially fewer parameters.

Synthetic dataset with noise
In this section, we train the sparse NAED method (see Section 2.5) and show the robustness of this method on a noisy dataset. To generate the noisy data, we contaminate the forced harmonic oscillator input/forcing x(t) from Section 3.2.1 with noise: where η(t) is a Gaussian process, mutually independent for different t. Note that noise is added on the original data x(t) ∈ [−5.8, 5.6] for t ∈ [0, T ].
For the noisy data, we apply the sparse NAED method within a cross-validation loop to select λ. For each value of λ ∈ {0.01, 0.03, 0.05, 0.1, 0.5, 1}, and within each fold of 5-fold cross-validation, we train with Algorithm 1 until convergence. We then choose λ to minimize the cross-validation test error. The last column of the Table 6 records the number of non-zero entries in the trained β. As shown in each block of Table 6, the performance of the sparse NAED method tends to be slightly better than competing methods. In particular, the sparse Fourier (2-2) method achieves the best test error with substantially fewer parameters than competing RNN methods. Note that in this synthetic example, the underlying dynamical system does possess a sparse representation.

UCR archive datasets
In this section, we compare the proposed NAED method with other algorithms on a few univariate time series datasets from the UCR archive [6]. For the NAED method, we present the most accurate model by varying the candidate functions of dictionary and cut-off values for sparsity. For the Neural CDE method, we used either 32 or 64 hidden channels; in this method, the vector field is represented using a feedforward neural network with one hidden layer with either 64 or 128 units. The total number of parameters is reported in 7. The results are summarized in Table 7. In Figure 2, we show an example trajectory for each class and use colored partitions to denote the classification regions and decision boundaries.
The Two Pattern dataset is synthetically generated and has 1000 training and 4000 test samples. There are four balanced classes and the sequence length for all samples is 128. As shown in Table 7, the best accuracy is obtained with a Fourier   dictionary. Compared with other methods, the NAED method provides slightly lower accuracy but is still close to 100% on both train and test data. The Plane dataset contains outlines of airplanes measured by a sensor. The classification problem is to distinguish the type of airplane where there are seven airplane shape classes: Mirage, Eurofighter, F-14 wings closed, F-14 wings opened, Harrier, F-22 and F-15. There are 105 instances in both the training and test sets, each having length 144. As presented in Table 7, the NAED method (with Fourier 1 dictionary) surpasses the test accuracy of the RNN, LSTM, and CFN. It does this even with 100-200 times fewer parameters than these competing methods. The NCDE method is the best on this dataset; it has over 200 times more parameters than NAED. This dataset shows that NAED works well on a multiclass classification problem.
The Kitchen Appliance dataset is behavioral data recorded from 251 households and measured by a device in two-minute intervals over a month. Note that each series has length 720. This problems classifies how consumers use electricity within their home, so there are three classes: Kettle, Microwave and Toaster. This data contains 375 instances in the training and test sets. In Table 7, NAED with a sparse-Fourier1 dictionary returns the best accuracy on this dataset with only 29 parameters. Here, the cutoff value is set to λ = 0.03 and two entries of β are dropped to zero.
The Computer dataset consists of 250 train and test instances for a consumer's electricity usage behavior in a home. Each sample consists of recordings made every two minutes over a month so that total length is 720. There are two classes: Desktop and Laptop. According to Table 7, the best accuracy is obtained by the LSTM method. NAED with sparse-Fourier2 dictionary and cutoff value λ = 0.05 nearly matches the LSTM's accuracy. The imposed sparsity condition replaces 16 entries in β with zero; consequently, the trained vector field is relatively simple and interpretable.
The FordB dataset contains 3636 training and 810 test instances. Each instance consists of 500 measurements of engine noise together with a label. The classification problem is to diagnose the existence of certain symptoms in the automotive subsystem, so there are two classes. Note that the training data were collected in typical conditions while test data were collected under noisy conditions. Hence the FordB dataset forces the classifier to generalize from clean to noisy data. Here, the NAED method is learned using the sparse-Fourier1 dictionary and thresholded by λ = 0.05.
For the Kitchen Appliance dataset, NAED achieves the best test set results; for the remaining four datasets, NAED's parameter count is on average > 200 times less than that of the method with the best test set performance. For the first three datasets considered in Table 7, NAED is the only method that achieves competitive test set results with a small number of parameters. For Computer, the parameter counts for NAED and LSTM are similar. For FordB, RNN performs surprisingly well with a low parameter count. Based on the RNN results here, we conjecture that NAED underfits this dataset; a more scalable implementation of the NAED method would enable us to explore larger values of the dimension of h and the largest Fourier multiplier K.
As we described in Section 2, the NAED method learns a representation of the underlying vector field based on a prespecified dictionary. With polynomial or harmonic basis functions, these vector fields can be approximated using only a few terms. By promoting sparsity, Algorithm 1 can further enhance parsimony. As shown in experiments, competing methods require at least 2 times and up to 500 times the number of parameters required by NAED.

Discussion
In this paper, we developed a framework for analyzing time signals based on nonautonomous dynamical systems. A time signal, x(t), is interpreted as a forcing function for a dynamical system (6) that governs a time-evolving hidden variable, h(t). As in equation discovery, the dynamical system is represented using a dictionary of prespecified candidate functions and the coefficients are learned from data. We refer to the resulting model as non-autonomous equation discovery (NAED). This framework is applied to the time signal classification problem, where the hidden variable, at a final time, h(t = T ), is used to make a prediction via the composition of the softmax function and an affine function. Using a cross-entropy loss function, we train the NAED model using a gradient based optimization method, where the gradients are efficiently computed using the adjoint method; see Theorem 2.2.
Through a variety of experiments-on both synthetic and real datasets-we demonstrated that the NAED method achieves accuracy that is comparable to RNN, LSTM, CFN and NCDE methods on binary and multi-class classification problems; see Section 3. Note that [18] shows that NCDE itself outperforms other RNN architectures, including continuous-time/ODE-like GRU models [2,16] and a method that merges an RNN with a neural ODE [26]. The NAED method generally requires far fewer parameters than neural network-based methods and the number of parameters can further be reduced by using a sparse version of the algorithm; see Algorithm 1. We also show in Section 3.4 that sparsity improves the trainability of the method and its robustness to noise in the data. Finally, by construction, our method is interpretable using the theory of dynamical systems. For example, using phase plots, we can visualize the trajectories of the underlying dynamical system and how they navigate the decision boundaries between classes.
Since our model is built on dynamical systems, we can generate synthetic data from a dynamical system and then pose the inverse problem of trying to recover the ground-truth system from the data. For a synthetic dataset based on the forced harmonic oscillator (Section 3.2.1), we showed that the NAED method for classification is not generally identifiable, i.e., the method does not always recover the ground-truth parameters. However, in the case of a linear dictionary, we recover the ground-truth parameters up to conjugation by an orthogonal matrix.
There are a variety of natural future directions for this work. Since the NAED method is built on dynamical systems, we could use dynamical systems theory to further analyze a particular trained NAED model. For example, one could use stability theory to further sharpen and generalize the misclassification estimates in Theorem 2.3 and 2.4. To enhance the method's ability to deal with noisy time signals, one could combine the NAED method with filtering methods (e.g., the Kalman filter). Since we interpret time signals as continuous objects and discretize within the method (the optimize-then-discretize approach), multi-scale methods could be used in training. A slight generalization of the model would be to let B in (7) be a parameterized operator, B = K k=0 B k ∂ k t , where B k ∈ R m×n are unknown coefficients. In the forced Lotka-Volterra equations (Section 3.2.4), we considered using as forcing either x orẋ and this generalization would avoid this. Finally, the NAED framework developed here could be applied to other time signal analysis tasks, such as prediction and forecasting, classification, segmentation, and denoising.

A Proofs
Proof (Theorem 2.1.) Let Φ(h, x(t)) be rewritten as Γ (h, t) = Φ(h, x(t)). For some r > 0 and a > 0, define Br = { h − h 0 ≤ r} ⊂ K, Ia = {|t| ≤ a}. Since x is continuous in time, there exists a constant M > 0 such that Also for every t ∈ Ia, h → Γ (h, t) satisfies the local Lipschitz condition on K: for every h 1 , h 2 ∈ Br, From the existence/uniqueness theorem in ordinary differential equations (see, e.g., Here, we have used integration by parts to rewrite the Lagrangian. Taking the variation of the Lagrangian with respect to h k (t) gives where D h Φ = βD h Ξ(h) is the Jacobian of Φ with respect to the h. Setting the variation to zero, we find that λ k (t) satisfies the adjoint equation given in (11). The gradients of the objective in (10) are then obtained by taking the partial derivatives of the Lagrangian with respect to the unknown parameters, Θ = {β, B, A, b}. The gradient with respect to β and B are given by  For the cross-entropy loss function in (8), a short computation shows that (y i − σ(Ah i (T i ) + b)).
Combining these results concludes the proof. where C = B e LT β . The softmax prediction function in (6c) is Lipschitz continuous with constant that we denote by Lσ. We have where L = LσC, as desired. Let L denote the Lipschitz constant for the dictionary Ξ. We estimate Continuity of W t implies continuity of sup 0≤s≤t |Ws|. Note that sup 0≤s≤t |Ws| is non-decreasing. Hence Gronwall's inequality yields where C = B e LT β . We combine this with the Lipschitz bound on softmax: where L = LσC as before. The remaining estimates can be derived from the density computed in [17, §2.8A]; for clarity, we provide a self-contained treatment. Let B t denote the Wiener process in R, and let τz = min{t : B t = z}, a first passage time. Then note that By symmetry of B t , term I is 1/2; by continuity of B t , term II is 0. Hence P (τz ≤ t) = 2P (B t ≥ z) = erfc(z(2t) −1/2 ) where erfc is the complementary error function. Now using the reflection principle, we have Let W t,j denote the j-th coordinate of W t ; each W t,j is an independent onedimensional Wiener process. With |w|p denoting the p-norm of the vector w ∈ R d , we have |w| = |w| 2 ≤ d 1/2 |w|∞. Putting these facts together, we estimate Combining this with (18) yields the conclusion of the theorem.