Nonlinear dynamics with neural networks

Preliminary Examination

Anil Radhakrishnan

Nonlinear Artificial Intelligence Laboratory, North Carolina State University

May 8, 2024

Table of Contents

Motivation

Neural Networks and optimization

▹ Background

▹ Metalearning Activation functions

Control and Chaos

▹ Background

▹ Neural Network control of chaos

Dynamics and symmetries

▹ Background

▹ Symmetries of controllers

Conclusion

References

Acknowledgements

Motivation

Motivation: Why Neural Networks?

• The preeminent computational tool of the contemporary world.
• Differentiable, optimizable, and scalable.
• Emergent uses in and out of Physics¹.

• Neural networks have established themselves as a great tool for a variety of tasks.
• They are not an answer to everything, but their scalability makes them a great tool in the age of big compute, big data, and big models.
• They have also found their place in fundamental science, both speeding up research and enabling new discoveries.

Motivation: Why Nonlinear Dynamics?

• Captures complexity.
• Formal, well-established framework.
• Emergent uses in and out of Physics.

Avalanche activity cascades in a sandpile automaton | Vortex street formed by flow past a cylinder | Turing patterns in a reaction-diffusion model

Animation by William Gilpin

Nonlinear dynamics is a powerful paradigm for understanding complex systems. And more and more systems are becoming more and more complex as we move forward.

Avalanche activity cascades in a sandpile automaton; a vortex street formed by flow past a cylinder; and Turing patterns in a reaction-diffusion model.

Motivation: Why Nonlinear Dynamics for Neural Networks?

• Optimization is inescapably nonlinear².
• Neural networks are inherently nonlinear.
• Scope of nonlinearities in neural networks is underexplored.

Trainability fractal by Jascha Sohl-Dickstein

Each pixel corresponds to training the same neural network from the same initialization on the same data — but with different hyperparameters. Blue-green colors mean that training converged for those hyperparameters, and the network successfully trained. Red-yellow colors mean that training diverged for those hyperparameters. The paler the color the faster the convergence or divergence

The neural network consists of an input layer, a tanh nonlinearity, and an output layer. In the image, the x-coordinate changes the learning rate for the input layer’s parameters, and the y-coordinate changes the learning rate for the output layer’s parameters.

Motivation: Why Neural Networks for Nonlinear Dynamics?

• Nonlinear dynamical systems are computationally expensive to solve.
• Paradigm of solution as an element of a distribution translates naturally to neural networks.
• Data-driven methods accommodate realistic complexity.

Background Neural Networks and optimization

Background: Differentiable Computing

• Paradigm where programs can be differentiated end-to-end automatically, enabling optimization of parameters of the program³.
• Techniques to differentiate through complex programs is more than just deep learning.
• Can be interpreted probabilistically or as a dynamical system.

• For the purposes of this study, it is better to think of these algorithms through the lens of differentiability.
• The large parameter space afforded by network structures are ripe for gradient based optimization.

Background: Neural Networks

\[ \newcommand{\defeq}{\mathrel{\mathop:}=} \newcommand{\R}{\mathbb{R}} \newcommand{\Param}{\theta} \newcommand{\Loss}{\mathcal{L}} \]

\[\varphi: \R^d \rightarrow \R^N_L \qquad T^l: \R^{N_{l-1}}\rightarrow \R^{N_{l}} \qquad \sigma: \R\rightarrow \R\] \[x \in \R^d \qquad W^l \in \R^{N_l\times N_{l-1}} \qquad b^l \in \R^{N_l}\]

The neural network is a composition of functions. We can swap out any of these components and and as long as its composable and differentiable, we can optimize it and call it a neural network.

Background: Backpropagation

Consider a differentiable loss function \(\Loss\),

\[ \begin{eqnarray*} \delta^L =& \nabla_{\varphi}\Loss\odot\sigma'(T^L)\\ \delta^l =& ((w^{l+1})^\intercal \delta^{l+1} ))\odot\sigma'(T^l)\\ \dfrac{\partial \Loss}{\partial b^l_j} =& \delta_j^l \\ \dfrac{\partial \Loss}{\partial w^l_{jk}} =& T^{l-1}_k\delta_j^l \\ \end{eqnarray*} \]

