1 Introduction↩︎

Deep learning has emerged as a transformative paradigm in modern artificial intelligence, driven by its success in tasks such as computer vision [1]–[3], natural language processing [4]–[6], autonomous driving [7], [8] and medical measurement [9], [10]. These achievements are largely attributed to the confluence of three key factors: the proliferation of large-scale annotated datasets [11], the rapid advancement in computational hardware [12], [13] (e.g., GPUs and TPUs), and the continuous evolution of neural network architectures [14], [15]. Despite these empirical successes, deep learning still faces challenges in interpretability and generalization [16], [17], which limit its applicability in domains requiring rigorous safety and reliability guarantees [18], [19].

ODEs and SDEs have emerged as theoretical frameworks providing novel insights and mathematical tools for tackling the challenges of interpretability and generalization in deep learning. Pioneering work by E et al. [20] conceptualized deep neural networks as discrete dynamic systems, framing the training phase as an optimal control problem. This perspective catalyzed the development of network architectures grounded in deterministic differential equations, encompassing ODEs and partial differential equations (PDEs), and leading to substantial advancements [21]–[24]. Stability in ODE-based models is often achieved by enforcing Lipschitz or contraction conditions on the function \(f\), which ensures reliable forward and backward propagation and mitigates gradient explosion or vanishing. Overall, existing research can be broadly categorized into two main directions: 1) leveraging differential equation frameworks, including ODEs, PDEs, and SDEs, for the theoretical interpretation of deep networks [20], [24]–[31]; 2) guiding the design of network architecture with DEs [22], [23], [32]–[39].

The intrinsic interpretability of ODEs arises from their deterministic nature; given an initial state, systems evolve along a single, predictable trajectory. This transition from discrete layer stacking to continuous evolution enhances the model’s physical intuitiveness and interpretability. For example, ODE-based models like ResNets [40] and LM-ResNets [22] interpret inter-layer evolution as a discretization of continuous flows, whereas PDE-based models such as PDE-CNNs [23] encode geometric structures and local features through spatial derivatives. Additionally, deterministic DEs facilitate deterministic generative modeling, including Neural ODEs [21] and flow matching methods [41], [42]. These frameworks generate data by solving initial value problems or construct continuous trajectories between distributions by deriving velocity fields through optimal transport. Furthermore, this continuous perspective extends to scientific machine learning through Neural PDE solvers, which operate under specified equation constraints to find numerical solutions for complex physical systems. These approaches illustrate the role of deterministic DEs in the unification of neural architectures, generative modeling and equation solving.

SDEs extend ODEs by incorporating a stochastic noise term, transforming deterministic systems into stochastic ones. This stochasticity induces trajectory variability from identical initial states, capturing richer distributional variations to improve generalization and enable the generation of diverse, high-quality samples. Randomness is pervasive in deep neural networks, arising from weight initialization [43], [44], minibatch sampling [45], and stochastic regularization methods such as Dropout [46]–[49] and Gaussian noise injection [50]–[52]. These strategies influence optimization dynamics, robustness, and generalization. Recent studies employ SDEs to rigorously analyze this randomness, forming two primary research trajectories. One trajectory focuses on stochastic regularization: explicitly embedding randomness via SDE frameworks (e.g., Dropout [53], [54] and stochastic depth [55]) to provide coherent mathematical models for network evolution under perturbations [22], [39], [56]. Another trajectory employs SDEs to model network dynamics [57]–[66], comprising two categories: Neural stochastic differential equations (Neural SDEs) [62] and SDE-guided generative models, exemplified by Denoising Diffusion Probabilistic Models (DDPMs) [57].

Taking a broader perspective, ODEs, PDEs, and SDEs collectively span a systematic theoretical hierarchy, ranging from deterministic temporal evolution to complex stochastic spatiotemporal systems. As illustrated in Fig. 1, these differential equations serve as rigorous foundations for network design and regularization. Furthermore, Fig. 2 demonstrates how deep learning leverages these continuous dynamics to either drive data evolution for generative modeling or compute accurate numerical solutions. We defer a comprehensive review of the fundamental neural architectures and mathematical prerequisites to Section S1 of the online supplementary material, thereby ensuring the main narrative remains concise and focused.

This survey limits its scope to model-level paradigms whose architectural topologies, generative continuous dynamics, or algorithmic equation-solving procedures are explicitly governed by ODEs, PDEs, or SDEs. Consequently, layer-level innovations such as Structured State-Space Models (e.g., S4 [67], Mamba [68]) and Liquid Time-Constant Networks (LTCs) [69] fall outside our primary focus. We refer interested readers to [70] for a comprehensive review of these layer-centric architectures. Existing surveys predominantly categorize models by specific domains, techniques, or isolated equation types, as detailed in 1. For the systematic organization of this field, we introduce a comprehensive framework based on dual-perspective classification criteria. First, we classify the literature by mathematical foundations: ODEs, PDEs, and SDEs. Second, we categorize models by methodological roles into two main areas. The first regime, Networks Guided by Equations, encompasses DE-driven architectures and stochastic regularization. The second regime, Networks for Learning Equations, subdivides into neural solvers for finding numerical solutions and networks learning data-driven dynamics, such as Neural DEs and DE-guided generative models.

In summary, this survey provides a systematic taxonomy and a holistic review of deep learning paradigms underpinned by differential equations. Specifically, our methodological classification successfully organizes the literature into two primary functional regimes. The first regime, Networks Guided by Equations, encompasses differential equation-driven architectures and stochastic regularization techniques. The second regime, Networks for Learning Equations, is further subdivided into neural numerical solvers and data-driven continuous dynamics. Beyond this structural taxonomy, we perform a rigorous qualitative meta-analysis of these DE-centric models across various standardized evaluation benchmarks. Finally, we synthesize prevailing open challenges and chart prospective research trajectories to navigate this rapidly evolving domain.

2 From Deterministic DEs to Deep Neural Networks↩︎

Since the pioneering work by E et al. [20] that reinterpreted the essence of neural networks from a dynamical systems perspective, DE-guided neural networks design has emerged as an active research direction [21]–[24], [27], [80]–[82]. This perspective not only established a rigorous mathematical framework for deep learning but also led to new architectures that combine robustness and interpretability, substantially enhancing model adaptability in complex real-world applications.

To systematically incorporate the mathematical rigor of deterministic DEs into deep learning, current mainstream approaches can be categorized into two key paradigms:

• Deterministic DE-Driven Architectures (Networks Guided by Equations): Exemplified by ResNets [40] and LM-ResNets [22], this paradigm formulates neural network forward propagation as solving initial value problems. By treating layers as discrete deterministic DE approximations via Euler or Runge-Kutta schemes [83], this perspective elucidates the networks’ function approximation capabilities and guides novel architectural design.

• Networks for Learning Deterministic Equations: This paradigm deploys deep learning to directly learn or compute continuous dynamical systems. For data representation, deterministic dynamics models like Neural ODEs [21] and ODE-guided generative models [41] describe the continuous evolution of data in latent space, leveraging temporal continuity for flexible trajectory interpolation and precise gradient computation. Extending beyond data modeling, this paradigm also encompasses Neural PDE solvers [84], [85], which constrain networks with specific physical laws to provide efficient numerical solutions for complex physical systems.

2.1 Deterministic DE-Driven Network Architectures↩︎

We first present a structured classification of deterministic DE-driven approaches for deep learning architecture design. A detailed summary of representative works is provided in 2, and the corresponding architectural illustrations are presented in Fig. S1 of the online supplementary material.

