A Deep Second-Order Stochastic Residual Method for Fully Nonlinear Parabolic PDEs

We introduce the Deep Second-Order Stochastic Residual Method (D2SRM) for high-dimensional, Hessian-dependent fully nonlinear parabolic PDEs. A single scalar space--time network generates derivative-consistent approximations of the solution, gradient, and Hessian, which are trained jointly through second-order Brownian one-step residuals and terminal value and gradient penalties. For globally Lipschitz equations with identity diffusion and sufficiently weak Hessian coupling, we establish well-posedness in a Brownian occupation space and develop a population-level convergence theory. Under additional regularity, an a posteriori estimate bounds the squared full-jet occupation error of any admissible candidate by the time step and its population objective. For approximate population minimizers, the error bound separates time discretization, neural approximation, and population suboptimality; when the latter two terms are O(h), the full-jet occupation norm is O(h^{1/2}). Experiments on a 100-dimensional manufactured benchmark compare terminal treatments, probe Hessian couplings inside and outside the proved small-gain range, and show decreasing errors as the time step decreases. The code is available at https://github.com/ZZHPKU/D2SRM.