The tunable parameters can then be updated by: \(\Param\gets \Param -\eta \nabla_{\Param}\Loss\)

• Computational efficiency aside, backpropagation is essentially the chain rule.
• We take the derivative of the loss with respect to the output of the network and then backpropagate the gradient through the network.
• In the simplest case, we just subtract the gradient from the parameters to update them.

Background: Loss Functions

• Neural networks compute a probability distribution on the data space.
• Constructing a suitable differentiable loss function gives the path to optimization of a neural network. Common loss functions include:
  • Mean Squared Error \(\Loss(\theta) = \sum_{i=1}^N (y_i - f(x_i;\theta))^2\)
  • Cross Entropy \(\Loss(\theta) = -\sum_{i=1}^N y_i \log(f(x_i;\theta))\)
  • Kullback-Leibler Divergence \(\Loss(\theta) = \sum_{i=1}^N y_i \log\left(\dfrac{y_i}{f(x_i;\theta)}\right)\)

Loss functions are combined and regularized to balance the tradeoff between model complexity and data fit.

• The art of crafting the loss function can make or break a model.
• The loss needs to be computationally efficient, and smooth enough to allow for stable gradient based optimization.

Background: Optimization

• To find the best model parametrization, we minimize the loss function with respect to the model parameters.
• That is we compute, \(\mathcal{L^\star} \defeq \inf_{\theta \in \Theta} \Loss(\theta)\) assuming an infimum exists.
• To converge to a minima the optimizer needs an oracle \(\mathcal{O}\), i.e. evaluation of the Loss function, its gradients, or higher order derivatives. Then, for an algorithm \(\mathcal{A}\), \[ \theta ^ {t+1} \defeq \mathcal{A}(\theta ^ 0, \ldots, \theta ^ t, \mathcal{O}(\theta ^ 0), \ldots, \mathcal{O}(\theta ^ t), \lambda), \] where \(\lambda \in \Lambda\) is a hyperparameter.

Stochastic Gradient Descent, Adam, and RMSProp are common optimization algorithms for training neural networks.

• The optimizer is an important choice in the training process.
• Optimizer takes in loss related info and updates the parameters.

Adam stands for adaptive moment estimation and is a stochastic gradient descent optimization algorithm that uses an adaptive learning rate based on estimates of the first and second moments. It maintains exponential moving averages of the weights and gradients, which it uses to scale the learning rate. In other words, Adam uses estimates of the mean and variance of the gradients to adaptively scale the learning rate during training, which can improve the speed and stability of the optimization process.

When you want a fast and efficient optimization algorithm: Adam requires relatively little memory and computation, making it a fast and efficient choice for training deep learning models. When you have noisy or sparse gradients: Adam is well-suited for optimizing models with noisy or sparse gradients, as it can use the additional information provided by the gradients to adapt the learning rate on the fly.

The algorithm is an extension of the famous stochastic gradient descent (SGD) algorithm. The key idea behind RMSProp is to scale the gradient of each weight in the model by dividing it by the root mean square (RMS) of the gradients of that weight. This helps prevent weights with high gradients from learning too quickly while allowing weights with low gradients to continue learning faster. If you are training a model with many parameters and are experiencing issues with the model diverging or oscillating during training, RMSProp can help stabilize the training process by adjusting to the gradient.

Background: Meta-Learning

• Improve the learning algorithm itself given the experience of multiple learning episodes.
• Base learning: an inner learning algorithm solves a task defined by a dataset and objective.
• Meta-learning: an outer algorithm updates the inner learning algorithm.

Algorithms for meta-learning are still in a nascent stage with significant computational overhead.

Meta-learning allows the learning algorithm to optimize itself by adding a further layer of abstraction.

Figure by John Lindner

Background: Physics-Informed Neural Networks(PiNNs)

• Synthesizing data with differential equation constraints.
• Physics as a weak constraint in a composite loss function or a hard constraint with architectural choices.
• Symplectic constraints to the loss function gives Hamiltonian Neural Networks⁴.

One way to integrate physics into the neural network is to add the physics as a constraint to the loss function. This can be done by adding a term to the loss function that penalizes the network for not satisfying the physics. This is known as a weak constraint because the network is not required to satisfy the physics exactly, but rather to minimize the difference between the physics and the network’s predictions.