2.1.1 ResNets [40]↩︎

ResNets revolutionized deep learning by introducing skip connections to overcome the degradation problem in very deep DNNs. The fundamental building block implements a residual transformation through the update rule: \[\label{ResNet1} x_{n+1} = x_n + f(x_n),\tag{1}\] where \(x_n\) is the layer \(n\) feature representation and \(f\) is a learnable convolutional transformation. Such skip connections stabilize gradient propagation and facilitate the training of ultra-deep architectures like ResNet-50.

Beyond their empirical success, ResNets have been rigorously analyzed from a dynamical systems perspective. Haber et al. [82] and E et al. [20] established that the residual update in 1 can be viewed as a forward Euler discretization of an underlying continuous-time ODE: \[\label{ODERN} \frac{dx(t)}{dt} = v(x,t),\,\, x(0)=x_0,\,\, t\in [0,T],\tag{2}\] where \(v(x,t)\) denotes a velocity field satisfying \(v(x, t) = f(x) / \Delta t\). This connection enables the numerical analysis of network stability and inspires continuous-depth models using ODE solvers for adaptive computation, motivating ResNet variants and advancing our theoretical understanding.

Since the introduction of ResNets, several variants have been proposed to enhance its performance. Architecturally, Wide ResNet (WRN) [87] and ResNeXt [26] increase network width and cardinality, while ResNet-S and ResNet-D [32] focus on optimizing downsampling paths. These adaptations balance depth and width for more efficient training. Numerically, variants like Residual-in-Residual Network (RiR) [28], Polynomial Network (PolyNet) [81], and Fractal Network (FractalNet) [27] explore different formulations to improve stability and expressiveness. Methods such as leapfrog block (LF-Block) [33], Runge-Kutta-2 block (RK2-Block) [33], and higher-order RK4/RK8 [34] integrate time-stepping techniques from numerical analysis to enhance learning accuracy and convergence.

2.1.2 Linear Multistep ResNets [22]↩︎

LM-ResNets introduce multistep dependencies inspired by Linear Multistep Methods (LMMs). Unlike standard ResNets that update solely from the current state \(x_n\), LM-ResNets compute \(x_{n+1}\) as a neural-parameterized nonlinear function of \(x_n\) and past states (e.g., \(x_{n-1}\)). The layer update follows the general form of LMMs: \[x_{n+1} = (1-k_n)x_n+k_n x_{n-1}+f_n(x_n),\] where \(k_n\in\mathbb{R}\) is a trainable parameter. LM-ResNets can be interpreted as a high-order ODE and written as: \[f(x_n,t)=(1+k_n)\dot{x}_n + (1-k_n)\frac{\Delta t}{2}\ddot{x}_n.\] Subsequently, HNO [35] generalizes deep unfolding networks as ODE solvers, employing linear multistep methods to accelerate convergence and improve accuracy. Concurrently, DenseNet [29] introduces a densely connected topology that concatenates features from all preceding layers, enhancing both feature reuse and gradient flow.

2.1.3 Reversible Residual Networks (RevNets) [24]↩︎

Reversible architectures, such as RevNets [24], invertible Residual Networks (i-ResNets) [30], and invertible Reversible Networks (i-RevNets) [88], address the memory bottleneck of traditional ResNets by reconstructing activations during backpropagation. This eliminates the need to store intermediate states within residual blocks, which can be formally expressed as: \[\label{Eq:RevNet} y_{n+1} = y_n + f(x_n),\, x_{n+1} = x_n + g(y_{n+1}),\tag{3}\] where \(f\) and \(g\) are learnable residual functions. The inputs are reconstructed as: \[x_n = x_{n+1} - g(y_{n+1}),\, y_n = y_{n+1} - f(x_n).\] RevNets reduce memory usage with limited performance loss but incur higher computational cost. I-ResNets [30] achieve full-network invertibility for both generative and discriminative tasks, while i-RevNets [88] add nonlinear mappings to further improve memory efficiency. Hamiltonian Networks [36], [86] extend RevNets by modeling reversible dynamical systems, enabling stable, memory-efficient training and accurate physical simulations.

2.1.4 PDE-CNNs [23]↩︎

In 1 , the function \(f\) is treated as a black box without structural constraints. For tasks involving structural data (e.g., speech, image, or video processing) where local spatial interactions are crucial, layers are standardly instantiated via convolutions. Ruthotto et al. [23] proposed a PDE-based framework that interprets the forward propagation of residual CNNs as a continuous-time nonlinear diffusion process. Formally, given the initial state as the input, the underlying continuous dynamics are modeled via the forward Euler method: \[\frac{\partial x(\theta,t)}{\partial t} = f(\theta,x,t),\, x(\theta,0)=x_0,\,f(\theta,x,t) = K_2(\theta^{(3)})\,\sigma\Big( \mathcal{N}\big(K_1(\theta^{(1)})Y,\theta^{(2)} \big) \Big)\] where \(f(\theta,x,t)\) represents convolution, nonlinearity, and normalization, with specific parameterizations for each operation. Specifically, \(K_1(\theta^{(1)})\) and \(K_2(\theta^{(3)})\) denote convolutional operators parameterized by \(\theta^{(1)}\) and \(\theta^{(3)}\), respectively. \(\mathcal{N}\) represents a normalization layer parameterized by \(\theta^{(2)}\), and \(\sigma\) denotes a nonlinear activation function (e.g., ReLU). PDE-CNNs [23] encompass parabolic and hyperbolic models, including Hamiltonian and second-order variants. While PDE formulations improve stability and interpretability, deep architectures often incur high computational costs. To mitigate this, LeanConvNets [37] and Global Feature Layers [38] reduce structural complexity. Furthermore, extending PDE frameworks to ResNets and RNNs enhances parameter efficiency [31], and Convection-Diffusion Networks (COIN) [39] further demonstrate robust generalization against perturbations.

Deterministic DE-driven neural architectures can be understood through a unified trade-off among interpretability, stability, expressivity, and computational cost, governed by how the underlying dynamical system constrains state evolution. First-order ODE-inspired models use simple single-step updates, offering strong stability and efficiency but limited expressivity, requiring deep compositions for complex mappings. Higher-order ODEs introduce multi-step dependencies, increasing degrees of freedom and representation capacity at the cost of additional state and computation. DE systems impose structured dynamics (e.g., invertibility or conservation laws), improving interpretability and memory efficiency while potentially limiting flexibility. PDE-based architectures extend dynamics to the spatial domain, enabling effective spatio-temporal modeling but incurring higher computational cost due to spatio-temporal discretization. Overall, improving expressivity and interpretability generally requires relaxing structural simplicity and increasing computational complexity, highlighting the need to match the equation type to the target task.

2.2 Networks for Learning Deterministic Equations↩︎

As illustrated in Fig. 2, we partition deterministic architectures designed to learn or solve continuous equations into two primary domains. The first encompasses dynamics modeling via Neural DEs [21] and ODE-guided generative models [41], while the second focuses on Neural PDE solvers. Neural DEs [21] are particularly effective for modeling sequential and continuous-time data, as they learn the governing dynamics \(f_\theta(x,t)\) through neural network parameterization. In contrast, ODE-guided generative models [41] focus on learning the velocity field \(v_\theta(x, t)\) in the Liouville equation via neural networks, thereby enabling direct control over the evolution of probability densities. When solving equations, Neural PDE solvers [84], [120], [121] leverage deep learning to compute numerical solutions for complex physical systems by enforcing specific equation constraints. The details and taxonomy of Neural DEs, representative formulations and recent advances in ODE-guided generative models, and the core architectures of Neural PDE solvers are summarized in 3, [tab:flow95model] and 4, respectively.

