Deep-Picard Iteration for Space-time Fractional Diffusion PDEs

1. Problem Statement

The paper addresses the numerical solution of high-dimensional nonlinear space-time fractional partial differential equations (PDEs). The specific problem is the Dirichlet problem for the equation:
$$\begin{cases} \partial_t^\beta u(t, x) + (-\Delta)^{\alpha/2} u(t, x) = f(t, u(t, x), x), & (t, x) \in (0, T] \times \Omega \\ u(t, x) = g(t, x), & (t, x) \in [0, T] \times \Omega^c \\ u(0, x) = u_0(x), & x \in \Omega \end{cases}$$
where:

• $\Omega \subset \mathbb{R}^d$ is a bounded domain.
• $\partial_t^\beta$ is the Caputo fractional derivative ($\beta \in (0, 1)$), representing time-fractional memory and trapping effects.
• $(-\Delta)^{\alpha/2}$ is the integral fractional Laplacian ($\alpha \in (0, 2)$), representing nonlocal spatial transport driven by Lévy flights.
• $f$ is a nonlinear reaction term.

Challenges:

1. Curse of Dimensionality: Classical grid-based methods (FEM, FDM) suffer from exponential cost growth as dimension $d$ increases.
2. Nonlocality & Memory: The fractional Laplacian couples every spatial point to the entire exterior domain, creating dense algebraic systems. The Caputo derivative introduces history dependence, requiring substantial storage.
3. Nonlinearity: Standard linear Feynman-Kac representations do not apply directly to nonlinear reaction terms, preventing straightforward probabilistic solutions.

2. Methodology: Deep-Picard Iteration (DPI)

The authors propose a framework that combines probabilistic representations with deep learning to avoid direct discretization of fractional operators.

A. Theoretical Foundation: Nonlinear Feynman-Kac Fixed-Point

The core innovation is recasting the PDE as a stochastic fixed-point equation.

• Stochastic Processes: The solution is linked to a coupled process involving:
  • A $\beta$-stable subordinator ($Y_s^\beta$) to handle the time-fractional derivative (inverse stable time change).
  • A rotationally symmetric $\alpha$-stable Lévy process ($X_s^x$) to handle the fractional Laplacian.
• Fixed-Point Formulation: The solution $u$ satisfies:
  $$u(t, x) = \mathbb{E}_x \left[ \mathcal{S}_{t,x} + \int_0^{\tau(t,x)} f(t_s, u(t_s, X_s^x), X_s^x) \, ds \right]$$
  where $\mathcal{S}_{t,x}$ is a terminal payoff (initial or boundary data) and $\tau$ is a stopping time determined by the process exiting the domain or the stochastic clock reaching zero.
• Contraction: The authors prove that the associated nonlinear operator $\mathcal{T}$ is a strict contraction in a weighted Banach space, guaranteeing the existence and uniqueness of the fixed point.

B. Algorithmic Implementation

The method approximates the fixed point via Picard iteration ($u^{(k+1)} = \mathcal{T}u^{(k)}$), where each iteration is solved using a neural network.

1. Trajectory Simulation (Mesh-Free):

   • Time Component: The $\beta$-stable subordinator is discretized to simulate the backward physical time $t_s = t - Y_s^\beta$.
   • Space Component: The $\alpha$-stable Lévy process is simulated using a Walk-on-Spheres (WoS) algorithm adapted for fractional Laplacians. This allows sampling exit locations from balls without a spatial mesh.
   • Stopping: Trajectories stop when the spatial process exits $\Omega$ or the time clock hits zero.

2. Neural Regression:

   • At each Picard step $k$, training data is generated by sampling points $(t_i, x_i)$ and simulating $M$ independent trajectories.
   • Labels: For each trajectory, the integral term is approximated (using trapezoidal rule) using the previous iterate $u^{(k)}$ (approximated by a neural network $u_{\theta_k}$). The Monte Carlo average of these trajectories forms the regression target $\hat{y}^{(k)}$.
   • Training: A new neural network $u_{\theta_{k+1}}$ is trained via supervised learning to minimize the mean squared error between its output and the Monte Carlo labels.

3. Relaxed Iteration:

   • To stabilize convergence against Monte Carlo noise and nonlinearity, a relaxed update is used:
     $$u^{(k+1)} = (1-\eta)u^{(k)} + \eta \mathcal{T}u^{(k)}$$
     where $\eta \in (0, 1]$ is a relaxation parameter.

3. Key Contributions

• Novel Framework: First application of Deep Picard Iteration to nonlinear space-time fractional PDEs, successfully handling both Caputo time derivatives and fractional Laplacians simultaneously.
• Avoidance of Direct Discretization: The method completely bypasses the dense linear systems and memory overhead associated with discretizing fractional operators, replacing them with stochastic simulations.
• High-Dimensional Scalability: The trajectory-based approach is naturally parallelizable and independent of tensor-product grids, making it effective in high dimensions ($d \le 100$).
• Theoretical Justification: Provides a rigorous proof of the contraction property of the nonlinear Feynman-Kac operator in a weighted space, ensuring the convergence of the Picard scheme.
• Mesh-Free Boundary Handling: Utilizes the Walk-on-Spheres method to naturally handle exterior Dirichlet conditions without complex boundary meshing.

4. Numerical Results

The authors tested the method on various benchmarks with $\alpha=1.5$ and $\beta=0.6$:

• 2D Unit Disk:
  • Demonstrated convergence with respect to operational time step ($\Delta s$) and Monte Carlo paths ($M$).
  • Showed that relaxation ($\eta < 1$) significantly reduces error (by ~30%) and stabilizes the iteration, especially for cubic nonlinearities where label noise is amplified.
  • Achieved high accuracy with errors in the range of $10^{-3}$ to $10^{-2}$.
• 2D Unit Square (Non-smooth Boundary):
  • Successfully handled non-smooth boundaries and profiles without closed-form fractional Laplacians (using precomputed FFT-based forcing).
  • Errors concentrated at corners, consistent with boundary regularity issues.
• Spatial Heterogeneity (Double-Bump):
  • The method successfully resolved complex, multi-modal spatial profiles, capturing localized features and saddle regions.
• High Dimensions ($d = 20, 50, 100$):
  • The method remained stable and accurate up to $d=100$.
  • Errors did not deteriorate with increasing dimension, confirming the method's ability to overcome the curse of dimensionality.
  • The "concentration of measure" in high-dimensional balls (where most points are near the boundary) was noted as a factor aiding stability.

5. Significance and Future Work

Significance:
This work bridges the gap between probabilistic numerical methods and deep learning for a class of PDEs that are notoriously difficult to solve due to nonlocality, memory effects, and high dimensionality. It offers a mesh-free, scalable, and naturally parallelizable alternative to traditional deterministic solvers for fractional PDEs.

Future Directions:

• Extension to more complex geometries and general boundary conditions.
• Implementation of advanced fixed-point acceleration techniques (e.g., Anderson acceleration, adaptive relaxation).
• Variance reduction strategies (Multilevel Monte Carlo, control variates) to improve efficiency.
• Rigorous discretization error analysis quantifying the interplay between stochastic simulation, quadrature, and neural approximation.

In summary, the paper presents a robust, theoretically grounded, and computationally efficient framework for solving high-dimensional nonlinear fractional PDEs, demonstrating that deep learning combined with stochastic fixed-point iterations can effectively tackle problems previously considered intractable in high dimensions.