Background: Coordinates matter

Neural networks are coordinate dependent. The choice of coordinates can significantly affect the performance of the network.

Adapted from Heliocentrism and Geocentrism by Malin Christersson

As with any data driven modelling, the way you represent the data can have a significant impact on the performance of the model.

Metalearning Activation functions

(Published, US and International Patent Pending)

Foundation

• Most complex systems showcase diversity in populations.
• Artificial neural network architectures are traditionally layerwise-homogenous
• Will neurons diversify if given the opportunity?

Insight

• Activation functions can themselves be modeled as neural networks.
• Activation function subnetworks are optimizable via metalearning.
• Multiple subnetwork initializations allow the activations to diversify and evolve into different communities.

Figure from Publication

Methodology

• Developed a metalearning framework to optimize activation functions.
• Tested the algorithm on classification and regression tasks for conventional and physics-informed neural networks.
• Showed a regime where learned diverse activations are superior.
• Gave preliminary analysis to support diversity in activation functions improving performance.

Figure from Publication

Results: Scaling

Figure from Publication

Using 1DMNIST dataset, we showed that the learned diverse activations outperformed the traditional activations in the regime of low data few shot learning of around 5 epochs. These results are averaged over 100 initializations.

Results: Real World Data

Figure from Publication

We also showed that it can work for real world data, in this case a real pendulum. We used an HNN for this task to further illustrate the generalizability of the learned diverse activations.

Analysis: Participation Ratio

Estimate change in dimensionality of network activations

\[ Nr = \mathcal{R} = \frac{(\operatorname{tr}{\bf C})^2}{\operatorname{tr}{\bf C}^2} = \frac{\left(\sum_{n=1}^N\lambda_n \right)^2}{\sum_{n=1}^N \lambda_n^2} \]

where \(\lambda_n\) are the co-variance matrix \(\bf C\)’s eigenvalues for neuronal activity data matrix. The normalized participation ratio \(r = \mathcal{R} / N\)⁵.

Diverse activation functions use more of the network’s capacity.

Figure from Publication

If all the variance is in one dimension, say \(\lambda_n = \delta_{n1}\), then \(\mathcal{R} = 1\)

If the variance is evenly distributed across all dimensions, so \(\lambda_n = \lambda_1\), then \(\mathcal{R} = N\).

Typically, \(1 < \mathcal{R} < N\), and \(\mathcal{R}\) corresponds to the number of dimensions needed to explain most of the variance

Conclusions

• Learned Diversity Networks discover sets of activation functions that can outperform traditional activations.
• These networks excel in the regime of low data few shot learning.
• Due to the nature of metalearning, the memory overhead is significant concern for scaling.

Speculations: Achieving stable minima

 

Optimization of a neural network with shuffled data is a noisy descent.

This can be modeled with the Langevin equation:

\[ \DeclareMathOperator{\Tr}{Tr} d\theta_{t} = - \nabla \Loss(\theta_{t})\, dt + \sqrt{2\mathbf{D}} \cdot d\mathcal{W}_{t} \] with noise intensity \(\mathbf{D} = (\eta) \Loss(\theta) \mathbf{H}(\theta\star)\)⁶

Minima from diverse neurons is flatter than from homogeneous ones, as measured by both the trace \(\Tr\mathbf{H}\) of the Hessian and the fraction \(f\) of its eigenvalues near zero: \(\Tr \mathbf{H}_1 > \Tr \mathbf{H}_2 > \Tr \mathbf{H}_{12}\) and \(f_1 < f_2 < f_{12}\)

Since the noise term aligns with Hessian bear the minima, this suggests that the minima found by diverse neurons are flatter and more stable.

Speculations: Structure and universality to diverse activations

• Learned activation functions appear qualitatively independent of the base activation function.
• The odd and even nature of the learned functions suggest that the network is learning to span the space of possible activation functions.

The learned activations seem to form spanning sets for 1-d activation functions. There can be future work to explore the universality of this with multidimensional activations and large activation communities. But computational overhead is a significant concern.

Background: Control and Chaos

Background: Hamilton Jacobi Bellman(HJB) Equation

