1 Introduction↩︎

Understanding the surprising empirical effectiveness of gradient-based optimization on highly overparameterized, nonconvex objectives remains one of the central challenges in deep learning theory. Explanation of non-overfitting in such problems is in a short supply. A widely studied approach focuses on representing the network output as an integral over a distribution of hidden units with \[\label{eq-def-fN-intro} f_N(x):=\tfrac1N\textstyle\sum_{i=1}^N w_i h(\theta_i,x)\,, \quad \text{ where \((w_i,\theta_i)_{i=1}^N\) stand for the unknown parameters,}\tag{1}\] and studies the minimization of a risk functional over the space of measures, with gradient descent interpreted as a flow on that space. This mean-field perspective has revealed deep connections to interacting particle systems, optimal transport, and Wasserstein gradient flows [1]–[7], and has led to convergence guarantees, insights into feature learning, and PDE-based analyses of training dynamics. In parallel, the neural tangent kernel (NTK) framework [8] and related overparameterization results [9], [10] explain why gradient descent succeeds in a linearized infinite-width regime. More recent efforts extend mean-field PDE analyses to deeper, residual architectures [11], sharpen particle and Langevin approximations [12], [13].[14] study existence of minimizers in the finite-width landscape via a closure of the search space; [15] establish nearly optimal approximation rates for shallow \(\mathrm{ReLU}^k\) networks on Sobolev classes. Both works operate at finite networks, whereas we concentrate on the regularity of the continuum minimizer. Finally, [16] add a Sobolev penalty on the network output over the input domain \(D\) to improve optimization conditioning. In contrast, our regularization acts on the parameter density over \(\Omega\) to enforce variational well-posedness and regularity. This makes the two approaches complementary yet structurally distinct.

Despite this progress, two principal frameworks dominate the field, each with characteristic trade-offs. Mean-field/Wasserstein approaches support a particle interpretation of gradient descent but offer only displacement convexity and require non-trivial approximation to pass to finite-width models. The NTK framework achieves a truly convex objective but only by linearizing around initialization, limiting guarantees to the lazy-training regime. Both leave partly open a different set of questions about the variational problem itself: Does the infinite-width objective admit a unique, stable minimizer? What regularity do optimal parameter distributions possess? Can one obtain global convexity for a nonlinear continuum model without linearization?

In this paper we propose a third route that addresses these questions directly. We formulate shallow neural network training as a variational problem over parameter densities in a weighted Sobolev space \(\mathcal{W} = W^{1,2}(\Omega)\cap L^2_\omega(\Omega)\) forming a linear Hilbert-space setting that retains the infinite-width mean-field spirit while gaining direct access to convex analysis, elliptic PDE theory, and quantitative gradient flow estimates. The regularized objective on \({\mathcal{W}}\) which is given by \[{\mathscr{F}}^{(f)}_{\alpha,\beta}(u) = \mathcal{R}(f, u\,\mathrm{d}\theta) + \alpha\|u\|_{L^2_\omega}^2 + \beta\|\nabla u\|_{L^2}^2\] combines a squared-risk term measuring data fidelity with \(L^2_\omega\) and Sobolev regularization. The \(L^2_\omega\) term induces coercivity and \(\lambda\)-convexity for \(\lambda>0\); the Sobolev term promotes smoothness and plays a role analogous to entropy regularization in the well-established, yet computationally hard, Wasserstein formulations [17]–[19]. In turn, the minimizer \(u^*_f\) to \({\mathscr{F}}^{(f)}_{\alpha,\beta}\) can be interpreted as the implicit target of NN in the overparameterized regime, while the regularization parameters \(\alpha,\beta\) play the role of weight decay and smoothness bias.

Our Hilbert-space framework is built on the classical calculus of variations staying fully compatible with standard shallow network parametrizations. Although our analysis does not directly model stochastic gradient descent, the exponential convergence rate \(e^{-\lambda t}\) of the gradient flow associated to our variational problem and the structure of the Euler–Lagrange equation suggest that the \(L^2_\omega\) geometry may serve as a useful proxy for understanding the implicit bias in overparameterized networks. In particular, the \(\lambda\)-convex structure implies that any discretization of the gradient flow inherits stability consistent with the empirical robustness of SGD near overparameterized minima [1], [2]. In our study, the minimization of \({\mathscr{F}}_{\alpha,\beta}^{(f)}\) can be reduced to the optimization of a quadratic functional depending on the coordinates of the competitors in a chosen basis of \(L^2_\omega(\Omega)\). Since the projection of the minimizer onto a finite dimensional subspace of \(L^2_\omega(\Omega)\) can be obtained by merely solving a linear system of equations, the optimization is a fast and easy procedure, as it does not require any iteration or gradient descent arguments.

The variational formulation we propose yields a cluster of remarkable properties.

No Lavrentiev gap. A key consistency result shows that the infimum of the risk is the same whether one optimizes over atomic measures (finite-width networks), Sobolev densities, or smooth compactly supported functions. No approximation error is introduced by passing to the continuum, and finite-width networks of width \(N\) achieve the continuum optimum up to an \(O(1/N)\) error. This places the continuum theory as an exact relaxation of the finite discrete training problem. The absence of Lavrentiev’s phenomenon is a typical key step in inferring regularity of minimizers [20]–[23].

\(\lambda\)-convexity. The functional \({\mathscr{F}}_{\alpha,\beta}^{(f)}\) is \(\lambda\)-convex on \(L^2_\omega\) with \(\lambda = 2\alpha\geq 0\), which gives the exponential convergence of the gradient flow, and it is \(2\min(\alpha,\beta)\)-convex on \({\mathcal{W}}\), which leads to the existence and the uniqueness of the minimizer, as well as its stability with respect to the data. These are global properties of the model, not consequences of linearization, equipping with strong tools, which are typically not available in the Wasserstein approach, [24]. The \(L^2_\omega\)-gradient flow of \(\mathcal{F}^{(f)}_{\alpha,\beta}\), the continuous-time analogue of gradient descent, converges exponentially fast to \(u^*\), cf. [25]. A time-discrete implicit Euler scheme provides an approximation with explicit error, combining the exponential rate with a \(\sqrt{\tau}\) term, where \(\tau\) is the discretization step.

