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

Abstract: This paper studies the numerical approximation of divergence-free vector fields by linearized shallow neural networks, also referred to as random feature models or finite neuron spaces. Combining the stable potential lifting for divergence-free fields with the scalar Sobolev integral representation theory via ReLU$<sup>k$ networks, we derive a core integral representation of divergence-free Sobolev vector fields through antisymmetric potentials parameterized by linearized ReLU$<sup>k$ neural networks. This representation, together with a quasi-uniform distribution argument for the inner parameters, yields optimal approximation rates for such linearized ReLU$<sup>k$ neural networks under an exact divergence-free constraint. Numerical experiments in two and three spatial dimensions, including $L<sup>2$ projection and steady Stokes problems, confirm the theoretical rates and illustrate the effectiveness of exactly divergence-free conditions in computation.

• The paper constructs exactly divergence-free finite neuron spaces by combining stable potential liftings with antisymmetric ReLU^k features, eliminating penalties and Lagrange multipliers from the ansatz.
• The resulting spaces achieve the bound O(n^{-(r-s)/d}) in H^s for H^r-regular fields, with a degrees-of-freedom rate that includes a dimension-independent n^{-1/2} term and can outperform Taylor–Hood elements in tested cases.
• Numerical experiments in two and three dimensions confirm predicted convergence for projection and Stokes problems, while cavity-flow results show that corner singularities and conditioning—not divergence enforcement—limit performance.

Overview and motivation

This paper develops a structure-preserving, mesh-free approximation framework for divergence-free vector fields based on linearized shallow neural networks—also called finite neuron spaces (FNS) or random feature models. The central contribution is an integral representation theorem showing that every Sobolev-regular divergence-free field on a starlike Lipschitz domain can be written as a continuous superposition of antisymmetric ReLU$^k$ features, from which optimal approximation rates for the corresponding finite-dimensional trial space follow. The work addresses two recognized gaps: constrained Gaussian process models encode structure but require expensive covariance algebra with cubic cost, while conservation-aware neural solvers typically lack explicit approximation spaces with functional-analytic guarantees.

The practical motivation is well established in the literature: violating $\nabla\cdot u = 0$ at the discrete level produces spurious forces in MHD computations (2603.28638), degrades pressure robustness in incompressible flow solvers, and compromises stability in Maxwell discretizations. Classical remedies—mixed formulations requiring LBB inf–sup compatibility, Raviart–Thomas/Nédélec elements, grad-div penalization, or projection schemes—all depend on mesh topology and specialized element construction, which limits their use in mesh-free or high-dimensional settings.

Construction of the divergence-free finite neuron space

The framework rests on two ingredients. First, the stable potential lifting of Costabel and McIntosh (2603.28638): on a domain starlike with respect to a ball, every $u \in H^r_{\mathrm{div}}(\Omega)$ admits a $(d-2)$-form potential $\mu \in H^{r+1}$ with $u = \mathcal{D}\mu = *d\mu$ and the stability bound $\|\mu\|_{H^{r+1}} \le C_\Omega \|u\|_{H^r}$. In coordinates, $\mathcal{D}\mu$ is the row-wise divergence of a skew-symmetric matrix field; for $d=2$ this recovers the stream function and for $d=3$ the classical vector potential curl.

Second, the scalar Sobolev integral representation theory for ReLU$\nabla\cdot u = 0$0 features (2603.28638): $\nabla\cdot u = 0$1 is the RKHS of ReLU$\nabla\cdot u = 0$2 features and every member admits a representation as an $\nabla\cdot u = 0$3 superposition of $\nabla\cdot u = 0$4, with norm equivalence to the minimal density norm.

Combining these, the paper's main structural result establishes that each component of the potential is itself representable via ReLU$\nabla\cdot u = 0$5 features, so that

$\nabla\cdot u = 0$6

with basis functions $\nabla\cdot u = 0$7, and with the density norm equivalent to $\nabla\cdot u = 0$8 at regularity $\nabla\cdot u = 0$9. Antisymmetry of the potentials is enforced by taking $u \in H^r_{\mathrm{div}}(\Omega)$0. This representation encodes incompressibility exactly at the ansatz level rather than through penalties or Lagrange multipliers, and it serves as the analytic bridge transferring scalar FNS theory to the vector-valued setting.

Optimal approximation rates

The discrete space $u \in H^r_{\mathrm{div}}(\Omega)$1 is defined as the span of $u \in H^r_{\mathrm{div}}(\Omega)$2 over quasi-uniform inner parameters $u \in H^r_{\mathrm{div}}(\Omega)$3, giving $u \in H^r_{\mathrm{div}}(\Omega)$4 outer linear coefficients. Because only outer coefficients are optimized, fitting reduces to convex least squares.

The main convergence theorem states that for $u \in H^r_{\mathrm{div}}(\Omega)$5 with $u \in H^r_{\mathrm{div}}(\Omega)$6 and $u \in H^r_{\mathrm{div}}(\Omega)$7,

$u \in H^r_{\mathrm{div}}(\Omega)$8

