Christian Kuehn

John Guckenheimer, Christian Kuehn

Slow manifolds are important geometric structures in the state spaces of dynamical systems with multiple time scales. This paper introduces an algorithm for computing trajectories on slow manifolds that are normally hyperbolic with both stable and unstable fast manifolds. We present two examples of bifurcation problems where these manifolds play a key role and a third example in which saddle-type slow manifolds are part of a traveling wave profile of a partial differential equation. Initial value solvers are incapable of computing trajectories on saddle-type slow manifolds, so the slow manifold of saddle type (SMST) algorithm presented here is formulated as a boundary value method. We take an empirical approach here to assessing the accuracy and effectiveness of the algorithm.

John Guckenheimer, Christian Kuehn

This paper investigates travelling wave solutions of the FitzHugh-Nagumo equation from the viewpoint of fast-slow dynamical systems. These solutions are homoclinic orbits of a three dimensional vector field depending upon system parameters of the FitzHugh-Nagumo model and the wave speed. Champneys et al. [A.R. Champneys, V. Kirk, E. Knobloch, B.E. Oldeman, and J. Sneyd, When Shilnikov meets Hopf in excitable systems, SIAM Journal of Applied Dynamical Systems, 6(4), 2007] observed sharp turns in the curves of homoclinic bifurcations in a two dimensional parameter space. This paper demonstrates numerically that these turns are located close to the intersection of two curves in the parameter space that locate non-transversal intersections of invariant manifolds of the three dimensional vector field. The relevant invariant manifolds in phase space are visualized. A geometrical model inspired by the numerical studies displays the sharp turns of the homoclinic bifurcations curves and yields quantitative predictions about multi-pulse and homoclinic orbits and periodic orbits that have not been resolved in the FitzHugh-Nagumo model. Further observations address the existence of canard explosions and mixed-mode oscillations.

Christian Kuehn

Hopf bifurcations in fast-slow systems of ordinary differential equations can be associated with surprising rapid growth of periodic orbits. This process is referred to as canard explosion. The key step in locating a canard explosion is to calculate the location of a special trajectory, called a maximal canard, in parameter space. A first-order asymptotic expansion of this location was found by Krupa and Szmolyan in the framework of a "canard point"-normal-form for systems with one fast and one slow variable. We show how to compute the coefficient in this expansion using the first Lyapunov coefficient at the Hopf bifurcation thereby avoiding use of this normal form. Our results connect the theory of canard explosions with existing numerical software, enabling easier calculations of where canard explosions occur.

Christian Kuehn

Numerical continuation methods for deterministic dynamical systems have been one of the most successful tools in applied dynamical systems theory. Continuation techniques have been employed in all branches of the natural sciences as well as in engineering to analyze ordinary, partial and delay differential equations. Here we show that the deterministic continuation algorithm for equilibrium points can be extended to track information about metastable equilibrium points of stochastic differential equations (SDEs). We stress that we do not develop a new technical tool but that we combine results and methods from probability theory, dynamical systems, numerical analysis, optimization and control theory into an algorithm that augments classical equilibrium continuation methods. In particular, we use ellipsoids defining regions of high concentration of sample paths. It is shown that these ellipsoids and the distances between them can be efficiently calculated using iterative methods that take advantage of the numerical continuation framework. We apply our method to a bistable neural competition model and a classical predator-prey system. Furthermore, we show how global assumptions on the flow can be incorporated - if they are available - by relating numerical continuation, Kramers' formula and Rayleigh iteration.

Maxime Breden, Christian Kuehn

Differential equations with random parameters have gained significant prominence in recent years due to their importance in mathematical modelling and data assimilation. In many cases, random ordinary differential equations (RODEs) are studied by using Monte-Carlo methods or by direct numerical simulation techniques using polynomial chaos (PC), i.e., by a series expansion of the random parameters in combination with forward integration. Here we take a dynamical systems viewpoint and focus on the invariant sets of differential equations such as steady states, stable/unstable manifolds, periodic orbits, and heteroclinic orbits. We employ PC to compute representations of all these different types of invariant sets for RODEs. This allows us to obtain fast sampling, geometric visualization of distributional properties of invariants sets, and uncertainty quantification of dynamical output such as periods or locations of orbits. We apply our techniques to a predator-prey model, where we compute steady states and stable/unstable manifolds. We also include several benchmarks to illustrate the numerical efficiency of adaptively chosen PC depending upon the random input. Then we employ the methods for the Lorenz system, obtaining computational PC representations of periodic orbits, stable/unstable manifolds and heteroclinic orbits.

Christian Kuehn, Patrick Kuerschner