Near-\(C^3\) smoothness. The minimizer \(u^*\) of \({\mathscr{F}}_{\alpha,\beta}^{(f)}\) solves an explicit Euler–Lagrange equation. Moreover, it is \(C^{2,s}\) for every \(s<1\), which gives a major structural advantage: it prevents pathological concentrations of mass in parameter space and opens the door to higher-order numerical schemes via classical PDE methods. In contrast, in the Wasserstein approach, regularity of minimizers is inherited from the activation function, whereas here we get surprising smoothness even in the case of merely locally Lipschitz activation functions, see [1]. This suggests that the predictions of the neural network concentrate around graphs of almost \(C^3\)-functions, which illustrates the non-overfitting phenomenon, robustness of training, and suggests that the neural networks are attracted to the lower-dimensional manifolds, cf. [26]–[28]. The stability results we infer provide a concrete form of implicit bias toward smooth solutions and explain why the variational formulation yields stable and well-generalizing models.

Contributions. Our framework provides an entirely variational formulation for infinite-width shallow networks. Replacing the discrete training process with a regularized surrogate, we identify a class of globally well-posed and stable objectives. This approach follows three steps: (i) transforming the \(N\)-neuron optimization into a continuum problem over weighted Sobolev spaces; (ii) exploiting the structure to ensure a unique and almost \(C^3\) regular minimizer; (iii) bridging to finite-width networks via atomic measures with \(O(1/N)\) consistency. This offers novel insight into the implicit bias of neural networks, as detailed below.

• A variational theory of shallow network training. We formulate shallow neural network training as a Sobolev-space variational problem (Section 2) and concentrate on the properties of solutions. We prove that this continuum formulation is exact: there is no Lavrentiev gap, so atomic, Sobolev, and smooth minimization classes coincide (Theorem 1). We show that this infinite-width problem is quantitatively faithful to finite networks via an \(O(1/N)\) approximation bound (Proposition 2).

• Global convexity and well-posedness beyond existing frameworks. We introduce a regularized functional \(\mathcal{F}_{\alpha,\beta}^{(f)}\) that is globally \(2\min(\alpha, \beta)\)-convex on \({\mathcal{W}}\) yielding existence, uniqueness, and stability of minimizers. We establish that the associated \(L^2_\omega\) gradient flow converges exponentially fast to equilibrium at rate \(e^{-2\alpha t}\) (Corollary 3). In particular, this provides an explicit analytical mechanism for non-overfitting in the infinite-width regime.

• Implications to ML. The minimizers to the regularized functional admit an Euler–Lagrange characterization and exhibit near-\(C^3\) regularity via elliptic PDE arguments (Theorem 3), a level of smoothness not accessible in standard infinite-width analyses. Propositions 4 and 5, together with Corollary 2, show quantitative dependence of the minimizer under perturbations providing a link between regularization, stability, and generalization.

Shallowness. Despite the dominance of deep architectures, the single-hidden-layer setting already captures the core variational complexity of neural network training: a nonlinear, non-convex objective over an infinite-dimensional space. Resolving these obstructions in this setting isolates the essential structure and provides a foundation for multilayer extensions [11].

Limitations. The current formulation is for shallow (one-hidden-layer) networks. Moving from one to more layers complicates the structure of the optimization problem, which becomes strongly nonlinear, see Section 3. In fact, for two hidden layers the parameter density lives on \(\Omega_1 \times \Omega_2\) and the risk functional becomes a composition, making \({\mathscr{R}}\) possibly nonconvex even after Sobolev regularization. The Euler–Lagrange system becomes a coupled nonlinear PDE rather than a linear equation, so existence of smooth minimizers requires additional structural assumptions such as small initial condition or a perturbative regime around a known solution. The gradient flow no longer reduces to a linear semigroup and convergence is elusive.

Summary. Our approach is related to the mean-field literature [1]–[3], [12], [13], but operates in a different analytical regime. The connection to Sobolev-type objectives recently explored in [16] is suggestive, though the focus there is on conditioning of the optimization problem rather than on variational well-posedness or regularity of minimizers. Our results also complement convex and function-space formulations of infinite-width networks [29], [30], which concentrate on description of the learning process; we focus instead on the variational structure of the objective and the regularity of its minimizers. Finally, while the NTK framework [8] achieves convexity through linearization and is well-suited to the lazy-training regime [9], [10], our \(\lambda\)-convexity holds for the same nonlinear problem that arises in the infinite-width limit of shallow network training, without any linearization assumption.

2 Variational formulation↩︎

We develop a continuum approach to shallow networks based on parameter distributions.

2.1 Exact variational characterization of neural network training↩︎

Let \(D\subset {\mathbb{R}}^d\) be the input (a bounded Borel set equipped with the Lebesgue measure), let \(f\in L^2(D)\) and let \(\Omega\subset {\mathbb{R}}^{d+1}\) be the parameter space of hidden units, assumed to be open. Let \(\theta=(\theta_0,\theta')\in\Omega\) parametrize the hidden layer with \(\theta_0\in {\mathbb{R}}\) being interpreted as the bias, \(\sigma\) be an activation function that we assume to be continuous with the typical choices of ReLU or sigmoid functions. We also consider the function \(h: \Omega\times D\to {\mathbb{R}}\) given by \[\label{eq:h-def}h(\theta,x):=\sigma(\theta_0+\theta'\cdot x)\,\tag{2}\] and we require that there exists \(C_h>0\) such that for every \((\theta,x)\in \Omega\times D\), \[\label{eq:h-ass} |h(\theta,x)|\leq C_h(1+|x|)(1+|\theta|)\,.\tag{3}\] This condition holds true in the typical choices of ReLU and sigmoid functions. In this framework training a shallow neural network consists in approximating \(f\) with functions \(f_N\) that are given by 1 .

