Sobolev Approximation By ReLU Networks

Why Sobolev Approximation Tracks Derivatives

A function can be close to a target in value while having a noticeably different slope. In learning problems, that distinction matters when the target represents a response surface whose gradients are meaningful, or when a loss penalizes derivative error. Sobolev approximation makes the distinction precise: the target’s regularity is measured through its weak derivatives, and approximation error can be measured in the corresponding function and derivative norms.

For ReLU networks, there is an important boundary on this idea. Their outputs are continuous and piecewise affine, so their first weak derivatives are well defined as functions almost everywhere. But their second derivatives are generally not ordinary integrable functions. Thus, standard guarantees in this setting concern function values and first derivatives, not arbitrary higher derivative norms.

Definition: For an open set $\Omega\subset\mathbb R^d$, an integer $s\geq 0$, and $1\leq p<\infty$, the Sobolev space $W^{s,p}(\Omega)$ consists of functions $f\in L^p(\Omega)$ whose weak partial derivatives $D^\alpha f$ belong to $L^p(\Omega)$ for every multi-index $\alpha$ with $|\alpha|\leq s$. We use the norm $$ \|f\|_{W^{s,p}(\Omega)} = \left(\sum_{|\alpha|\leq s}\|D^\alpha f\|_{L^p(\Omega)}^p\right)^{1/p}. $$ For $s=0$, this is the $L^p$ norm. The $W^{1,p}$ norm measures both function error and error in all first weak partial derivatives.

A ReLU network here means a feedforward network whose hidden units apply $\operatorname{ReLU}(t)=\max\{t,0\}$, with an affine output layer. Its size can be measured by the number of nonzero weights and biases. Network size is not the same as grid resolution: a construction may use a grid with $M$ subdivisions per coordinate, so that the number of grid cells grows like $M^d$.

The Sobolev Approximation Rate

The following standard quantitative approximation theorem gives the rate on the unit cube. It separates the effect of smoothness from the effect of dimension. The constants may depend on the smoothness, dimension, integrability exponent, and domain, but not on the target function or the resolution $M$.

Theorem: Let $\Omega=(0,1)^d$, let $1\leq p<\infty$, and let $s\geq 2$ be an integer. There is a constant $C$ such that for every $f\in W^{s,p}(\Omega)$ and every integer $M\geq 2$, a ReLU network $\Phi_M$ can be chosen with at most $C M^d\log M$ nonzero parameters and depth at most $C\log M$, satisfying $$ \|f-\Phi_M\|_{L^p(\Omega)} \leq C\|f\|_{W^{s,p}(\Omega)}M^{-s}, \qquad \|f-\Phi_M\|_{W^{1,p}(\Omega)} \leq C\|f\|_{W^{s,p}(\Omega)}M^{-(s-1)}. $$ The $W^{1,p}$ bound includes both the function error and all first-derivative errors. The theorem is stated for a cube; analogous results hold on bounded Lipschitz domains when the relevant Sobolev extension and approximation conditions hold, with constants depending on the domain.

The rate is naturally read in two stages. A resolution $M$ gives the displayed error orders, while the number of parameters is on the order of $M^d\log M$. Ignoring logarithmic factors, a budget of $P$ parameters corresponds to $M$ of order $P^{1/d}$. The resulting rates are therefore approximately $P^{-s/d}$ in $L^p$ and $P^{-(s-1)/d}$ in $W^{1,p}$. Higher smoothness improves these powers; higher dimension weakens them.

The theorem is an existence guarantee, not a claim that every network trained by a particular algorithm will attain these errors. Its construction uses local approximations to the target on small regions and ReLU networks to assemble approximations with controlled complexity. The regularity assumption supplies the approximation order; the network construction supplies a representation with the stated size.

Worked Example: Reading the Rate for a Smooth Surface

Let $f(x,y)=x^2+y^2$ on $(0,1)^2$, and use $p=2$. Here $d=2$ and $s=2$. Its Sobolev norm can be calculated directly. First, $$ \|f\|_{L^2}^2 =\int_0^1\int_0^1(x^2+y^2)^2\,dx\,dy =\frac25+\frac29 =\frac{28}{45}. $$

The first derivatives are $2x$ and $2y$, so the sum of their squared $L^2$ norms is $4/3+4/3=8/3$. The nonzero second derivatives are $f_{xx}=2$ and $f_{yy}=2$, whose squared $L^2$ norms sum to $8$. Under the stated norm convention, $$ \|f\|_{W^{2,2}}^2 =\frac{28}{45}+\frac83+8 =\frac{508}{45}. $$

