Aggregation-diffusion equations skip the mesh: a neural velocity field keeps mass and positivity intact

The aggregation-diffusion equation is rewritten as an optimal control problem whose control is the transport velocity and whose constraint is the continuity equation, with density as the state, so mass conservation and density positivity follow from the constraint rather than from post-hoc correction. Most earlier structure-preserving numerical methods relied on mesh-based spatial discretization, whose computational cost and memory requirements grow rapidly with spatial dimension, making high-dimensional problems hard to handle. The equivalence is stated as a theorem: the constrained problem attains global minimum value zero, and its minimizers correspond to solutions of the aggregation-diffusion equation, with density and optimal velocity uniquely determined in the relevant function spaces.

The velocity field is parameterized by a neural network and the variational objective is evaluated along particle trajectories, with density and generalized score equations coupled into one system advanced by the Euler method, so the high-dimensional density equation need not be solved explicitly. Neural solvers that approximate the density directly generally do not automatically preserve positivity, mass conservation and free-energy dissipation, and small violations can accumulate over long-time simulations. Numerical experiments cover the one-dimensional Doi-Onsager equation, the two-dimensional Patlak-Keller-Segel system, the porous medium equation, and three- and ten-dimensional Ornstein-Uhlenbeck problems; in the three- and ten-dimensional cases the time per iteration stays between 1.2 and 1.5 seconds and peak memory grows approximately linearly in dimension and particle number.

Total error is decomposed into neural approximation, optimization, Monte Carlo sampling and time-discretization components with high-probability bounds, separating the contributions of network capacity, optimization residual and sampling error. Comparable neural approaches did not provide such a unified error decomposition, and convergence discussions usually covered only some of these sources. The bounds hold under strictly positive, sufficiently smooth initial density and bounded, Lipschitz interaction forces; for the Newtonian potential singularity of the two-dimensional Patlak-Keller-Segel model the theory covers a mollified kernel, while the numerics use a blob-regularized kernel with self-interaction set to zero.

The framework targets aggregation-diffusion equations with free-energy dissipation structure, including chemotaxis, biological aggregation and rod-like polymer dynamics; numerical validation covers the one-dimensional Doi-Onsager equation, the two-dimensional Patlak-Keller-Segel system, the porous medium equation, and three- and ten-dimensional Ornstein-Uhlenbeck problems. The error analysis holds under strictly positive, sufficiently smooth initial density and bounded, Lipschitz interaction forces, and the Newtonian potential singularity of the two-dimensional Patlak-Keller-Segel model is handled through a mollified kernel. The authors propose extensions to coupled PDE systems, more general mobility functions, anisotropic transport, non-Euclidean variational structures, and adaptive strategies for efficient training over long-time dynamics.

The error bound is a conditional decomposition rather than an unconditional convergence guarantee, and converting it into an asymptotic convergence result requires explicit hypotheses on neural network expressivity and optimization performance. In the porous-medium case the KL-based stability argument does not directly apply, and a separate analysis with an appropriate discrepancy measure would be needed. In the high-dimensional experiments the relative residual on training particles stays at a few percent, while the error on fresh particles is up to five times larger at ten dimensions, indicating that finite-sample generalization rather than optimization or time discretization dominates. Direct pairwise evaluation of the nonlocal interaction costs quadratically in particle number per time step, manageable at the particle counts used here but requiring fast summation methods at larger scale.