The learning process amounts to the minimization of \(\int_D (f(x)-f_N(x))^2{\,{\rm d}x}\) over \((w_i,\theta_i)_{i=1}^N\). Each function \(f_N\) can be rewritten as \[\label{eq:fN-muN} f_N(x)=\int_\Omega h(\theta,x){\,{\rm d}\mu}_N(\theta)\quad\text{for }\;\mu_N = \frac{1}{N} \sum_{i=1}^N w_i \delta_{\theta_i}\in{{\mathscr{M}}^{\rm at}}(\Omega)\,.\tag{4}\] Observe that \(\mu_N\) is a finite atomic signed measure (with finite number of atoms). Denote by \({\mathscr{M}}_a(\Omega)\) the set of those finite signed Borel measures \(m\) such that \(\int_{\Omega}|\theta|^a{\,{\rm d}}|m|(\theta)<\infty.\) Then, we generalize 4 for \(m\in{\mathscr{M}}_1(\Omega)\) by introducing the function \[\label{eq:fm} f_m(x)=\int_\Omega h(\theta,x){\,{\rm d}}m(\theta)\,.\tag{5}\] In view of 3 and the fact that \(m\in {\mathscr{M}}_1(\Omega)\), the function \(f_m\) is well-defined and belongs to \(L^{2}(D)\). We can thus consider the population risk functional \[\label{eq:risk} {\mathscr{R}}(f,m):=\int_D (f(x)-f_m(x))^2{\,{\rm d}x}\,.\tag{6}\] The training problem becomes the minimization of \({\mathscr{R}}(f,m)\) over \(m\in{\mathscr{M}}_1(\Omega)\), as suggested in [3]. To obtain a Hilbert-space formulation, we restrict attention to measures that admit a density, writing \({\,{\rm d}}m(\theta)=u(\theta){\,{\rm d}\theta}\) for some Sobolev function \(u\). Since \(\Omega\) is typically unbounded, we work in a weighted space that controls integrability, so that all quantities under consideration are well-posed. Let \(\omega:\Omega\to[1,\infty)\) be a smooth weight such that \[\label{eq:weight-assumptions} c_\omega:=\textstyle\int_\Omega {(1+|\theta|)^4}\tfrac{1}{\omega(\theta)}{\,{\rm d}\theta}<\infty\,.\tag{7}\] Moreover, we define \[L^2_\omega(\Omega) = \left\{ u:\textstyle\int_\Omega u(\theta)^2\omega(\theta){\,{\rm d}\theta}<\infty \right\} \qquad\text{and}\qquad {\mathcal{W}}:= W^{1,2}(\Omega)\cap L^2_\omega(\Omega)\,.\] Condition 7 and the Schwarz inequality imply that for every \(v\in L^{2}_\omega(\Omega)\), the measure \(v{\,{\rm d}x}\) belongs to \({\mathscr{M}}_1(\Omega)\). Our main result establishes a striking equivalence: neural network training admits an exact variational characterization. In particular, the continuum formulation is not merely an approximation, but an exact relaxation of the finite discrete problem, which reveals its intrinsic variational structure.

From a variational perspective, this shows that no Lavrentiev gap occurs for \({\mathscr{R}}(f,\cdot)\) across measures, Sobolev maps, and smooth functions. In particular, minimizing sequences can be transferred between these classes without loss of optimality, which is a fundamental structural property for regularity theory. This bridges the discrete and continuum regimes at a deeper level, ensuring that analytical and numerical approximations faithfully capture the true variational behavior of the problem, and places neural network training firmly within the scope of variational analysis.

Let us point out the rate of convergence of the risk with respect to the width of the neural network in the spirit of [1]. Observe that \(L^{2}_\omega(\Omega)\) is a natural space to consider \({\mathscr{R}}\) and further regularity results will be provided for the problem on the right-hand side below. In the next proposition, we denote by \({\mathscr{M}}^{\rm at}_N(\Omega)\) the set of those purely atomic signed measures on \(\Omega\) that involve at most \(N\) atoms. Observe that such measures can be written as in 4 for suitable coefficients \(w_i\).

2.2 Towards justification of non-overfitting: stability of training↩︎

By expanding the risk term and using the Fubini theorem, we obtain for \(f\in L^{2}(D)\) and \(m\in {\mathscr{M}}_1(\Omega)\), that \[\label{eq:quadratic-decomposition} {\mathscr{R}}(f,m) = \|f\|_{L^2(D)}^2 -2\int_\Omega Q(\theta){\,{\rm d}}m(\theta) +\int_\Omega\int_\Omega K(\theta,{\vartheta}){\,{\rm d}}m(\theta){\,{\rm d}}m({\vartheta})\,,\tag{9}\] where \[\label{eq:bK} Q(\theta):=\int_D f(x)h(\theta,x){\,{\rm d}x}\qquad\text{and}\qquad K(\theta,{\vartheta}):=\int_D h(\theta,x)h({\vartheta},x){\,{\rm d}x}\,.\tag{10}\] In view of 3 for every \(\theta, {\vartheta}\in \Omega\), for \(C_{h}':=C_{h}\|1+|x|\|_{L^{2}(D)}\), one has: \[\label{def-CKCQ} |Q(\theta)|\leq C_{h}\,(1+|\theta|)\int_{D}(1+|x|)|f(x)|{\,{\rm d}x}\qquad\text{and} \qquad |K(\theta,{\vartheta})|\leq (C_{h}')^2\,(1+|\theta|)\,(1+| {\vartheta}|)\,.\tag{11}\] The functional \({\mathscr{R}}(f,\cdot)\) is convex on \({\mathscr{M}}_1(\Omega)\) (see Lemma 3 below) but not necessarily strictly convex. We endow the Hilbert space \({\mathcal{W}}\) with the norm \(\|u\|_{{\mathcal{W}}}^2 := \|u\|_{L^{2}_\omega(\Omega)}^2 + \|\nabla u\|_{L^{2}(\Omega)}^2\). We are thus led to introduce for every \(\alpha, \beta\geq 0\) the natural regularized objective \[\label{eq:regularized-functional} {\mathscr{F}}^{(f)}_{\alpha,\beta}(u) := \begin{cases} {\mathscr{R}}(f,u {\,{\rm d}x}) +\alpha\|u\|_{L^2_\omega}^2 +\beta\|\nabla u\|_{L^2}^2\,, & u\in \mathcal{W}\,,\\ +\infty\,, & u\in L^{2}_\omega(\Omega)\setminus {\mathcal{W}}\,. \end{cases}\tag{12}\] The objective combines a term measuring data fidelity with Sobolev and \(L^2_\omega\) regularization terms that provide coercivity and enforce regularity in the parameter space. In particular, the Sobolev term \(\|\nabla u\|_{L^2}^2\) induces a diffusive effect in the associated gradient flow, analogous to entropy regularization in Wasserstein formulations well-designed for the study of the learning process, cf. [17]–[19]. We avoid it as, in contrast to the computationally challenging Wasserstein approach, the \(L^2_\omega\) setting is naturally compatible with the study of the solutions to standard shallow neural network parametrizations. This enables the use of convex-analytic techniques together with classical regularity theory.