2.2.1 Neural Deterministic Differential Equations↩︎

Integrating deterministic differential equations as continuous-depth architectural backbones provides a powerful deep learning paradigm for enhancing spatiotemporal physical simulations and complex sequential representation learning.

Neural ODEs [21] transcend the rigid, discrete layers of traditional ResNets. By modeling continuous-time dynamics, they enable adaptive computation and formulate hidden state evolution via an ODE: \[\frac{dx(t)}{dt} =f_\theta(x,t),\] where \(f_\theta\) parameterizes the continuous vector field. This formulation replaces discrete layers and solves the hidden state trajectory via an ODE, guaranteeing memory-efficient training through the adjoint method [122]. The continuous-time paradigm offers advantages such as constant memory cost regardless of model depth, suitability for irregularly sampled time series [21], [123], and adaptive computational depth. However, Neural ODEs come with drawbacks, including high computational cost during training and inference, as well as challenges with numerical instability and sensitivity to integration dynamics [21], [89], [124], [125].

Second-Order Neural ODEs (SONODEs) [90], [126] capture second-order dynamics by incorporating initial velocity and acceleration. This physics-inspired inductive bias yields more faithful representations of physical systems, enhancing both interpretability and numerical stability. Specifically, we consider SONODEs whose initial position \(x(t_0)\), initial velocity \(\dot{x}(t_0)\), and acceleration \(\ddot{x}\) are given by: \[\ddot{x}=f_\theta(x,\dot{x},t),\,\dot{x}(t_0)=g_\theta(x(t_0)),\,x(t_0)=X_0,\] where \(g_\theta\) and \(f_\theta\) are neural networks predicting the initial velocity and acceleration. While convertible to a first-order ODE for training, this second-order formulation demonstrates superior physical interpretability, robustness, and extrapolation over ANODEs. High computational cost during training and inference remains a key challenge for Neural ODEs. Consequently, a range of methods have been proposed to improve their efficiency [89], [90], [124], [127]–[136]. In parallel, to enhance the stability and robustness of Neural ODEs under long-term integration and noisy inputs, several approaches have been developed [62], [137]–[145].

Building upon the success of SONODEs, continuous deep architectures have naturally extended to spatiotemporal systems governed by PDEs. Shifting from traditional black-box dynamics fitting to white-box PDE discovery, NeuPDE [92] incorporates strong physical inductive biases. Combined with sparse regularization, it explicitly identifies the underlying governing equations: \[\frac{\partial u(t, x)}{\partial t} = \mathcal{N}_\theta (\mathcal{D}(u)) \approx \sum_{k=1}^{K} c_k \phi_k(u, \nabla u, \nabla^2 u, \dots),\] where \(c_k\) represent learnable sparse coefficients and \(\phi_k\) denotes a library of candidate functions, comprising spatial derivatives and their nonlinear products. This formulation rigorously preserves physical consistency and enables direct mathematical extraction. Driven by a similar philosophy, PDE-Net [146], [147] establishes mathematical mappings via constrained filter moment matrices to learn differential operators, achieving robust long-term forecasting and discovery of hidden PDE structures.

Standard Neural ODEs are limited in modeling irregularly sampled time series and in capturing complex temporal dynamics. Thus, researchers have developed several important extensions. Among these, Neural Controlled Differential Equations (Neural CDEs) [83] formulate the hidden state evolution as a controlled differential equation, providing: \[z_t = z_{t_0} + \int_{t_0}^t f_\theta (z_s)\,dX_s, \,t\in(t_0,t_n),\] where \(z(t)\) denotes the hidden state, \(X_s\) is an interpolated input path and \(f_\theta\) controls how the state responds to input variations. Neural CDEs enable learning from irregularly sampled data by treating inputs as continuous trajectories. Neural Rough Differential Equations (Neural RDEs) [93] extend these methods to rough or non-differentiable paths using rough path theory [148], capturing high-frequency or noisy dynamics through higher-order integration terms. Neural Fractional Differential Equations (Neural FDEs) [94] incorporate long-range memory and non-local dependencies by modeling \[D^\alpha z(t) = f_\theta(z(t),t),\,0<\alpha<1,\] where \(D^\alpha\) is a Caputo or Riemann–Liouville fractional derivative. Neural FDEs consider history-dependent dynamics, making them suitable for non-Markovian systems.

As summarized in 3, continuous architectures evolved to overcome the theoretical and numerical limits of early deterministic models. While standard Neural ODEs [21] pioneered continuous-depth modeling, their topology-preserving nature and Markovian assumptions restrict physical realism and induce severe numerical stiffness. Architectures resolve this along two trajectories. The first injects physical inductive biases to create interpretable models, utilizing augmented state spaces in ANODEs [89] and second-order fields in SONODEs [90]. The second relaxes autonomous constraints to capture history-dependent dynamics via control paths in Neural CDEs [83] or fractional derivatives in Neural FDEs [94]. Ultimately, enhancing physical interpretability and memory shifts the computational bottleneck toward complex numerical integration, demanding strict alignment between the DE framework and computational resources.

2.2.2 ODE-Guided Generative Models↩︎

In generative modeling, deterministic differential equations serve as optimal transport frameworks that smoothly map simple base distributions to complex target data manifolds.

ODE-guided Generative Models [41], [42] deterministically transport a simple base distribution to a complex target distribution. This transformation is achieved by learning a time-dependent velocity field \(v_\theta(x,t)\) that governs the continuous evolution of data points via an ODE: \[\frac{dx}{dt} = v_\theta(x,t),\,x(0)\sim p(0),\,x(1)\sim p_{\text{target}} (x).\] Under this deterministic flow, probability density evolution follows the Liouville equation. Flow Matching [41] introduces a simulation-free objective to overcome the expensive ODE simulations traditionally required to track this density. It constructs a target velocity field via conditional path interpolation and trains \(v_\theta\) to regress it directly, thereby completely bypassing explicit density computation and enabling scalable high-dimensional training.

Standard flow matching is restricted to continuous Euclidean spaces [97] and linear trajectories [105], [149]. To overcome these limitations, subsequent research advances along two primary directions: extending the framework to non-Euclidean topologies, such as discrete domains [98], [150], graphs [151], and manifolds [97], [152], and optimizing transport paths [104], [105] to accelerate sampling, ultimately enabling one-step generation [109].

Generalized Flow Matching: A series of generalized flow matching approaches have emerged, each systematically extending the framework toward broader classes of data and structural priors. Notable examples include:

• Conditional Flow Matching (CFM) [41] is a method that learns a velocity field \(v_\theta(x, t, y)\), which is conditioned on an auxiliary variable \(y\), allowing for controllable generation. The system is governed by the following ODE: \[\frac{dx}{dt} = v_\theta(x, t),\, x(0) \sim p(0), \, x(1) \sim p_{\text{target}}(x|y),\] which facilitates the direct generation of samples with specified attributes to ensure alignment with conditioning information. This capability is particularly beneficial for tasks such as class-conditional image synthesis.

• Discrete Flow Matching (DFM) [98], [150] generalizes flow matching to discrete and categorical data domains by replacing the continuous flow with a continuous-time Markov jump process. The dynamics are governed by a time-dependent rate matrix \(Q_\theta(t)\): \[\frac{dp(t)}{dt} = Q_\theta(t)^\top p(t),\] where \(p(t)\) denotes the time-evolving probability mass over discrete states. This approach enables the efficient, non-autoregressive generation of sequence, graph, and molecular data.