With grid resolution $M=16$, the theorem gives a network with at most $C\,16^2\log(16)$ nonzero parameters and depth at most $C\log(16)$. Its errors satisfy $$ \|f-\Phi_{16}\|_{L^2} \leq C\sqrt{\frac{508}{45}}\,16^{-2}, \qquad \|f-\Phi_{16}\|_{W^{1,2}} \leq C\sqrt{\frac{508}{45}}\,16^{-1}. $$ The constants are not specified numerically by the theorem. The calculation shows the distinct orders: the function error decreases as $M^{-2}$, while the joint function-and-gradient error decreases as $M^{-1}$.

A Fully Worked One-Dimensional Derivation

In one dimension, linear interpolation makes the role of ReLU particularly transparent. Divide $[0,1]$ into $M$ equal intervals with nodes $x_i=i/M$, and let $I_M f$ be the continuous function that agrees with $f$ at each node and is linear on each interval. The following derivation proves a useful Sobolev error bound for this interpolant.

Identify the interval slope.
On an interval $I=[a,b]$ of length $h=1/M$, the slope of $I_M f$ is $(f(b)-f(a))/h$. For $f\in W^{2,p}(0,1)$, the one-dimensional absolutely continuous representatives of $f$ and $f'$ give $f(b)-f(a)=\int_a^b f'(t)\,dt$. Thus the slope is the average of $f'$ on $I$.

Bound the derivative error.
The derivative error on $I$ is $f'(x)-h^{-1}\int_I f'(t)\,dt$. The mean-subtraction Poincaré inequality on an interval gives its $L^p(I)$ norm at most $h\|f''\|_{L^p(I)}$. Summing the $p$-th powers across the disjoint intervals gives $\|f'-(I_M f)'\|_{L^p(0,1)}\leq h\|f''\|_{L^p(0,1)}$.

Bound the function error.
The error $e=f-I_M f$ vanishes at both endpoints of every interval. For $x\in I$, $e(x)=\int_a^x e'(t)\,dt$. Hölder’s inequality gives $\|e\|_{L^p(I)}\leq h\|e'\|_{L^p(I)}$. Summing across intervals and using the derivative bound yields $\|f-I_M f\|_{L^p(0,1)}\leq h^2\|f''\|_{L^p(0,1)}$.

Together these inequalities prove $$ \|f-I_M f\|_{L^p(0,1)} \leq M^{-2}\|f''\|_{L^p(0,1)}, \qquad \|f'-(I_M f)'\|_{L^p(0,1)} \leq M^{-1}\|f''\|_{L^p(0,1)}. $$ They apply for $1\leq p<\infty$; the same interval argument also yields the corresponding essential-supremum bounds for $p=\infty$. This one-dimensional result displays the two powers directly: one additional integration of the derivative error gives the extra factor of $M^{-1}$ for function values. $\square$

Worked Example: Exact Errors for a Quadratic

Take $f(x)=x^2$ on $[0,1]$. On an interval $[a,a+h]$, the interpolating line has slope $$ \frac{(a+h)^2-a^2}{h}=2a+h. $$ Writing $x=a+t$, where $0\leq t\leq h$, gives $$ (I_M f)(x)-f(x)=t(h-t), \qquad (I_M f)'(x)-f'(x)=h-2t. $$

The first expression is nonnegative, has maximum $h^2/4$ at $t=h/2$, and the absolute derivative error has maximum $h$ at either endpoint. Thus the $L^\infty$ errors are exactly $h^2/4$ for function values and $h$ for the derivative.

For $p=2$, the squared function error on one interval is $$ \int_0^h t^2(h-t)^2\,dt=\frac{h^5}{30}. $$ There are $M$ intervals and $h=1/M$, so the total squared error is $Mh^5/30=h^4/30$, giving $L^2$ error $h^2/\sqrt{30}$. The squared derivative error on one interval is $$ \int_0^h(h-2t)^2\,dt=\frac{h^3}{3}. $$ Summing gives $Mh^3/3=h^2/3$, so the derivative $L^2$ error is $h/\sqrt{3}$. Both calculations verify the respective second- and first-order rates.

Why These Splines Are ReLU Networks