We notice that, when restricted to its domain \({\mathcal{W}}\), the regularized functional from 12 is Fréchet differentiable and \(2\min(\alpha, \beta)\)-convex, i.e., the functional \(u\mapsto{\mathscr{F}}^{(f)}_{\alpha,\beta}(u)-\min(\alpha, \beta)\|u\|_{{\mathcal{W}}}^2\) is convex on \({\mathcal{W}}\). Moreover, on \(L^{2}_\omega(\Omega)\), the functional \({\mathscr{F}}^{(f)}_{\alpha,\beta}\) is \(2\alpha\)-convex.

Some of our results require to assume that \(h\) is locally Lipschitz with respect to the first variable, in the sense that for every \(R>0\), there exists \(c_{h,R}>0\) such that \[\label{eq-Phi-loc-Lip} |h(\theta, x) - h({\vartheta},x)|\leq c_{h,R}|\theta-{\vartheta}|(1+|x|)\quad \text{for every \(\theta, {\vartheta}\in \Omega\cap B_R\), for every \(x\in D\)}\,.\tag{13}\] The regularized functional enjoys strong properties due to the convex structure of \({\mathscr{F}}^{(f)}_{\alpha,\beta}\). It ensures global existence and uniqueness of the minimizer \(u_f^*\) together with an explicit Euler–Lagrange PDE. The almost \(C^{3}\) regularity of \(u_f^*\) is a major advantage: it lifts the optimal density far beyond generic solutions, opening the way to prove sharp asymptotic analysis, stability estimates, and a justification of particle approximation tools that are essential yet scarce in mean-field neural network theory.

We stress that due to convexity, our variational problem is stable with respect to the target.

As a consequence, given \(f_1, f_2\in L^{2}(D)\), the corresponding minimizers \(u_{f_1}^*, u_{f_2}^{*}\) satisfy \[\|u_{f_1}^* - u_{f_2}^*\|_{L^2_\omega(\Omega)} \le \tfrac{C}{\alpha} \|f_1 - f_2\|_{L^2(D)} \, \quad \text{and} \quad \|\nabla u_{f_1}^*-\nabla u_{f_2}^*\|_{L^{2}(\Omega)} \leq \tfrac{C}{\sqrt{\alpha\beta}} \|f_1-f_2\|_{L^{2}(D)}\,.\] The continuous dependence of the minimizers with respect to the targets \(f\) entails a similar property for the minimal values of the functionals \({\mathscr{F}}^{(f)}_{\alpha, \beta}\):

Let us also emphasize the local uniform stability property. We say that \(u \in C^2(\overline{\Omega})\) if \(u,\; \nabla u,\; D^2 u\) admit continuous extensions to \(\overline{\Omega}\) and \(\|u\|_{C^2(\overline{\Omega})} := \|u\|_{L^\infty(\Omega)} + \|\nabla u\|_{L^\infty(\Omega)} + \|D^2u\|_{L^\infty(\Omega)}<\infty.\)

It follows that for every \(f_1, f_2\in L^{2}(D)\) and every \(\Omega'\Subset \Omega\), the corresponding minimizers \(u^*_{f_1}\) and \(u^*_{f_2}\) satisfy \[\|u^*_{f_1}-u^*_{f_2}\|_{C^{2}(\overline{\Omega'})} \leq C \|f_1-f_2\|_{L^{2}(D)}\,.\]

As a consequence of Proposition 4, Proposition 5 and the Fubini theorem, we get:

Propositions 4 and 5, together with Corollary 2, establish Lipschitz dependence of the minimizer on the target, with generalization error scaling linearly in the noise level \(\varsigma\) and controlled by \(1/\alpha\). The near-\(C^3\) regularity ensures that the local geometry of the learned representation varies smoothly with the data, providing a concrete form of implicit bias toward smooth, well-generalizing solutions.

In order to describe the target of the neural network related to \({\mathscr{F}}_{\alpha,\beta}^{(f)}\), we consider the gradient flow of this functional in the Hilbert space \(L^2_\omega(\Omega)\), following [25]. This approach is natural in the analysis of overparametrized neural networks because the gradient flow corresponds to the continuous-time limit of gradient descent (the steepest descent curve of \({\mathscr{F}}^{(f)}_{\alpha,\beta}\) in the \(L^2_\omega\)-geometry). It provides a powerful analytic framework for our novel study of convergence, implicit bias, and generalization properties that are often difficult to capture in discrete time.

Since \({\mathscr{F}}^{(f)}_{\alpha, \beta}\) is quadratic, the existence of the gradient flow follows from the Hille–Yosida theory. We introduce the unbounded linear operator \(A:D(A)\subset L^{2}_\omega(\Omega)\to L^{2}_\omega(\Omega)\) defined by its domain \[D(A):=\left\lbrace u\in {\mathcal{W}}: \exists p=p(u)\in L^{2}_{\omega}(\Omega) \textrm{ such that } \forall v\in {\mathcal{W}}, \int_{\Omega}\nabla u \cdot \nabla v {\,{\rm d}\theta}= -\int_{\Omega}pv\omega{\,{\rm d}\theta}\right\rbrace\,\] and \(A u:=2\alpha u -2\beta p(u)+\frac{2}{\omega}\int_{\Omega}K(\theta, \cdot)u(\theta){\,{\rm d}\theta}\). Note that here the function \(p(u)\) is \(\Delta u/\omega\).