These developments mark a shift in flow matching research toward modeling data in non-Euclidean [97], [152] and structured spaces [99], [100], [102]. Riemannian Flow Matching (RFM) [97] adapts flow matching to non-Euclidean domains by respecting manifold constraints, enabling generative modeling on curved data spaces. Generator Matching [101] generalizes flow matching to arbitrary modalities and Continuous Time Markov Processes (CTMPs), supporting non-autoregressive generation for discrete sequences and graphs. Methods like Continuous-State Discrete Flow Matching (CS-DFM) [102], Fisher-Flow [99], and Statistical Flow Matching (SFM) [100] bring in information geometry, aligning probability transport with natural statistical metrics such as the Fisher–Rao distance. Diffeomorphism Conditional Flow Matching (DiffeoCFM) [103] further extends flow matching to positive-definite matrix manifolds via pullback geometry. Altogether, these innovations substantially broaden the applicability and expressivity of flow matching, providing a robust geometric foundation for ODE-guided generative modeling in complex domains.

Trajectory Optimization: Trajectory Optimization explores advanced strategies to optimize both the generative trajectory and its numerical integration, aiming to improve sampling efficiency, stability, and accuracy [153]–[155].

• Rectified Flow [104], [105] simplifies the generative process by constraining the transport path between noise and data to a straight line in data space, allowing for highly efficient and deterministic sampling. The underlying ODE is \[\frac{dx}{dt} = x(1) -x(0),\] where \(x(0)\) and \(x(1)\) are the endpoints of the trajectory, allowing the model to learn a velocity field that maps base samples directly along a rectified path.

• MeanFlow [109], [156] adopts a one-step generation approach by learning a time-averaged velocity field \(\overline{v}_\theta(x, t)\). This enables the model to directly map a base sample to the target in a single step: \[x(1) = x(0) + \int_0^1 v_\theta(x, t)\, dt \approx x(0) + \overline{v}(x(0)).\]

Within this trajectory-optimization perspective, subsequent studies mainly revisit three closely related design issues: how to train the transport dynamics more robustly, how to choose more expressive transport paths, and how to reduce the cost of numerical sampling. Work on training objectives and optimization dynamics addresses issues such as velocity variance and trajectory crossings, which can affect training stability and sampling efficiency [157]–[161]. Path-oriented studies move beyond a single global straight-line transport by introducing hierarchical, piecewise, or more flexible trajectories, allowing the generative path to better adapt to complex data geometries [107], [162], [163]. Efficiency-oriented methods, including InstaFlow [106], Optimal Flow Matching [164], and SlimFlow [110], seek to reduce the number of function evaluations through distillation, shortcut mappings, or optimized transport objectives. Taken together, these studies show that trajectory optimization in flow-matching models has progressed from fixed straight-line transports toward more stable, adaptive, and computationally efficient transport formulations [110], [165], [166].

2.2.3 Neural PDE solvers↩︎

PDEs are fundamental for modeling physical and engineering systems. Although traditional numerical methods such as FEM and FVM are reliable, they often suffer from high computational cost in high-dimensional or parametric settings. Neural PDE solvers address this limitation by using DNNs to approximate specific PDE solutions or the underlying solution operators.

Physics-Informed Neural Networks (PINNs). PINNs [120], [167]–[169] incorporate governing equations into network training to approximate either specific PDE solutions. For a PDE \(\mathcal{N}[u](x,t)=0\), the surrogate solution \(u_{\theta}\) is optimized by minimizing \[\label{eq:loss} \mathcal{L} = \omega_{\text{pde}}\mathcal{L}_{\text{pde}} + \omega_{\text{bound}}\mathcal{L}_{\text{bound}} + \omega_{\text{init}}\mathcal{L}_{\text{init}} \quadwith~~ \mathcal{L}_{\text{pde}} = \frac{1}{N_r}\sum_{i=1}^{N_r} \left| \mathcal{N}[u_{\theta}](x_i,t_i) \right|^2.\tag{4}\] The main challenge lies in optimization instability caused by PDE stiffness and imbalance among loss terms. Current improvements mainly follow two directions. (1) Loss balancing strategies: adaptive weights [170] are introduced using NTK theory [171], self-adaptive mechanisms [172], [173], and gradient-enhanced methods [174], [175]. (2) Variational and energy-based formulations: weak-form PINNs such as VPINNs [176] minimize \[\label{eq:variation95loss} \mathcal{L}_{\text{pde}} = \frac{1}{K}\sum_{i=1}^{K} \left| \int_{\Omega} \varphi_i(x)\mathcal{N}[u_{\theta}](x)\,dx \right|^2,\tag{5}\] and are further extended by local test functions [177], adversarial weak forms [178], and energy minimization [179].

Operator Learning. Unlike PINNs, which approximate solutions for specific instances, operator learning aims to learn the mapping \(\mathcal{G}: \mathcal{F} \to \mathcal{U}\) between infinite-dimensional function spaces \(\mathcal{F}\) and \(\mathcal{U}\) [121]. Given \(N\) input-output function pairs and \(P\) query points per sample, the model is typically optimized via the empirical supervised loss: \[\mathcal{L}_{\text{data}} = \frac{1}{NP} \sum_{i=1}^{N}\sum_{j=1}^{P} \left| \mathcal{G}_{\theta}(f^{(i)})(y_j^{(i)}) - \mathcal{G}(f^{(i)})(y_j^{(i)}) \right|^2.\]

DeepONet [84] leverages the universal approximation theorem for operators by decomposing the mapping into branch and trunk networks. The operator is represented as a generalized basis expansion: \[\label{eq:basis95epansion} \mathcal{G}_{\theta}(f)(y) = \sum_{k=1}^{p} b_k\big(f(x_1),\dots,f(x_m)\big)t_k(y)+b_0.\tag{6}\] Notable extensions include POD-DeepONet [180] for reduced-order modeling, GreenONet [181] based on Green’s functions, and variants tailored for geometry or sequential data such as Geom-DeepONet [182] and Sequential DeepONet [183].

The Fourier Neural Operator (FNO) [85] parameterizes the integral kernel in the spectral domain, enabling the capture of global dependencies. The architecture updates the hidden state \(v_k\) through a combination of a linear transform and a spectral convolution: \[v_{k+1}(y) = \sigma \left( W v_k(y) + \mathcal{F}^{-1} \big( R_k \cdot \mathcal{F}(v_k) \big) (y) \right),\] where \(\mathcal{F}\) denotes the Fourier transform and \(R_k\) is a learnable weight matrix in the frequency domain. This framework has been generalized to other functional bases, including wavelets, polynomials, and Laplace transforms [184]–[187].

Physics-Informed Operator Learning [188] integrates operator learning with the PDE residual as a regularization ter,m to enforce physical constraints \[\label{eq:pino95loss} \mathcal{L} = \mathcal{L}_{\text{data}} + \mathcal{L}_{\text{pde}},\tag{7}\] where \(\mathcal{L}_{\text{pde}}\) penalizes violations of the underlying governing equations. This paradigm enhances physical consistency and generalization in data-scarce regimes. Subsequent developments include PINO [189], PIWNO [190], and multi-resolution frameworks [191], [192].