For a control system with state \(x(t)\) and control \(u(t)\), \[ \dfrac{dx}{dt} = f(x(t),u(t))dt \\ \]

\[ \mathcal{H}(x,u,t_0,t_f) = Q(x(t_f), t_f) + \int_{t_0}^{t_f} \Loss(x(\tau),u(\tau))d\tau\\ V(x(t_0),t_0,t_f) = \min_{u(t)} \mathcal{H}(x(t_0),u(t),t_0,t_f)\\ -\dfrac{\partial V}{\partial t} = \min_{u(t)} \left[ \Loss(x(t),u(t)) + \dfrac{\partial V}{\partial x}^T f(x(t),u(t)) \right] \]

V(x,t) is the value function, the minimum cost-to-go from state x at time t. If we can solve for V then we can find from it a control u that achieves the minimum cost.

Extension of Hamilton Jacobi equation/ Continuous time analog of Dynamic Programming.

Background: Model Predictive Control

Control scheme where a model is used for predicting the future behavior of the system over finite time window⁷.

Animations from do-mpc documentation

Based on these predictions and the current measured/estimated state of the system, the optimal control inputs with respect to a defined control objective and subject to system constraints is computed. After a certain time interval, the measurement, estimation and computation process is repeated with a shifted horizon

Proactive control action: The controller is anticipating future disturbances, set-points etc.

Non-linear control: MPC can explicitly consider non-linear systems without linearization

Arbitrary control objective: Traditional set-point tracking and regulation or economic MPC

constrained formulation: Explicitly consider physical, safety or operational system constraints The dotted line indicates the current prediction and the solid line represents the realized values.

This simulation uses moving horizon estimation with the problem discretized by orthogonal collocation.

Background: Chaos

Let \(X\) be a metric space. A continuous map \(f: X \rightarrow X\) is said to be chaotic on \(X\) if:

\(f\) has sensitive dependence on initial conditions is topologically transitive and has dense periodic orbits8

Devaney’s definition of chaos is a general useful definition for our purposes. Any two of the rules give the other in ℝ.

There is an equivalent definition where f is indecomposable and the periodic points are dense in X.

topologically transitive dynamical system has points which eventually move under iteration from one arbitrarily small open set to any other. Consequently, such a dynamical system cannot be decomposed into two disjoint sets with nonempty interiors which do not interact under the transformation.

Background: Traditional Chaos Control

Relies on the ergodicity of chaotic systems^9,10.

OGY relies on stabilizing unstable periodic orbits in the chaotic system. When a trajectory is near the UPO, the OGY method uses a perturbation to move the trajectory off the unstable manifold of the UPO and onto the stable manifold.

In the neighborhood of the UPO, the next point of the system’s Poincar´e section can be approximated using the Jacobian. And then the shift can be approximated by linear approximation

Shown attractor of the duffing oscillator with the UPOs marked in purple.

Background: Chaotic Pendulum Array¹¹

\[ l_n^2\ddot\theta_n =-\gamma\dot\theta_n\nonumber-l_n\sin{\theta_n}\nonumber +\tau_{0}+\tau_{1}\sin{\omega t}\nonumber +\kappa(\theta_{n-1}+\theta_{n+1}-2\theta_n) \]

This damped driven pendulum array is a classic example of a chaotic system that was deeply explored by Y. Braiman, John and Bill in the 90s exploring disorder taming chaos.

If we remove the coupling we get a simple driven damped pendulum. This system can itself be chaotic as well.

Background: Kuramoto Oscillator¹²

\[\dot{\theta}_i = \omega_i + \lambda\sum_{j=1}^N \sin(\theta_j - \theta_i)\]

Synchronizing fireflies video from Robin Meier

Kuramoto can be thought of as a first order mean field theory of the complex Ginzburg-Landau equation.

This model was initially proposed to capture synchronization in large weakly-coupled oscillator systems. It found a lot of applications like josephson junction arrays, power grids, and even modelling fireflies.

Neural Network control of chaos

Insight

• Optimal control of network dynamics involves minimizing a cost functional.
• Traditional control approaches like Pontryagin’s (maximum) principle or Hamilton-Jacobi-Bellman equations are analytically and computationally intractable for complex systems.
• Neural networks based on neural ODES can approximate the optimal control policy¹³.