The interpolant just used is not merely piecewise linear; it has an exact ReLU representation. Let $s_i=(f(x_{i+1})-f(x_i))/(x_{i+1}-x_i)$ be the slope on interval $[x_i,x_{i+1}]$. Then $$ I_M f(x) =f(x_0)+s_0(x-x_0) +\sum_{i=1}^{M-1}(s_i-s_{i-1})\operatorname{ReLU}(x-x_i). $$

To verify the formula, take $x$ in interval $[x_k,x_{k+1}]$. The active hinge terms are exactly those with $i\leq k$. The slope of the right-hand side there is $$ s_0+\sum_{i=1}^k(s_i-s_{i-1})=s_k, $$ because the sum telescopes. At $x=x_0$, the expression equals $f(x_0)$. It is continuous at every knot, and on each interval it has the same slope and starting value as the interpolant. It therefore agrees with $I_M f$ everywhere. The affine term can also be represented using ReLU units because $x=\operatorname{ReLU}(x)-\operatorname{ReLU}(-x)$; an affine output bias handles the constant. This proves exact representation by a one-hidden-layer ReLU network with a number of units proportional to $M$.

Worked Example: An Exact ReLU Representation of a Cusp

Consider $f(x)=|x-\tfrac12|$ on $[0,1]$. The identity $|z|=\operatorname{ReLU}(z)+\operatorname{ReLU}(-z)$ gives $$ f(x)=\operatorname{ReLU}(x-\tfrac12)+\operatorname{ReLU}(\tfrac12-x). $$ Thus one small ReLU network represents this function exactly, with zero error in both $L^p$ and $W^{1,p}$ for every finite $p$. Indeed, its weak derivative is $-1$ on $(0,\tfrac12)$ and $1$ on $(\tfrac12,1)$, apart from the irrelevant value at the knot.

The cusp also shows why the smoothness assumption matters. The first derivative jumps at $1/2$, so it has no absolutely continuous representative there and is not in $W^{1,p}(0,1)$. Consequently $f\notin W^{2,p}(0,1)$. Exact representation of this particular nonsmooth function does not contradict a theorem requiring $W^{2,p}$; a smoothness-based rate is a uniform guarantee over a function class, not a statement that every individual target needs that much regularity.

Domain, Norm, and Smoothness Conditions

The rates depend on more than the activation function. The domain in the theorem is a cube, where regular subdivisions are available. On a bounded Lipschitz domain, extension results can transfer Sobolev regularity to a surrounding region, but the constants then depend on the domain and the construction must handle its boundary. Highly irregular domains do not automatically inherit the same statement with the same constants.

The target must also have the indicated weak derivatives in $L^p$. In the theorem, $s\geq2$ gives enough regularity for both the function and first-derivative bounds. The stronger $M^{-s}$ function rate and the $M^{-(s-1)}$ first-derivative rate reflect that a derivative consumes one order of smoothness. They should not be interchanged: a small function error alone does not imply a small gradient error.

Finally, the derivative norm is restricted to order one for a structural reason. A ReLU network is piecewise affine. Its first derivatives are piecewise constant and are weak $L^p$ functions, but slope changes across cell boundaries create singular contributions to second derivatives. In general, a ReLU network is not in $W^{2,p}$. A claim of $W^{2,p}$ approximation by such networks would therefore require special qualifications or a different approximation notion.

Common Mistake

Wrong reasoning: “The network approximates $f$ very closely in $L^p$, so it must also approximate the derivatives closely.” These are different norms; small function error alone gives no such implication.
Correction: State the target smoothness, the domain, and the norm in the claimed result. Use a theorem that explicitly controls $W^{1,p}$ when gradients matter. For ReLU, do not casually extend the claim to second derivatives, which generally are not $L^p$ functions.

Key takeaway: On a regular $d$-dimensional domain, Sobolev smoothness $s$ leads to function error of order $M^{-s}$ and first-derivative error of order $M^{-(s-1)}$ at grid resolution $M$, with network size growing roughly like $M^d$ up to logarithmic factors. The derivative order, domain conditions, and target regularity are part of the guarantee.

Check Your Understanding

Use the definitions and rates above to answer these questions.

For $s=3$, what are the stated powers of $M$ for $L^p$ error and $W^{1,p}$ error?
Why does converting grid resolution $M$ to a parameter budget introduce dimension into the rate?
For linear interpolation of $x^2$, what are the exact $L^2$ errors in function value and derivative when the interval width is $h$?
What feature of a ReLU spline gives rise to the telescoping-slope representation?
Why does exact ReLU representation of $|x-\tfrac12|$ not establish a general $W^{2,p}$ approximation guarantee?