In this work, we study the numerical approximation of local fluctuations of certain classes of parabolic stochastic partial differential equations (SPDEs). Our focus is on effects for small spatially-correlated noise on a time scale before large deviation effects have occurred. In particular, we are interested in the local directions of the noise described by a covariance operator. We introduce a new strategy and prove a Combined ERror EStimate (CERES) for the five main errors: the spatial discretization error, the local linearization error, the noise truncation error, the local relaxation error to steady state, and the approximation error via an iterative low-rank matrix algorithm. In summary, we obtain one CERES describing, apart from modelling of the original equations and standard round-off, all sources of error for a local fluctuation analysis of an SPDE in one estimate. To prove our results, we rely on a combination of methods from optimal Galerkin approximation of SPDEs, covariance moment estimates, analytical techniques for Lyapunov equations, iterative numerical schemes for low-rank solution of Lyapunov equations, and working with related spectral norms for different classes of operators.

Christian Kuehn

The parameter space of dynamical systems arising in applications is often found to be high-dimensional and difficult to explore. We construct a fast algorithm to numerically analyze "quantitative features" of dynamical systems depending on parameters. Using a classical problem from mathematical ecology as an example, we demonstrate how to apply the algorithm to investigate the amplitude of a limit cycle depending on seven parameters. We stress the practical value of the algorithm but we also provide a rigorous error analysis to justify the overall strategy. Our approach turns out to be particularly useful in the case of comparing experimental data to a model defined by differential equations and to investigate whether the equations can approximate the modeled system.

Christian Kuehn, Sara-Viola Kuntz, Tobias Wöhrer

We analyze the universal approximation constraints of narrow Residual Neural Networks (ResNets) both theoretically and numerically. For deep neural networks without input space augmentation, a central constraint is the inability to represent critical points of the input-output map. We prove that this has global consequences for target function approximations and show that the manifestation of this defect is typically a shift of the critical point to infinity, which we call the ``tunnel effect'' in the context of classification tasks. While ResNets offer greater expressivity than standard multilayer perceptrons (MLPs), their capability strongly depends on the signal ratio between the skip and residual channels. We establish quantitative approximation bounds for both the residual-dominant (close to MLP) and skip-dominant (close to neural ODE) regimes. These estimates depend explicitly on the channel ratio and uniform network weight bounds. Low-dimensional examples further provide a detailed analysis of the different ResNet regimes and how architecture-target incompatibility influences the approximation error.

Christian Kuehn

Convergence properties of binary stationary subdivision schemes for curves have been analyzed using the techniques of z-transforms and eigenanalysis. Eigenanalysis provides a way to determine derivative continuity at specific points based on the eigenvalues of a finite matrix. None of the well-known subdivision schemes for curves have complex eigenvalues. We prove when a convergent scheme with palindromic mask can have complex eigenvalues and that a lower limit for the size of the mask exists in this case. We find a scheme with complex eigenvalues achieving this lower bound. Furthermore we investigate this scheme numerically and explain from a geometric viewpoint why such a scheme has not yet been used in computer-aided geometric design.

Christian Kuehn, Sara-Viola Kuntz

Neural Ordinary Differential Equations (Neural ODEs), which are the continuous-time analog of Residual Neural Networks (ResNets), have gained significant attention in recent years. Similarly, Neural Delay Differential Equations (Neural DDEs) can be interpreted as an infinite depth limit of Densely Connected Residual Neural Networks (DenseResNets). In contrast to traditional ResNet architectures, DenseResNets are feed-forward networks that allow for shortcut connections across all layers. These additional connections introduce memory in the network architecture, as typical in many modern architectures. In this work, we explore how the memory capacity in neural DDEs influences the universal approximation property. The key parameter for studying the memory capacity is the product $K τ$ of the Lipschitz constant and the delay of the DDE. In the case of non-augmented architectures, where the network width is not larger than the input and output dimensions, neural ODEs and classical feed-forward neural networks cannot have the universal approximation property. We show that if the memory capacity $Kτ$ is sufficiently small, the dynamics of the neural DDE can be approximated by a neural ODE. Consequently, non-augmented neural DDEs with a small memory capacity also lack the universal approximation property. In contrast, if the memory capacity $Kτ$ is sufficiently large, we can establish the universal approximation property of neural DDEs for continuous functions. If the neural DDE architecture is augmented, we can expand the parameter regions in which universal approximation is possible. Overall, our results show that by increasing the memory capacity $Kτ$, the infinite-dimensional phase space of DDEs with positive delay $τ>0$ is not sufficient to guarantee a direct jump transition to universal approximation, but only after a certain memory threshold, universal approximation holds.

Tobias Wöhrer, Christian Kuehn

We investigate finite-time Lyapunov exponents (FTLEs), a measure for exponential separation of input perturbations, of deep neural networks within the framework of continuous-depth neural ODEs. We demonstrate that FTLEs are powerful organizers for input-output dynamics, allowing for better interpretability and the comparison of distinct model architectures. We establish a direct connection between Lyapunov exponents and adversarial vulnerability, and propose a novel training algorithm that improves robustness by FTLE regularization. The key idea is to suppress exponents far from zero in the early stage of the input dynamics. This approach enhances robustness and reduces computational cost compared to full-interval regularization, as it avoids a full ``double'' backpropagation.

