# On the limit set of Neural Network functions trainable by Machine Learning

Bruno Després

<figure><img src="image_1.png" /></figure>

## To cite this version:

Bruno Després. On the limit set of Neural Network functions trainable by Machine Learning. 2026. [hal-05497591v2](https://hal.sorbonne-universite.fr/hal-05497591v2)

HAL Id: hal-05497591

https://hal.sorbonne-universite.fr/hal-05497591v2

Preprint submitted on 13 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# On the limit set of Neural Network functions trainable by Machine Learning*

Bruno Després†‡

April 13, 2025

## Abstract

We define a functional space together with a convenient norm for the set of functions generated by modern Neural Network structures [4] which can be trained by Machine Learning algorithms. The key is to concentrate on an effective version of the chain rule for functions which are slightly more regular than Lipschitz continuous. Two applications are detailed with one numerical illustration in Pytorch.

## 1 Introduction

We define and study a functional framework for functions generated by modern Neural Network structures [4] which can be trained by Machine Learning algorithms. The mathematical nature of the framework is twofold. Firstly we establish some differential rules for long chains or layers of composition of elementary functions of minimal regularity. Secondly we provide a functional space such that one can pass to the limit with respect to various parameters such as the number of layers or any similar parameter. This space is the subset of Lipschitz-continuous functions with an associated gradient compatible with the chain rule in the sense of Murat-Trombetti [32, 18].

We concentrate on the mathematical structure of these functions, independently of any application of Neural Networks. Such applications are numerous in present days. Suffice to mention applications to images processing, NLP (natural Neural language), supervised learning, transformers, SciML (Scientific Machine Learning), modeling of meteorology, game theory, IA (Artificial Intelligence), computational protein design and prediction, optimization and numerical optimization, and many other fields.

We start with a vague definition. We will say $f: \mathbb{R}^{a_0} \to \mathbb{R}^{a_{\ell+1}}$ is a Neural Network function if it can be written as the composition of finite number of finite

<footnote>*This work received funding by Agence Nationale de la Recherche, program France 2030, reference ANR-23- PEIA-0004.</footnote>

<footnote>†Sorbonne Université, Université Paris Citée, CNRS, LJLL, F-75005 Paris, France</footnote>

<footnote>‡MEGAVOLT, INRIA, Paris, France</footnote>
dimensional prescribed functions

$$f = f_{\ell} \circ S_{\ell} \circ f_{\ell-1} \circ S_{\ell-1} \circ \cdots \circ f_1 \circ S_1 \circ f_0, \quad (1)$$

where the prescribed finite dimensional functions $f_r$ and $S_r$ belong to a given dictionary and have parameters that can be optimized for a given task. The parameter $\ell$ is identified with the number of hidden layers. The functions $f_r$ for $r = 0, \dots, \ell$ are usually affine functions with varying input and output dimensions $(a_0, a_1, \dots, a_{\ell+1}) \in \mathbb{N}^{\ell+2}$ with $a_r > 0$, for $r = 0, \dots, \ell$

$$f_r(x_r) = W_r x_r + b_r \in \mathbb{R}^{a_{r+1}} \text{ for all } x_i \in \mathbb{R}^{a_r}.$$

The parameters [4] are the matrices $W_r$ for $r = 1, \dots, \ell$ and the vectors $b_r$ for $r = 1, \dots, \ell$. The hyper-parameters [4] $a_r$ for $r = 1, \dots, \ell$ are identified as the number of neurons on a given hidden layer. The intermediate functions $S_r$ for $r = 1, \dots, \ell$ belong to the dictionary of nonlinear activation functions [39] which play a key role in Neural Networks. Activation functions are defined for single scalar input then applied component wise to any vector of arbitrary size. Activation functions are normalized $0 \le S'_r(x) \le 1$ for almost all $x \in \mathbb{R}$. Among standard activation functions, one distinguishes sigmoidal functions which are monotone increasing from 0 to 1 and regular $S_i \in C^{\infty}(\mathbb{R})$. A popular sigmoidal function is $S_r(x) = 1/(1 + e^{-x})$. Another popular activation function [4] is the ReLU function $S_r = R$

$$R(x) = \max(0, x).$$

Clearly the ReLU function is not differentiable at $x = 0$. The point-wise non differentiability of ReLU is shared with other basic functions in Neural Networks. For example maxpooling used in convolutional Neural Networks [4] is based on the maximum function $(a, b) \mapsto \max(a, b)$ which is also not uniquely differentiable for $a = b$. In practice maxpooling is activated on tensors of arbitrary dimensions. A function $f: \mathbb{R}^a \to \mathbb{R}^b$ is Lipschitz-continuous if there exists $L \ge 0$ such that $|f(x) - f(x')| \le L|x - x'|$ for all $x, x' \in \mathbb{R}^a$. Since all activation elementary functions $f_i$ and $S_i$ are Lipschitz-continuous by assumption, then the resulting function $f$ in (1) is also Lipschitz continuous, that is

$$f \in \mathrm{Lip}(\mathbb{R}^{a_0} : \mathbb{R}^{a_{\ell+1}}). \quad (2)$$

The Rädemacher Theorem [36, 31] states that Lipschitz functions are differentiable almost everywhere, that is up to sets of zero Lebesgue measure. The gradient, interpreted in the sense of distributions, is a bounded matrix valued function. One has $\nabla f_r \in L^\infty(\mathbb{R}^{a_0} : \mathcal{M}_{a_{r+1}, a_0}(\mathbb{R}))$, $\nabla S_r \in L^\infty(\mathbb{R}^{a_0} : \mathcal{M}_{a_r, a_0}(\mathbb{R}))$ and

$$\nabla f \in L^\infty(\mathbb{R}^{a_0} : \mathcal{M}_{a_{\ell+1}, a_0}(\mathbb{R})) \quad (3)$$

where real matrices with $a$ lines and $b$ columns are denoted as $\mathcal{M}_{a, b}(\mathbb{R})$. Function $f_r$ is linear by assumption so $\nabla f_r = W_r$ is a constant matrix which is of course bounded as a function.

We will assume that all internal functions $f_r$ and $S_r$ in the composition formula (1) are at least Lipschitz continuous. However this class of regularity
reveals not sufficient for our purposes. The needed class of regularity is explain below.

*We will say that a Neural Network function $f$ is trainable by Machine Learning if one can write a chain rule formula for the derivative $\nabla f$ where all intermediate matrices are defined everywhere (that is not only up to sets of vanishing Lebesgue measure).*

This functional constraint comes from the fact that Machine Learning is a computational procedure which itself has its own constraints. As noticed in [29] the derivative of a function implemented in Machine Learning softwares [4] needs to have a generalized derivative (denoted as an intentional derivative in [29]) defined for value of the input. The reason is that Machine Learning is universally based on the chain rule formula to evaluate a discrete gradient needed for gradient descent techniques [4]. The chain rule is also referred to as backpropagation, automatic differentiation or autodiff in Machine Learning literature [4]. Automatic differentiation is also used for the implementation of PINNs [28, 37], DeepRitz methods [44] and Variational PINNs [25, 6]. Without a generalized derivative compatible with the chain rule, any Machine Learning computational tool is useless. In terms of mathematical regularity nevertheless, there is a clear difficulty if the different functions are only Lipschitz continuous as in (2-3). This is known as the problem of $R'(0)$ in [7], also referred to as non smooth autodiff in [10, 12, 8, 9, 14]. In [12], the problem is addressed within path differentiable functions. As noticed in [18], similar mathematical difficulties have been identified in the partial differential community in the past [40, 3, 26, 32, 5].

*We will say that a function is in the limit set of Neural Network functions trainable by Machine Learning if there is a Banach space such that the function is in the closure of the set of Neural Network functions trainable by Machine Learning.*

Indeed discrete Neural Networks models are more and more seen as discrete instances of continuous limit models. This issue is the topic of a recent work [34] and references therein. The continuous limit is natural for Neural Networks interpreted as the transport of empirical measures [2]. Transformers [42] are nowadays influential architectures in Neural Networks and Machine Learning research. See also [20] and references therein. For transformers the limit of an infinite number of tokens raises questions which are close to the ones of this work. The theory of deep networks with a very large or infinite number of layers $\ell \to \infty$ is similar. Another parameter which can tend to $\infty$ is also the number of neurons per layer $a_r \to \infty$. In our setting, we will show that a large number of layers or a large number of neurons per layer induce a large number of pieces for the representation of functions and their associated gradients.

The organization of the work is as follows. In Section 2 we review some definitions and basic properties which were proposed in the seminal work [32], then adapted to Neural Networks functions in [18]. The notion of an associated derivative is introduced and the chain rule is justified in Theorem 1 which is the Murat-Trombetti Theorem. The increased complexity of the representation is explained in Remarks 1-3 which come from [18]. A Banach space is constructed
in Section 3, which is the main theoretical contribution of this work. This Banach space is made of functions which admit an associated derivative as a Borel function. In Section 4, we detail two examples. The first example is the smooth case of a Takagi's function where the limit associated derivative can be calculated with an explicit formula. In this case the limit associated derivative is the sum of an affine polynomial plus a bounded Borel function which is zero everywhere except at a countable set of points. A numerical illustration calculated with Pytorch confirms the theoretical analysis. The other example is about the differentiability of the solution of ODEs. Simple generalizations are finally evoked.

This article tries to be self contained and the notations are kept to the minimum. The properties of Borel functions that we use can be found in [30, 38].

## 2 Definitions

A computational heuristic needed in Machine Learning is that the generalized derivative must be defined everywhere in compatibility with the chain rule formula. This computational heuristic is now elevated at the level of a mathematical design principle, independently of any other considerations.

But before, let us consider the example [32, 18] which stressed that Lipschitz continuity is not enough to address the issue. Take two functions $f_0$ and $S_1$. The first function $f_0: \mathbb{R} \to \mathbb{R}^2$ is linear $f_0(x) = (x, x) = W_0x$ with $W_0 = (1, 1)$. The second function $S_1$ called maxpooling in Neural Networks [4] is Lipschitz continuous

$$S_1(y) = \max(y_1, y_2), \text{ where } y = (y_1, y_2) \in \mathbb{R}^2. \qquad (4)$$

Then $f = S_1 \circ f_0$ is the identity function: $f(x) = x$ for all $x \in \mathbb{R}$, so that $f' \equiv 1$ is of course bounded. The gradient of $f_0$ is $\nabla f_0 = W_0$. The gradient of $S_1$ is defined almost everywhere. There are three cases: a) if $y_1 > y_2$ then $\nabla S_1(y) = (1, 0)^T$, b) if $y_1 < y_2$ then $\nabla S_1(y) = (0, 1)^T$, c) if $y_1 = y_2$ which is a set of zero measure, then $\nabla S_1(y)$ is not defined in compatibility with the fact that a Lipschitz continuous function has a bounded derivative almost everywhere. Now since $f(x) = x$, the natural chain rule formula becomes

$$1 = \nabla S_1(f_0(x))W_0. \qquad (5)$$

But $f_0(x) = (x, x)$ so only the case c) matters. Therefore $\nabla S_1(f_0(x))$ is nowhere defined as a function. The class of Lipschitz continuous functions is not regular enough to give meaning to the chain rule (5). For this example the difficulty is solved in Appendix A with a convenient representation of $S_1$ as a $C^1$-piecewise function.

### 2.1 $C^1$-piecewise functions

To solve the issue in (5), we start with the notion of a $C^1$-piecewise function, as it was introduced in the seminal work [32][Murat-Trombetti]. We use the notations from [18].
We take a decomposition of $\mathbb{R}^a$ in Borel pieces $P^\alpha \subset \mathbb{R}^a$

$$ \mathbb{R}^a = \bigcup_{\alpha} P^\alpha, \quad P^\alpha \cap P^\beta = \emptyset \text{ for } \alpha \neq \beta. $$

Following [32], the number of pieces is finite $1 \le \alpha \le \alpha_{\max} < \infty$. This constraint will be relaxed in next Sections which contain the original part of this work. For simple functions, the pieces $P^\alpha$ are usually piecewise affine with Lipschitz regularity. The appropriate notations in our context are taken from [18].

**Definition 1.** We say that a Lipschitz-continuous function $f: \mathbb{R}^a \to \mathbb{R}^b$ is $C^1$-piecewise if there exists finite number $1 \le \alpha \le \alpha_{\max}$ of Borel pieces $P^\alpha$ and functions $f^\alpha \in Lip(\mathbb{R}^a : \mathbb{R}^b) \cap C^1(\mathbb{R}^a : \mathbb{R}^b)$ with the representation

$$ f(x) = \sum_{\alpha} \mathbf{1}_{P^{\alpha}}(x) f^{\alpha}(x) \quad \forall x \in \mathbb{R}^{a}. \qquad (6) $$

Here the notation $\mathbf{1}_\omega$ denotes the indicatrix function of a set $\omega$, that is $\mathbf{1}_\omega(x) = 1$ for $x \in \omega$ and $\mathbf{1}_\omega(x) = 0$ for $x \notin \omega$. By definition $f^\alpha(x) = f(x)$ for all $x \in P^\alpha$ and $\nabla f^\alpha \in L^\infty(\mathbb{R}^a : \mathcal{M}_{b,a}(\mathbb{R})) \cap C^0(\mathbb{R}^a : \mathcal{M}_{b,a}(\mathbb{R}))$.

The set of all such $C^1$-piecewise functions will be denoted as $C_{pw}^1(\mathbb{R}^a : \mathbb{R}^b)$, with the embedding $C_{pw}^1(\mathbb{R}^a : \mathbb{R}^b) \subset Lip(\mathbb{R}^a : \mathbb{R}^b)$.

One defines the gradient associated to the above representation, which is in essence similar to the intentional derivative [29]. This is an a priori definition since it is defined independently of the chain rule that will be examined later.

**Definition 2.** We say that the gradient associated to the representation (6) is

$$ \tilde{\nabla} f(x) = \sum_{\alpha} \mathbf{1}_{P^{\alpha}}(x) \nabla f^{\alpha}(x) \in \mathcal{M}_{b,a}(\mathbb{R}) \quad \forall x \in \mathbb{R}^{a}. \qquad (7) $$

We refer to it as an associated gradient. The associated gradient is a real $b \times a$ matrix defined everywhere, that is for all $x \in \mathbb{R}^a$. It is not unique since it depends on the representation (6) which is not unique.

The main virtues of the associated gradient are that it is defined everywhere and that it is compatible with the chain rule.

**Theorem 1** (Murat-Trombetti [32]). Consider two functions. The first one $u \in Lip(\mathbb{R}^a : \mathbb{R}^b)$ is Lipschitz-continuous. The second one $v \in Lip(\mathbb{R}^b : \mathbb{R}^c)$ is Lipschitz-continuous and $C^1$-piecewise with a representation (6) and with an associated gradient (7). Then the chain rule identity holds in $L^\infty(\mathbb{R}^a : \mathcal{M}_{c,a}(\mathbb{R}))$

$$ \nabla(v \circ u) = \tilde{\nabla}v \circ u \nabla u $$

where $\tilde{\nabla}v \circ u$ denotes the function such that $\tilde{\nabla}v \circ u(x) = \tilde{\nabla}v(u(x))$ for all $x \in \mathbb{R}^a$.

The proof ultimately relies on a Stampacchia formula [40]. It must be noted that a similar framework was anticipated in [13] two years before [32]. The proof with notations adapted to the central definition of an associated gradient comes from [18].
Proof of Theorem 1. It is decomposed in three steps.

• Since $v$ is $C^1$-piecewise, it admits a piecewise smooth approximations denote as $v^\alpha \in C^1(\mathbb{R}^b)$ for all $\alpha$. Since $\nabla v^\alpha \circ u \in C^0(\mathbb{R}^a : \mathbb{R}^c)$ is a continuous function, the chain rule is applied without difficulty

$$ \nabla(v^\alpha \circ u) = \nabla v^\alpha \circ u \nabla u. \qquad (8) $$

This identity holds almost everywhere (a.e.) with respect to $x \in \mathbb{R}^b$. Let $U^\alpha$ be the set $U^\alpha = \{x \in \mathbb{R}^a : u(x) \in P^\alpha\}$. Since $U^\alpha$ is the pre-image of the Borel set $P^\alpha$ by the continuous function $u$, that is $U^\alpha = u^{-1}(P^\alpha)$, then $U^\alpha$ is also a Borel set so it is a measurable set.

• One notes that $v(u(x)) = v^\alpha(u(x))$ a.e. $x \in U^\alpha$. To use this identity one notes $w = v \circ u - v^\alpha \circ u$. Then one uses an intuitive but non trivial important property from Stampacchia [40, 26]. It writes

$$ \nabla w(x) = 0 \quad \text{a.e. } x \in \{x \in \mathbb{R}^a : w(x) = 0\}. \qquad (9) $$

Since $U^\alpha \subset \{x \in \mathbb{R}^a : w(x) = 0\}$, it yields using (8)

$$ \nabla(v \circ u)(x) = \nabla v^\alpha \circ u(x) \nabla u(x) \quad \text{a.e. } x \in U^\alpha. \qquad (10) $$

One also has

$$ \mathbf{1}_{U^\alpha}(x) = \mathbf{1}_{P^\alpha}(u(x)) \quad \text{a.e. } x \in \mathbb{R}^a. \qquad (11) $$

• Consider the difference of the two terms in the claim $A(x) = \nabla(v \circ u)(x) - \tilde{\nabla}v \circ u \nabla u$. It writes also

$$ \begin{align*} A(x) &= \nabla(v \circ u) - \left(\sum_{\alpha} \mathbf{1}_{P^\alpha}(u(x)) \nabla v^\alpha(u(x))\right) \nabla u(x) \\ &= \left(\sum_{\alpha} \mathbf{1}_{P^\alpha}(u(x))\right) \nabla(v \circ u)(x) - \left(\sum_{\alpha} \mathbf{1}_{P^\alpha}(u(x)) \nabla v^\alpha(u(x))\right) \nabla u(x) \\ &= \sum_{\alpha} \mathbf{1}_{P^\alpha}(u(x)) (\nabla(v \circ u)(x) - \nabla v^\alpha(u(x)) \nabla u(x)) \quad \text{(use (11))} \\ &= \sum_{\alpha} \mathbf{1}_{U^\alpha}(x) (\nabla(v \circ u)(x) - \nabla v^\alpha(u(x)) \nabla u(x)) \quad \text{(use (10))} \\ &= \sum_{\alpha} \mathbf{1}_{U^\alpha}(x) (0) \quad \text{(use (10))} \end{align*} $$

where all manipulations holds a.e. with respect to $x$. Therefore $A(x) = 0$ a.e. which is the claim. □

Taking $u(x) = x$ in Theorem 1, one obtains $\tilde{\nabla}v = \nabla v$ in $L^\infty(\mathbb{R}^a : \mathcal{M}_{c,a}(\mathbb{R}))$. That is an associated gradient is a representative of the gradient in the sense of distribution. Since it is defined everywhere from elements which are all bounded Borel functions, one has naturally $\tilde{\nabla}v \in l^\infty(\mathbb{R}^a : \mathcal{M}_{b,a}(\mathbb{R}))$, where the space of matrix-valued bounded Borel functions will be denoted as

$$ l^\infty(\mathbb{R}^a : \mathcal{M}_{b,a}(\mathbb{R})) \subset \mathcal{B}(\mathbb{R}^a : \mathcal{M}_{b,a}(\mathbb{R})). \qquad (12) $$

The space $l^\infty(\mathbb{R}^a : \mathcal{M}_{b,a}(\mathbb{R}))$ is closed in the vectorial space of matrix-valued bounded Borel functions equipped with the supremum norm [38][Corollary a) of Theorem 1.14] or [30][Theorem 2.5.1]. One has the embedding $l^\infty(\mathbb{R}^a : \mathcal{M}_{b,a}(\mathbb{R})) \subset L^\infty(\mathbb{R}^a : \mathcal{M}_{b,a}(\mathbb{R}))$ because bounded Borel measurable functions are also bounded Lebesgue measurable functions.
## 2.2 Properties of $C^1$-piecewise functions