It follows from the Hille–Yosida theorem complemented by the Brézis–Komura theorem that:

In the case of shallow networks, one can approximate the minimizer to \({\mathscr{F}}^{(f)}_{\alpha,\beta}\) by considering the restriction of the functional to a finite dimensional space \(E_\ell := \mathrm{span}\{e_0,\dots,e_\ell\} \subset \mathcal{W}\), where \(\{e_i\}_{i\ge 0}\) forms a complete basis of the Hilbert space \(L^2_\omega(\Omega)\). The problem thus reduces to solving a linear system, see Section 3. In order to pass to the multi-layer case, the problem starts to be strongly nonlinear and analytically infeasible. Hence, another learning algorithm is necessary. One of the ideas could be to adapt the implicit Euler scheme associated to the gradient flow, see e.g. [25]. Given \(\tau>0\) and an initial condition \(u\in {\mathcal{W}}\), we define \[\label{eq:jko} \tilde{\varrho}^0 = u, \qquad \tilde{\varrho}^{k\tau} = \arg\min_{v \in L^2_\omega(\Omega)} \left\{ {\mathscr{F}}_{\alpha,\beta}^{(f)}(v) + \tfrac{1}{2\tau} \|v - \tilde{\varrho}^{(k-1)\tau}\|_{L^2_\omega}^2 \right\}\,.\tag{15}\] The scheme admits minimizers due to the convexity and the lower semicontinuity of \({\mathscr{F}}_{\alpha,\beta}^{(f)}\). Defining the piecewise constant interpolation \(\tilde{\varrho}_\tau(t)\) by 37 , one obtains convergence toward the gradient flow as \(\tau \to 0\), cf. Theorem 11. Combining this with Corollary 3 yields the error estimate:

3 Numerical experiments↩︎

The experiments below illustrate the key properties established in the theoretical analysis: no Lavrentiev gap, \(\lambda\)-convexity, regularity of the optimal parameter density. We note that our numerical findings show that estimations of \({\mathscr{F}}^{(f)}_{\alpha,\beta}\) are both accurate and robust to noise and outliers. All simulations are done on standard personal computer and can be found in the supplementary material. Reduction to a linear system. We approximate the parameter density \(u \in {\mathcal{W}}\) by a polynomial ansatz \(u(\theta) \approx \sum_{i=1}^M a_i\,\widehat{u}_i(\theta).\) We have used three different types of basis functions \(\widehat{u}_i\): polynomials, orthonormal cosine basis and Legendre polynomials. We minimize \({\mathscr{F}}_{\alpha,\beta}^{(f)}\) over the coefficient vector \(\vec{a}\in{\mathbb{R}}^M\). Since the functional \({\mathscr{F}}_{\alpha,\beta}^{(f)}\) involves expectations with respect to the (unknown) data distribution, it cannot be evaluated directly. We therefore work with its empirical version, obtained by replacing the population quantities with averages over the observed data. Given data points \((x_k, f(x_k))_{k=1}^{N_D}\) , the functional reduces to the sum of a quadratic form and an affine term: \[\widehat{{\mathscr{F}}}_{\alpha,\beta}^{\;(f)}(a) = C_D\left( \vec{f}\,^\top \vec{f} - 2\vec{f}\,^\top U\vec{a} + \vec{a}^\top U^\top U \,\vec{a}\right)+{\vec{a}}^\top(\alpha V + \beta W)\,\vec{a}\,,\] where \(C_D=\frac{\mathcal{L}^d(D)}{N_D}\) (here, \(\mathcal{L}^d(D)\) is the Lebesgue measure of \(D\) ), \(\vec{f} = (f(x_1),\ldots,f(x_{N_D}))^\top\), \(\vec{a}=(a_1,\ldots,a_M)^\top\), and the matrices \[U_{ki} = \int_\Omega h(\theta, x_k)\,\widehat{u}_i(\theta){\,{\rm d}}\theta, \quad V_{ij} = \langle \widehat{u}_i , \widehat u_j \rangle_{L^2_\omega (\Omega)}, \quad W_{ij} = \langle \nabla \widehat{u}_i, \nabla \widehat{u}_j \rangle _{L^2 (\Omega)}\] are determined by the given data, the parameter domain \(\Omega\), and the weight \(\omega(\theta)=1+|\theta|^{2d+4}\). Crucially, no iterative optimization or gradient descent is required. We can find the minimizer directly by solving the linear equation: \[\nabla \widehat{{\mathscr{F}}}^{\;(f)}_{\alpha, \beta} (\vec{a}) = 0 \implies \vec{a} = \left(U^\top U +\frac{\alpha}{C_D} V+\frac{\beta}{C_D} W\right)^{-1}\left(U^\top \vec{f}\right) \, .\] The matrix \(U^\top U + \frac{\alpha}{C_D} V + \frac{\beta}{C_D} W\) is positive semidefinite (and moreover positive definite for \(\alpha>0\) in the case the basis functions \(\{\widehat{u}_i\}_{i=1}^{M}\) are linearly independent in \(L^2_\omega(\Omega)\) or \(\beta>0\) and the gradients of the basis functions are linearly independent in \(L^{2}(\Omega)\)) and the minimizer is found via pseudoinverse (or inverse); the problem reduces to ridge regression.

Let us compare: (i) the unregularized case \(\alpha=\beta=0\), which corresponds to the least-squares projection onto the polynomial ansatz; and (ii) a single-hidden-layer network \(f_N\) with \(N=10{\,}000\) ReLU neurons trained by Adam optimizer [31], representing the finite-width benchmark. In Example 4 NN was trained via SGD due to the instability of Adam.

3.0.0.1 Example 1: Sinus function for \(\boldsymbol{d=1}\).