Overall, Neural PDE solvers have evolved from pointwise residual learning to operator-level approximation, becoming a core direction in scientific machine learning. PINNs require little labeled data and preserve physical interpretability, but often suffer from optimization difficulties and poor scalability. Operator learning enables rapid inference and strong generalization across parameter spaces, but usually depends on large high-quality datasets. Physics-informed operator learning provides a promising compromise by combining data efficiency with physical consistency. A detailed summary is provided in 4. Furthermore, while the above methods assume a fully known governing PDE, Universal Differential Equations (UDEs) [193] have been developed for scenarios where the PDE is only partially known, learning the unknown terms by embedding trainable sub-networks directly into the solution procedure. The rapid development of these solvers is driven by the growing demand for fast, scalable, and reliable surrogate models in computational science, especially for high-dimensional PDEs and multi-query simulations.

3 From SDEs to Deep Neural Networks↩︎

Incorporating stochasticity into ODEs yields SDEs, providing a more expressive framework for modeling uncertainty. These equations typically feature additive noise to account for external disturbances, or multiplicative noise to capture state-dependent randomness like thermal variations. The integration of SDEs with neural networks has established a new machine learning paradigm, enabling powerful regularization techniques and the principled modeling of stochastic dynamics. This intersection primarily manifests in two domains:

• Stochastic Network Regularization (Networks Guided by Equations): Controlled noise is injected at various network levels to enhance robustness and generalization. This includes neuron-level Dropout [46], [54], weight perturbations [195], [196], activation noise [50], [197], [198], and dynamic depth adjustment [55], [199]. By mimicking stochastic perturbations in the optimization landscape, these methods help avoid poor local minima.

• Stochastic Dynamics Modeling: Forward propagation in neural networks can be mathematically characterized as a dynamical system driven by SDEs. Neural SDEs [62] parameterize both drift and diffusion coefficients to model complex stochastic dynamics. In contrast, SDE-guided generative models (e.g., DDPMs) [112], [200] define a fixed forward noise process and generate data by learning the drift term of the reverse SDE. Fundamentally, they differ in their learning objectives: Neural SDEs capture forward-time stochastic dynamics, whereas DDPMs learn reverse-time denoising dynamics.

3.1 Stochastic Network regularization↩︎

We begin by reviewing stochastic regularization methods that introduce randomness into different components of neural networks, such as neurons, weights, activations, and architectural paths. Representative methods are summarized in 5, along with their key design choices and regularization effects.

3.1.1 Dropout [53], [54]↩︎

Dropout is a fundamental regularization technique that randomly deactivates hidden units during training to mitigate overfitting. When integrated with ResNets, the forward propagation through a dropout-augmented residual block is formally expressed as: \[x_{n+1} = x_n + f(x_n,\omega_n)\odot \frac{z_n}{p},~~z_n\sim\text{Bernoulli(p)},\] where \(p\) is the retention probability, \(\odot\) denotes element-wise multiplication. Dropout has been theoretically interpreted as an approximation to Bayesian inference [214], with further formalization accomplished through a SDE [39], [215], given by: \[dx(t) = f(x)dt + \sqrt{\frac{1-p}{p}} \,f(x)\cdot d\mathcal{W}(t),\] where \(\mathcal{W}(t)\) has the same meaning as in SDEs. This perspective enables stability analysis and supports adaptive optimization.

Since the introduction of standard dropout, numerous neuron-level stochastic regularization variants have emerged. Fast Dropout [201] leverages Gaussian approximation for efficiency, while Annealed [202] and Curriculum Dropout [204] adjust noise levels over training. Variational Dropout [208] treats weights as Gaussian variables for adaptive regularization and sparsity. Other extensions include structure-aware methods such as SpatialDropout [203], Max-pooling and Convolutional Dropout [216], as well as selective masking strategies like Targeted Dropout [205], DropBlock [206], and EDropout [217].

3.1.2 DropConnect [207]↩︎

DropConnect regularizes neural networks by randomly masking weights rather than neurons. It applies a binary mask \(M\), sampled as \(M_{i,j} \sim \text{Bernoulli}(1-p)\), to the weight matrix \(W\), producing an output: \[y = f_a((W \odot M)\,x),\,x_{n+1}^i = \sum_j (W_{ij}\, x_{n}^j) M_{ij}\] where \(\odot\) denotes element-wise multiplication and \(f_a\) is a non-linear activation function. Theoretically, DropConnect reduces to Dropout in shallow networks [218] and can be represented by the following SDE: \[d x(t) = pW x\, dt + \sqrt{p(1 - p)}\, \left(W \odot W\right)^{1/2} \left(x \odot x\right)^{1/2}d\mathcal{W}(t),\] which has inspired numerous extensions, including Sparsity-aware masking [195], Output uncertainty modeling via Dirichlet distributions [219], Gradient-based dynamic drop rates [220], Adaptive drop rates using generalization gap and complexity measures [221].

3.1.3 Gaussian Noise Injection (GNI) [197]↩︎

GNI enhances model generalization by perturbing intermediate representations during training. In a standard ResNet with \(L\) residual blocks, the output of the \(n\)-th block with injected noise is given by: \[x_{n+1} = x_n + f(x_n) + \alpha\,\mathcal{N}(0,I),\] where \(\alpha\) controls the noise intensity and \(\mathcal{N}(0,I)\) is standard Gaussian noise. By setting \(\alpha = \sigma \sqrt{\Delta t}\) and interpreting \(f(x_n)\) as a discretized velocity field \(v(x_n, t_n)\), this update rule corresponds to the Euler–Maruyama discretization of SDEs: \[dx(t) = v(x,t)\,dt + \sigma\,d\mathcal{W}(t).\] Whiteout [198] injects adaptive Gaussian noise into activations to promote structured sparsity, an approach later grounded in a variational interpretation [222]. Theoretical analyses reveal that such activation-level noise penalizes high-frequency components [51], acts as ridge regularization in the infinite-sample limit [223], and benefits overparameterized regimes [224]. To enhance adversarial robustness, Parametric Noise Injection (PNI) [50] applies learnable noise to activations or weights. Ultimately, Gaussian noise injection serves as a principled regularization mechanism that smooths the loss landscape and bridges discrete training with continuous-time stochastic dynamics [224], [225].

3.1.4 Stochastic Depth [55]↩︎

Stochastic Depth mitigates overfitting by randomly bypassing residual blocks during training. For an \(L\)-layer residual network, the forward propagation at the \(n\)-th block is formulated as: \[x_{n+1} = x_n + z_n \,f(x_n),\quad z_n\sim\text{Bernoulli(p)},\] where \(p\) denotes the layer-wise survival probability. The network with stochastic drop can be written as an SDE \[dx(t) = p f(x) dt + \sqrt{p(1 - p)} F(x) d\mathcal{W}(t).\]

Recent studies have explored the theoretical and empirical properties of stochastic depth [226]. Building on foundational block-level stochasticity, subsequent research extends dropping mechanisms through fine-grained structural perturbations and dynamic probability scaling. Structurally, PyramidSepDrop [210] independently drops channels, DropPath [27] discards specific residual paths, and Swapout [209] perturbs both identity and residual branches. To balance regularization with training stability, survival probabilities are dynamically calibrated: Scheduled DropPath [212] modulates rates over training time, while LayerPath [213] scales them across network depth. Furthermore, ShakeDrop [211] introduces stochastic path blending to stabilize single-branch architectures. Collectively, these progressive refinements underscore the critical role of structured noise in robust network optimization and generalization.

3.2 Stochastic Dynamics Modeling↩︎