Methodology

Instead of training a neural network with a predefined dataset, we train the network in situ with the network serving as the controller that is integrated in the dynamical system at each epoch.

The time horizon for the integration is changed from the first half and second half of training by 5x to allow extrapolation of the controller.

Results: Kuramoto Grid

The simplest test for this method was a grid of kuramoto oscillators. The network was trained to control the phase of the oscillators to a desired state. This works very well. Interestingly, this method also extends naturally to noisy kuramoto, which is not guaranteed by HJB.

Results: Pendulum Array

We also tested the method on the chaotic pendulum array into synchrony while maintaining the mean of the lengths. The network was able to control this system as well navigating out of the sliver of chaos.

Future Work: Exotic Dynamics¹⁴

\[\dot{\theta}_i = \omega_i + \lambda\sum_{j=1}^N A_{ij}\sin(\theta_j - \theta_i-\alpha)\]

This is an extension of kuramoto system where the coupling is nonlocal and the phase difference is shifted by a phase \(\alpha\).

It cannot be ascribed to a supercritical instability of the spatially uniform oscillation, because it occurs even if the uniform state is stable. Furthermore, it has nothing to do with the partially locked/partially incoherent states seen in populations of non-identical oscillators with distributed frequencies. There, the splitting of the population stems from the inhomogeneity of the oscillators themselves; the desynchronized oscillators are the intrinsically fastest or slowest ones. Here, all the oscillators are the same.

It would be a good test of the model to see if it can move a system into and out of chimeric states.

Future Work: Recurrence matrices in \(\mathcal{L}\)¹⁵

Recurrence matrices are are way to capture the dynamics of a system by looking at the recurrences of states. This is similar to the notion of delay coordinate embeddings used for approximating attractors. As you can see here the recurrence matrices are invariant to the delays.

This can be used in the loss function to prescribe the exact dynamics to optimize towards,

Computed by subtracting the transpose of the timeseries from itself and then thresholding the absolute value.

Background: Dynamics and symmetries

Background: Group theory for Dynamical Systems

Let \(\Gamma\) act on \(\mathbb{R}^n\) and \(f: \mathbb{R}^n \rightarrow \mathbb{R}^m\) . Then \(f\) is said to be \(\Gamma\)-equivariant if \(f(\gamma x) = \gamma f(x)\) for all \(x \in \mathbb{R}^n\) and \(\gamma \in \Gamma\)¹⁶.

For dynamical systems, for a fixed point \(f(x\star)=0\), \(\gamma x\star\) is also a fixed point.

Isotropy subgroup: Let \(v\in \mathbb{R}^n: \Sigma_v = \{\gamma \in \Gamma: \gamma v = v\}\)

Fixed pt subspace: Let \(\sigma \in \Sigma \subseteq \Gamma\). \(\text{Fix}(\Sigma)=\{v \in \mathbb{R}^n:\sigma v=v\}\)

Thm: Let \(f: \mathbb{R}^n \rightarrow \mathbb{R}^m\) be a \(\Gamma\)-equivariant map and \(\Sigma \subseteq \Gamma\). Then \(f(\text{Fix}(\Sigma))\subseteq\text{Fix}(\Sigma)\).

Cyclic group Dihedral group = Cyclic group + Flips Symmetric group = m! S^1 circle group O(n) Measure preserving transformations

For ODEs that means if \(f(x)\) is a solution, then so is \(f(\gamma x)\) for all \(gamma \in \Gamma\). So we can classify solutions upto their symmetry.

If a system starts in a symmetric state, it will remain symmetric.

Isotropy subgroup: The set of all group elements which leave the solution unchanged. i.e Fixed-point subspaces are invariant under the flow generated by f

Background: Group Equivariant Neural Networks

Group equivariant NN work through lifting convolutions to the group space, and then projecting back to the original space¹⁷.

Figure adapted from Bart Smets et al.

One was thats groups ar integrated into neural networks are as an extension of the traditional convolution neural network setup. Instead of just lifting to some real space which preserves translations as in traditional CNNs, we lift to a group space which preserves the group symmetries.

Symmetries of controllers

Insight

• Dynamical systems often obey symmetries.
• Controllers are essentially dynamical systems.
• Is there a computationally viable mapping between the symmetries of the system and the controller?