The proof approximates each potential component by the scalar FNS in $u \in H^r_{\mathrm{div}}(\Omega)$9 and applies continuity of $(d-2)$0. A notable consequence highlighted by the authors: measured in degrees of freedom, the guaranteed rate always carries a dimension-independent leading term of $(d-2)$1, so convergence is never worse than $(d-2)$2. For example, in $(d-2)$3 with $(d-2)$4, the $(d-2)$5 rate $(d-2)$6 exceeds the Taylor–Hood $(d-2)$7 rate of $(d-2)$8 at matched degrees of freedom.

Two assumptions deserve emphasis. The potential lifting requires $(d-2)$9 to be starlike (hence topologically trivial); on general Lipschitz domains a Hodge decomposition adds a finite-dimensional harmonic field, which the authors note can be handled by adjoining harmonic bases but do not analyze here. Additionally, quasi-uniformity of the inner parameters is essential to the rate and is achieved empirically via Riesz energy minimization on $\mu \in H^{r+1}$0, not proven constructively.

Numerical methodology

All experiments use $\mu \in H^{r+1}$1, double precision, piecewise tensor-product Gauss–Legendre quadrature, and SVD-based least-squares solvers (scipy.linalg.lstsq with driver gelsd). Two solver strategies are compared: normal equations, which yield small $\mu \in H^{r+1}$2 systems but square the condition number, and direct least squares on the tall matrix, which avoids condition squaring at higher memory cost. An ablation at $\mu \in H^{r+1}$3, $\mu \in H^{r+1}$4 shows error saturation of the normal-equation approach at large $\mu \in H^{r+1}$5 while direct least squares continues to track the predicted rate—indicating that beyond moderate width, mass-matrix accuracy is limited by conditioning rather than approximation power. Neurons whose hyperplanes miss the domain are discarded.

For Stokes problems, the paper uses the reduced velocity formulation: restricting test functions to the exactly divergence-free space eliminates the pressure term, and no inf–sup condition arises. Homogeneous Dirichlet data are imposed weakly via a boundary penalty; sensitivity tests across penalty parameters spanning four orders of magnitude show essentially unchanged convergence rates.

Numerical results

Manufactured solutions. For divergence-free $\mu \in H^{r+1}$6 projection with targets built from stream functions ($\mu \in H^{r+1}$7) and curls of vector potentials ($\mu \in H^{r+1}$8), empirical log–log slopes match the theoretical bounds $\mu \in H^{r+1}$9: e.g., in $u = \mathcal{D}\mu = *d\mu$0, observed rates reach approximately 0.76, 1.30, 1.85, and 2.33 for $u = \mathcal{D}\mu = *d\mu$1 against predictions 0.75, 1.25, 1.75, 2.25. Stokes experiments with bubble-function manufactured solutions show analogous agreement in both $u = \mathcal{D}\mu = *d\mu$2 and $u = \mathcal{D}\mu = *d\mu$3 norms across $u = \mathcal{D}\mu = *d\mu$4 in both dimensions.

Lid-driven cavity. On the classical constant-lid benchmark, where corner singularities reduce solution regularity, global convergence degrades substantially—the global $u = \mathcal{D}\mu = *d\mu$5 rate falls below 0.1 while the interior-subdomain rate improves to roughly 0.3–0.44, localizing the difficulty near the top corners. With the regularized lid profile $u = \mathcal{D}\mu = *d\mu$6, high-order trends recover (rates up to about 1.67 in $u = \mathcal{D}\mu = *d\mu$7 for $u = \mathcal{D}\mu = *d\mu$8 and 2.10 for $u = \mathcal{D}\mu = *d\mu$9). At matched degrees of freedom within the tested range, the divergence-free FNS slightly outperforms a Taylor–Hood $\|\mu\|_{H^{r+1}} \le C_\Omega \|u\|_{H^r}$0 dolfinx reference on the smooth benchmark. The authors attribute the classical-cavity degradation to boundary-data regularity rather than to the discretization itself, though they do not provide a rigorous regularity analysis confirming this interpretation.

Limitations and open questions

Several restrictions are stated plainly in the paper. The theory requires domains starlike with respect to a ball; extension to general topologies via harmonic-field augmentation is proposed but unproven. Quasi-uniformity of inner parameters is established only empirically through Riesz energy refinement. The boundary condition is enforced by penalty rather than exactly, and while sensitivity studies show rate robustness, the boundary residual itself is not driven to zero. The comparison with Taylor–Hood elements is limited to one smooth benchmark and a specific DOF range. Finally, the analysis covers linearized (fixed inner parameter) networks; whether comparable guarantees extend to adaptively trained nonlinear dictionaries remains open, as does a theoretical explanation of the corner-singularity-induced rate degradation observed in the classical cavity problem.

Conclusion

The paper provides a complete pipeline—from stable Bogovskiĭ-type potential lifting through scalar ReLU$\|\mu\|_{H^{r+1}} \le C_\Omega \|u\|_{H^r}$1 RKHS theory to a provably optimal, exactly divergence-free linearized network ansatz—with numerical evidence in two and three dimensions confirming the predicted rates. Its main value lies in demonstrating that exact constraint enforcement and sharp functional-analytic guarantees can coexist with the scalability of convex linear least squares, while leaving domain generality, rigorous sampling conditions, and nonlinear training as open directions.