We generate \(N_D = 50\) observations of the function \(\sin(7x)\) plus small Gaussian noise on the interval \((-1,1)\). The minimum of \({\mathscr{F}}^{(f)}_{\alpha, \beta}\) is approximated by a polynomial 38 of order \(s=15\), thus \(M = 136\). We have chosen \(\Omega=(-R,R)\times (-L,L)\) with \(R=L=7\). See the results on Figure 1 (a) and compare the prediction for the single-layer neural network with \(N=10\,000\) neurons (NN), our prediction for \(\alpha=\beta=0\) and then with small penalties \(\alpha=8.8\times10^{-12}\), \(\beta=8.8\times 10^{-10}\).

Figure 1 (a) shows that the regularized solution tracks the target accurately and smoothly, while the unregularized ansatz overfits and the neural network baseline is noisier. The Sobolev penalty makes the optimal density far from oscillatory, producing a function concentrated near a smooth manifold reflecting near-\(C^3\) regularity of Theorem 3 and the implicit bias toward smooth solutions. In addition, surprisingly, the regularized version is able to recover almost fully the next period of the \(\sin\) function.

3.0.0.2 Example 2: Discontinuous target — sign function (\(d=1\))

We generate \(N_D = 50\) observations of \(\mathrm{sign}(x)\) on \((-1,1)\) with small Gaussian noise, including one outlier, and use trigonometric polynomials \(\{\widehat{u}_i\}_{i=1}^{M}\) 39 up to the coefficient \(s=60\) (\(M=1{,}891\)). We have chosen \(\Omega=(-R,R)\times (-L,L)\) with \(R=5\), \(L = 5.1\). For the regularized functional, we have taken \(\alpha = 4\times 10^{-8}\) and \(\beta = 2\times 10^{-7}\).

Figure 1 (b) illustrates the stability stated in Propositions 4–5: the minimizer varies continuously with the data (Corollary 2), avoiding the large excursions of the unregularized case, while the no-Lavrentiev-gap property (Theorem 1) guarantees that the polynomial approximation incurs no hidden error. Additional pictures are given in Appendix, see Figures 2 and 3.

3.0.0.3 Examples 3–4: Benchmark datasets (\(d>1\))

We validate the approach on two standard regression benchmarks, comparing against the single-hidden-layer network with \(\mathrm{ReLU}\) activation function with \(N=10{\,}000\) neurons. In both examples we have split the dataset into the train and test with ratio of test to be \(0.2\).

In Example 3, we have used Diabetes dataset [32] containing \(d=10\) features with \(N_D=442\) observations. We have taken \(\Omega=B_{1}^{2}\times (-1,1)\) with \(B_{1}^2\) the unit ball in \({\mathbb{R}}^2\), and \(s=5\) implying \(M=4{,}368\). The sparsity of the matrix \(U\) is \(64\,\%\). We have chosen very small coefficients \(\alpha/C_D = \beta/C_D =10^{-10}\).

In Example 4, we have selected California housing dataset [33]. We take \(\Omega=B_{1}^{8}\times (-1,1)\) (where \(B_{1}^{8}\) is the unit ball in \({\mathbb{R}}^8\)), \(\alpha/C_D = \beta/C_D =10^{-10}\), polynomials of order \(s = 6\), so that \(M = 5{,}005\). The sparsity of the matrix \(U\) is \(34\,\%\). We have selected single category of ocean proximity (near bay) leading to \(N_D = 1{,}860\) observations with \(d=8\) features.

The results presented in Table 1 confirm that the regularized functional achieves competitive or superior accuracy, aligning with Theorem 1. This yields strong finite-sample performance.

4 Conclusion↩︎

The right geometry unlocks the right tools. Casting shallow network training as a variational problem over parameter densities in a weighted Sobolev space gives simultaneous access to global \(\lambda\)-convexity, elliptic regularity, and gradient flow theory, which is a combination unavailable in either the NTK or Wasserstein frameworks. The outcome is a continuum model that is analytically controlled, dynamically stable, and exactly tied to finite-width networks via an explicit \(O(1/N)\) bound: no approximation error, no linearization, no weak regularity.

The near-\(C^3\) smoothness of optimal parameter densities is perhaps the most surprising consequence. It suggests that overparameterized networks are implicitly attracted to low-dimensional smooth structures in parameter space. This is a concrete, quantitative form of implicit bias that goes beyond what existing analyses can see. To be stressed out, it is not an artifact of the limit: the absence of a Lavrentiev gap ensures that the continuum picture faithfully reflects the discrete one. Concretely, the minimizer is found via a single ridge regression, with generalization error bounded by \(C/\alpha\), making \(\alpha\) a directly interpretable robustness parameter rather than a black-box regularizer.

The shallow setting was the necessary first step. Resolving convexity, regularity, and convergence cleanly here opens a concrete path toward multilayer architectures, sharper implicit bias descriptions, and new regularization principles grounded in the underlying variational framework.

Acknowledgments↩︎

MB was supported by internal grant for specific research No. FSI-S-26-8958. IC and BM were supported by National Science Centre Grant 2024/55/B/ST6/016.

5 Technical appendices and supplementary material↩︎

In the case it is clear from the context, we simplify the notation setting \({\mathscr{F}}_{\alpha,\beta}={\mathscr{F}}^{(f)}_{\alpha,\beta}\).

5.1 Proof of Theorem 1↩︎

Theorem 1 is a consequence of two approximation results. In the first one (Proposition 7), we regularize any \(m\in {\mathscr{M}}_1(\Omega)\) by relying on standard convolution and truncation arguments. We then prove that the functional \({\mathscr{R}}(f,\cdot)\) converges along the regularized sequence to \({\mathscr{R}}(f,m)\) by exploiting the growth assumptions satisfied by \(h\). In the second approximation result (Proposition 8), we approximate \(m\) by a sequence of purely atomic measures \(\mu_i\). In order to establish the convergence of \({\mathscr{R}}(f,\mu_i)\) to \({\mathscr{R}}(f,m)\), we strongly rely on the quadratic structure of \({\mathscr{R}}(f,\cdot)\).

We now have all the ingredients to present the proof of Theorem 1.

5.2 Proof of Proposition 2↩︎

