{'fraction_non_alphanumeric': 0.04117829457364341, 'fraction_numerical': 0.015007751937984495, 'mean_word_length': 4.537603021978022, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Differential Equation Units: Learning Functional Forms of Activation Functions from Data

Most deep neural networks use simple, fixed activation functions, such as sigmoids or rectified linear units, regardless of domain or network structure. We introduce differential equation units (DEUs), an improvement to modern neural networks, which enables each neuron to learn a particular nonlinear activation function from a family of solutions to an ordinary differential equation. Specifically, each neuron may change its functional form during training based on the behavior of the other parts of the network. We show that using neurons with DEU activation functions results in a more compact network capable of achieving comparable, if not superior, performance when is compared to much larger networks.


Introduction
Driven in large part by advancements in storage, processing, and parallel computing, deep neural networks (DNNs) have become capable of outperforming other methods across a wide range of highly complex tasks.
The advent of new activation functions such as rectified linear units (ReLU) Nair & Hinton (2010), exponential linear units (ELU) Clevert et al. (2015), and scaled exponential linear units (SELU) Klambauer et al. (2017) address a network's ability to effectively learn complicated functions, thereby allowing them to perform better on complicated tasks. The choice of an activation function is typically determined empirically by tuning, or due to necessity. For example, in modern deep networks, ReLU activation functions are often favored over sigmoid functions, which used to be a popular choice in the earlier days of neural networks. A reason for this preference is that the ReLU function is non-saturating and does not have the vanishing gradient problem when used in deep structures Hochreiter (1998).
In all of the aforementioned activation functions, the functional form of the activation function is fixed. However, depending on data, different forms of activation functions may be more suitable to describe data. In this paper, we introduce differential equation units (DEUs), where the activation function of each neuron is the nonlinear, possibly periodic solution of a parameterized second order, linear, ordinary differential equation. By learning the parameters of the differential equations using gradient descent on the loss function, the solutions to the differential equations change their functional form, thus adapting their forms to extract more complex features from data.
From an applicability perspective, our approach is similar to Maxout networks Goodfellow et al. (2013) and adaptive piece-wise linear units (PLUs) Agostinelli et al. (2014); Ramachandran et al. (2017). The problem of fitting a challenging function (y = (sin(t) − cos(2t) 2 )/2 + 4 * (1 + arccos(sin(t/2)))/3). A neural network with 25 DEUs in one hidden layer fits the function much more accurately than neural networks of 250 neurons with ReLU and Maxout activation function. Figure 1 shows nested noisy circles for which we trained neural networks with Maxout, ReLU, and our proposed DEU activation functions. For using ReLU, we needed at least two layers of four hidden units to separate the circles, while one layer of either two Maxout or two DEUs is enough to learn an appropriate decision boundary. While, the number of parameters learned by Maxout and PLU is proportional to the number of input weights to a neuron, and the number of linear units in that neuron, for each DEU, we learn only five additional parameters that give us highly flexible activation functions in return; altogether, resulting in more compact representations. Moreover, DEUs can represent a broad range of functions including harmonic functions, while Maxout is limited to piecewise approximation of convex functions (see Figure 2).
DEUs can also transform their forms during training in response to the behavior of other neurons in the network in order to describe data with interaction to the other part of the network, thus the neurons in a neural network may adopt different forms as their activation functions. The variety of activation function forms throughout a network enables it to encode more information, thus requiring less neurons for achieving the same performance comparing to the networks with fixed activation functions.
The main contribution of this paper is to explore learning the form of activation functions through introducing differential equation units (DEUs). We also propose a learning algorithm to train the parameters of DEUs, and we empirically show that neural networks with DEUs can achieve high performance with more compact representations and are effective for solving real-world problems.

Differential Equation Units
Inspired by functional analysis and calculus of variations, instead of using a fixed activation function for each layer in a network, we propose a novel solution for learning an activation function for each neuron. Experimentally, we show that by allowing each neuron to learn its own activation function, the network as a whole can perform on par with (or even outperform) much larger baseline networks. In the supplementary material, we have briefly outlined how in Calculus of Variations, the Euler-Lagrange formula can be used to derive a differential equation from an optimization problem with function variables Gelfand & Fomin (1963);Gelfand et al. (2000).
The main idea behind DEUs is to find the parameters of an ordinary differential equation (ODE) for each neuron in the network, whose solution would be used as the activation function of the corresponding neuron. As a result, each neuron learns a personalized activation function. We select (learn) the parameters of the differential equation from a low dimensional space (i.e. five dimensions). By minimizing the networks loss function, our learning algorithm smoothly updates the parameters of the ODE, which results in an uncountably 1 extensive range of possible activation functions.
We parameterize the activation function of each neuron using a linear, second order ordinary differential equation ay (t) + by (t) + cy(t) = g(t), parameterized by five coefficients (a, b, c, c 1 , c 2 ), where a, b, and c are scalars used to parameterize the ODE, c 1 and c 2 represent the initial conditions of the ODE's solution, and g(t) is a regulatory function that we call the core activation function. The coefficients are the only additional parameters that we learn for each neuron and are trained through the backpropagation algorithm. To simplify the math and because it is a standard practice in control theory, we have set g(t) = u(t), where u(t) is the Heaviside (unit) step function. Therefore, the ODE that we have chosen has the following form: This model is often used to denote the exchange of energy between mass and stiffness elements in a mechanical system or between capacitors and inductors in an electrical system Ogata & Yang (2002). Interestingly, by using the solutions of this formulation as activation functions, we can gain a few key properties: approximation or reduction to some of the standard activation functions such as sigmoid or ReLU; the ability to capture oscillatory forms; and, exponential decay or growth.

The learning algorithm
For fixed a, b and c, the solution of the differential equation is: for some functions f , f 1 , f 2 . Here, y(t) lies on an affine space parameterized by scalars c 1 and c 2 that represent the initial conditions of the solution (for desired values of y 0 and y 0 at some t = t 0 such that y(t 0 ) = y 0 and ∂y(t) ∂t | t=t0 = y 0 ).

Closed-form solutions
First, we solve the differential equations parametrically and take derivatives of the closed-form solutions: ∂y ∂t with respect to its input t, and ∂y ∂a , ∂y ∂b , and ∂y ∂c with respect to parameter a, b, and c. Moreover, the derivatives with respect to c 1 and c 2 are f 1 and f 2 , respectively. This is done once. We have solved the equations and taken their derivatives using the software package Maple (2018). Maple also generates optimized code for the solutions, by breaking down equations in order to reuse computations. 2 Although we used Maple here, this task could have been done by pen and paper.

Training the parameters
The parameters of the activation function (the ODE parameters and the appropriate initial conditions of its solution) are jointly trained with the neural networks' parameters using back-propagation. We adopt regular backpropagation to update the values of parameters a, b, c, c 1 and c 2 for each neuron, along with using ∂y ∂t for updating network parameters w (input weights to the neuron), and for backpropagating the error to the lower layers.
We initialize network parameters using current best practices with respect to the layer type (e.g. linear layer, convolutional layer, etc.). We initialize parameters a, b, c for all neurons with a random positive number less than one, and strictly greater than zero. We initialize c 1 = c 2 = 0.0. To learn the parameters θ = [a, b, c, c 1 , c 2 ] T along with the weights w on input values to each neuron, we deploy a gradient descent algorithm. Both the weights w, as well as θ are learned using the conventional backpropagation algorithm with Adam updates Kingma & Ba (2014). During training, we treat a, b, c, c 1 and c 2 like biases to the neuron (i.e., with input weight of 1.0) and update their values based on the direction of the corresponding gradients in each mini-batch.
Four forms based on the sign of ∆ = b 2 − 4ac, a and c

Singularity of solutions
If one or two of the coefficients a, b or c become zero, then the solution of the differential equation falls into a singularity subspace that is different than the affine function space of neighboring positive or negative values for those coefficients. For example, for b = 0 and a * c > 0 , the solution will be In this example, we observe that moving from c > 0 to c = 0 changes the resulting activation function from a pure oscillatory form to a parametric (leaky) rectified quadratic activation function. More formally, for a parametric variable quantity p ∈ {a, b, c, b 2 − 4ac}, if p = 0, then the solution of the differential equation may be different than the solution of the neighboring differential equations with p = 0 (see Table 1 for the complete set of singular subspaces). Moreover, with p → 0, then we may have y(t; p) → ∞ for certain values of t. In particular, when exponential forms are involved in the solution, this phenomenon can cause an extreme increase in the magnitude of the output value of the DEU. Therefore, we introduce a hyperparameter which is used as a threshold to clamp small values to zero. This simple trick avoids numerical computation errors as well as exponentiating large positive numbers when the denominator of the exponentiated fraction is very small.
Our learning algorithm allows the activation functions to jump over singular subspaces. However, if the absolute value of a,b, orc falls below , then the derivative with respect to that parameter becomes zero. The value of the parameter remains zero during the training if we use the regular derivative. In order to allow an activation function to "escape" a singular subspace, we introduce the concept of "outward gravitation" in the next subsection.
We do not allow a = b = c = 0, and in this rare case, we set b to . During the learning process at most two of a, b, and c can be zero. The sign of a, b, c, and b 2 − 4ac might also change the solution of the ODE, which create "subspaces" that are individually solved in closed-form.
When b 2 − 4ac is close to zero and ac > 0, the generic solution may become exponentially large. Therefore if − < b 2 − 4ac < and sign(a) == sign(c), we explicitly assume a = b 2 sign(a) c = b 2 sign(c) in our implementation to stabilize the solution and to avoid large function values (i.e., we force b 2 − 4ac = 0 in the solution of the ODE.).

Approximation of Dirac's delta function
The derivatives of the activation function with respect to t include Dirac's delta function δ(t), which is the derivative of the Heaviside function. In the parametric derivatives, we substitute the delta function with its approximation s e −s * t (1+e −s * t ) 2 , which is the derivative of σ(s * t) = 1 (1+e −s * t ) . This approximation is a commonly used in practice for the delta function Zahedi & Tornberg (2010).
The larger values of s result in more accurate approximation of the delta function. In all of our experiments, we set = .01, and s = 100 although further tuning might improve the results.

Outward gravitation
We need to let a DEU to jump out of a singular subspace when the network benefits from it. This way, no matter what the initialization has been and how the order of training samples has changed the geodesic path that an activation function follows to get the appropriate functional form, the activation function can recover from accidental falling into singularity points. To allow this, for a singular solution, we introduce a hypothetical non-singular differential equation that has the same initial conditions as the corresponding solution function at a reference point t * . 3 For a parametric variable p ∈ {a, b, c}, is the closest value outside of the singular subspace (p = p if the parameters are not in a singularity region). If the value of any of the parameters a, b, or c are in the singular subspace, we form a hypothetical differential equationãy +by +cy = u(t), whose solution has the same initial condition at the reference point t * as the current activation function (which is parameterized by a, b, c, c 1 and c 2 ). We use the derivative of the solution of this hypothetical differential equation for any of a, b, and c that are in the singular subspace. It is important to certify that the initial conditions of the hypothetical ODE (i.e., the value of the function and its derivative with respect to input t) are the same as the initial conditions of the actual activation function (otherwise, the derivative with respect to parameters a, b, and c would be from a function that behaves differently than the actual activation function). Let y(t; a, b, c, c 1 , c 2 ) be the solution of an activation function, then our goal is to choosẽ c 1 andc 2 such that forỹ(t;ã,b,c,c 1 ,c 2 ), we have: 3 The reference point can be t * = 0, or t * = t0 for some other choice of t0. In our implementation, we set ti that is the center of values in the training batch.
To findc 1 andc 2 , we assign the values of t * , a, b, c, c 1 , c 2 ,ã,b andc into equations (3) and (4), and sinceỹ is affine with respect toc 1 andc 2 , the results forms a two-variable, two-equation system of linear equations that is solved in constant time (see appendix for the details).

Implementation
Due to the individual activation functions associated parameters per neuron, each neuron can have a different activation. For example, one neuron could have a sinusoidal activation function while another has a quadratic form of ReLU. The direct way to implement our method would involve iterating over all neurons and computing the corresponding activation values individually, and would result in significantly higher latency than common activation functions such as ReLU. However, considering the ODE's closed form solutions, such computations are parallelizable on modern computing units like GPUs. We first compute a mask to assign each neuron to the subspace in which its current activation function resides. Next, we iterate over all possible functional subspaces, and compute in parallel the output of each neuron assuming that its activation function lies in the current subspace while ignoring the validity of the solutions or parameters. Finally, we use the subspace masks computed earlier to select the correct output of the neuron. A similar trick can be applied during the backward pass to compute the gradients in parallel using the masks computed during the forward pass. The pseudo-code is detailed in Algorithm 1. Return output 10: end procedure 2.2 Neural networks with DEUs are universal approximators Feedforward neural networks with monotonically-increasing activation functions are universal approximators Hornik et al. (1989); Barron (1993). Networks with radial basis activation functions that are bounded, increasing and then decreasing are also shown to be universal approximators Park & Sandberg (1991). A neural network with DEU is also a universal approximator (see Appendix for the proof).

A geometric interpretation
The solution set of a differential equation forms a functional manifold that is affine with respect to c 1 and c 2 , but is nonlinear in a, b, and c. Clearly, this manifold has a trivially low dimensional representation in R 5 (i.e., {a, b, c, c 1 , c 2 }). Gradient descent changes the functionals in this low dimensional space, and the corresponding functional on the solution manifold is used as the learned activation function. Figure 4 attempts to visually explain how, for example, a ReLU activation function transforms to a cosine activation function.
In appendix (toy datasets), we empirically show that only one DEU neuron that is initialized with ReLU (a = c = c 1 = c 2 = 0, b = 1) can transform itself to a sine function. In this scenario, the learned differential equation coefficients {a, b, c} and initial condition coefficients {c 1 , c2} perfectly represented a sine function after training. In contrast, the learned model by ordinary fixed activation FFNNs were much less accurate with significantly larger networks. Figure 7 (in the appendix) visually shows how changing a parameter changes the activation function's behavior. The red arrow shows a path that an initialized function takes to gradually transform itself to a different one.

Reduction to common activation functions
Two common activation functions are the sigmoid σ(t) = 1 1+e −t , and the rectified linear unit ReLU(t) = max(0, t). The sigmoid function is a smooth approximation of the Heavyside step function, and ReLU can be approximated by integrating sigmoid of s * t for a large enough s: 1+e −sz dz = log(1+e st )/s. Equivalently, y(t) = log(1+e t )+c 1 ≈ ReLU(t)+c 1 will be a solution of the following first order linear differential equation: y (t) = 1 1+e −st ≈ u(t) We can set the core activation function g(t) to σ(t) = 1 1+e −st , or to the step function u(t). For g(t) = σ(t) when a = 0, the solutions of the differential equation will involve the Gauss hypergeometric and Li 2 functions, which are expensive to evaluate. Fortunately, if we set the right hand side of the ODE to u(t), then the particular solutions will only involve common functions such as linear, exponential, quadratic, trigonometric, and hyperbolic functions.
In practice, if the learning algorithm decides that a and b should be zero, we use g(t) = 1 1+e −st (i.e. y(t) = 1 c * (1+e −st ) ). Otherwise, we use the step function to avoid complex-valued solutions that involve special mathematical functions. With these conditions in place, if a = 0, b = 0, and c = 1, we recover the sigmoid function; if a = 0, b = 1, and c = 0, we recover the ReLU function; if a = 1, b = 0, and c = 0, we obtain a parametric rectified quadratic form y = ReLU(t) 2 + c 1 t + c 2 (similar to parametric ReLU He et al. (2015); Xu et al. (2015)), which is the solution of y (t) = u(t). When b 2 − 4ac < 0, we observe an oscillatory behaviour. Depending on the sign of b, this can be decaying or exploding, but when b = 0, we observe a purely oscillating activation function.
The above-mentioned cases are only a few examples of solutions that could be chosen. The point to emphasize here is that an extensive range of functions can be generated by simply varying these few parameters (Figure 3 illustrates several examples).

Related Work
Recently, Chen et al. (2018) propose using ODE solvers in neural networks. However, the nature of our work is very different from Neural ODE. While that work maps iterative parts of the network to an ODE and uses external solver to find the solution and estimating the gradient, here we use the flexible functions that are the solutions of ODE as activation function. At the training time there is no ODE solver involved, but only functional forms of the solutions.
Oscillatory neurons are generally believed to be important for information processing in animal brains. Efforts have been made to explore the usefulness of periodic oscillations in neural networks since the 1990s, especially for medical applications Minami et al. (1999). However overall their applicability has not yet been appreciated Sopena et al. (1999). In recent times, however, researchers have begun re-exploring the potential of periodic functions as activations Parascandolo et al. (2016), and demonstrated their potentials for probabilistic time-series forecasting Hatalis & Kishore (2017) and lane departure prediction Tan et al. (2017). Furthermore recent work has demonstrated how pattern recognition can be achieved on physical oscillatory circuits Velichko et al. (2019).  For vision, speech and other applications applications on mobile, or other resource-constrained devices, research has been ongoing to make compact networks. ProjectionNet Ravi (2017) andMobileNet Howard et al. (2017) are both examples of methods that use compact DNN representations with the goal of on-device applications. In ProjectionNet, a compact projection network is trained in parallel to the primary network, and is used for the on-device network tasks. MobileNet, on the other hand, proposes a streamlined architecture in order to achieve network compactness. In these approaches, the network compactness is achieved at the expense of performance. We propose a different method for learning compact, powerful, stand-alone networks: we allow each neuron to learn its individual activation function enabling a compact neural network to achieve higher performance.

Results and Discussion
We have conducted several experiments to evaluate the performance and compactness of DEU networks.

Classification
We evaluate DEU on different models considering the classification performance and model size. We first use MNIST and Fashion-MNIST as our datasets to assess the behavior of DEUs with respect to the commonly used ReLU activation function, as well as Maxout and SELU. The neural network is a 2-layer MLP with 1024 and 512 dimensional hidden layers. While the CNN used is a 2-layer model made by stacking 32 and 16 dimensional convolutional filters atop one another followed by average pooling. DEUs are competitive or better than normal networks for these tasks while having substantially smaller number of parameters (see Table 2). For these experiments, we keep the network architecture fixed to ResNet-18 He et al. (2016b) and use the hyperparameter settings as in He et al. (2016b). We further show that this improvement persists across other model designs. First we use a preactivation ResNet He et al. (2016a), which is a ResNet-like architecture with a slightly smaller size. Second, to assess suitability for reducing the model size, we experiment with a stunted ResNet-18. The stunted model is created by removing all 256 and 512 dimensional filters from the standard ResNet-18 model. The result of this comparison (presented in Table 3) indicates that DEUs constantly work better than ReLU, PReLU, SELU, and Swish. Although Maxout is slightly better than DEU, it is using much more parameters (13272k vs 11174k for Resnet, 13267k vs 11174k for Preact-ResNet, and 801k vs 678k for ResNet-Stunted), which makes the comparison biased toward Maxout. Figure 5 shows the classification error on MNIST across different activations as the training progresses. It is clear that DEUs are better at almost all steps. Also as one might expect, they are better in the initial epochs due to the greater expressiveness of the activation function.

Computational cost
Prima-facie our method seems extremely compute-intensive. However as described earlier, with judicious design one can parallelize the overall model to be quite fast. In the worst case, our method will be |S| times slower than the same architecture with standard activation like ReLU, where |S| is the number of singular subspace solution of the ODE. Nevertheless, in practice all subspaces are unlikely to occur together. In most of our experiments on real data, we observed that three to five different solutions appear. Furthermore, we evaluate computation time of Resnet-18 models with DEUs and ReLUs on CIFAR-10 using a batch size of 128 on a Tesla K40c GPU. The model with DEUs takes 640ms total for combined forward and backward pass per batch, while the one with ReLUs requires 180ms per step. During prediction on test data, the time taken per batch is 111ms for the DEU-based model and 65ms for ReLU-based model.

Regression
We compare the performance of neural networks with one hidden-layer of one, two, four, eight, and 16 neurons with DEU, ReLU, LeakyReLU (LReLU), SELU, Swish and Maxout activation functions applied to a standard diabetes regression dataset. 4 We use 3-fold cross validation and report the average performance. Figure 6 shows that the neural networks with DEUs achieve specifically better performance with more compact networks. We see that other activation functions do not surpass the single DEU neuron performance until they are eight or more neurons.  (2015)(2016) hourly data from the California Department of Transportation (6 and 12 hours look ahead prediction); Solar-Energy. Solar power production records in the year of 2006 sampled every 10 minutes from 137 PV plants in Alabama State (6 and 12 hours look ahead prediction); Electricity. The electricity consumption in kWh recorded every 15 minutes from 2012 to 2014(12 hours look ahead prediction) Lai et al. (2018)  We test our activation function on multiple regression tasks involving timeseries forecasting. For these tasks, we use the LSTNet model Lai et al. (2018) and evaluate the functions based on the root relative squared error (RSE, lower is better) and correlation (CORR, higher is better) metrics. The results are presented in Table 4. DEU gives improvement in all cases and substantial improvements in some datasets.

Conclusion
In this paper we introduce differential equation units (DEUs), as novel activation functions based on the solutions of second-order ordinary differential equations (ODEs). DEUs can adapt their function form based on the features of data during training by learning the parameters of ODEs using gradient descent. We have showcased the ability of neural networks with DEUs to learn complicated concepts with a compact network representation. We have demonstrated DEUs' potential to outperform conventional activation function across a number of tasks, and with a reduced network size. Modern DNNs achieve performance gains in large by increasing the size of the network, which is not a sustainable trend. In response, we believe that this line of research can open future directions to explore more complex activation functions such as using the solutions of partial differential equations in order to compactly represent complex functions.
Proof. If b 2 − 4ac < 0, the roots of the characteristic equation of the ODE will be By substituting the solutions in Euler's formula the resulting sine and cosine functions will have frequency ω = Particularly, the solution of the ODE ay (t) + cy(t) = u(t) for ac > 0 is y(t) = sin( c a t)c 2 + cos( c a t)c 1 − u(t) c cos( c a t) − 1 . We will use this in the proof of the following theorem: Theorem 2. A neural network with DEU is a universal approximator.
Proof. By Lemma 1, for b = 0 and c = 1, the oscillation frequency of the activation function of a neuron will be ω = 1/ √ a. Therefore, by varying a ∈ R + , a neuron activation function can generate all oscillation frequencies. Now, consider a neural network with one hidden layer of differential equation neurons and one linear output layer. For a target function h(t), we take the Fourier transformation of h(t), and find the set of present frequencies Ω. For each present frequency ω i ∈ Ω, we add two differential equation neuron i 1 and i 2 to the hidden layer with c = 1, b = 0 and a = 1/ω 2 i . Let β i and γ i be the corresponding coefficient of the sine and cosine functions from the Fourier transformation. Then, we let c 1i1 = γ i , c 2i1 = β i and c 1i2 = c 2i2 = 0, and corresponding weights w i1 and w i2 from neurons i 1 and i 2 to the output linear layer to be equal to 1.0 and −1.0, respectively. This way, the two neurons will cancel the particular solution term of the function, and we'll have: y i1 (t) − y i2 (t) = β i sin(ω i t) + γ i cos(ω i t). By construction, after adding all frequencies from the Fourier transformation of h(t), the neural network will output There are multiple other ways to prove this theorem. For example, since the Sigmoid, ReLU, and ReLU-like functions are among the solutions of differential equations, we can directly apply results of the Stone-Weierstrass approximation theorem De Branges (1959).

B The Euler-Lagrange Equation in Calculus of Variations
As practiced in calculus of variations, a functional J is a function of functions, and it is usually expressed as a definite integral over a mapping L from the function of interest y and its derivatives at point t to a real number: J[y] = t∈T L(t, y, y , y , . . . , y (n) )dt. y(t) is the function of interest, and y , y , . . . , y (n) are its derivatives up to the nth order.
Intuitively, L is a combining function, which reflects the cost (or the reward) of selecting the function y at point t, and the integral sums up the overall costs for all t ∈ T . Therefore, the functional J : (R → R) → R can represent the cost of choosing a specific input function. The minimum of the functional is a function y * (t) = arg min y J(y) that incurs the least cost. Among other methods, y * (t) can be found by solving a corresponding differential equation obtained by the Euler-Lagrange equation Gelfand & Fomin (1963);Gelfand et al. (2000). In particular, the extrema to J[y] = t2 t1 L(t, y, y , y )dt are the same as the solutions of the following differential equation: C Outward Gravitation, Findingc 1 andc 2 Let y(t; a, b, c, c 1 , c 2 ) be the current activation function solution, andã,b,c be values outside the singularity subspace. We need to selectc 1 andc 2 such that: ỹ(t * ;ã,b,c,c 1 ,c 2 ) = y(t * ; a, b, c, c 1 , c 2 ) ∂ỹ(t;ã,b,c,c 1 ,c 2 ) ∂t t=t * = ∂y(t; a, b, c, c 1 , c 2 ) ∂t t=t * Since y = f (t; a, b, c) + c 1 f 1 (t; a, b, c) + c 2 f 2 (t; a, b, c) andỹ =f (t;ã,b,c) +c 1f1 (t;ã,b,c) + c 2f2 (t;ã,b,c), we will need:   Figure 7: The spectra of functions generated by varying one of a, b, c, and fixing the other two with c 1 = c 2 = 0. Figure 7 shows how changing a coefficient in the low dimensional differential equation space representation will affect the resulting functional on the manifold.

E Toy Datasets
We explore the behavior of a DEU neuron on sine function toy data by training it numerous times with different initialization, each time the differential transformed itself to have a sine-like solution. This supports that DEU activation functions transform its functional form during the training process. We also demonstrated the capability of DEU to learning complex concepts, and with a significantly reduced network size.
In these tests we compared fixed activation networks using ReLU, LeakyReLU and SELU activations to a significantly smaller neural network with DEUs.
We trained on two periods of a sine function and extrapolated half of one period of unseen values (expect for DEU that we test for one full period.) As seen in Figure 8, a single DEU (initialized with ReLU) can learn the sine function almost perfectly, while the fixed activation baseline networks fail to fit comparable to the training data, even with a 10x larger network.  Figure 8: Learning the sine function. The differential equation neuron learns the function significantly better than larger baselines. Also, the DEU was initialized with ReLU. It can easily transform itself to an oscillating form.