Divergence-free Linearized Neural Networks: Integral Representation and Optimal Approximation Rates

1. Problem Statement

Many physical systems, particularly in fluid dynamics (incompressible flows), magnetohydrodynamics (MHD), and electromagnetism, are governed by the divergence-free constraint ($\nabla \cdot \mathbf{u} = 0$).

• The Challenge: In numerical simulations, violating this constraint at the discretization level can lead to instability, spurious forces, and a loss of physical fidelity.
• Limitations of Existing Methods:
  • Mesh-based methods (e.g., Finite Elements) often require complex saddle-point formulations, specific mesh topologies (like staggered grids), or strict inf-sup compatibility conditions (LBB condition) to enforce incompressibility.
  • Physics-Informed Neural Networks (PINNs) typically enforce constraints via soft penalties in the loss function, which introduces optimization difficulties and does not guarantee exact satisfaction of the constraint.
  • Constrained Gaussian Processes (GPs) offer exact structure preservation but suffer from cubic computational complexity ($O(N^3)$) due to dense covariance algebra, making them unscalable for large datasets.

The paper aims to develop a mesh-free, linearized neural network framework that enforces the divergence-free constraint exactly at the function space level while providing provable optimal approximation rates.

2. Methodology

The proposed method constructs a Divergence-Free Finite Neuron Space (FNS), denoted as $V^k_n$, using linearized shallow neural networks with ReLU$^k$ activation functions.

A. Theoretical Foundation: Integral Representation

The core theoretical contribution is an integral representation theorem for divergence-free Sobolev vector fields.

• Potential Lifting: The authors utilize the fact that on star-shaped domains, any divergence-free field $\mathbf{u} \in H^r_{div}(\Omega)$ can be represented as the generalized curl (or exterior derivative) of an antisymmetric potential tensor $\mu$:
  $$\mathbf{u} = D\mu = \nabla \times \mu \quad (\text{generalized to } d \text{ dimensions})$$
• Scalar-to-Vector Bridge: By combining the stable potential lifting with the known integral representation of scalar Sobolev spaces via ReLU$^k$ features, they derive that the potential $\mu$ can be represented as an integral over the unit sphere $S^d$ involving antisymmetric coefficients.
• Resulting Representation:
  $$\mathbf{u}(x) = \int_{S^d} \sum_{1 \le i < j \le d} \Phi_{\theta, ij}(x) \psi_{ij}(\theta) d\theta$$
  where $\Phi_{\theta, ij}(x) = \partial_j \sigma_k(\theta \cdot \tilde{x})\mathbf{e}_i - \partial_i \sigma_k(\theta \cdot \tilde{x})\mathbf{e}_j$. This ensures $\nabla \cdot \mathbf{u} = 0$ by construction.

B. Discrete Approximation Space ($V^k_n$)

• Construction: The continuous integral is discretized by fixing the inner parameters (neuron weights and biases) $\{\theta_\ell\}_{\ell=1}^n$ to be quasi-uniformly distributed on the sphere $S^d$.
• Basis Functions: The space is spanned by vector-valued basis functions derived from the antisymmetric potential:
  $$\Phi_{\ell, ij}(x) = \partial_j \phi_\ell(x)\mathbf{e}_i - \partial_i \phi_\ell(x)\mathbf{e}_j$$
  where $\phi_\ell(x) = \sigma_k(\theta_\ell \cdot \tilde{x})$.
• Linearity: The outer coefficients are learned via convex linear least squares. This avoids the non-convex optimization issues of standard deep learning training.

C. Numerical Implementation

• Quadrature: Deterministic Gauss-Legendre quadrature is used to approximate inner products and boundary integrals.
• Solvers: The resulting linear systems are solved using either Normal Equations or Direct Least Squares (SVD/QR). The paper notes that for high dimensions or large $n$, direct least squares is more numerically stable as it avoids squaring the condition number.
• Boundary Conditions: Homogeneous Dirichlet conditions are enforced weakly via a penalty term, while the divergence-free constraint is hard-coded into the trial space.

3. Key Contributions

1. Exact Divergence-Free Ansatz: The method constructs a finite-dimensional function space where every linear combination of basis functions is strictly divergence-free, eliminating the need for pressure variables or penalty terms for incompressibility.
2. Integral Representation Theorem: The paper establishes a rigorous integral representation for $H^r_{div}(\Omega)$ using antisymmetric potentials parameterized by ReLU$^k$ networks. This bridges scalar approximation theory with vector-valued divergence-free fields.
3. Optimal Approximation Rates: The authors prove that under quasi-uniform sampling of inner parameters, the method achieves optimal convergence rates consistent with the underlying scalar theory:
   $$\inf_{\mathbf{u}_n \in V^k_n} \|\mathbf{u} - \mathbf{u}_n\|_{H^s} \lesssim n^{-(r-s)/d} \|\mathbf{u}\|_{H^r}$$
   This confirms that the "curse of dimensionality" is handled similarly to scalar problems, and the divergence-free constraint does not degrade the convergence order.
4. Scalability: Unlike constrained GPs, this approach replaces infinite-dimensional covariance algebra with scalable, convex linear algebra in a fixed feature space.

4. Experimental Results

The authors validate the theory through numerical experiments in 2D and 3D:

• L2 Approximation of Divergence-Free Fields:
  • Tested on synthetic vector fields with varying smoothness and activation powers ($k=1,2,3,4$).
  • Result: Observed error decay rates matched the theoretical predictions ($O(n^{-(r-s)/d})$).
• Steady Stokes Equations (Homogeneous BCs):
  • Solved the Stokes problem using the reduced formulation (no pressure variable).
  • Result: The method achieved high-order convergence in both $L^2$ and $H^1$ seminorms, outperforming standard finite element baselines (Taylor-Hood P2/P1) in terms of error per degree of freedom for specific activation orders.
• Lid-Driven Cavity Flow (Non-Manufactured Solution):
  • Tested on a benchmark with non-homogeneous boundary conditions.
  • Findings:
    • For the classical case (discontinuous lid velocity), convergence was limited by the low regularity of the solution near the corners, not the method itself.
    • For the regularized case (smooth lid velocity), the method recovered high-order convergence, demonstrating its effectiveness on realistic problems.
• Numerical Stability: Experiments comparing Normal Equations vs. Direct Least Squares showed that for large networks, direct solvers are necessary to avoid accuracy loss due to ill-conditioning.

5. Significance

This work represents a significant step forward in Scientific Machine Learning (SciML) for fluid dynamics and vector field problems:

• Structure Preservation: It provides a mathematically rigorous way to embed physical conservation laws (divergence-free) directly into the neural network architecture, ensuring stability and physical consistency without soft penalties.
• Theoretical Guarantees: It moves beyond empirical success by providing provable optimal approximation rates, a rarity in neural network-based PDE solvers.
• Efficiency: By linearizing the problem (fixing inner parameters), it offers a computationally efficient alternative to both traditional mesh-based methods (which struggle with complex geometries) and constrained Gaussian Processes (which are computationally expensive).
• Generalizability: The framework is dimension-independent ($d \ge 2$) and applicable to various vector field problems beyond fluid dynamics, such as Maxwell's equations.

In summary, the paper successfully unifies differential geometry (exterior calculus), approximation theory (Sobolev integral representations), and machine learning (linearized neural networks) to create a powerful, exact, and theoretically grounded solver for divergence-free vector fields.