Inspired by the proof of [1] written for probability measures, we provide related result for the more general setting of finite Borel measures. In all this section, we fix \(f\in L^{2}(D)\). For every \(N\geq 1\), for every \(\theta_i\in \Omega\) and \(w_i\in {\mathbb{R}}\) with \(1\leq i \leq N\), the measure \(\mu_N\) defined in 4 satisfies: \[\label{eq-R0-rhoN} {\mathscr{R}}(f,\mu_N)=\|f\|_{L^2(D)}^2+\frac{1}{N^2}\sum_{i,j=1}^N w_i w_j K(\theta_i, \theta_j)-\frac{2}{N}\sum_{i=1}^N w_i Q(\theta_i)\,.\tag{23}\]

We begin with the following observation:

Remember that \({\mathscr{M}}_2(\Omega)\) denotes the set of all those finite Borel measures on \(\Omega\) such that \[\int_{\Omega}|\tau|^2{\,{\rm d}}|m|(\tau)<+\infty\,.\]

We conclude this section with the proof of Proposition 2.

5.3 Proof of Theorem 3↩︎

In order to establish Theorem 3, and more specifically the existence and uniqueness of the minimizer of \({\mathscr{F}}^{(f)}_{\alpha, \beta}\), we rely on the uniform convexity of that functional. In turn, the latter is a simple consequence of the observation already formulated in 24 .

For later use, we introduce for every \(g\in L^{2}_{\omega}(\Omega)\) and every \(\alpha, \beta \geq 0\), the following variant of the functional \({\mathscr{F}}_{\alpha, \beta}\): \[\begin{gather} {\mathscr{J}}_{\alpha, \beta, g}(u):=\alpha \|u\|_{L^{2}_\omega(\Omega)}^2+\beta \|\nabla u\|_{L^{2}(\Omega)}^2 + \int_{\Omega\times \Omega}K(\theta, {\vartheta})u(\theta)u({\vartheta}){\,{\rm d}\theta}{\,{\rm d}}{\vartheta}- \int_{\Omega}g(\theta)u(\theta)\omega(\theta){\,{\rm d}\theta}\,. \end{gather}\]

The regularity of the minimizer \(u^*\) stated in Theorem 3 relies on standard elliptic estimates satisfied by the Euler equation associated to the functional \({\mathscr{F}}^{(f)}_{\alpha, \beta}\). In order to fully exploit this elliptic structure, we need to establish some regularity properties of the right-hand side of the Euler equation.

The existence and regularity of the minimizer in Theorem 3 can now be obtained through classical tools from convexity theory and bootstrap arguments often used in elliptic regularity.

We proceed with the stability results, first in \({\mathcal{W}}\) (Proposition 4) and then in \(C^{2}_{loc}(\Omega)\) (Proposition 5). We also justify the continuity of \(\inf_{L^{2}_{\omega}(\Omega)}{\mathscr{F}}^{(f)}_{\alpha, \beta}\) with respect to \(f\) (Corollary 1).

5.4 On the gradient flow in \(L^{2}_\omega(\Omega)\)↩︎

5.4.1 The Hille–Yosida approach↩︎

To study the gradient flow associated to the minimization of \({\mathscr{F}}_{\alpha, \beta}\), we first rely on the Hille–Yosida approach, see [35].

Remember that \(A:D(A)\subset L^{2}_\omega(\Omega)\to L^{2}_\omega(\Omega)\) is the unbounded linear operator defined by \[D(A)=\left\lbrace u\in {\mathcal{W}}: \exists p=p(u)\in L^{2}_{\omega}(\Omega) \textrm{ such that } \forall v\in {\mathcal{W}}, \int_{\Omega}\nabla u \cdot \nabla v = -\int_{\Omega}pv\omega\right\rbrace\,,\] \[Au=2\alpha u -2\beta p(u)+\frac{2}{\omega}\int_{\Omega}K(\theta, \cdot)u(\theta){\,{\rm d}\theta}\,.\] Observe that for every \(u\in D(A)\), the function \(p(u)\) is \(\Delta u/\omega\). Hence, \(\Delta u\in L^{2}_{loc}(\Omega)\), so that by [34], the function \(u\) is in \(W^{2,2}_{loc}(\Omega)\).

In the next lemma, we check that \(A\) satisfies all the required properties to apply the Hille–Yosida theorem.

In the proof of Proposition 6, we exploit two important results related to gradient flows in Hilbert spaces: the Hille–Yosida theorem on the one hand, and the Brézis–Komura theorem on the other hand. The latter is well-adapted to lower semicontinuous and \(\lambda\)-convex functionals. We have already checked that \({\mathscr{F}}_{\alpha, \beta}\) is \(2\alpha\)-convex on \(L^{2}_\omega(\Omega)\). We now verify that it is lower semicontinuous.

Now that all their assumptions have been verified, it remains to apply the Hille–Yosida theorem and the Brézis–Komura theorem to obtain Proposition 6.

We conclude this section by an interesting estimate on the approximation of the solution of a gradient flow by the solution of the corresponding implicit Euler scheme.

5.5 Numerical simulations↩︎

In this section, we provide more detailed explanations and derivations accompanying Section 3.

Derivation of quadratic form↩︎

