Do neural networks with similar weights necessarily compute similar functions? Let $f_1(\cdot|\Theta_1)$ and $f_2(\cdot|\Theta_2)$ be two Feedforward Neural Networks with the same architectures (number of layers, width of each layer, activation function...) and parameter arrays $\Theta_1$ and $\Theta_2$ respectively. Since $f_1(\cdot|\Theta_1)$ and $f_2(\cdot|\Theta_2)$ are entirely parametrized by their respective parameters, I am interested in knowing whether a uniform bound on the distance between $\Theta_1$ and $\Theta_2$ implies a uniform bound on the distance between $f_1(\cdot|\Theta_1)$ and $f_2(\cdot|\Theta_2)$ as well. More precisely, my question can be formulated as follows :

Does $\|\Theta_1-\Theta_2\|_\infty\le\epsilon\implies\|f_1(\cdot|\Theta_1)-f_2(\cdot|\Theta_2)\|_\infty\le g(\epsilon)$ ? If this is not true in general, under what assumptions can this be verified ?

(Where $g(\epsilon)\to0$ when $\epsilon\to0$).
Intuitively, I'd expect this relation to hold as long as the activation functions are well-behaved enough, but the answer might not be so simple, as suggested by this paper by Petersen, Raslan and Voigtlaender, in which the authors show that the converse of my desired identity is false.
In fact, I would be happy with the probabilistic version of the implication, i.e.

Does $\mathbb P\left(\|\Theta_1-\Theta_2\|_\infty>\epsilon\right)\le\delta\implies\mathbb P\left(\|f_1(\cdot|\Theta_1)-f_2(\cdot|\Theta_2)\|_\infty>g(\epsilon)\right)\le h(\delta)$ ? If this is not true in general, under what assumptions can this be verified ?

(Where $g(\epsilon)\to0$ when $\epsilon\to0$ and $h(\delta)\to0$ when $\delta\to0$).
I would also be grateful for references which might be helpful in proving the desired results.

My thoughts :
Feedforward Neural Networks can be mathematically written down as
$$f(\cdot|\Theta) : x\mapsto W_L\sigma_{V_L}W_{L-1}\sigma_{V_{L-1}}\ldots W_1\sigma_{V_1}W_0 x,\quad x\in\mathbb R^d $$
Where the $W_i$ are the weight matrices, $\sigma$ are the bias activation functions between each layers and $V_i$ are biases such that for $x'\in\mathbb R^J$, $\sigma_{V_i}(x') =  (\sigma(x'^1-V_i^1,\ldots,\sigma(x'^J - V_i^J))^T$. Here $\Theta$ would be the collection of all the $W_i$ and $V_i$.
The thing is that this expression is really awkward to work with in practice, and I don't see a straightforward way to translate the bound on $\|\Theta_1-\Theta_2\|_\infty$ to a bound on $\|f_1(\cdot|\Theta_1)-f_2(\cdot|\Theta_2)\|_\infty$, so I haven't made much progress so far.

Yes, under some boundedness conditions on inputs.
The output is differentiable with respect to all the weights, because that's how we get back-propagation to work. So, the question reduces to whether the derivative is bounded.
All the common activation functions $\sigma$ are Lipschitz, so $|\sigma(a)-\sigma(b)|<C|a-b|$ for some $C$, and in fact for $C=1$.  So for any single neuron $g(z;w)$ with inputs $z$ and weights $w$, the output satisfies
$$|g(z;w_1)-g(z;w_2)|\leq |z\cdot (w_1-w_2)|$$
which is bounded if $z$ is bounded. But the bound depends on the maximum value of $z$ and on the number of input connections ($m$,say). If the difference in the weights is less than $\epsilon$, the difference in the output is less than $(\max z-\min z)\times m\times\epsilon$.
Now, each $z$ itself is the output of a previous neuron, so there is a similar bound $(\max z'-\min z')m'\epsilon$ for its change, where the $z'$ and $m'$ are the same things for that input neuron. So, you end up  with a bound that's a sum of over input connections of a product of the bound for this layer and the previous layer. More layers, more sum-and-products.
For any fixed architecture, this is all still a finite series of sums and products, so if all the data inputs are bounded and all the parameters are bounded the Lipschitz constant will still be bounded everywhere and you get a uniform bound out.  You do need boundedness for both the raw data inputs and the parameters.  For any fixed $f_1$ with finite parameter values, all the $f_2$ whose parameters are close will also have bounded parameter values, so you don't need additional assumptions on $\Theta$ if you only want a uniform bound for $f_1-f_2$ for a fixed $f_1$ and nearby $f_2$, but you do if you want it uniform for all close pairs $f_1$, $f_2$.