Dennis Chemnitz, Maximilian Engel, Christian Kuehn et al.

In this chapter, we utilize dynamical systems to analyze several aspects of machine learning algorithms. As an expository contribution we demonstrate how to re-formulate a wide variety of challenges from deep neural networks, (stochastic) gradient descent, and related topics into dynamical statements. We also tackle three concrete challenges. First, we consider the process of information propagation through a neural network, i.e., we study the input-output map for different architectures. We explain the universal embedding property for augmented neural ODEs representing arbitrary functions of given regularity, the classification of multilayer perceptrons and neural ODEs in terms of suitable function classes, and the memory-dependence in neural delay equations. Second, we consider the training aspect of neural networks dynamically. We describe a dynamical systems perspective on gradient descent and study stability for overdetermined problems. We then extend this analysis to the overparameterized setting and describe the edge of stability phenomenon, also in the context of possible explanations for implicit bias. For stochastic gradient descent, we present stability results for the overparameterized setting via Lyapunov exponents of interpolation solutions. Third, we explain several results regarding mean-field limits of neural networks. We describe a result that extends existing techniques to heterogeneous neural networks involving graph limits via digraph measures. This shows how large classes of neural networks naturally fall within the framework of Kuramoto-type models on graphs and their large-graph limits. Finally, we point out that similar strategies to use dynamics to study explainable and reliable AI can also be applied to settings such as generative models or fundamental issues in gradient training methods, such as backpropagation or vanishing/exploding gradients.

Boumediene Hamzi, Christian Kuehn, Sameh Mohamed

We study the maximum mean discrepancy (MMD) in the context of critical transitions modelled by fast-slow stochastic dynamical systems. We establish a new link between the dynamical theory of critical transitions with the statistical aspects of the MMD. In particular, we show that a formal approximation of the MMD near fast subsystem bifurcation points can be computed to leading-order. In particular, this leading order approximation shows that the MMD depends intricately on the fast-slow systems parameters and one can only expect to extract warning signs under rather stringent conditions. However, the MMD turns out to be an excellent binary classifier to detect the change point induced by the critical transition. We cross-validate our results by numerical simulations for a van der Pol-type model.

Christian Kuehn

We study the Allen-Cahn equation with a cubic-quintic nonlinear term and a stochastic $Q$-trace-class stochastic forcing in two spatial dimensions. This stochastic partial differential equation (SPDE) is used as a test case to understand, how numerical continuation methods can be carried over to the SPDE setting. First, we compute the deterministic bifurcation diagram for the PDE, i.e. without stochastic forcing. In this case, two locally asymptotically stable steady state solution branches exist upon variation of the linear damping term. Then we consider the Lyapunov operator equation for the locally linearized system around steady states for the SPDE. We discretize the full SPDE using a combination of finite-differences and spectral noise approximation obtaining a finite-dimensional system of stochastic ordinary differential equations (SODEs). The large system of SODEs is used to approximate the Lyapunov operator equation via covariance matrices. The covariance matrices are numerically continued along the two bifurcation branches. We show that we can quantify the stochastic fluctuations along the branches. We also demonstrate scaling laws near branch and fold bifurcation points. Furthermore, we perform computational tests to show that, even with a sub-optimal computational setup, we can quantify the subexponential-timescale fluctuations near the deterministic steady states upon stochastic forcing on a standard desktop computer setup. Hence, the proposed method for numerical continuation of SPDEs has the potential to allow for rapid parametric uncertainty quantification of spatio-temporal stochastic systems.

Christian Kuehn

Numerical continuation calculations for ordinary differential equations (ODEs) are, by now, an established tool for bifurcation analysis in dynamical systems theory as well as across almost all natural and engineering sciences. Although several excellent standard software packages are available for ODEs, there are - for good reasons - no standard numerical continuation toolboxes available for partial differential equations (PDEs), which cover a broad range of different classes of PDEs automatically. A natural approach to this problem is to look for efficient gluing computation approaches, with independent components developed by researchers in numerical analysis, dynamical systems, scientific computing and mathematical modelling. In this paper, we shall study several elliptic PDEs (Lane-Emden-Fowler, Lane-Emden-Fowler with microscopic force, Caginalp) via the numerical continuation software pde2path and develop a gluing component to determine a set of starting solutions for the continuation by exploting the variational structures of the PDEs. In particular, we solve the initialization problem of numerical continuation for PDEs via a minimax algorithm to find multiple unstable solution. Furthermore, for the Caginalp system, we illustrate the efficient gluing link of pde2path to the underlying mesh generation and the FEM MatLab pdetoolbox. Even though the approach works efficiently due to the high-level programming language and without developing any new algorithms, we still obtain interesting bifurcation diagrams and directly applicable conclusions about the three elliptic PDEs we study, in particular with respect to symmetry-breaking. In particular, we show for a modified Lane-Emden-Fowler equation with an asymmetric microscopic force, how a fully connected bifurcation diagram splits up into C-shaped isolas on which localized pattern deformation appears towards two different regimes.