In this subsection we derive the approximating functional \(\widehat{\mathscr{F}}_{\alpha, \beta}\). We denote \(u \approx \widehat{u} = \sum_{i=1}^{M} a_i \widehat{u}_i\). We consider a given dataset \(\{(x_i, f(x_i))\}_{i=1}^{N_D}\). We start with the approximating risk function: \[\begin{align} \mathscr{R} (\vec{f}, \vec{a}) &= C_D\sum_{j=1}^{N_D} \Big(f(x_j) - \sum_{i=1}^{M} a_i \int_\Omega h(\theta, x_j) \,\widehat{u}_i(\theta) {\,{\rm d}}\theta\Big)^2 = C_D\lvert \vec{f} - U \vec{a} \rvert^2\\ &= C_D\left(\vert \vec{f}\,\rvert^2 - 2\vec{f} \,^{\top} U \vec{a} + \lvert U \vec{a}\rvert^2\right) \, , \end{align}\] where \(C_D=\mathcal{L}^d(D)/N_D\) (here, \(\mathcal{L}^d(D)\) is the Lebesgue measure of \(D\)) and \[U_{ki} = \int_\Omega h(\theta, x_k) \,\widehat{u}_i(\theta) {\,{\rm d}}\theta \, .\] Now we derive the weighted Lebesgue norm \(\lVert\widehat u\rVert^2_{L^2_\omega(\Omega)}\): \[\begin{align} \lVert\widehat u&\rVert^2_{L_\omega(\Omega)} = \Big\lVert \sum_{i=1}^{M} a_u \widehat{u}_i\Big\rVert^2_{L^2_\omega(\Omega)} = \sum_{i,j}^{M} a_i a_j \langle \widehat{u}_i, \widehat{u}_j \rangle_{L^2_\omega (\Omega)} = \vec{a}^{\top} V \, \vec{a} \, ,\quad V_{ik} = \langle \widehat{u}_i, \widehat{u}_j \rangle _{L^2_\omega(\Omega)}\, . \end{align}\] And finally the gradient norm \(\|\nabla \widehat{u}\|_{L^2(\Omega)}^2\): \[\begin{align} \|\nabla \widehat{u}\|_{L^2(\Omega)}^2 &= \sum_{\ell=1}^{d+1} \| \partial_\ell u \|_{L^2(\Omega)}^2 = \sum_{\ell=1}^{d+1} \Big\lVert \sum_{i=1}^{M} a_i \partial_\ell \widehat{ u}_i \Big\rVert_{L^2(\Omega)}^2 = \sum_{d=1}^{\ell+1} \sum_{i,j = 1}^{M} a_i a_j \langle \partial_\ell\widehat{u}_i, \partial_\ell\widehat{u}_j \rangle_{L^2 (\Omega)} \\ & = \vec{a} \cdot W \vec{a} \, , \quad W_{ij} = \langle \nabla \widehat{u}_i, \nabla \widehat{u}_j \rangle_{L^{2}(\Omega)} \, . \end{align}\]

Basis functions↩︎

We have used the following type of basis functions \(\{\widehat {u}_i \}_{i=1}^{M}\):

1. Polynomials on \(\Omega = B_{R}^d \times (-L,L)\) (where for \(d \geq 1\) the notation \(B_{R}^d\) refers to the ball of radius \(R\) and center \(0\) in \({\mathbb{R}}^d\)): \[\label{eq:polynomial} \widehat{u}_i(\theta) = \prod_{j=1}^{d+1} \theta_j^{p^{i}_j}\tag{38}\]

2. Trigonometric orthonormal basis on \(\Omega = (-R,R)\times (-L,L)\) (for \(d=1\)): \[\label{eq:harmonic} \widehat{u}_i (\theta) = c_1 c_2 \cos \Big(\tfrac{p_1^i \pi (\theta_0 + L)}{2L} \Big)\cos \Big(\tfrac{p_2^i \pi (\theta_1 + R)}{2R} \Big)\tag{39}\] with \[c_1 = \textstyle\begin{cases} 1/\sqrt{2L} & p_1^i >0 \,,\\ 1/\sqrt{L} & p_1^i =0\,, \end{cases} \qquad \qquad c_2 = \textstyle\begin{cases} 1/\sqrt{2R} & p_2^i >0\,, \\ 1/\sqrt{R} & p_2^i =0\,. \end{cases}\] These functions are the solutions of the eigenvalue problem for the Laplacian: \(\Delta \widehat{u}_i = \lambda\widehat{u}_i\) with zero Neumann boundary condition.

Number of basis functions \(M\)↩︎

Given \(s\in \mathbb{N}_0\), we consider the set of polynomials of degree not larger than \(s\) on the space \(\Omega \subseteq {\mathbb{R}}^{d+1}\). The exponents of the monomials involve all possible tuples \(p^i=(p^i_j)_{1\leq j\leq d+1}\) such that \(\sum_{j=1}^{d+1} p^i_j \leq s\). The number \(M\) of such tuples can be computed as follows.

We define an additional element of each tuple variable \(p^i_{d+2} := s - \sum_{j=1}^{d+1} p^i_j\) and thus the condition \(\sum_{j=1}^{d+1} p^i_j \leq s\) is equivalent to \(p^i_{d+2} \ge 0\). We obtain \(\sum_{j=1}^{d+2} p^i_j = s\) with all variables satisfying \(p^i_j \in \mathbb{N}_{0}\). So the number of all such possible tuples is given by the stars and bars formula: \[\textstyle \binom{s + (d+2) - 1}{(d+2) - 1} = \binom{s + d + 1}{d + 1} \,.\]

Additional examples↩︎

Example 2’: Discontinuous target — sign function (\(d=1\))

We consider the same dataset as in Example 2. We approximate by polynomial functions \(\{u_i\}_{i=1}^{M}\). The minimum of \({\mathscr{F}}^{(f)}_{\alpha, \beta}\) is approximated by a polynomial 38 of order \(s=60\), thus \(M = 1{,}891\). We have chosen \(\Omega=(-R,R)\times (-L,L)\) with \(R=L=1.5\). For the regularized functional, we have taken \(\alpha = 3.2\times 10^{-6}\) and \(\beta = 2 \times 10^{-5}\). See the results on Figure 2. It is worth to emphasize that the smooth minimizer of the regularized functional does not reflect the outlier.

Example 5: Sinus of norm value for (\(d=2\))

We generate \(N_D = 100\) observations on the domain \((-1,1)^2\) of the function \(f=\sin(3\lvert x\vert)\). We approximate the minimum by polynomials 38 of order \(s=5\), so that \(M=112\). We have chosen \(R=L=5\) and \(\alpha = 4\times 10^{-9}\) and \(\beta=4\times 10^{-8}\). We note that the sparsity of the matrix \(U\) is \(\approx10\,\%\). See the result of approximation by the regularized functional \(\widehat{\mathscr{F}}_{\alpha,\beta}^{\;(f)}\) on Figure 3.

References↩︎