The fact that the previous construction uses only a finite number of pieces lacks generality on mathematical grounds. It has important limitations for some of the applications discussed in the Introduction and for theoretical reasons which are detailed in this Section. A first evident property is the following.

**Lemma 1.** The set $C_{pw}^1(\mathbb{R}^a : \mathbb{R}^b)$ is a vectorial space.

*Proof.* Take $f, g \in C_{pw}^1(\mathbb{R}^a : \mathbb{R}^b)$ with finite representations. One has $f(x) = \sum_{\alpha} \mathbf{1}_{P^\alpha}(x)f^\alpha(x)$ for all $x \in \mathbb{R}^a$ together with $g(x) = \sum_{\beta} \mathbf{1}_{Q^\beta}(x)g^\beta(x)$ for all $x \in \mathbb{R}^a$. Define $R^{\alpha, \beta} = P^\alpha \cap Q^\beta$ which is a Borel set by intersection. Define $h^{\alpha, \beta} = f^\alpha + g^\beta \in \text{Lip}(\mathbb{R}^a : \mathbb{R}^b) \cap C^1(\mathbb{R}^a : \mathbb{R}^b)$. Then $h = f + g$ admits the representation formula

$$h(x) = \sum_{\alpha, \beta} \mathbf{1}_{R^{\alpha, \beta}}(x)h^{\alpha, \beta}(x) \text{ for all } \forall x \in \mathbb{R}^a. \quad (13)$$