SDEs provide a powerful mathematical framework for modeling random dynamical systems in deep learning. Two major classes of SDE-driven approaches have emerged. The first, Neural SDEs [62] utilize neural networks to parameterize both drift and diffusion terms, providing a flexible framework for modeling complex stochastic dynamics and effectively capturing temporal uncertainty in time-series analysis. A detailed categorization and summary of these methods can be found in 3. The second class consists of SDE-guide generative models, represented DDPMs [57]. The evolution of probability densities is described by the FPE, and the sampling process can be expressed using the SDE. DDPMs learn the score or noise component to guide the evolution of probability densities, avoiding explicit density computation. Detailed comparisons are summarized in [tab:flow95model].

3.2.1 Neural SDEs [62]↩︎

Neural SDEs extend Neural ODEs [21] by introducing stochasticity through a learnable diffusion term, enabling the modeling of systems with intrinsic randomness. Unlike the deterministic trajectories of Neural ODEs, Neural SDEs capture uncertainty via stochastic processes parameterized by neural networks. The latent state \(x_t\) evolves according to the Itô SDE \[\label{SDE1} dx(t) = f_\theta(x,t)\,dt+g_\theta(x,t)\,d\mathcal{W}(t),\tag{8}\] where \(f_\theta\) and \(g_\theta\) denote the drift and stochastic driving networks, respectively. Eq. 8 is typically simulated via numerical methods like Euler–Maruyama [227], with gradients computed through reparameterization or adjoint methods [122]. While Neural SDEs excel at modeling noisy or irregular data across time series and generative tasks, they suffer from significant computational overhead and limited interpretability.

Recent advancements in Neural SDEs can be categorized into four main areas. Firstly, they have been developed as expressive generative models: Tzen et al. [228] link them to diffusion-limit approximations of latent Gaussian models, Li et al. [64] extend latent ODEs by integrating stochastic priors, and Kidger et al. [65] adapt both the generator and discriminator of GANs as Neural SDEs. Secondly, in terms of efficiency and stability, efforts have concentrated on enhancing solvers and objectives. Reversible solvers [229] and deterministic approximations [230] have been employed to mitigate integration overhead, while alternative objectives such as signature kernel losses [231] and marginal-matching criteria [232] bolster training robustness. Complementary improvements include cubature methods on Wiener space [233], finite-dimensional distribution matching [234], stochastic optimal control formulations [235], and high-dimensional gradient estimators [236], all contributing to enhanced scalability and optimization efficiency. Thirdly, Neural SDEs have been applied across various fields such as robust financial modeling [237], spatiotemporal dynamics [238], and unified sequence modeling [239]. Finally, recent studies have improved the modeling of rare and irregular dynamics: Neural Jump SDEs [63] address discontinuities through Poisson jumps, while Stable SDEs [66] introduce variants like Langevin-type, linear-noise, and geometric models to enhance the analysis of irregular time series. Collectively, these developments position Neural SDEs as a versatile and robust framework for modeling complex stochastic systems in continuous time. A detailed overview is presented in 3 for Neural SDEs and their variants.

The evolution of Neural SDEs embodies a critical trade-off between uncertainty quantification and optimization tractability. Theoretically, Neural SDEs succeed exactly where deterministic ODEs fail: their learnable diffusion terms break rigid topological constraints to naturally capture aleatoric uncertainty and robustly model highly irregular, noisy data. However, this stochastic expressivity introduces severe computational failures. The inherent non-differentiability of Brownian motion forces stochastic integrators to use vanishingly small steps, while computing gradients via the stochastic adjoint method frequently yields catastrophic variance, destabilizing optimization. Recent advancements succeed by mathematically circumventing these numerical bottlenecks. Techniques such as signature kernels, cubature methods, and deterministic approximations shift the objective from tracking high-variance individual paths to matching robust topological signatures or finite-dimensional distributions. Ultimately, deploying Neural SDEs requires strictly weighing the necessity of modeling intrinsic randomness against the steep costs of numerical integration and gradient stabilization.

3.2.2 SDE-Guided Generative Models [112], [240]↩︎

SDE-guided generative models generate samples by first adding noise to data in a forward process, then reconstructing data from noise in a reverse denoising process. The forward process can be formulated as a SDE: \[\label{SDE2} dx = f(x,t)\,dt+g(x,t)\,d\mathcal{W}(t),\tag{9}\] where \(f(x,t)=-\frac{1}{2}\beta(t)x\) and \(g(t)=\sqrt{\beta(t)}\) in the standard DDPMs setting. In DDPMs, the generative process is realized by solving the reverse SDE, which simulates the dynamic evolution from the noise distribution to the target data distribution. The evolution of the system’s probability density is described by the FPE, which accurately characterizes the statistical properties of the underlying stochastic dynamics. The corresponding reverse SDE can be expressed as: \[\label{RSDE1} \begin{align} dx \!=\! \Big[ f(x,t)\!-\!g^2(t)\nabla_x \log p(x,t)\Big] dt\!+\! g(t) d\overline{\mathcal{W}}(t), \end{align}\tag{10}\] where the score function \(\nabla_x \log p(x,t)\) is typically approximated by a parameterized neural network \(s_\theta(x,t)\) and trained using score matching methods. Therefore, the generative process of DDPMs essentially involves the discrete numerical solution of the aforementioned reverse SDE.

Guided by the DDPMs paradigm, the development of SDE-guided generative models can be systematically divided into three major directions: the design of forward SDEs, the design of reverse SDEs, and the development of solvers for reverse SDEs.

a) Forward SDE: In diffusion models, the forward process describes the transformation of the data distribution into a prior noise distribution. This process is typically represented by a SDE, as shown in Equation 9 . By flexibly selecting specific forms for the drift term \(f(x, t)\) and the diffusion term \(g(t)\), or by extending the diffusion term to a state-dependent function \(g(x, t)\) to include multiplicative noise, various diffusion model variants can be developed. Several representative works will be introduced in the following section.

Mean-Reverting SDE (MRSDE) [59] extends the diffusion framework to model processes with a natural tendency to revert to a long-term mean, making it suitable for mean-reverting scenarios. The process is described by: \[dx(t)=\theta(t)(\mu-x)dt+g(t)d\mathcal{W}(t),\] which captures mean-reverting trajectories and stabilizes long-term generative dynamics. Beyond classical additive noise SDEs, various forward processes inject Gaussian noise using distinct schedules and parameterizations. Notable examples include Noise Conditional Score Networks (NCSNs) [111], Multi-Speed Diffusion (MSDiff) [58], Bayesian Flow Networks (BFNs) [241], and Ornstein-Uhlenbeck diffusion [60].

When the diffusion term \(g(x, t)\) depends on \(x\), the additive SDE framework extends to multiplicative noise SDEs, enabling the modeling of state-dependent stochasticity. This extension addresses more complex phenomena such as variable uncertainty or structural constraints in data [114], [115].

Flexible Diffusion Processes (FP-Diffusion) [115] enhances model expressiveness by simultaneously learning the drift \(f(x)\) and data-dependent diffusion \(R(x)\) through a multiplicative noise SDE, surpassing traditional additive noise methods: \[dx = f(x)\beta'(t)dt+\sqrt{2\beta'(t)R(x)}d\mathcal{W}(t),\] which allows the model to adapt flexibly to heterogeneous data distributions and to capture more intricate generative dynamics by learning both the drift and diffusion structures.

Recent works have generalized the framework by considering non-Gaussian noise [242]–[246], which broadens the class of data distributions that can be modeled. Further extensions introduce physical constraints to the diffusion process, for example, by enforcing boundary conditions or reflecting barriers [113], [247]–[250].