Methodology

• Analyze the controllers from the previous sections for symmetries using group equivariant autoencoders
• Construct controllers that respect the symmetries of the system
• Compare the performance of symmetric to conventional controllers.

Expectation

• Viability test for group equivariant neural networks in control systems.
• Mapping between the symmetries of the system and the controller.
• Performance analysis of symmetric controllers.

Conclusion

• Non-traditional treatments of neural networks let us better capture nonlinearity.
• Standard paradigms of Geometry, Statistics, and Algebra for understanding nonlinear systems are augmented by neural networks.
• The interplay of Physics, Mathematics, and Computer Science gives us the best shot at understanding complex systems

Deliverables

● Publication and Patent on Metalearning Activation functions.

◕ Publishable article on Neural Network Control of Chaos.

◔ Thorough study of the symmetries of controllers.

◑ Codebase for the above studies.

◔ Dissertation detailing the discoveries and pitfalls found throughout.

Estimated time of completion: August 2025

Acknowledgements

Background image from Jacqueline Doan

References

1.

Vicentini, F. et al. NetKet 3: Machine Learning Toolbox for Many-Body Quantum Systems. SciPost Physics Codebases 007 (2022) doi:10.21468/SciPostPhysCodeb.7.

2.

Sohl-Dickstein, J. The boundary of neural network trainability is fractal. (2024) doi:10.48550/arXiv.2402.06184.

3.

Blondel, M. & Roulet, V. The Elements of Differentiable Programming. (2024) doi:10.48550/arXiv.2403.14606.

4.

Greydanus, S., Dzamba, M. & Yosinski, J. Hamiltonian Neural Networks. (2019) doi:10.48550/arXiv.1906.01563.

5.

Gao, P. et al. A theory of multineuronal dimensionality, dynamics and measurement. (2017) doi:10.1101/214262.

6.

Mori, T., Ziyin, L., Liu, K. & Ueda, M. Power-law escape rate of SGD. (2022) doi:10.48550/arXiv.2105.09557.

7.

Fiedler, F. et al. Do-mpc: Towards FAIR nonlinear and robust model predictive control. Control Engineering Practice 140, 105676 (2023).

8.

Devaney, R. L. An Introduction to Chaotic Dynamical Systems. (Westview Press, Boulder, Colo, 2003).

9.

Ott, E., Grebogi, C. & Yorke, J. A. Controlling chaos. Physical Review Letters 64, 1196–1199 (1990).

10.

Ditto, W. L., Rauseo, S. N. & Spano, M. L. Experimental control of chaos. Physical Review Letters 65, 3211–3214 (1990).

11.

Braiman, Y., Lindner, J. F. & Ditto, W. L. Taming spatiotemporal chaos with disorder. Nature 378, 465–467 (1995).

12.

Kuramoto, Y. In International Symposium on Mathematical Problems in Theoretical Physics. in (1975).

13.

Böttcher, L., Antulov-Fantulin, N. & Asikis, T. AI Pontryagin or how artificial neural networks learn to control dynamical systems. Nature Communications 13, 333 (2022).

14.

Abrams, D. M. & Strogatz, S. H. Chimera States for Coupled Oscillators. Physical Review Letters 93, 174102 (2004).

15.

Kennel, M. B., Brown, R. & Abarbanel, H. D. I. Determining embedding dimension for phase-space reconstruction using a geometrical construction. Physical Review A 45, 3403–3411 (1992).

16.

Zee, A. Group Theory in a Nutshell for Physicists. (Princeton University Press, Princeton, 2016).

17.

Cohen, T. & Welling, M. Group Equivariant Convolutional Networks. in Proceedings of The 33rd International Conference on Machine Learning 2990–2999 (PMLR, 2016).

18.

Leshno, M., Lin, V. Ya., Pinkus, A. & Schocken, S. Multilayer feedforward networks with a nonpolynomial activation function can approximate any function. Neural Networks 6, 861–867 (1993).

19.

Cybenko, G. Approximation by superpositions of a sigmoidal function. Mathematics of Control, Signals and Systems 2, 303–314 (1989).

20.