So $C_{pw}^1(\mathbb{R}^a : \mathbb{R}^b)$ is invariant by additivity. It is evidently invariant by multiplication by constants. $\square$

An evident consequence is the additivity of associated gradients.

**Lemma 2.** Take $f, g \in C_{pw}^1(\mathbb{R}^a : \mathbb{R}^b)$ and $h = f + g \in C_{pw}^1(\mathbb{R}^a : \mathbb{R}^b)$. Then $\tilde{\nabla}h = \tilde{\nabla}f + \tilde{\nabla}g$ in $L^\infty(\mathbb{R}^a : M_{c,a}(\mathbb{R}))$.

*Proof.* It is sufficient to determine the associated gradient of $h = f + g$ directly from the representation (13). $\square$

**Remark 1.** A consequence of the two previous lemmas is that any time one sums up such functions, the number of pieces may increase multiplicatively.

The next result shows that $C^1$-piecewise functions are an algebra for the operation of composition. This property is not evident since the definition of an associated gradient is based on additive representations such as (6). There is in principle no reason that additivity is compatible with composition. So for general functions, such a result would clearly be wrong. However the key is that indicatrix functions establish some compatibility between additivity and composition. This was proved for the first time under this form in [18].

**Proposition 1.** Consider two functions $f \in C_{pw}^1(\mathbb{R}^a : \mathbb{R}^b)$ and $g \in C_{pw}^1(\mathbb{R}^b : \mathbb{R}^c)$. Let $h = g \circ f$ be the composition. Then $h \in C_{pw}^1(\mathbb{R}^a : \mathbb{R}^c)$.

Moreover a representation and an associated gradient for $h$ can be determined from the representation and the associated gradient of $f$ and $g$ and one has a chain rule formula for associated gradients

$$\tilde{\nabla}h(x) = \tilde{\nabla}g \circ f(x)\tilde{\nabla}f(x) \text{ for all } x \in \mathbb{R}^a. \quad (14)$$

**Remark 2.** In brief one has the embedding

$$C_{pw}^1(\mathbb{R}^b : \mathbb{R}^c) \circ C_{pw}^1(\mathbb{R}^a : \mathbb{R}^b) \subset C_{pw}^1(\mathbb{R}^a : \mathbb{R}^c).$$
Proof. The proof is similar to [18][Lemma 1]. Take $f \in C_{\text{pw}}^1(\mathbb{R}^a : \mathbb{R}^b)$ and $g \in C_{\text{pw}}^1(\mathbb{R}^b : \mathbb{R}^c)$. The representation of $f$ with finite number of pieces is

$$f(x) = \sum_{\alpha} \mathbf{1}_{P^{\alpha}}(x) f^{\alpha}(x) \text{ for all } x \in \mathbb{R}^{\alpha}.$$

The representation of $g$ with finite number of pieces is

$$g(y) = \sum_{\beta} \mathbf{1}_{Q^{\beta}}(y) g^{\beta}(y) \text{ for all } y \in \mathbb{R}^{\beta}.$$

Consider the composition of these functions $h = g \circ f$.

One defines $R^{\alpha, \beta} = \{x \in P^{\alpha} | f(x) \in Q^{\beta}\} = P^{\alpha} \cap f^{-1}(Q^{\beta})$. Since $f$ is a continuous function and $Q^{\beta}$ is a Borel piece, then $f^{-1}(Q^{\beta})$ is Borel. By intersection, $R^{\alpha, \beta}$ is a Borel set. By construction one has $R^{\alpha, \beta} \cap R^{\alpha', \beta'} = \emptyset$ for $(\alpha, \beta) \neq (\alpha', \beta')$. One also has $\bigcup_{\alpha, \beta} R^{\alpha, \beta} = \mathbb{R}^{\alpha}$. So the finite number of Borel pieces $R^{\alpha, \beta}$ realizes a decomposition of the space. One also defines the functions $h^{\alpha, \beta} = g^{\beta} \circ f^{\alpha} \in \text{Lip}(\mathbb{R}^{\alpha} : \mathbb{R}^{c}) \cap C^{1}(\mathbb{R}^{\alpha} : \mathbb{R}^{c})$.

Now to show that

$$h(x) = \sum_{\alpha, \beta} \mathbf{1}_{R^{\alpha, \beta}}(x) h^{\alpha, \beta}(x), \qquad (15)$$

one assembles $A(x) = h(x) - \sum_{\alpha, \beta} \mathbf{1}_{R^{\alpha, \beta}}(x) h^{\alpha, \beta}(x)$. Take $x \in R^{\alpha, \beta}$. Then $A(x) = g(f(x)) - h^{\alpha, \beta}(x) = g(\sum_{\alpha} \mathbf{1}_{P^{\alpha}}(x) f^{\alpha}(x)) - h^{\alpha, \beta}(x)$. Since $R^{\alpha, \beta} \subset P^{\alpha}$, one gets $A(x) = g(f^{\alpha}(x)) - h^{\alpha, \beta}(x)$. Since $f^{\alpha}(x) \in Q^{\beta}$, one gets $A(x) = g^{\beta}(f^{\alpha}(x)) - h^{\alpha, \beta}(x) = 0$ for all $x \in R^{\alpha, \beta}$. Since it is true for all $\alpha, \beta$, one gets $A \equiv 0$. We have just proved that $h \in C_{\text{pw}}^1(\mathbb{R}^a : \mathbb{R}^c)$. From (15), one immediately deduces an associated gradient for $h$ under the form

$$\tilde{\nabla}h(x) = \sum_{\alpha, \beta} \mathbf{1}_{R^{\alpha, \beta}}(x) \tilde{\nabla}h^{\alpha, \beta}(x) \text{ for all } x \in \mathbb{R}^{\alpha}.$$