b) Reverse SDE: While the reverse SDE 10 provides a principled trajectory, its strict dependence on the forward process restricts modeling flexibility. Consequently, recent studies mathematically redesign the reverse SDE by learning more expressive drift functions or jointly optimizing the drift and diffusion terms. Unified Framework [61] introduces a tunable parameter \(\gamma \in [0, 1]\) into the reverse SDE, allowing smooth interpolation between deterministic ODEs and stochastic SDEs: \[dx=[f(x,t)-\gamma\cdot\frac{1}{2} g^2(t) s_\theta(x,t)]dt+\sqrt{1-\gamma^2}g(t)d\overline{\mathcal{W}}(t),\] which can flexibly balance between sample diversity (stochasticity) and sample fidelity (determinism), and unifies a wide range of generative modeling dynamics within a single formulation.

This perspective connects reverse-SDE formulations with ODE-based flow models. Flow matching [41] and rectified flow [104] rely on deterministic velocity fields to construct efficient transport dynamics, whereas reverse-SDE formulations retain diffusion terms to model stochasticity and improve distributional coverage. Hybrid ODE/SDE dynamics may thus be understood as a continuum between deterministic transport and stochastic diffusion, in which the drift or velocity term characterizes the transport geometry and the diffusion term modulates the stochastic component. Other recent advances, such as DiffFlow [116], DEIS [117] and ER SDE [118], further explore hybrid drift formulations and independent tuning of drift and diffusion terms, respectively. These extensions provide additional flexibility for balancing sample quality, diversity, and robustness in high-variance regions.

c) Solvers for reverse SDE: Beyond redesigning the reverse dynamics, sampling inherently relies on integrating learned ODEs or SDEs. Consequently, advanced numerical solvers can directly accelerate generation and enhance precision without altering the underlying model architecture. A variety of solvers have been proposed to approximate the solution of the reverse SDEs or its deterministic ODEs counterpart. For example, the reverse process in Probability Flow ODE (PF ODE) [57] can be formulated as the following ODE: \[dx = [f(x, t) - \frac{1}{2} g^2(t) s_\theta(x, t)]dt.\] Several representative solvers have been developed to improve sampling efficiency and quality in practice. Denoising Diffusion Implicit Models (DDIMs) [119] can be interpreted as a deterministic numerical solver for the underlying probability flow ODE, accelerating sampling. Generalized denoising diffusion implicit models (gDDIM) [251] interpolates between deterministic and stochastic sampling. Diffusion Probabilistic Model Solver (DPM-Solver) [252], [253], Pseudo Numerical methods for Diffusion models (PNDMs) [254], Approximate Mean-Direction Solver (AMED-Solver) [255], Stochastic Runge-Kutta Methods (SRK) [256] and Stochastic Adaptive Solver (SA-Solver) [257] apply high-order integration methods to improve sampling efficiency and quality. These solvers enable fast and diverse sample generation with fewer steps. Beyond generation, integrating DDPMs’ stochastic denoising into training acts as a effective regularizer, enhancing stability and generalization while bypassing iterative inference in downstream tasks like image inpainting [258] and segmentation [259].

SDE-guided generative models succeed over deterministic continuous flows by utilizing stochastic noise as a smoothing operator to break topology-preserving constraints, guaranteeing robust coverage of complex data distributions. However, this stochastic success introduces a inference bottleneck: the reverse denoising SDE generates highly curved transport trajectories that cause massive truncation errors during discretization, forcing basic solvers into thousands of costly evaluation steps. To overcome these integration bottlenecks, research progresses along two complementary trajectories. First, advanced models like FP-Diffusion [115] and the Unified Framework [61] redesign the governing dynamics to balance deterministic transport with stochastic exploration, thereby reducing reverse-path curvature while preserving the distributional coverage benefits of diffusion. Second, high-order numerical solvers including DPM-Solver [252] and SA-Solver [257] are deployed to directly combat curvature-induced errors. Ultimately, this reveals a fundamental computational trade-off: optimizing diffusion models requires mitigating reverse-path geometric curvature to maximize sampling efficiency without sacrificing the vital stochastic smoothing benefits of the forward process.

4 Applications and Performance↩︎

This section reviews the experimental performance of a range of DE-driven neural networks on CIFAR-10 and CIFAR-100 benchmarks. While Neural ODEs and Neural SDEs particularly excel in modeling sequential and time-series data, we refer the reader to  [79], [260] for detailed results in that domain. Here, our focus is on architectures and generative models informed by deterministic DEs and SDEs. We structure our discussion as follows: We first review the experimental setup used in the surveyed CIFAR-based studies, and then conduct a qualitative analysis of deterministic DE-driven neural architectures evaluated primarily on CIFAR-10 and CIFAR-100, as summarized in [tab:de-nn], alongside ODE-guided and SDE-guided generative models evaluated mainly on CIFAR-10, as detailed in 6. Section S3 of the online supplementary material complements this focused performance discussion by summarizing representative tasks, application domains, and benchmark datasets for the three model families in Tab. S1, Tab. S2, and Tab. S3. These appendix tables serve as applicability-oriented references rather than additional quantitative comparisons. Through a comprehensive analysis of empirical results across these categories, we clarify the advantages and drawbacks of different models, providing actionable insights for designing equation-inspired architectures with improved training stability and task performance.

4.1 Experiment setup↩︎

For CIFAR experiments, a standard data augmentation (\(4\)-pixel padding followed by \(32\times32\) random cropping or horizontal flipping) is uniformly applied, and pixel values are rescaled to \([0, 1]\). The experimental configurations in the surveyed literature exhibit significant diversity in optimization trajectories.

Standard classification tasks [22], [261], [262] typically utilize a \(300\)-epoch regimen with an initial learning rate (LR) of \(0.1\) and a mini-batch size ranging from \(32\) to \(128\). In contrast, alternative strategies have been explored to enhance convergence or stability: a staged decay over 80k steps is employed in [86], whereas the Adam optimizer with a fixed learning rate of \(0.001\) is used for \(1000\) epochs in [263]. Furthermore, piecewise constant schedules and momentum-based constant LR strategies are adopted in [23] and [264], spanning \(100\) to \(200\) epochs to accommodate specific architectural constraints.

Generative frameworks [250], [265] follow the specialized training protocols established by [111]. These models are typically optimized for \(1.3\)M iterations, with specific SDE variants requiring only \(50000\), using the Adam optimizer with a mini-batch size of \(128\) and a constant learning rate of \(0.001\). To capture multi-scale score functions, a geometric sequence of 10 noise scales \(\{\sigma_i\}_{i=1}^{10}\) is configured, ranging from \(1\) to \(0.01\). During the sampling phase via annealed Langevin dynamics, \(T=100\) steps and a step size of \(\epsilon=2\times10^{-5}\) are adopted. Recent advancements in consistency and flow-matching models further optimize these trajectories. Liu et al. [266] and Ma et al. [267] employ the RAdam optimizer with mini-batch sizes ranging from \(128\) to \(512\) over \(100\)k iterations. Zhou et al. [255] train a lightweight predictor on \(10\)k images using a polynomial time schedule, while Zhang et al. [162] adopt an effective batch size of \(256\) for \(1000\) epochs (approximately \(391\)k iterations) with a polynomial decay schedule and a \(45\)k-step warmup phase.

4.2 Deterministic DE-Driven Network Architectures↩︎

We provide the experimental results of networks guided by first-order ODEs, higher-order ODEs, as well as systems of equations and PDEs in [tab:de-nn]. Consistent with convergence order theory in numerical analysis, networks employing higher-order or more sophisticated schemes generally achieve superior performance.