Choudhary, A., Radhakrishnan, A., Lindner, J. F., Sinha, S. & Ditto, W. L. Neuronal diversity can improve machine learning for physics and beyond. Scientific Reports 13, 13962 (2023).

21.

Smets, B., Portegies, J., Bekkers, E. & Duits, R. PDE-based Group Equivariant Convolutional Neural Networks. Journal of Mathematical Imaging and Vision 65, 209–239 (2023).

Backup slides

Hyperparameters, integrators, and tolerances

Metalearning Activation functions:

MNIST1D: 1 hidden layer of 100 neurons, activation of 50 hidden neurons. 100 initializations averaged. RMSprop, Real Pendulum: same other than 50 initializations. Pytorch and Jax frameworks used.

Neural Network control of chaos:

Control and dynamics integration via diffrax Jax library for neural ODEs. TSit5 algorithm for ODE integration with PID controller with rtol 1e-7 and atol 1e-9. Stratanovich Milstein solver for SDE integration with PID controller with rtol 1e-7 and atol 1e-9. 1000 epochs Controller training with only 1/5 data for the first half. Implemented in Jax.

RMSProp addresses the issue of a global learning rate by maintaining a moving average of the squares of gradients for each weight and dividing the learning rate by this average. This ensures that the learning rate is adapted for each weight in the model, allowing for more nuanced updates. The general idea is to dampen the oscillations in directions with steep gradients while allowing for faster movement in flat regions of the loss landscape.

Tsit5 - Tsitouras 5/4 Runge-Kutta method. (free 4th order interpolant). Basically Dormand Prince with better coefficients.

Stratanovich Milstein - order 1 strong Taylor scheme

Universal approximation theorem

\[M(\sigma) = \text{span}\{\sigma(wx-\theta):w\in \mathbb{R}^n, \theta\in\mathbb{R}\}\]

Theorem 1: Let \(\sigma \in C(\mathbb{R})\). Then \(M(\sigma)\) is dense in \(C(\mathbb{R}^n)\), in the topology of uniform convergence on compacta, if and only if \(\sigma\) is not a polynomial¹⁸.

For Neural Networks,

Theorem 2: Let \(\sigma\) be any continuous sigmoidal function. Then finite sums of the form \(G(x) = \sum_{i=1}^N \alpha_i \sigma(w_i^T x + b_i)\) are dense in \(C(I_n)\). i.e. \(\forall f\in C(I_n)\) and \(\epsilon>0\), there is a sum \(G(x)\) such that \(\max_{x\in I_n}|f(x)-G(x)|<\epsilon\)¹⁹.

Informally, at least one neural network exists that can approximate any continuous function on \(I_n=[0,1]^n\) with arbitrary precision.

SGD via gradient flows

Gradient descent: \(x_{n+1} = x_n - \gamma\nabla f(x_n)\) where \(\gamma > 0\) is the step-size. Gradient flow is the limit where \(\gamma\rightarrow 0\).

There are two ways to treat SGD in this regime,

Consider fixed times \(t=n\gamma\) and \(s=m\gamma\).

• Convergence to gradient flow

Given the recursion \(x_{n+1} = x_n-\gamma\nabla f(x_n) - \gamma\epsilon_n\),

applying this \(m\) times, we get:

\[ \begin{align*} X(t+s) - X(t) &= x_{n+m}-x_n\\ &= -\gamma \sum_{k=0}^{m-1}\nabla f(X(t+\frac{sk}{m})) - \gamma\sum_{k=0}^{m-1}ε_{k+n} \end{align*} \] in the limit, \[ \begin{align*} X(t+s) - X(t) &= -\int_t^{t+s}\nabla f(X(u))du + 0 \end{align*} \] which is just the gradient flow equation.

• Convergence to Langevin diffusion

Given the recursion \(x_{n+1} = x_n-\gamma\nabla f(x_n) - \sqrt{2\gamma}\epsilon_n\),

applying this \(m\) times, we get:

\[ \begin{align*} X(t+s) - X(t) &= x_{n+m}-x_n\\ &= -\gamma \sum_{k=0}^{m-1}\nabla f(X(t+\frac{sk}{m}))- \sqrt{2\gamma}\sum_{k=0}^{m-1}ε_{k+n} \end{align*} \] The second term has finite variance \(\propto 2s\). When \(m\) tends to infinity,