It remains to prove the last part of the claim. One assembles $A(x)$ where

$$A(x) = \tilde{\nabla}h(x) - \tilde{\nabla}g \circ f(x) \tilde{\nabla}f(x)$$

$$= \sum_{\alpha, \beta} \mathbf{1}_{R^{\alpha, \beta}}(x) \nabla h^{\alpha, \beta}(x) - \left( \sum_{\beta} \mathbf{1}_{Q^{\beta}}(x) \nabla g^{\beta}(f(x)) \right) \left( \sum_{\alpha} \mathbf{1}_{P^{\alpha}}(x) \nabla f^{\alpha}(x) \right).$$

Take $x \in R^{\alpha, \beta}$, which by definition implies $f(x) = f^{\alpha}(x) \in Q^{\beta}$. Then

$$A(x) = \nabla h^{\alpha, \beta}(x) - \nabla g^{\beta}(f^{\alpha}(x)) \nabla f^{\alpha}(x) = 0$$

because $h^{\alpha, \beta} = g^{\beta} \circ f^{\alpha}$ by definition. Since it holds for all values of $\alpha$ and $\beta$, then $A(x) = 0$ for all $x$. $\square$

**Remark 3.** So any time one composes such functions, the number of pieces may also increase multiplicatively.
Remarks 1-3 have consequences for the theoretical understanding of Neural Networks functions trainable by Machine Learning. Indeed the basic operations of addition and composition generate a number of pieces which soon may become extremely large in the deep regime $l \gg 1$ (a similar observation is in [15], but without discussion of the differentiability properties). In this context it is needed to address the closure in a convenient space when the number of pieces $\alpha_{max}$ tends to infinity. This is the objective of the next Section.

## 3 Closure

We now come to the central theoretical contribution of this work, which is the definition of a Banach space such that one can pass to the limit of a countable number of pieces.

For a given $v \in \text{Lip}(\mathbb{R}^a : \mathbb{R}^b)$, we define

$$G(v) = \{g \in l^{\infty}(\mathbb{R}^a : \mathcal{M}_{b,a}(\mathbb{R})) \mid g = \nabla v \text{ in } L^{\infty}(\mathbb{R}^a : \mathcal{M}_{b,a}(\mathbb{R})).\}$$

It is known that a bounded Lebesgue measurable function is equal almost everywhere (up to a set of zero Lebesgue measure) to a bounded Borel measurable function. So $G(v) \neq \emptyset$. We define the functional space

$$X_{ab} = \{(v, g) \text{ where } v \in \text{Lip}(\mathbb{R}^a : \mathbb{R}^b), g \in G(v)\}. \quad (16)$$

For $v, w \in X_{ab}$, one has the additivity $G(v + w) = G(v) + G(w)$. Therefore $X_{ab}$ is a vectorial space. We equip $X_{ab}$ with the norm

$$\|(v, g)\|_{X_{ab}} = |v(0)| + \sup_{x \in \mathbb{R}^a} |g(x)|, \quad v \in X_{ab}. \quad (17)$$

Since $g = \nabla v$ in $L^{\infty}(\mathbb{R}^a : \mathcal{M}_{b,a}(\mathbb{R}))$, one has

$$\|\nabla v\|_{L^{\infty}(\mathbb{R}^a; \mathcal{M}_{b,a}(\mathbb{R}))} \le \sup_{x \in \mathbb{R}^a} |g(x)|, \quad v \in X_{ab}. \quad (18)$$

**Proposition 2.** The space $X_{ab}$ is a Banach space.

*Proof.* It is essentially immediate as a consequence of the closeness of bounded Borel functions with respect to the supremum norm, [38][Rudin, Corollary a) of Theorem 1.14] or [30][Malliavin, Theorem 2.5.1].

Take a Cauchy-sequence $(v_n, g_n) \in X_{ab}$, $n = 0, 1, \dots$. Using (18) $v_n$ is also a Cauchy sequence in $\text{Lip}(\mathbb{R}^a : \mathbb{R}^b)$ and it has a limit $v \in \text{Lip}(\mathbb{R}^a : \mathbb{R}^b)$. Also for all $x \in \mathbb{R}^a$, the sequence $n \mapsto g_n(x)$ is a Cauchy-sequence in $\mathbb{R}^b$. So it has a limit which defines a function $g(x) = \lim_{n \to \infty} g_n(x)$ for all $x$. The main point is that the pointwise limit of bounded Borel-measurable function is as well a bounded Borel-measurable function. Therefore the limit $g$ is itself a bounded Borel function. Clearly $g = \nabla v$ in the sense of distribution, and $(v, g) \in X_{ab}$. The rest of the proof is immediate. $\square$
The vectorial space of $C^1$-piecewise functions accompanied with their associated gradients is denoted as

$$ Y_{ab} = \{(v, g) \text{ where } v \in C_{\text{pw}}^1(\mathbb{R}^a : \mathbb{R}^b), g = \tilde{\nabla} v \}. \quad (19) $$

Clearly $Y_{ab} \subset X_{ab}$. The closure of $Y_{ab}$ in $X_{ab}$ is

$$ Z_{ab} = \overline{Y_{ab}} \text{ where the closure is wrt to the norm in } X_{ab}. \quad (20) $$

The set $Z_{ab}$ is a vectorial space since $Y_{a,b}$ and $X_{ab}$ are vectorial spaces.

**Definition 3.** We define the limit vectorial space of Neural Network functions trainable by Machine Learning as

$$ \overline{C_{\text{pw}}^1(\mathbb{R}^a : \mathbb{R}^b)} = \{v \in \text{Lip}(\mathbb{R}^a : \mathbb{R}^b) \mid \exists(\alpha, \beta) \in Z_{ab} \text{ such that } v = \alpha\}. $$

The associated gradient of a function $v \in \overline{C_{\text{pw}}^1(\mathbb{R}^a : \mathbb{R}^b)}$ is any matrix valued Borel function $\beta$ such that $(v, \beta) \in Z_{ab}$. We will write $\tilde{\nabla} v = \beta$.

Once again an associated gradient $\tilde{\nabla} v$ of a function $v \in \overline{C_{\text{pw}}^1(\mathbb{R}^a : \mathbb{R}^b)}$ cannot be unique. Firstly because it was already non unique in previous definition (7). Secondly because of the closure procedure (20). This closure procedure allows for limit functions which may an infinite countable number of Borel pieces.

**Remark 4.** A convenient notation would be to identify $C_{\text{pw}}^1(\mathbb{R}^a : \mathbb{R}^b)$ and its limit $\overline{C_{\text{pw}}^1(\mathbb{R}^a : \mathbb{R}^b)}$. It results in shorter notations. However for the sake of clarity we keep on using the notation $\overline{C_{\text{pw}}^1(\mathbb{R}^a : \mathbb{R}^b)}$ throughout this work. One has the embedding by construction $\overline{C_{\text{pw}}^1(\mathbb{R}^a : \mathbb{R}^b)} \subset \text{Lip}(\mathbb{R}^a : \mathbb{R}^b)$.

A remaining question is to determine if $\overline{C_{\text{pw}}^1(\mathbb{R}^a : \mathbb{R}^b)}$ is equal to $\text{Lip}(\mathbb{R}^a : \mathbb{R}^b)$. Considering the way the limit space is constructed and pathological cases of Lipschitz functions discussed in [43], in particular example 7, some care has to be taken before answering this question. This question which is not treated in this text.

The Murat-Trombetti Theorem admits a generalization in this setting.

**Theorem 2.** Consider two functions $u \in \text{Lip}(\mathbb{R}^a : \mathbb{R}^b)$ and $v \in \overline{C_{\text{pw}}^1(\mathbb{R}^a : \mathbb{R}^b)}$. Then the chain rule identity holds in $L^\infty(\mathbb{R}^a : \mathcal{M}_{c,a}(\mathbb{R}))$.

$$ \nabla(v \circ u) = \tilde{\nabla}v \circ u \nabla u \quad (21) $$

where $\tilde{\nabla}v$ is any gradient associated to $v$ and $\tilde{\nabla}v \circ u$ denotes the function such that $\tilde{\nabla}v \circ u(x) = \nabla v(u(x))$ for all $x \in \mathbb{R}^a$.

*Proof.* The function $v \circ u$ is Lipschitz continuous, so it has a bounded gradient in consequence of the Rädemacher Theorem. An associated gradient $\tilde{\nabla}v$ is Borel measurable. Since $u$ is a continuous function, then $\tilde{\nabla}v \circ u$ is Borel measurable. It is therefore Lebesgue measurable. The norm (17) controls the supremum
(with respect to the space variable) of associated gradients. Therefore $\tilde{\nabla}v \circ u \in L^\infty(\mathbb{R}^a : \mathbb{R}^c)$, so the right hand side and the left hand side of the statement (21) both belong $L^\infty(\mathbb{R}^a : \mathcal{M}_{c,a}(\mathbb{R}))$.