The experimental results highlight a clear progression, where innovations in DEs modeling translate into advances in both network architectures and empirical performance. ResNets exemplify the explicit Euler discretization of first-order ODEs, where the skip connection directly mirrors the explicit integration step, achieving a strong baseline accuracy of \(94.07\%\) on CIFAR-10 and \(74.84\%\) on CIFAR-100. As variants of ResNet, RiR and ResNeXt further improve accuracy by exploring different formulations. With further advances in mathematical formulation, FractalNet employs a multi-path architecture inspired by Runge-Kutta integration, while RKCNN adopts the fourth-order Runge-Kutta method. These higher-order discretizations enable more accurate modeling of feature dynamics, resulting in notable improvements, such as RKCNN’s \(95.52\%\) on CIFAR-10 and \(77.31\%\) on CIFAR-100. Furthermore, LM-ResNets and LM-ResNeXt employ a linear multi-step strategy to significantly enhance the accuracy of first-order networks under the same configurations.

Extending beyond ODEs, PDEs introduce spatiotemporal priors whose benefits scale with model capacity, as detailed in [tab:de-nn]. Lightweight architectures such as Hamiltonian CNNs [23] utilize sub-million parameter budgets to prioritize stability and interpretability over absolute accuracy. Conversely, at standard capacities, PDE networks yield explicit empirical gains. Specifically, \(11.17\)M-parameter Parabolic and Hyperbolic PDE models [264] achieve up to \(95.70\%\) on CIFAR-10, outperforming the \(94.47\%\) ResNet baseline. Thus, PDE modeling provides mathematical stability for resource-constrained scenarios and superior accuracy for scaled networks.

Overall, DE-driven architectures elevate network design from heuristic discrete mappings to theoretically grounded continuous dynamical systems. This paradigm provides a principled architectural toolkit: first-order ODEs yield computational efficiency, higher-order dynamics capture complex non-linearities, and PDEs enforce rigorous mathematical stability and spatial interactions. Collectively, these continuous formulations significantly enhance network expressiveness and accuracy, offering profound insights for modern model design.

4.3 DE-guided generative models↩︎

6 presents a comprehensive taxonomy of representative generative models, systematically evaluating them based on their governing dynamical frameworks (i.e., ODE- versus SDE-guided paradigms). Over time, advancements in mathematical modeling and solver design have led to consistent improvements in both generative quality and sampling efficiency on CIFAR-10, as reflected by the trends in Fréchet Inception Distance (FID) and Number of Function Evaluations (NFE) metrics.

Benchmark results on CIFAR-10 illustrate a shift in generative modeling from stochastic diffusion to highly efficient deterministic trajectories. Early models like DDPMs require substantial computational budgets exceeding \(1000\) steps and experience severe fidelity degradation under sparse sampling, with FID escalating from \(3.17\) to \(137.7\) at \(20\) steps. This failure primarily stems from the large truncation errors caused by highly curved stochastic paths. Specialized continuous-time formulations and numerical solvers help address this limitation. Score-based SDEs stabilize generation quality, while the EDM framework utilizes refined preconditioning to achieve a state-of-the-art FID of \(1.96\) in only \(35\) steps. Specialized ODE solvers, including DPM-Solver and AMED-Solver, further compress the required budget to under 10 steps while maintaining robust fidelity. Ultimately, methods such as Rectified Flow redefine efficiency limits by straightening transport paths, enabling near-linear velocity fields and high-quality one-step generation.

Despite these algorithmic accelerations, achieving ultra-fast inference introduces a strict trade-off among sample fidelity, inference latency, and training overhead. The transition to single-iteration models represents a redistribution rather than an elimination of computational burden. While classical SDEs maximize sample quality at the cost of high latency, advanced ODE solvers provide a highly practical balance by accelerating inference without requiring supplementary optimization. Conversely, single-step generation relies on computationally expensive offline distillation. Therefore, practical deployment dictates model selection: ODE solvers remain optimal for environments with constrained training resources that prioritize generation quality, whereas distillation approaches are suitable for latency-critical applications capable of absorbing large upfront training costs.

5 Discussions and Conclusion↩︎

5.1 Discussions↩︎

The integration of DEs with deep learning has opened up several promising avenues for advancing the field, but significant opportunities and challenges remain [273], [274]. Building upon the empirical synthesis and fundamental limitations identified in our review, we argue that moving beyond empirical trial-and-error requires a deeper mathematical synthesis. As the relationship between DEs and deep learning evolves, future research should focus on the following key directions to address the limitations of current models.

Overcoming Markovian Limitations. Standard ODEs and CDEs often fail to capture long-range dependencies due to inherent Markovian assumptions. To mathematically embed historical memory and bypass these failures, future architectures should integrate higher-order [275], delay [276], and fractional-order equations [277], [278], which offer theoretically grounded solutions to memory decay.

Breaking Gaussian Assumptions. While SDEs successfully capture aleatoric uncertainty, they fail under extreme discontinuous shifts due to Gaussian assumptions. Exploring Lévy processes [279] naturally accommodates heavy-tailed distributions and sudden jumps, providing rigorous regularization for highly volatile dynamics where standard continuous SDEs fail.

Synthesizing ODE and SDE Paradigms. Generative modeling highlights a trade-off between stochastic coverage and deterministic efficiency. SDE-based diffusion models improve distributional coverage and sample diversity through stochastic perturbations, but typically require iterative reverse-time integration. In contrast, ODE-based approaches, such as flow matching and rectified flow, learn deterministic velocity fields for efficient transport, but provide less explicit stochastic exploration. Recent hybrid formulations address this trade-off at the level of governing dynamics. For example, the Unified Framework [61] interpolates between probability-flow ODEs and reverse SDEs, while DEIS [117] and ER SDE [118] adjust drift and diffusion terms to balance fidelity, diversity, and computational cost. This view also links naturally to flow-matching-based models, where velocity fields characterize transport geometry and diffusion terms modulate stochasticity. Future work may further develop principled ODE/SDE hybrid dynamics for efficient and robust generative modeling.

Aligning with Domain Constraints. The practical failure of DE-driven models often stems from the misalignment neural architectures with domain requirements [280], [281]. Because solving physical PDEs demands exact conservation laws [282] while computer vision prioritizes robustness [283], future research should aim to explicitly embed mathematical properties (e.g., Lyapunov stability) into network designs to ensure theoretically rigorous and safely deployable models.

To sum up, these emerging directions underscore the immense potential of DEs in driving innovation within deep learning, guiding the development of models that are both mathematically principled and highly effective across a wide range of applications.

5.2 Conclusion↩︎

This survey provides a systematic review of the intersection between deep learning and DEs. Subsequently, we establish a unified dual taxonomy that organizes the field by its mathematical formulations and methodological roles. Our analysis demonstrates that model success depends on the synergy between neural architectures and numerical schemes. We clarify the fundamental trade-offs among interpretability, stability, and computational cost to provide principled guidelines for framework selection. Furthermore, we elucidate how SDEs bridge discrete training with continuous dynamics for advanced uncertainty quantification and generative modeling.

Future research should focus on overcoming Markovian limitations by integrating delay and fractional-order equations, while relaxing restrictive Gaussian assumptions through Lévy processes to model more volatile dynamics. Another promising direction is to design hybrid dynamics and align neural architectures with domain-specific constraints, thereby supporting the development of mathematically grounded, reliable, and adaptable intelligent systems. This structured perspective provides researchers with a clearer foundation for future innovation in this evolving field.

References↩︎