\[ X(t+s) - X(t) =\\ -\int_t^{t+s} \nabla f(X(u))du + \sqrt{2}[B(t+s)-B(t)] \]

This limiting distribution is the Gibbs distribution with density \(\exp{(-f(x))}\)

γ∑m−1k=0εk+n has zero expectation and variance equal to γ2m=γs times the variance of each εk+n and thus tends to 0

Argument from Francis Bach

1DMNIST

Neural ODEs

• Forward propagation: \[ x(t_1) = \text{ODESolve}(f(x(t),t,\theta), x(t_0), t_0, t_1)\\ \text{Compute loss: } \Loss(x(t_1))\\ a(t_1) = \frac{\partial \Loss}{\partial x(t_1)}\\ \]

• Back propagation: \[ \begin{bmatrix} x(t_0) \\ \frac{\partial \Loss}{\partial x(t_0)}\\ \frac{\partial \Loss}{\partial \theta} \end{bmatrix} = \text{ODESolve}\left(\begin{bmatrix} f(x(t),t,\theta) \\ -a(t)^T \frac{\partial f(x(t),t,\theta)}{\partial x} \\ -a(t)^T \frac{\partial f(x(t),t,\theta)}{\partial \theta} \end{bmatrix} , \begin{bmatrix} x(t_1) \\ \frac{\partial \Loss}{\partial x(t_1)}\\ 0_{|\theta|} \end{bmatrix} , t_1, t_0 \right) \]

HJB Derivation Sketch

Bellman optimality equation: \[ V(x(t_0),t_0,t_f) = V(x(t_0),t_0,t) + V(x(t),t,t_f) \] \[ \begin{align*} \dfrac{dV(x(t),t,t_f)}{dt} &= \dfrac{\partial V}{\partial t} + \dfrac{\partial V}{\partial x}^T \frac{dx}{dt}\\ &= \min_{u(t)} \dfrac{d}{dt} \left[ \int_0^{t_f}\Loss(x(\tau),u(\tau))d\tau + Q(x(t_f), t_f) \right]\\ &= \min_{u(t)} \left[ \dfrac{d}{dt}\int_0^{t_f}\Loss(x(\tau),u(\tau))d\tau \right] \\ \end{align*} \]

\[ \implies \boxed{-\dfrac{\partial V}{\partial t} = \min_{u(t)} \left[ \Loss(x(t),u(t)) + \dfrac{\partial V}{\partial x}^T f(x(t),u(t)) \right]} \]

Pontryagin’s Principle

Let \((x\star,u\star)\) be an optimal trajectory-control pair. Then there exists a \(\lambda\star\) such that, \[ \lambda\star =-\dfrac{\partial \mathcal{H}}{\partial x} \] and \[ \mathcal{H}(x\star,u\star,\lambda\star, t) = \min_{u} \mathcal{H}(x\star,u,\lambda\star, t) \] and by definition, \[ x\star = \dfrac{\partial \mathcal{H}}{\partial \lambda} \]

Stochasticity and noise

We can add intrinsic noise to it by adding a noise term \(\xi\) with a strength \(\sigma\): \[ \frac{dx}{dt}=f(x,t)+ σ(x,t)ξ \]

Then the rectangular reimann construction is equivalent to: \(f(t)≈\sum_{i=1}^n f(\hat{x_i})χ_{Δx_i}\)

where \(f(\hat{x_i})\) are constants. This approximation is exact in the limit of \(n\to\infty\).

For the stochastic case, the rectangular construction is equivalent to: \(\phi(t)≈\sum_{i=1}^n \hat{e}_iχ_{Δx_i}\) Then the integral is: \[ \int_0^T \phi(t) dW_t=\sum_{i=1}^n \hat{e}_i ΔW_i \]

Here, the coefficients \(\hat{e}_i\) are not constants but random variables since we are allowed to integrate \(σ\) which depends on \(X_t\).

Recurrence plots caveat emptor

Roessler lyap

Roessler rec

Lyapunov exponent and Recurrence plot from Recurrence Plot Pitfalls

Neural network architectures

Autoencoder

GRU

CNN

All images from associated wikipedia pages (Autoencoder, GRU, CNN)