It remains to show that the equality holds. By construction the pair $(v, \tilde{\nabla}v)$ belongs to $Z_{ab}$ which is the closure of $Y_{ab}$. Therefore there exists a sequence $(v_n, \tilde{\nabla}v_n) \in Y_{ab}$ such that $\lim_{n \to \infty} \|(v - v_n, \tilde{\nabla}v - \tilde{\nabla}v_n)\|_{X_{ab}} = 0$. The Murat-Trombetti Theorem yields that $\nabla(v_n \circ u) = \tilde{\nabla}v_n \circ u \nabla u$ or all $n \ge 0$ and $u \in \text{Lip}(\mathbb{R}^a : \mathbb{R}^b)$. One has

$$ \nabla(v \circ u) - \tilde{\nabla}v \circ u \nabla u = A_n - B_n $$

where $A_n = \nabla(v \circ u) - \nabla(v_n \circ u)$ and $B_n = \tilde{\nabla}v \circ u \nabla u - \tilde{\nabla}v_n \circ u \nabla u$. Since $v_n$ is bounded in $W^{1,\infty}(\mathbb{R}^b : \mathbb{R}^c)$ as a Lipschitz function uniformly with respect to $n$, then $v_n \circ u$ is Lipschitz uniformly with respect to $n$. Also $v_n \to v$ is in $L^\infty(\mathbb{R}^b : \mathbb{R}^c)$ uniformly with respect to $n$. So one has

$$ A_n \to 0 \text{ weak} $$

where the weak convergence can be taken in the sense of distributions, or in the $L^\infty_{w,\star}$ sense (it has no importance). One also has

$$ \|B_n\|_{L^\infty(\mathbb{R}^a; \mathbb{R}^c)} \le \|\tilde{\nabla}v \circ u - \tilde{\nabla}v_n \circ u\|_{L^\infty(\mathbb{R}^a; \mathbb{R}^c)} \|\nabla u\|_{L^\infty(\mathbb{R}^a; \mathbb{R}^c)} \\ \le \left( \sup_{x \in \mathbb{R}^b} \left| \tilde{\nabla}v(x) - \tilde{\nabla}v_n(x) \right| \right) \|\nabla u\|_{L^\infty(\mathbb{R}^a; \mathbb{R}^c)} $$

So $\|B_n\|_{L^\infty(\mathbb{R}^a; \mathbb{R}^c)} \to 0$. Passing to the limit $n \to \infty$ (weak for $A_n \to 0$, strong for $B_n \to 0$), one obtains $\nabla(v \circ u) - \tilde{\nabla}v \circ u \nabla u = 0$ in $L^\infty(\mathbb{R}^a : \mathbb{R}^c)$ which is the claim. $\square$

**Lemma 3.** Associated gradient are bounded Borel representative of gradients in the sense of distribution.

*Proof.* Take $u(x) = x$ in (21). $\square$

# 4 Examples and final considerations

We review two cases where the notions presented before can be of some use, then we will review some natural extensions of the framework proposed in this work. Our first example details a known case of a Takagi's function applied to Neural Networks [15]. We explain why the derivative of this function calculated with Pytorch is indeed a non trivial bounded Borel function. So it makes a connection between the above theory and some practical calculations. This example corresponds to algorithmic issues which are similar to the ones encountered for in implicit autodifferentiation [27, 21, 11].
## 4.1 Takagi's functions

Long time after the seminal work of Takagi [41] in 1900, the theory of such functions has been developed in [23, 22]. See the recent review [1]. Recently in 2017 [45], Yarotsky noticed that it can be an entry to develop an approximation theory of Neural Networks with ReLU as the main activation function. Since then, the numerical analysis community started some works that establish the connection with the Finite Element theory [33, 24, 19], with special mention to [15]. More connections with representation of polynomial functions as Neural Networks can be found in [16, 17] with an extension to a bivariate polynomial in [35].

We present Takagi's functions on one example with notations directly adapted to associated gradients. The interval of study is $I = [0, 1]$. One starts with the hat function

$$g(x) = 2x \text{ for } 0 \le x \le \frac{1}{2}, \quad g(x) = 2 - 2x \text{ for } \frac{1}{2} \le x \le 1.$$

This function is continuous and piecewise affine. It can be constructed from the ReLU function $R(x) = \max(0, x)$ in at least two different ways:

$$\begin{array}{lll} \text{either} & g_1(x) = 2R(x) - 4R(x - 1/2), & x \in I, \\ \text{or} & g_2(x) = 2R(1-x) - 4R(1/2-x), & x \in I. \end{array} \qquad (22)$$

Other possibilities exist. Of course $g_1 = g_2$ are the same function for $x \in [0, 1]$. The point is that different representations for ReLU induce different associated derivatives (i.e. associated gradients with respect to one variable) for $g_1$ and $g_2$. A standard choice [7] that we adopt is $\tilde{R}'(0) = 0$ which corresponds to the representation with two pieces which are $R^1(x) = 0$ for $x \in P^1 = (-\infty, 0]$ and $R^2(x) = x$ for $x \in P^2 = (0, \infty)$. In the two pieces, the functions are just affine so very simple. Then $\tilde{g}'_1(x) = \tilde{g}'_2(x) = 2$ for $0 \le x < \frac{1}{2}$, $\tilde{g}'_1(x) = \tilde{g}'_2(x) = -2$ for $\frac{1}{2} < x \le 1$, $\tilde{g}'_1(1/2) = 2$ and $\tilde{g}'_2(1/2) = -2$. For the sake of generality, we will note

$$\tilde{g}'_1(1/2) = z_1 \text{ and } \tilde{g}'_2(1/2) = z_2$$

where $z_1, z_2 \in \mathbb{R}$ are arbitrary numbers depending on the construction (22) of $g_1$ and $g_2$. The associated derivatives are different $z_1 \neq z_2 \implies \tilde{g}'_1 \neq \tilde{g}'_2$.

Now take a number $\alpha \in \mathbb{R}$ and write a functional equation where the unknown is the function $u$

$$u(x) = g_1(x) + \alpha u \circ g_2(x), \quad x \in I. \qquad (23)$$

We use on purpose the two versions of function $g$. This functional equation is a natural way to generate functions by composition which can easily be implemented as Neural Networks [45, 15, 17]. A discussion of different possibilities in terms of $g_1, g_2$ and other parameters is in [23, 22, 1].
**Lemma 4.** Let $|\alpha| < 1$. The generating equation (23) has a unique solution $u \in C^0(I)$ expressed as a series $u(x) = \sum_{n \ge 0} \alpha^n g_1 \circ g_2^{(n)}(x)$ where the functions $g_2^{(n)}$ are constructed by iterations

$$g_2^{(0)}(x) = x, \quad g_2^{(n+1)}(x) = g_2^{(n)} \circ g(x).$$

*Proof.* Since $|\alpha| < 1$, the fixed point series converges in the maximum norm. The limit $u$ is continuous because it is a uniform limit of continuous functions. □

The converging fixed point [17] which assembles the general series is

$$\left\{ \begin{array}{l} u_0(x) = 0, \\ u_{k+1}(x) = g_1(x) + \alpha u_k(g_2(x)), \quad k \in \mathbb{N}. \end{array} \right. \qquad (24)$$

The truncated series can be written as is

$$u_n = \sum_{k=0}^{n} \alpha^k g_1 \circ \underbrace{g_2 \circ \cdots \circ g_2}_{\text{k times}}. \qquad (25)$$

The implementation of this truncated series in a Neural Network with $\approx n$ hidden layers and 3 neurons per layers is very natural and is explained in [45, 15, 33, 24, 16, 17]. We refer to functions generated by such series (24) or (25) as Takagi's functions. Takagi's functions are a simple example of Neural Network functions with an arbitrary large number of hidden layers.

Two important cases are $\alpha = \frac{1}{2}$ which is the original Takagi's function. The Takagi's function is continuous but nowhere differentiable [23, 22, 1]. A independent proof that it is not a BV function can be found in [19]. The other case is $\alpha = \frac{1}{4}$ because the solution is simple and explicit. Actually it is a polynomial function. These two cases show that Takagi's functions are able to capture in the same setting very irregular functions ($\alpha = \frac{1}{2}$) and smooth functions ($\alpha = \frac{1}{4}$).

**Lemma 5.** Take $\alpha = \frac{1}{4}$. Then $u(x) = 4x(1-x)$.

*Proof.* A proof by direct summation of the series can be found in [19]. Another simpler proof is as follows. The key observation is that the proposed solution is a polynomial of degree two so it is easy to show it is the solution of (23). Indeed one has

$$0 \le x \le \frac{1}{2}: \quad g_1(x) + \frac{1}{4}u(g_2(x)) = 2x + \frac{1}{4} \times 4 \times 2x(1-2x) = 4x - 4x^2 = u(x).$$

The proposed solution $u$ is right-left invariant $u(1-x) = u(x)$ and the functions $g_1 = g_2$ are also right-left invariant. So the equality also holds on $\frac{1}{2} \le x \le 1$. □

A standard mathematical question is to establish in what sense the derivative of the series is the series of the derivatives. To answer to this question, we
calculate the associated derivatives. Starting with a zero associated gradient, the algebra of associated derivatives deduced from (24) is

$$ \left\{ \begin{array}{l} \tilde{u}'_0(x) = 0, \\ \tilde{u}'_{n+1}(x) = \tilde{g}'_1(x) + \alpha \tilde{g}'_2(x) \tilde{u}'_n(g_2(x)), \quad n \in \mathbb{N}. \end{array} \right. \qquad (26) $$

The result at step n can be written directly as

$$ \tilde{u}'_n = \sum_{k=0}^{n} \alpha^k \left( g_1 \circ \underbrace{g_2 \circ \dots \circ g_2}_{k \text{ times}} \right)' $$

where the associated derivative of the long chain of composition are calculated in accordance of the algebra of associated derivatives. Adapting the notation (16) to the interval $I = [0, 1]$, one notes that $(u_n, \tilde{u}'_n) \in X_{11}(I)$.

**Theorem 3.** Assume $|\alpha| \max(2, |z_2|) < 1$. There exists a Borel function $\beta \in l^\infty(I)$ such that

$$ (u_n, \tilde{u}'_n) \xrightarrow{X_{11}} (u, \beta) \text{ as } n \to \infty. \qquad (27) $$

*Proof.* The convergence of the function in $W^{1,\infty}(I)$ is immediate of course. The convergence of the associated derivative is also immediate because the iterative loop is contractive in the set of Borel function equipped with the sup norm (17). The multiplicative factor is upper bounded by

$$ \sup_{0 \le x \le 1} |\alpha \tilde{g}'_2(x)| \le |\alpha| \max(2, |z_2|) < 1 $$

for our definition of the associate gradient $\tilde{g}'_2$.

<figure><img src="image_3.png" /></figure>

Next we focus on the case $\alpha = \frac{1}{4}$ and we explicitly calculate the limit associated derivative $\beta$. We will use the Borel function $\psi$

$$ \psi(x) = 0 \text{ for } x \in \left[0, \frac{1}{2}\right) \cup \left(\frac{1}{2}, 1\right], \quad \psi\left(\frac{1}{2}\right) = 1. $$

**Proposition 3.** Take $\alpha = \frac{1}{4}$. Then the limit associated derivative in (27) is the sum of two functions, where the first one is a polynomial function and the second one is a bounded Borel function in $l^\infty(\mathbb{R} : \mathbb{R})$ assembled as a series

$$ \beta(x) = u'(x) + (z_1 - z_2) \left( \psi(x) + \sum_{n \ge 1} \frac{\Pi_{k=0}^{n-1} \tilde{g}'_2(g_2^{(k)}(x))}{4^n} \psi(g_2^{(n)}(x)) \right) \quad (28) $$

Moreover the sum of series is non zero only at the countable number of points $S = \{\frac{2p+1}{2k+1} : \text{ for } k \in \mathbb{N} \text{ and } 0 \le p \le 2^k - 1\}$.
*Proof.* Write for convenience the classical derivative as $v(x) = u'(x) = 4 - 8x$. Our first objective is to determine in what sense $v$ is a fixed point of the iterative process (26). To do so, one assembles the Borel function

$$b(x) = v(x) - \tilde{g}'_1(x) - \frac{1}{4}\tilde{g}'_2(x)v(g_2(x))$$

which is just the residual of the equation which is the fixed point of the iterative process (26). The guess is equal to $v$.

For $0 \le x < \frac{1}{2}$, one has $\tilde{g}'_1(x) = \tilde{g}'_2(x) = g'(x) = 2$. So

$$b(x) = 4 - 8x - 2 - \frac{1}{4} \times 2 \times (4 - 8(2x)) = 0, \quad 0 \le x < \frac{1}{2}.$$

For $\frac{1}{2} < x \le \frac{1}{2}$, one has $\tilde{g}'_1(x) = \tilde{g}'_2(x) = g'(x) = -2$. So

$$b(x) = 4 - 8x + 2 + \frac{1}{4} \times 2 \times (4 - 8(2 - 2x)) = 0, \quad 0 \le x < \frac{1}{2}.$$

For $x = \frac{1}{2}$, one has $\tilde{g}'_1(x) = z_1$ and $\tilde{g}'_2(x) = z_2$. So $b(x) = 0 - z_1 - \frac{1}{4}z_2v(1) = z_2 - z_1$. So one has

$$b(x) = \mu\psi(x) \qquad \mu = z_2 - z_1$$

where $\psi$ is the Borel function null everywhere except at the midpoint

$$\psi(x) = 0 \text{ for } x \ne \frac{1}{2}, \quad \psi(\frac{1}{2}) = 1.$$

By subtraction, the Borel function $w_n = \tilde{u}'_n - v$ satisfies the iterative loop

$$w_{n+1}(x) = \tilde{u}'_{n+1}(x) - v(x) = \begin{aligned}[t] & \left(\tilde{g}'_1(x) + \frac{1}{4}\tilde{g}'_2(x)\tilde{u}'_n(g_2(x))\right) \\ & - \left(b(x) + \tilde{g}'_1(x) + \frac{1}{4}\tilde{g}'_2(x)v(g_2(x))\right), \end{aligned}$$

that is

$$w_{n+1}(x) = \mu\psi(x) + \frac{1}{4}\tilde{g}'_2(x)w_n(g_2(x)) \quad \mu = z_1 - z_2, \ n \ge 0. \quad (29)$$

This is an iteration loop similar to (26), but the constant term $\tilde{g}'_1$ is replaced by $\mu\psi$ which support is restricted to one point. The limit $w$ is the convergent series which is easy to identify and which does not depend on the initial condition $w_0$

$$w(x) = \mu\psi(x) + \mu \sum_{n \ge 1} \frac{\prod_{k=0}^{n-1} \tilde{g}'_2(g_2^{(k)}(x))}{4^n} \psi(g_2^{(n)}(x)). \quad (30)$$

Indeed (30) implies

$$\frac{1}{4}\tilde{g}'_2(x)w(g_2(x)) = \frac{\mu\tilde{g}'_2(x)}{4}w(g_2(x)) + \frac{\mu\tilde{g}'_2(x)}{4} \sum_{n \ge 1} \frac{\prod_{k=0}^{n-1}\tilde{g}'_2(g_2^{(k+1)}(x))}{4^n} \psi(g_2^{(n+1)}(x))$$
<figure><img src="image_4.png" /><figcaption>Figure 1: Evaluation with Pytorch of the associated derivative $\beta = \tilde{u}'$ of $u(x) = 4(x - x^2)$. The associated derivative $\beta$ is the sum of two parts as in formula (28). One part is the classical derivative $u'(x) = 4 - 8x$. The other part is a bounded Borel function in $l^\infty(\mathbb{R} : \mathbb{R})$ with a countable number of pikes. The value of the pikes is in complete accordance with formula (28).</figcaption></figure>

$$ = \mu \sum_{n \ge 1} \frac{\prod_{k=0}^{n-1} \tilde{g}'_2(g_2^{(k)}(x))}{4^n} \psi(g_2^{(n)}(x)) = w(x) - \mu\psi(x) $$

from which one deduces that $w(x) = \mu\psi(x) + \frac{1}{4}\tilde{g}'_2(x)w(g_2(x))$. Since the iterative loop (29) is contractive in $l^\infty(I)$, then $w_n \to w$ is the sup norm. So $\beta(x) = v(x) + w(x)$.

Next all functions $\psi \circ g_2^{(k)}$ (for $k \ge 0$) in (30) are non zero only at the finite number of points. The first set of points are

$$ \begin{aligned} n = 0: & \quad g(x) = \frac{1}{2} & \text{that is } x \in S_0 = \left\{ \frac{1}{2} \right\}, \\ n = 1: & \quad g(x) \in S_0 & \text{that is } x \in S_1 = \left\{ \frac{1}{4}, \frac{3}{4} \right\}, \\ n = 2: & \quad g(x) \in S_1 & \text{that is } x \in S_2 = \left\{ \frac{1}{8}, \frac{3}{8}, \frac{5}{8}, \frac{7}{8} \right\} \end{aligned} $$
By recurrence

$$ \psi \circ g_2^{(k)}(x) \neq 0 \iff x \in S_k = \left\{ \frac{2p+1}{2^{k+1}} \right\}_{0 \le p \le 2^{k-1}} \quad (31) $$

All terms in the series (30) vanish except the countable set $S = \cup S_k$. $\square$

A numerical illustration is proposed in Figure 1. It has been computed with Pytorch. We assemble the truncated series $u_n$ with (24), then we differentiate the function with the autodiff capability of Pytorch. We took $n = 9$ in this numerical test. One has to take care that all points in $S_k$ for $k \le n$ are used in the visualisation. If not the illustration misses all the pikes visible in the Figure. Clearly, the function represented in the Figure is the sum of a smooth polynomial $4 - 8x$ plus a bounded Borel function. One can check that the amplitude of the pikes is the exact value given in formula (28).

Numerical tests with different values of $\alpha \in (0, 1/2)$ show the same behavior, except that the bulk part is no more an affine polynomial but a discontinuous function. For all values of $\alpha \in (0, 1/2)$, a Borel part

$$ \beta(x) = b_0\psi(x) + \sum_{n \ge 1} b_n \prod_{k=0}^{n-1} \tilde{g}_2' \circ g_2^{(k)}(x) \psi(g_2^{(n)}(x)) \quad \text{for some numbers } b_n \in \mathbb{R} $$

is visible in the numerical tests. In this formula all associated derivatives $\tilde{g}_2' \circ g_2^{(k)}$ can be replaced by $g_2' \circ g_2^{(k)}$ because of (31).

## 4.2 ODEs

We discuss the application of the previous tools to a simple ordinary differential equation (ODE)

$$ \left\{ \begin{array}{ll} \frac{d}{dt}x(t, x_0) = f(x(t, x_0)), & t > 0, \\ x(0, x_0) = x_0, & t = 0. \end{array} \right. \qquad (32) $$

We assume that $f \in \text{Lip}(\mathbb{R})$ is a Lipschitz function so that the solution of (35) exists on a maximal interval $I = [0, T)$ with $T > 0$. One can say if $f$ is $C^1$-piecewise, that is in a space of functions slightly more regular than Lipschitz?

A simple example is $f(x) = |x|$ which can be represented with two Borel pieces. Then the solution of (35) is

$$ \left\{ \begin{array}{ll} x(t, x_0) = e^t x_0, & x_0 \ge 0, \\ x(t, x_0) = e^{-t} x_0, & x_0 < 0. \end{array} \right. \qquad (33) $$

For a given $t > 0$, the regularity is $\{x_0 \mapsto x(t, x_0)\} \in C_{\text{pw}}^1(\mathbb{R}, \mathbb{R})$. In this case the function $f$ can be represented with two Borel pieces and the solution can be represented with the same number of pieces.

It is then natural to think that other features of $C^1$-piecewise functions can be transmitted from $f$ to the solution at a given positive time $\{x_0 \mapsto x(t, x_0)\}$. An evident result is the following, where one can check that the formula (34) describes the associated derivative of the analytical solution (33).
**Lemma 6.** Take $f \in \overline{C}_{\text{pw}}^{1}(\mathbb{R}^a : \mathbb{R}^a)$. Then the solution at time $t > 0$ of the ODE (35) has the regularity $\{x_0 \mapsto x(t, x_0)\} \in \overline{C}_{\text{pw}}^{1}(\mathbb{R}^a : \mathbb{R}^a)$. It particular this function has an associated gradient under the form

$$ \tilde{\nabla}_{x_0} x(t, x_0) = e^{\int_{s=0}^{t} \tilde{\nabla}_x f(x(s, x_0)) ds}. \quad (34) $$

*Proof.* The formula is well known for classical derivatives in case $f \in C^1(\mathbb{R}, \mathbb{R})$ is a regular function. The only point to check is that $\tilde{\nabla}_{x_0}x(t, x_0)$ is Borel as a function of $x_0$. It is actually evident. Indeed $(s, x_0) \mapsto x(s, x_0)$ is continuous with respect to both variables $(s, x_0)$. Since the associated gradient $\tilde{\nabla}_x(x)$ is Borel, then the function $(s, x_0) \mapsto \tilde{\nabla}_x f(x(s, x_0))$ is Borel. Then the integral in time yields a Borel function with respect to $x_0$. Fundamentally [38] it is because the Borel measure in higher dimensions is the tensorization of the Borel measure in smaller dimensions (this is something that is wrong with the Lebesgue measure). The exponentiation is a continuous function, so the result is a Borel function. $\square$

Such formulas may be useful in the context of hybrid modeling and implicit differentiation [27, 21, 11] which is a modern use of Machine Learning coupled with the calculation of solutions of ODEs or PDEs. Here we describe the problematic on a practical scenario for the ODE (33). We assume that $f$ is the main unknown. It is a function parameterized as feed-forward Neural Network as in (1). Once the hyperparameters are chosen such as the number of hidden layers, the number of neurons per layer and the activation functions, the function $f$ depends only on the parameters $(W_r, b_r)$ for $r = 0, \dots, \ell$. The collection of parameters is written as $W = (W_\ell, b_\ell, \dots, W_0, b_0)$. We will note the function as $x \mapsto f_W(x)$ to have remembrance that it is a parametric function.

Here we use a vocabulary which comes from data sciences. Assume one observes initial solutions which are in a bounded restricted part of the whole domain, that is $x_0 \in \Omega \subset \mathbb{R}^a$. Assume one also observes final solutions $x_1 \approx x(T)$. We write $x_1 = \varphi(x_0)$ for $x_0 \in \Omega$ where the function $\varphi$ is ultimately the unknown.

A practical question in Machine Learning is to "learn" the function $\varphi$. To do so one solves the parametric ODE

$$ \left\{ \begin{array}{ll} \frac{d}{dt}x(t, x_0, W) = f_W(x(t, x_0, W)), & t > 0, \\ x(0, x_0) = x_0, & t = 0, \end{array} \right. \quad (35) $$

where the function $f_W$ is typically assembled with a compositional structure (1). The task is to find parameters $W_{\text{opt}}$ such that the final solution $x(T, x_0, W_{\text{opt}})$ is as close as possible to the observed solution $x_1 = \varphi(x_0)$.

It is usually done by assembling a cost function and minimizing it. This is what is called training [4]. For the sake of generality, we include in the cost function a regularization term with a parameter $\varepsilon > 0$

$$ J(W) = \int_{x_0 \in \Omega} (|x(T, x_0, W) - \varphi(x_0)|^2 + \varepsilon |\nabla_{x_0} x(T, x_0, W) - \nabla_{x_0} \varphi(x_0)|^2) dx_0. \quad (36) $$
The regularization term is the $H^1$ semi-norm of the error. Many other terms of various regularity can be incorporated in the functional.

In Machine Learning softwares, such a functional is necessarily assembled with an automatic evaluation of the derivative. So a more meaningful notation consists to replace $\nabla_{x_0}x(T, x_0, W)$ with the associated derivative $\tilde{\nabla}_{x_0}x(T, x_0, W)$

$$J(W) = \int_{x_0 \in \Omega} \left( |x(T, x_0, W) - \varphi(x_0)|^2 + \varepsilon |\tilde{\nabla}_{x_0}x(T, x_0, W) - \nabla_{x_0}\varphi(x_0)|^2 \right) dx_0. \quad (37)$$

We keep the standard gradient notation for $\nabla_{x_0}\varphi(x_0)$ because it depends only on the data and does not need the chain rule. On the contrary the effective calculation (34) of $\tilde{\nabla}_{x_0}x(T, x_0, W)$ involves the chain rule. Then one can invoke Lemma 3 and Lemma 6 to justify that the integral with the associated gradient is correctly evaluated, that is (37) is equal to (36).

Neural Network approximation of solutions of PDES calculated with PINNs [28, 37], DeepRitz methods [44] or Variational PINNs [25, 6] contain similar integrals. These integrals can be justified as well with the same methodology.

The correctness of the evaluation of the gradient with respect to the parameters $W$ also yields similar questions which are important for the minimization procedure (the training in Machine Learning language).

# 5 Final considerations

The remark 4 concerns the mathematical nature of the theoretical framework and has already been formulated in Section 3. Additional remarks can be made.

It is clear that the method presented in Section 2 has a non void intersection with the theory of path differentiable functions [12]. But the clear connections between the two approaches are still to be established.

One can impose some constraints in the representation formula (6). For example one can add the constraint $\|\nabla f^\alpha\|_{L^\infty(P^\alpha; \mathcal{M}_{ba})} \le \|\nabla f\|_{L^\infty(\mathbb{R}^a; \mathcal{M}_{ba})}$. Then the norm (17) satisfies $\|(v, g)\|_{X_{ab}} = |v(0)| + \|\nabla v\|_{L^\infty(\mathbb{R}^a; \mathcal{M}_{ba})}$. By density it is also satisfies for functions in $C_{pw}^1(\mathbb{R}^a; \mathbb{R}^b)$.

Smooth activation functions of class $C^1$ are easy to characterize in this framework. It is sufficient to ask that the representation formula (6) contains just one piece, that is $\alpha_{max} = 1$. Then $\alpha_{max} = 1$ is of course propagated by addition and composition, as well as the $C^1$ regularity. By density with respect to the sup-norm of the gradients, $C^1$ regularity is also preserved in the limit space.

The set $C_{pw}^1(\mathbb{R}^a; \mathbb{R}^b)$ is an algebra for addition (for all $a$ and $b$) and composition (for $b = a$). One can develop a framework such that is also an algebra for the multiplication of two functions (for all $a$ and $b$). The multiplication of functions is needed to described some applications of Neural Networks. For doing so, it is sufficient to replace everywhere Lipschitz continuity $f \in \text{Lip}(\mathbb{R}^a; \mathbb{R}^b)$ with local Lipschitz continuity $f \in \text{Lip}_{\text{loc}}(\mathbb{R}^a; \mathbb{R}^b)$. All properties are preserved up to this change.
## A Associate gradient of maxpooling

To construct an associated gradient for the maxpooling function $S_1$ (4), one can distinguish three Borel pieces which are $P^1 = \{(y_1, y_2) \in \mathbb{R}^2 : y_1 < y_2\}$, $P^2 = \{(y_1, y_2) \in \mathbb{R}^2 : y_1 > y_2\}$ and $P^3 = \{(y_1, y_2) \in \mathbb{R}^2 : y_1 = y_2\}$. Then the smooth functions are $f^1$, $f^2$ and $f^3$. Clearly $f^1(y_1, y_2) = y_1$ in $P^1$, so that $\nabla f^1(y_1, y_2) = (1, 0)$ in $P^1$. For similar reasons, $\nabla f^2(y_1, y_2) = (0, 1)$ in $P^2$.

The only issue concerns $P^3$. The construction principle (6) yields that $y_1 = y_2 = f^3(y_1, y_2)$ in $P^3$. It can be written $y = f^3(y, y)$ for all $y \in \mathbb{R}$. The function $f^3$ being continuously differentiable, one has necessarily $1 = \partial_{y_1} f^3(y, y) + \partial_{y_2} f^3(y, y)$ for all $y$, that is

$$1 = \partial_{y_1} f^3(y_1, y_2) + \partial_{y_2} f^3(y_1, y_2) \text{ on } P^3. \qquad (38)$$

The modified chain rule formula with associated gradients writes

$$1 = \tilde{\nabla} S_1(f_0(x)) W_0 \text{ for all } x \in \mathbb{R}. \qquad (39)$$

The key observations are that $W_0 = (1, 1)$ and that $\tilde{\nabla} S_1(f_0(x)) = \nabla f^3(f_0(x))$ since $f_0(x) \in P^3$. Then (39) reduces to the identity (38) which holds by definition. Now the new Borel function $\tilde{\nabla} S_1 \circ f_0$ is defined everywhere and the chain rule is rigorous in the set of Borel functions.

## References

[1] Pieter Allaart and Kiko Kawamura. The Takagi function: a survey. *Real Analysis Exchange*, 37, 10 2011.

[2] L. Ambrosio. Gradient flows in metric spaces and in the spaces of probability measures, and applications to fokker-planck equations with respect to log-concave measures. *Bollettino dell'Unione Matematica Italiana*, 1(1):223-240, 2008.

[3] Luigi Ambrosio and Gianni Dal Maso. A general chain rule for distributional derivatives. *Proceedings of the American Mathematical Society*, 108(3):691-702, 1990.

[4] Yoshua Bengio, Ian Goodfellow, and Aaron Courville. Deep learning. MIT press, 2017.

[5] Julius Berner, Dennis Elbrächter, Philipp Grohs, and Arnulf Jentzen. Towards a regularity theory for relu networks-chain rule and global error estimates. In *2019 13th International conference on Sampling Theory and Applications (SampTA)*, pages 1-5. IEEE, 2019.

[6] Stefano Berrone, Claudio Canuto, and Moreno Pintore. Variational physics informed neural networks: the role of quadratures and test functions. *Journal of Scientific Computing*, 92(3), 2022.
[7] David Bertoin, Jérôme Bolte, Sébastien Gerchinovitz, and Edouard Pauwels. Numerical influence of relu'(0) on backpropagation. *Advances in Neural Information Processing Systems*, 34:468–479, 2021.

[8] David Bertoin, Jérôme Bolte, Sébastien Gerchinovitz, and Edouard Pauwels. Erratum: Numerical influence of relu'(0) on backpropagation. Technical report, 2023.

[9] Louis Béthune. *Deep learning with Lipschitz constraints*. PhD thesis, Université de Toulouse, 2024.

[10] Mathieu Blondel and Vincent Roulet. The elements of differentiable programming, 2025.

[11] Jerome Bolte, Tam Le, Edouard Pauwels, and Tony Silveti-Falls. Nonsmooth implicit differentiation for machine learning and optimization. *Advances in neural information processing systems*, 34:13537–3549, 2021.

[12] Jérôme Bolte and Edouard Pauwels. Conservative set valued fields, automatic differentiation, stochastic gradient methods and deep learning. *Mathematical Programming*, 188:19–51, 2021.

[13] Francois Bouchut. Renormalized solutions to the Vlasov equation with coefficients of bounded variation. *Archive for Rational Mechanics and Analysis*, 157:75–90, 02 2001.

[14] Ryan Boustany. On the numerical reliability of nonsmooth autodiff: a maxpool case study. *Transactions on Machine Learning Research*, 2024.

[15] Ingrid Daubechies, Ronald A. DeVore, Simon Foucart, Boris Hanin, and Guergana Petrova. Nonlinear approximation and (deep) relu neural networks. *Constructive Approximation*, 55:127–172, 2019.

[16] B. Després and M. Ancellin. A functional equation with polynomial solutions and application to neural networks. *Comptes Rendus Mathématique*, 358(9-10):1059–1072, 2020.

[17] Bruno Després. A convergent deep learning algorithm for approximation of polynomials. *Comptes Rendus. Mathématique*, 361:1029–1040, 09 2023.

[18] Bruno Després. A functional framework for nonsmooth autodiff with max-pooling functions. *Transactions on Machine Learning Research*, 2025.

[19] Bruno Després. *Neural Networks and Numerical Analysis*. De Gruyter, 2022.

[20] Borjan Geshkovski, Cyril Letrouit, Yury Polyanskiy, and Philippe Rigollet. A mathematical perspective on transformers. *Bulletin of the AMS*, 62:427–479, 04 2025.
[21] Riccardo Grazzi, Massimiliano Pontil, and Saverio Salzo. Nonsmooth implicit differentiation: Deterministic and stochastic convergence rates. In *Proceedings of the 41st International Conf. on Machine Learning*, volume 235, pages 16250–16274. PMLR, 2024.

[22] M. Hata and M. Yamaguti. Weierstrass’s function and chaos. *Hokkaido Mathematical Journal*, 12:333–342, 1983.

[23] M. Hata and M. Yamaguti. The Takagi function and its generalization. *Japan J. AppL Math.*, 1:83–199, 1984.

[24] Juncai He, L. Li, J. Xu, and Chunyue Zheng. Relu deep neural networks and linear finite elements. *Journal of Computational Mathematics*, 38(3):502–527, 2020.

[25] Ehsan Kharazmi, Zhongqiang Zhang, and George E.M. Karniadakis. hp-vpinns: Variational physics-informed neural networks with domain decomposition. *Computer Methods in Applied Mechanics and Engineering*, 374:113547, 2021.

[26] David Kinderlehrer and Guido Stampacchia. *An introduction to variational inequalities and their applications*. SIAM, 2000.

[27] Dmitrii Kochkov et al. Neural general circulation models for weather and climate. *Nature*, 632:1–7, 07 2024.

[28] I.E. Lagaris, A. Likas, and D.I. Fotiadis. Artificial neural networks for solving ordinary and partial differential equations. *IEEE Transactions on Neural Networks*, 9(5):987–1000, 1998.

[29] Wonyeol Lee, Hangyeol Yu, Xavier Rival, and Hongseok Yang. On correctness of automatic differentiation for non-differentiable functions. In *Advances in Neural Information Processing Systems*, volume 33, pages 6719–6730, 2020.

[30] Paul Malliavin. *Integration and probability*, volume 157. Springer Science & Business Media, 2012.

[31] C. B. Morrey Jr. *Multiple integrals in the calculus of variations*. Springer, 2009.

[32] François Murat and Cristina Trombetti. A chain rule formula for the composition of a vector-valued function by a piecewise smooth function. *Bollettino dell'Unione Matematica Italiana*, 6(3):581–595, 2003.

[33] Joost A. A. Opschoor, Philipp C. Petersen, and Christoph Schwab. Deep ReLU networks and high-order finite element methods. *Anal. Appl. (Sin gap.)*, 18(5):715–770, 2020.

[34] Gabriel Peyré. Optimal and diffusion transports in machine learning, 2025.
[35] Moreno Pintore and Bruno Després. Computable lipschitz bounds for deep neural networks. Technical report, 2024.

[36] Hans Rademacher. Über partielle und totale Differenzierbarkeit von Funktionen mehrerer Variabeln und über die Transformation der Doppelintegrale. *Mathematische Annalen.*, 79 (4):340–359, 1919.

[37] M. Raissi, P. Perdikaris, and G.E. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. *Journal of Computational Physics*, 378:686–707, 2019.

[38] Walter Rudin. *Real and complex analysis*. McGraw-Hill, Inc., 1987.

[39] Sagar Sharma, Simone Sharma, and Anidhya Athaiya. Activation functions in neural networks. *Towards Data Sci*, 6(12):310–316, 2017.

[40] Guido Stampacchia. Equations elliptiques du second ordre à coefficients discontinus. *Séminaire Jean Leray*, (3):1–77, 1963.

[41] Teiji Takagi. *A simple example of the continuous function without derivative*, pages 5–6. Springer Japan, Tokyo, 1990.

[42] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and I. IPolosukhin. Attention is all you need. *CoRR*, abs/1706.03762, 2017.

[43] Xianfu Wang. *Fine and pathological properties of subdifferentials*. National Library of Canada, Ottawa, 2001.

[44] E. Weinan and T. Yu. The Deep Ritz method: A deep learning-based numerical algorithm for solving variational problems. *Comm. in Mathematics and Statistics*, 6:1–12, 2017.

[45] D. Yarotsky. Error bounds for approximations with deep ReLU networks. *Neural Networks*, 97:103–114, 2017.