Exact Upper and Lower Bounds for the Output Distribution of Neural Networks with Random Inputs
Authors: Andrey Kofnov
Daniel Kapla
Ezio Bartocci
Efstathia Bura
Venue: ICML 2025 poster (2025)
Keywords: Neural Networks
Uncertainty Propagation
Predictive Output Distribution
Guaranteed Bounds
Primary area: theory->probabilistic_methods

Abstract: We derive exact upper and lower bounds for the cumulative distribution function (cdf)  of the output of a neural network (NN) over its entire support subject to noisy (stochastic) inputs. The upper and lower bounds converge to the true cdf over its domain as the resolution increases. Our method applies to any feedforward NN using continuous monotonic piecewise twice continuously differentiable activation functions (e.g.,  ReLU, tanh and softmax) and convolutional NNs, which were beyond the scope of competing approaches. The novelty and instrumental tool of our approach is to bound general NNs with ReLU NNs. The ReLU NN-based bounds are then used to derive the upper and lower bounds of the cdf of the NN output. Experiments demonstrate that our method delivers guaranteed bounds of the predictive output distribution over its support, thus providing exact error guarantees, in contrast to competing approaches.

Exact Upper and Lower Bounds for the Output Distribution of Neural Networks with Random Inputs
Andrey Kofnov! Daniel Kapla
Ezio Bartocci Efstathia Bura
Abstract
where Y is the output and X the input (typically multivari- ate), and f is known function modeling the relationship be- tween X and Y parametrized by O. Model (T incorporates uncertainty neither in Y nor in X and NN fitting is numer- ical algorithm for minimizing a loss function. Lack un- certainty quantification, such assessment mechanisms for prediction accuracy beyond the training data, prevents neu- ral networks, despite their potential, from being deployed in safety-critical applications ranging from medical diagnostic systems (e.g., (Hafiz & Bhat 2020)) to cyber-physical sys- tems such as autonomous vehicles, robots or drones (e.g.. (Yurtsever et al., (0200). Also, the deterministic nature of NÑs renders them highly vulnerable to not only adversarial noise but also to even small perturbations in inputs ((ibbi et al.
We derive exact upper and lower bounds for the cumulative distribution function (edf) of the out- put of a neural network (NN) over its entire sup- port subject to noisy (stochastic) inputs. The up- per and lower bounds converge to the true odf over its domain as the resolution increases. Our method applies to any feedforward NN using con- tinuous monotonic piecewise twice continuously differentiable activation functions (e.g., ReLU, tanh and softmax) and convolutional NNs, which were beyond the scope of competing approaches. The novelty and instrumental tool of our approach is to bound general NNs with ReLU NNs. The ReLU NN-based bounds are then used derive the upper and lower bounds of the edf of the NN output. Experiments demonstrate that method delivers guaranteed bounds of the predictive out- put distribution over its support, thus providing exact error guarantees, in contrast to competing approaches.
2018 Fawzi et al. 2018, Goodfellow et al. 2014
2015
Hosseini et al. 2017))
Uncertainty in modeling is typically classified as epistemic or systematic, which derives from lack of knowledge of the model, and aleatoric or statistical, which reflects the inherent randomness in the underlying process being mod- eled (see, e.g., (Hällermeier & Waegeman) 2021)). The universal approximation theorem (UAT) (Cybenko) 1989 Hornik et al. 1989) states that NN with one hidden layer can approximate any continuous function for inputs within a specific range by increasing the number of neurons. In the context of NNs, epistemic uncertainty is of secondary importance to aleatoric uncertainty. Herein, we focus on studying the effect of random inputs on the output distribu- tion of NNs and derive uniform upper and lower bounds for the cdf of the outputs of NN subject to noisy (stochastic) input data.
1. Introduction
Increased computational power, availability of large datasets, and the rapid development of new NN architec- tures contribute to the ongoing success of NN based learning in image recognition, natural language processing, speech recognition, robotics, strategic games, etc. A limitation of NN machine learning (ML) approaches is that they lack built-in mechanism assess the uncertainty or trustwor- thiness of their predictions, especially on unseen out-of- distribution data. A NN is model of the form:
We evaluate our proposed framework on four benchmark datasets (Iris (Fisher 1936). Wine (Aeberhard & Forina 1992). Diabetes (Efron et al. 2004). and Banana (Jacchan daran 2017)) and demonstrate the efficacy of our approach to bound the cdf of the NN output subject to Gaussian and Gaussian mixture inputs. We demonstrate that our bounds cover the true underlying edf over its entire support. In contrast, the similar but approximate approach of Krapf) et al. (2024), as well as high-sample Monte-Carlo simula- tions, produce estimates outside the bounds over areas of the output range where the bounds are tight.
y f(x +),
'Faculty of Mathematics and Geoinformation, TU Wien, Vienna, Austria 'Faculty of Informatics. TU Wien, Vi- enna, Austria. Correspondence to: Andrey Kofnov <an- drey.kofnov tuwicn.ac.at>.
Proceedings of the 42nd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Exact Upper and Lower Bounds for the Output Distribution of NNs with Random Inputs
2. Statement of the Problem
Our contributions:
A NN is a mathematical model that produces outputs from inputs. The input is typically vector of predictor variables, X R", and the output Y, is univariate or multivariate, continuous or categorical.
1.
We develop a method to compute the exact cdf of the output of ReLU NNs with random input pdf, which is a piecewise polynomial over a compact hyperrectan- gle. This result, which can be viewed as a stochastic analog to the Stone- Weierstrass theorem I significantly contributes to the characterization of the distribution of the output of NNs with piecewise- heeeerse-llnner activation functions under any input continuous pdf.
A seedforward NN with I, layers from R + R is a composition of L functions,
fx () = f(l) o o o (x),
(2)
where the 1-th layer is given by
2.
We derive guaranteed upper and lower bounds of the NN output distribution resulting from random input perturbations on a fixed support. This provides exact upper and lower bounds for the output edf provided the input values fall within the specified support. No prior knowledge about the true output edf is required to guarantee the validity of our bounds.
/ (x W(), (1) ()) (1) (W x + +b(),,
with weights W R bias terms ( E R and non-constant, continuous activation function o! R - R that is applied component-wise. The NN parameters are collected in e- = [vec(W; W1), bi...., vec(WL). L). by) E (0
The first layer that receives the input
3. We show the convergence of our bounds to the true edf; that is, our bounds can be refined to arbitrary accuracy.
is called the input layer, and the last layer is the outpur layer. All other layers are called hidden. For categorical outputs, the class label is assigned by a final application of a decision function, such as arg max.
4. We provide a constructive proof that any feedforward NN with continuous monotonic piecewise swice con- tinuously differentiable activation functions can be approximated from above and from below by a fully connected ReLU network, achieving any desired level of accuracy. Moreover, we enable the incorporation of multivariate operations such as max, product and softmax, as well as some non-monotonic functions such as |ar| and x", n E N.
Despite not being typically acknowledged, the training data in NNs are drawn from larger populations, and hence they contain only limited information about the corresponding population. We incorporate the uncertainty associated with the observed data assuming that they are random draws from an unknown distribution of bounded support. That is, the data are comprised of m draws from the joint distribution of (X, Y), and the network is trained on observed (x,, y:). x; = (Ira1, F22. (Fnng)) and ys. i 1.... .m 3 A NN with I. layers and ny neurons each layer, 1 ..... T., is trained on the observed (x,, M.). i 1,..., m, m,, to pro- duce m outputs Dt and the vector of the NN parameters, C = (vec( W1). bi...., vee( W by). is obtained. θ uniquely identifies the trained NN. Given o we aim to quantify the robustness of the corresponding NN, to pertur- bations in the input variables. For this, we let
5.
We prove a new universal distribution approximation theorem (UDAT), which states that we can estimate the cdf of the output of any continuous function of a random variable (or vector) that has a continuous distribution supported on compact hyperrectangle, achieving any desired level of accuracy.
3. Our Approximation Approach
We aim to estimate the cdf (v) of the output = 3) of the NN in @ under (3) i.e., subject to random per- turbations of the input X. We do so by computing upper and lower bounds of Fe: that is, we compute F Fz. E' such that
Fx,
(3)
where X E Ro stands for the randomly perturbed input variables with edf Fx and probability density function (pdf) o(x) that is piecewise continuous and bounded on a compact support. We study the propagation of uncertainty (effect of the random perturbation) in the NN by deriving upper and lower bounds of the edf F (v) = P(Y < y) of the random output, Y = fx }
Fv(v) S F+(1) S FT((), Vy
(4)
We refer the NN in 2 as prediction NN when needed for clarity. We estimate the functions F Fy F' on a "super- set" of the output domain of the prediction NN via an
'The operation vec R 1 stacks the columns of a matrix one after another.
'A significant corollary to Stone Weierstrass theorem that any continuous function defined on a compact set can be uniformly approximated as closely as desired by a polynomial. 'A wide class of the most common continuous activation func- tions, including ReLUU, tanh and logistic function.
"e use the convention of denoting random quantities with capital letters and their realizations (ohserved) by lowercase letters.
'The notation fz. (X +) signifies that Q, equivalently the NN, is fixed and only X varies.
2

---

Exact Upper and Lower Bounds for the Output Distribution of NNs with Random Inputs
Definition 3.6 (Piecewise polynomial). A function p:: K + R is piecewise polynomial, if there exists finite set of nn-simplices such that K = Ü k and the function p 2-1
16
Exact cdf Empiriaa adf
0.6
constrained to the interior k? of k2 a polynomial; that is, pla. : kg + IR is polynomial for all 1..... .-
02
01
Remark 3.7. We do not require piecewise polynomials to be continuous everywhere on the hyperrectangle. Specifi- cally, we allow discontinuities the borders of simplices. However, the existence of left and right limits of the func- tion at every point the bounded support is guaranteed by properties of polynomials.
-123
-015
-233
4010
Outcomn of the lirst element of the NN
Figure 1: Exact edf of the ReLU NN outcome for class Setosa in
the Iris data, assuming Beta-distributed inputs.
integration procedure. The edf of Y is given by
RRghu et al. 2017) showed that ReLU deep networks di- vide the input domain into activation patterns (see (Sud) jianto et al. 2020)) that disjoint convex polytopes (P,} over which the output function is locally represented as the affine transformation (x.(x) NN (x) = e² + V for x {Pj}. the number of which grows at the order )m.o)) (Sudjianto et al. 2020) outline an algorithm for extracting the full set of polytopes and determining local affine transformations, including the co- efficients e² V for all (Py}. by propagating through the layers. For our computations, we utilize recent GPU- accelerated algorithm from (Berzins) 2023).
Pv())= = P(Ý u) = 1 (8su) d(x)dxx
(5)
where p(x) is the pdf of X. To bound Fv. we bound ø by its upper σ and lower estimates the bounded support of お as described in Section 3.3 If Φ is a piecewise polynomial, then 5 can be computed exactly for a ReLU prediction Once o and o are
network, as we show in Section 3.1 estimated, then
(PSW Ev(1) S o(x)do S PP(1)) s
(6)
We aim to derive a superset of the range of the network output. For this, we exploit the technique of Interval Bound Propagation (IBP), based on ideas from (Gowal et al. 2019 Wang et al., 2022 Gehr et al., 2018). Propagating the no- dimensional box through the network leads to the superset of the range of the network output. We compute the edf of the network's output at each point of a grid of the superset of the output range.
1 o(x)dr= = Fxx)d==F&(0
Remark 3.1. Ev( and (v) are not always true cdfs since we allow the lower estimator not to achieve I, while the upper bound is allowed take the smallest value greater than 0.
Theorem 3.8 (Exact edf of ReLU NN w.I.t. piecewise polynomial pdf). Let Y R R be a feedforward ReLU NN, which splits the input space into set of almost disjoint polytopes (P,) yay with local affine transformations = = ( for x E Pa. Let d(()) denote the pdf of the random vector X that is a piecewise polynomial with local polynomials, Φ(x) = dx(xx) for all E & over an almost disjoint set of simplices
3.1. Exact edf evaluation for fully connected NN with ReLU activation function
Definition 3.2 (Almost disjoint sets). We say that sets A and B are almost disjoint with respect to measure a, if a(An B3- = O.
Definition 3.3 (Closed halfspace). A no-dimenniional closed halfspace is set H fx E I2"0 [vTT S c} for E R and some v E IR called the normal of the halfspace.
9e and a compact
hyperrectangle support K C TR.o, Then, the cdf of at point y E Ra.. is
It is known that convex polytope can be represented as intersection of halfspaces, called H H-representation (Ziegler 1995).
qa qy Σ
py = Y
Y Sy
I
Definition 3.4 (H-polyyope). A no-dimensiooal H- polytope p n H; is the intersection of finitely many
2o qy S
I (x); Tr..
closed halfspaces.
Definition 3.5 (Simplex). A no -i-dimensional simplex is a no dimeinensional polytope with no + 1 vertices.
where I (x) is the integral of the polynomial (x(x) over the simplex Teg. such that the reduced poly-
3

---

Exact Upper and Lower Bounds for the Output Distribution of NNs with Random Inputs
lope
The proof is presented in Appendix A.2
S. P :Nnnn(xNNNN( (x) sy) sy)= U Tajs
Theorem cannot be directly applied in practice. We develop an approach to approximate the activation functions of a NN provided they are non non-decreasing piecewise twice continuously differentiable.
is defined by the intersection of polytopes Pj and ka and the intersection of halfspaces
Definition 3.12. Let Ω = |a, a| C R be closed interval and g 52 - R is well-defined, continuous on Q. We will say that is piecewise twice continuously differentiable on Ω2 that is E p.u. (2), if there exists a finite partition of Ω into closed subintervals U [a,, a[+1] = S where
RL
fx: N.:NN' (x) S y)
x: (x) - ve
a = #1 < <a² < ... < an+1 7 and E N, such that [11.0.11]
Theorem 3.8 is shown in Appendix A.1
The proof relies on the algorithm for evaluating the integral of a polynomial over a simplex, as described in (Lasserre 202I). The right- hand side of (7) results from the Delaunay triangulation (Delaunay 1934), dividing the reduced polytope Pj,, into
twice continuously differentiable.
Our method provides a constructive proof for the restricted case of Theorem 3.11 for non-decreasing activation func- tions in CP p.u. (S2, at each node i of NN, where Ω2 is the input domain of node i, and applies to any fully connected or convolutional (CNN) feedforward NN with such activation functions analyzed on a hyperrectangle. The key features of our approach are:
Si, almost disjoint simplices (Ti,,=))
Remark 3.9. Theorem 3.8 is a tool for approximating the output edf of any feedforward NN with piecewise linear activation functions on compact domain, given random inputs of arbitrary continuous distribution any desired degree of accuracy.
Local adaptability: The algorithm adapts to the cur- vature of the activation function, providing an adaptive approximation scheme depending on whether the func- tion is locally convex or concave.
Example 3.10. We compute the output cdf of 3-layer, 12-neuron fully connected ReLU NN with the last (be- fore softmax) linear 3-neuron layer trained on the Iris dataset (Fisher, 1936). The Iris dataset consists of 150 samples of iris flowers from three different species: Setosa, Versicolor, and Virginica. Each sample includes four fea- tures: Sepal Length, Sepal Width, Petal length, and Petal width. We focus on two features, Sepal Length and Sepal Width (scaled to [0, 1), to classify flowers into three classes. Specifically, we recover the distribution of the first compo- nent (class Setosa) before applying the softmax function, assuming Beta-distributed inputs with parameters (2, 2) and (3. 2). The exact edf is plotted in purple in Figure [ with additional details provided in Appendix D The agreement with the empirical cdf is almost perfect.
Streamlining: By approximating the network with piecewise linear functions, the complexity of analyzing the network output is significantly reduced.
How it works:
Input/Output range evaluation: Using IBP (Gowal et al., 2019), we compute supersets of the input and output ranges of the activation function for every neuron and every layer.
Segment Splitting: First, input intervals are divided into macro-areas based on inflection points (different curvature areas) and points of discontinuity in the first second derivative (e.g., 0 for ReLU). Next, these macro-areas are subdivided into intervals based on user-specified points or their predefined number within each range. The algorithm utilizes knowledge about the behavior of the activation func- tion and differentiates between concave and convex regions of the activation function, which impacts how the approx- imations are constructed and how choose the points of segment splitting. A user defines the number of splitting segments and the algorithm ensures the resulting disjoint sub-intervals are properly ordered and on each sub-interval the function is either concave or convex. If the function is linear in a given area, it remains unchanged, with the upper and lower approximations equal to the function itself.
3.2. Algorithm for Upper and Lower Approximation of the Neural Network using ReLU activation functions.
We start by showing a general result in Theorem 3.11
Theorem 3.11. Let Y be a feedforward NN with L lay- ers of arbitrary width with continuous activation functions. There exist sequences fully connected ReLU NNs {Y..), v.). (Y,)) which are monotonically decreasing and increasing, respectively, such that for any > O and any compact hy- perrectangle K C Ro there exists N E N such that for all n 2 N
0 Y (x) YX, <E, 0 O - F(x) < for all x e K.
4

---

Exact Upper and Lower Bounds for the Output Distribution of NNs with Random Inputs
The method con- structs tighter upper and lower bounds through the specific choice of points and subsequent linear interpolation. It cal- culates new points within each interval (one per interval) and uses them to refine the approximation, ensuring the lin- car segments closely follow the curvature of the activation function.
Upper and Lower Approximations:
For monotonically increasing functions, this procedure guar- antees that the constructed approximators exact upper and lower approximations. Moreover, decreasing the step size (increasing the number of segments) reduces the error at each point, meaning that the sequences of approxima- tors for the activation functions at each node are monotonic: f S (1), Jn+1 (x) 2 [x)),. for all æ in the IBP domain. This procedure of piecewise linear approximation of each activation function at each neuron is equivalent to forming a one-layer ReLU approximation network for the given neuron.
Our method differentiates between upper and lower ap- proximations over concave and convex segments. For up- per (lower) approximations on convex (concave) segments [ax, ak+1]. we employ local linear interpolation by intro- ducing an intermediate point ax between the segment end- points. That is, we let
By the UAT, for any continuous activation function Iti at each neuron and any positive 61,; we can always find a ReLU network-approximato N NN,,, such that [N [x) 01,((z)l < C1, for all x in the input domain, defined by the IBP. To find such an approximating network we need to choose the corresponding number of splitting segments of the IBP input region. Uniform convergence is preserved by Dini's theorem, which states that a monotonic sequence of continuous functions that converges pointwise on compact domain to continuous function also converges uniformly. Moreover, the approximator always stays within the range of the limit function, ensuring that the domain in the next layer remains unchanged and preserves its uniform convergence.
iin. ag dg ag E <lax, ak+1]
J (x) for x E [ag; dx+1]
Conversely, for upper (lower) approximations on concave (convex) segments [ak: (I+1]. we construct a piecewise tan- gent approximation by inserting a linking point ag between the function tangent at the segment boundaries. That is,
aw a pieston aK E [ax, ak++1] [(x)= joie tan (x) for x E [ux, a841]
The method guarantees that the piecewise linear approxi- mation of the activation function for each neuron (a) is a non-decreasing function, and (b) the output domain remains the same.
Algorithm Piecewise Linear U/L Approximation
1: INPUT: Interval bounds (a., T.y) and activation func- tion f for neuron 1) OUTPUT: Sets {lax, 11+1]), (/)(())), {[l(())
To see this, consider layer neuron. For upper (lower) ap- proximation on a convex (concave) segment, we choose midpoint ax' = (ax + ax+1)/2 for each subinterval (ax, +1+1] and compute a linear interpolation, as follows:
3: PROCEDURE:
4: Split [a., a.p] into segments {lax. ax+1]} such that the A! [ak, ax+1] + I with (1) f() is either
linear or twice continuously differentiable with constant
sign of flK
R1 = [(xx) - -(ax)
R2 = Man+1) 1)- - [(ax)
for each segment [ag, ax+1] do
ak' 0kk
ak+1 02*'
6:
if fle is linear then
(()) = f(ax) + (a-a))m11
T E (ak: avl
7:
flx(x) fle(x) = fle(x)
j(r) = s(ax) + (T an) )x2)
E [aw, 01++1]
8:
else if f 2 0 then
Flx jlsm (xx)
For upper (lower) approximation on a concave (convex) segment, compute derivatives and look for tangent lines at border points of the aub-interval [as, dx+1]. We choose point as' 04 I ag* I a to be the intersection of the tangent lines. The original function is approximated by the following two tangent line segments:
10:
flx = fpin_tam la (z)
11:
else if f S then
12:
Tl fpie_tan (x)
13:
flx(x) flim | ((())
14: end if
15:
end for
f(as) - (ax+1)) - (f4(un (ax))a )ax - 1alla=---
av-
f_(44)1)= fz (ax)
Transformation into ReLU-Equivalent Form: The ap- proximating piecewise linear upper-lower functions are con- verted into form that mimics the ReLU function's 's behavior; i.e., f(x) = (2) ()LU(W(()) ((W())++bl) +)())++()+++)))))))))))++)) +b(2) This involves creating new weighting coefficients and intercepts that repli- cate the ReL.U's activation pattern across the approximation
f(r) = f(ax) + (T=) - ax),
TE [ak: arl ī(T) f(ax+1) + 1'' (ax+1)(+ -x++), (ax, ax+1);
where f-,,,, are left and right derivatives, respec- tively.
5

---

Exact Upper and Lower Bounds for the Output Distribution of NNs with Random Inputs
propagation of approximators, to create an upper approxi- mation, we connect the upper approximation of the external layer with the upper approximation of the internal subnet- work if the internal subnetwork has a positive coefficient, or with the lower approximation if it has a negatíve coefficient. The reverse applies to the lower approximation. Detailed ex- planation provided Lemma B. [nd
1.5
andano
0.5
and B. B.12]j in Appendix B.1 This leads to a composition of uniformly convergent sequences, guaranteeing the overall uniform convergence of the final estimator to the original neural network, we show in Theorem 3.14 The final output is a set of piecewise linear approximations that bound the output of the original neural network, which can then be used for further analysis or verification.
20
0.0
Feature 2 (Sepal Width)
0.4
0.6
0.8
1.0
Feature 1 (Sepal Length)
0.8
0.4
Example 3.13. Using the same setup as in Example 3.10 we train a fully connected neural network with 3 layers of 12 neurons each with tanh activation followed by 1 layer of output neurons with linear activation on the Iris dataset, focusing on the re-scaled features Sepal Length and Sepal Width. We construct upper and lower approximations of the network's output by ReLU neural networks with linear output layer. Two approximations are performed: with 5 and 10 bounding segments per convex section at each node (upper and lower panels of Figure 2 respectively). Notably, with 10 segments, the original network and its approximations are nearly indistinguishable.
1.5
3.0
3nogno
0.5
Feature 2 (Sepal Width)
0.2 0.4
0.6
0.8
1.0
1.0
Feature 1 (Sepal Length)
0.8
0.6
0.4
0.2
0.0
Figure 2: Tonk NN output for the Setosa class (blue) in the Iris dataset, and its ReLU NN upper (red) and lower (green) approx- imations with (upper panel) and 10 (lower panel) segments of bounding per convex section.
This procedure applies to various non mon-monttoiic functions, such as monomials x², n E N. These functions can be attained by a sequence transformations that ensure mono- tonicity all intermediate steps. These transformations can be represented as a subnetwork. Furthermore, multivari- ate functions like softmax and the product operation also have equivalent subnetworks with continuous monotonic transformations, as shown in Appendix C
intervals.
Specifically, we receive a set of intervals {[ aml with the corresponding set of parameters of affine transfor- mations {(v,, ::)
The main result of this section is Theorem 3. 14 that shows that the upper and lower ReLU bounds converge uniformly and monotonically to the target NN.
f(r) + 0T, T T
Theorem 3.14 (Uniform Monotonic Convergence of the ReLU Bounds) The sequences of estimating functions that are ReLU NNs generated by the method in Section 3.2 es- tablish upper and lower bounds for the target NN. These se- quences are monotonically decreasing for the upper bounds and monotonically increasing for the lower bounds, and converge uniformly to the target network.
and set 0g- 0. Then, the corresponding ReLU ReLU-equivvlent definition of approximation ff on [zo, xn] is
F(T) = un +c + Σ RReLU(( - -u-1((T x)),)),
& = signnn - 1)-
The proof is given in Appendix B.2 steps, each outlined in different lemmas in Appendix B.I
It involves multiple
Full Neural Network Approximation: The entire NN is approximated by applying the above techniques to each neuron, layer by layer, and then merging all intermediate weights and biases. For each neuron, both upper and lower approximations are generated, capturing the range of pos- sible outputs under different inputs. To ensure the correct
Remark 3.15. Theorem 3.11 follows from the UAT (Hornik et al. 1989). Theorem 3.14 is a special case of Theorem 3.11 but our proof is constructive and explicitly outlines the sequences that establish the bounds.
6

---

Exact Upper and Lower Bounds for the Output Distribution of NNs with Random Inputs
3.3. Application to an arbitrary function on compact domain
In all experiments, we compare our guaranteed bounds for the output edf with edf estimates obtained via a Monte Carlo (MC) simulation (100 million samples), and with the "Piecewise Linear Transformation" (PLT) method of Krapi et al. (2024) following their default subdivision setup. For each simplex an edgewise subdivision is performed with the number of subdivisions depending on the input dimension ng. For input dimension no from subdivisions was set to 250, 100 and 30, respectively. In contrast, our implementation uses a bound (arbitrarily set to 50 000) on the total number of vertices based on the iteratively refined triangulation.
We present the universal distribution approximation theo- rem, which may serve as a starting point for further research in the field of stochastic behavior of functions and the NNs that describe them.
Theorem 3.16 (Universal distribution approximation theo-
3. the number of
rem). Let X be a random vector with continuous pdf 6(x) supported over a compact hyperrectangle K C IR", Let Y=W(x) E R be a continuous function X with domain K, and let F(y) denote its cdf. Then, there exist sequences of cdf bounds {Fa). {Fn). ,. n 1,2,, ..., which can be constructed by bounding the distributions of sequences of ReLU NNs and such that
Since the true edf is contained within the limits computed by our method, the experiments assess the tightness of our bounds compared both MC and PLT by tallying the num- ber of "out of bounds" instances. Essentially, achieving very tight bounds makes it challenging to stay within those limits, whereas even imprecise estimates can fall within broader bounds.
En(() S F()) S (v)))
for all y E (W(x): x E }} and
En(y) + Fly), Fv( - Fly)
Our experimental setup is based on small pre-trained (fixed) neural networks. For the Banana data, we trained a ReLU network with 2 hidden layers of 16 neurons each. For the Diabetes dataset, we trained a ReLU network with 3 hidden layers with 32, 16, and 8 neurons, respectively. The uni- variate output has no activation. Training was performed on 70% of the data, which were randomly selected. The re- maining 30% comprise the test set. We added normal noise to the or randomly selected features (see no in Table [ of every observation with variance the sample variance of the selected feature(s) in the test set (different features in dif- ferent observations). This is simplified simulation setting of that of (Krapf et al. 2024).
for all y where Fly)+ is continuous. Moreover, if WX) is almost surely nowhere locally constant, that is
o(x)dx=0 =0
(8)
[WiX)=y)
for all y E (W(x): x E K}, then both bounds Eng F.. converge uniformly to the true edf F.
The proof is presented in Appendix A.3 It leverages the UAT (Hornik et al. 1989) to approximate the function W(x) and input pdf (xx) with ReLU NNs to arbitrary accuracy. The edf bounds are then computed over polytope intersec- tions. Increasing the NN complexity results in more sim-
For both Iris and Wine, we replicated the exact experimen- tal setup from (Krapf et al. 2024) using the same test sets, Gaussian mixtures as randomness as well as the same pre- trained networks The 1 to 3 dimensional Gaussian mix- tures were computed by first deleting 25% or 50% of the test data. Subsequently, 50 new observations are imputed using MICE (van Buuren & Sonohhui-OOdssoom) 201I)) The new imputed observations where used as the dataset for a Gaussian kernel density estimate providing the Gaus- sian mixture for each of the original deleted observations. The only difference in the experimental setup is the MC estimate, which we recomputed with a higher sample count (108 as opposed to 8. 10a till 1.024 - 107 with a median of 6.4 - 10 deduced by an incremental doubled sample size convergence criteria). Here, the "ground truth" are the exact cdf bounds our method computes as opposed to the MC estimate in the (Krapf et al. 2024) experiments, which was
plices and thus leads to finer local affine approximations of the pdf.
To bound the edf of given NN with respect to a specified input pdf, we construct upper and lower bounding ReLU NNs to approximate the target NN. Next, we subdivide the resulting polytopes of the bounding ReLU NNs into simplices as much as needed achieve the desired accuracy and locally bound the input pdf with constant values (the simplest polynomial form) these simplices, transforming the problem into the one described in Section 3.1
4. Experiments
We consider four datasets, Banana, Diabetes (Efron et al. 2004 Iris (Fisher 1936). and Wine (Aeberhard & Forina 1992
wttpuuu.ccommmtt cOm/URNI2/ / Pi Siecewise-Lieend Srtevssseeiieeer-Taataattttttto Transformation] accessed Jan 25, 2025
"GitHub page
Tris, Wine and Diahetes are provided in the Python package scik boikit kit weit-lear 1ll-earr Tearn rn and Banana Kaggle
7

---

Exact Upper and Lower Bounds for the Output Distribution of NNs with Random Inputs
an approximation to the true edf.
robustness certificates for neural networks under Gaussian and Sub- Gabb-Ga input perturbations with bounded support with a given probability. It employs CROWN (Zhang et al. 2018 2020) to compute deterministic affine bounds and subsequently leverages straightforward probabilistic tech- niques based on Hoeffding's 's inequality (Hoetiding) [963). PROVEN provides a probabilistically sound solution to en- suring the output of NN is the same for small input per- turbations with given probability, its effectiveness hinges on the activation functions used. It cannot refine bounds or handle various input distributions, which may limit its ability to capture all adversarial attacks or perturbations in practical scenarios.
We summarize our results in Table When the input dimen- sion is I, our bounds are very tight and we thus observe a high ratio of "out of bounds" samples for both MC and PLT compared to the number of grid points the edf was evaluated at. This has two components: (a) In regions where the cdf is very flat, we obtain very tight bounds leading to small errors in bucketed estimation approach easily falling outside of these tight bounds; (b) due either pure random effect in the case of MC or numerical estimation inaccuracies in case of PLT, MC and PLT estimates are outside the bounds. Al- though the PLT estimator directly targets the pdf and would be expected to achieve greater precision than ours (as we target bounding the edf instead), we note that in these exam- ples, especially as regards PLT, "coarse grid" can cause inaccuracies in areas where the pdf fluctuates significantly. Nevertheless, this behavior diminishes as the input dimen- sion increases. Due to our hard bound on the maximum number of vertices in this simulation (namely 50 000), our estimated edf bounds are wider in higher dimensions as a consequence of the curse of dimensionality.
The most relevant published work to ours we could find in the literature is Krapf et al. (2024). They propagate in- put densities through NNs with piecewise linear activations like ReL.U, without needing sampling or specific assump- tions beyond bounded support. Their method calculates the propagated pdf in the output space using the piecewise linearity of the ReLU activation. (Krapf et al. 2024) esti- mate the output pdf and show experimentally that is very close to the Monte Carlo based pdf. Despite its originality, the approach has drawbacks, as they compare histograms rather than the actual pdfs in their experiments. Theorem 5 (App. C) in (Krapf et al.) (024) suggests approximating the distribution with fine bin grids and input subdivisions, but this is difficult to implement in practice. Without knowl- edge of the actual distribution, it is challenging to define a sufficiently "fine" grid. In contrast, we compute exact bounds of the true output edf over its entire support (at any point, grid required) that depicts the maximum error over its support, and show convergence to the true cdf. Krapp ct al.] (2024) use a piecewise constant approximation for input pdfs, which they motivate by their Lemma 3 (App. C) to deduce that exact propagation of piecewise polynomi- als through neural network is infeasible. We show that, in effect, it is possible and provide method for exact in- tegration over polytopes. Additionally, their approach is limited to networks with piecewise linear activations, ex- cluding locally nonlincar functions. In contrast, our method adapts to NNs with continuous, monotonic piecewise twice continuously differentiable activations.
5. Related Work
The literature on NN verification is not directly related to ours as it has been devoted to standard non-stochastic input NNs, where the focus is on establishing guarantees of local robustness. This line of work develops testing algorithms for whether the output of NN stays the same within a specified neighborhood of the deterministic input (see, e.g.,
Gowal et al. (2019); Xu et al. 2020); Zhang al. (2018
2020); Shi et al. (2025) Bunel et al. (2020); Ferrari et al. (2022): Katz et al. (2017: 2019): Wu et al. (2024)).
To handle noisy data or aleatoric uncertainty (random in- put) in NNs, two main approaches have been proposed: sampling- based and probability density function (pdf) approximation-busee Sampling-based methods use Monte Carlo simulations to propagate random samples through the NN (see, e.g.. Abdelaziz et al. (2015): Ji et al. (2020). but the required replications to achieve similar accuracy to theo- retical approaches such as ours, as can be seen in Table [ can be massive. Pdf approximation- based methods assume specific distributions for inputs or hidden layers, such as Gaussian (Abdelaziz al., 2015) or Gaussian Mixture Mod- els (Zhagg& & Shin 2021), but these methods often suffer from significant approximation errors and fail to accurately quantify predictive uncertainty. Comprehensive summaries and reviews of these approaches can be found in sources like Siiciing et al. (2022) and Gawlikowski et al. 20033). )
6. Conclusion
We develop novel method to analyze the probabilistic behavior of the output of neural network subject to noisy (stochastic) inputs. We formulate an algorithm to compute bounds (upper and lower) for the edf of a neural network's 's output and prove that the bounds are guaranteed and that they converge uniformly to the true edf.
In the context of verifying neural network properties within probabilistic framework, (Weng et al.] 2019) proposed PROVEN, general probabilistic framework that ensures
Our approach enhances deterministic local robustness ver- ification using non-random function approximation. By
8

---

Exact Upper and Lower Bounds for the Output Distribution of NNs with Random Inputs
Table Comparison our approach U/L-Dist (guaranteed upper and lower bounds) with pointwise estimators from Monte-Carlo simulations and PLT (Krapf et al. 2024). Under column headed by no are the numbers of imput variables, # Tests is test set size grouped by under U/L-dist are the mean distance (standard deviation) between our upper and lower bounds, under OOBMC and OOBrLT are the minimum, median and maximum number of points outside our hounds relative to Grid Size many odf evaluation points for Monte-Carlo simulations and PI.T, respectively. Runtime gives the computation time for PLT and our approach in seconds (MC took ahout 20 all cases).
Dataset
Y4g
Tests
U/L-Dist
Out of Bounds
Grid
Runtime [sec]
MC
PLT
Size
PLT
OUR
Menn (Sud)
MelnMMeeinnMMaa
(Min
Mean
Mean
banana
1591
0.0001 (0.0000)
0/ 103 /
896
77 565
999
1000
1.4
1.0
diabetes
133
0.0001 (0.0000)
2 265 /
815
107
678
993
1000
1.2
1.0
diabetes
2
133
0.0015 (0.0009)
0/ 1/
92
0/
207
992
1000
504
4.0
iris25
1
25
0.0000
(0.0000)
44 2822
868
1510
2245
2594
8000
0.23
5.9
iris25
2
10
0.0059 (0.0040)
0// 9
127
0/
0/
74
8000
38
20
iris25
3
irisso
1
20
0.0000 (0.0000)
115 321 /
828
1794
2311
2976
8000
0.41
5.9
iris50
2
23
0.0064 (0.0064)
0/ 37
413
0/
9
119
8000
43
20
irisso
3
0.2019 (0.0414)
0/ 18
128
0/
31 /
178
8000
270
100
wine225
32
0.0000 (0.0000)
23 /29
853
1864
2294
/ 3026
8000
0.36
5.8
wine225
0.0043 (0.0016)
11 34 /
341
0/
4
22
8000
44
19
wine225
3
0.1777 (0.0068)
300 / 449
598
60
63
8000
340
99
wine250
27
0.0000 (0.0000)
169 / 491 2129
1892
2191
2955
8000
0.41
5.9
wine250
2
21
0.0053 (0.0028)
0/ 7
615
0/
4
46
8000
42
20
wine250
3
13
0.1859 (0.0560)
13 / 110/
489
1/
39
93
8000
300
99
bounding intermediate neurons with piecewise affine trans- formations and known ranges of activation functions evalu- ated with IBP (Gowal et al. 2019), we achieve more precise functional bounds. These bounds converge the func- tions of input variables as local linear units increase.
Acknowledgments
We would like to thank three reviewers for their helpful feedback and suggestions, which improved our work.
The research in this paper has been partially funded by the Vienna Science and Technology Fund (WWTF) grant [10.47379/CC19008] (ProbInG), the TU Wien Doctoral College (SecInt), the FWF research project P 30690-N35, and WWTF project ICT22-023 (TAIGER).
Our method targets neural networks with continuous mono- tonic piecewise twice continuously differentiable activation functions using tools like Marabou (Wu et al. 2024), orig- inally designed for piecewise linear functions. While the current approach analyzes the behavior of NNs on a compact hyperrectangle, can easily extend our theory unions of bounded polytopes. In future research, plan to bound the edf of a neural network where the input admits arbitrary distributions with bounded piecewise continuous pdf sup- ported on arbitrary compact sets. Moreover, we intend to improve the algorithmic performance so that our method applies to larger networks.
Impact Statement
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none which we feel must specifically highlighted here.
9

---

Exact Upper and Lower Bounds for the Output Distribution of NNs with Random Inputs
References
Gawlikowski, J., Tassi, C. R. N., Ali, M., Lee, J., Humt, M., Feng, J., Kruspe, A., Triebel, R., Jung, P., Roscher, R., Shahzad, M., Yang, W., Bamler, R., and Zhu, X. X. A survey of uncertainty in deep neural networks. Artificial Intelligence Review,
Abdelaziz, A. H., Watanabe, S., Hershey, J. R., Vincent, E., and Kolossa, D. Uncertainty propagation through deep neural networks. In Interspeech 2015, pp. 3561-3565, 2015.
2023.
Gehr, T., Mirman, M., Drachsler- D., Tsankov, P., Chaudhuri, S., and Vechev, M. Aí2: Safety and robustness certification of neural networks with abstract interpreta- tion. In 2018 IEEE Symposium on Security and Privacy (SP), pp. 3-18, 2018.
Aeberhard, and Forina, M. Wine. UCI Machine Learning Repository, 1992.
Berzins, A. Polyhedral complex extraction from ReLU net- works using edge subdivision. In Krause, A., Brunskill, E., Cho, K., Engelhardt, B., Sabato, S., and Scarlett, (eds.), Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Ma- PMLR, 23-29
Goodfellow, I. J., Pouget-Abadiis J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., and Bengio, Y. Generative adversarial nets. In Proceedings of the Neural Information
chine Learning Research, pp.
28th International Conference Processing Systems Volume 2, NIPS' 14, pp. 2672-2680, Cambridge, MA, USA, 2014. MIT Press.
Jul 2023.
Bibi, A., Alfadly, M., and Ghanem, B. Analytic expressions for probabilistic moments of pl-dnn with gaussian input. In 2018 JEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 9099-9107, 2018.
Goodfellow, I. J., Shlens, J., and Szegedy, C. Explaining and harnessing adversarial examples. In Bengio, and LeCun, Y. (eds.), 3rd International Conference on Learn- ing Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, 2015.
Bunel, R., Lu, J., Turkaslan, I., Torr, P. H., Kohli, P., and Kumar, M. P. Branch and bound for piecewise linear neural network verification. Journal of Machine Learning Research, 21(42):1-39, 2020.
Gowal, S., Dvijotham, K. D., Stanforth, R., Bunel, R., Qin, C., Uesato, J., Arandjelovic, R., Mann, T., and Kohli, P. Scalable verified training for provably robust image clas- sification. In Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV), October 2019.
Cybenko, G. Approximation by superpositions of a sig- moidal function. Mathematics of Control, Signals and Systems, 2(4):303-314, 1989.
Delaunay, B. Sur la sphère vide. Bulletin de l'Académie des Sciences de I'URSS. Classe des sciences mathématiques et na, 193466):93-8000, 1934.
Hafiz, A. M. and Bhat, G. M. A survey of deep learning techniques for medical diagnosis. In Tuba, M., Akashe, S., and Joshi, A. (eds.), Information and Communication Technology for Sustainable Development, pp. 161-170, Singapore, 2020. Springer Singapore.
Driscoll, T. A. and Braun, R. J. Fundamentals of numerical computation. Society for Industrial and Applied Mathe- matics, Philadelphia, 2018.
Hoeffding, W. Probability inequalities for sums of bounded random variables. Journal of the American Statistical Association, 58(301): 1963.
Efron, B., Hastie, T., Johnstone, I,, and Tibshirani, R. Least angle regression. The Annals of Statistics, pp. 407-451, 2004.
Hornik, K., Stinchcombe, M., and White, H. Multilayer feedforward networks are universal approximators. Neu- ral Networks, 255):359-366, 1989.
Fawzi, A., Fawzi, H., and Fawzi, O. Adversarial vulner- ability for any classifier. In Proceedings of the 32nd International Conference on Neural Information Process- ing Systems, NIPS' 18, pp. 1186-1195, Red Hook, NY, USA, 2018. Curran Associates Inc
Hosseini, H., Xiao, B., and Poovendran, R. Google's cloud vision api is not robust to noise. In 2017 16th IEEE International Conference on Machine Learning and Ap- plications (ICMLA), pp. 101-105, 2017.
Ferrari, C., Mueller, M. N., Jovanović, N., and Vechev, M. Complete verification via melti-neuron relaxation guided branch-and-bount In International Conference on Learning Representations, 2022.
Hüllermeier, E. and Waegeman, W. Aleatoric and epistemic uncertainty in machine learning: an introduction to con- cepts and methods. Machine Learning, 11003)4577-06. 2021.
Fisher, R. A. The use of multiple measurements taxo- nomic problems. Annals of Human Genetics, 7:179-188, 1936.
Jaichandaran, S. Standard classification with banana dataset, 2017.
10

---

Exact Upper and Lower Bounds for the Output Distribution of NNs with Random Inputs
Ji, W., Ren, Z., and Law, C. K. Uncertainty Propagation in Deep Neural Network Using Active Subspace, 2020.
van Buuren, S. and Groothuis-Oudshone K. mice: Multi- variate imputation by chained equations in R. Journal of Statistical Software, 45(3):1-67, 2011.
Katz, G., Barrett, C., Dill, D. L., Julian, K., and Kochender- fer, M. J. Reluplex: An efficient SMT solver for verifying deep neural networks. In Majumdar, R. and Kunčak, V. (eds.), Computer Aided Verification, pp. 97-117, Cham, 2017. Springer International Publishing.
Wang. Z., Albarghouthi, A., Prakriya, G., and Jha, S. Inter- val universal approximation for neural networks. In Proc. ACM Program. Lang., volume 6, New York, NY, USA, Jan 2022. Association for Computing Machinery.
Katz, G., Huang, D. A., Ibeling, D., Julian, K., Lazarus, C., Lim, R., Shah, P., Thakoor, S., Wu, H., Zeljić, A., Dill, D. L.,, Kochenderfer, J., and Barrett, C. The Marabou framework for verification and analysis of deep neural networks. In Dillig, and Tasiran, S. (eds.), Computer Aided Verification, pp. 443-452, Cham, 2019. Springer International Publishing.
Weng, L., Chen, P.-Y., Nguyen, L., Squillante, M., Boopa- thy, A., Oseledets, I., and Daniel, L. PROVEN: Verifying robustness of neural networks with probabilistic ap- proach. In Chaudhuri, K. and Salakhutdinov, R. (eds.), Proceedings of the 36th International Conference on Mu- chine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 6727-6736 6736. PMLR, 09-15 Jun 2019.
Krapf, T., Hagn, M., Miethaner, P., Schiller, A., Luttner, L., and Heinrich, B. Piecewise linear transformation propagating aleatoric uncertainty in neural networks. In Proceedings of the AAAI Conference on Artificial Intelli- Mar 2024.
Wu, H., Isac, O., Zeljie, A., Tagomori, T., Daggitt, M. L.,, Kokke, W., Refaeli, I., Amir, G., Julian, K., Bassan, S., Huang, P., Lahav, O., Wu, M., Zhang, M., Komen- dantskaya, E., Katz, G., and Barrett, C. W. Marabou 2.0: A versatile formal analyzer of neural networks. In Gurfinkel, A. and Ganesh, V. (eds.), Computer Aided Verification 36th International Conference, 2024, Montreal, QC, Canada, July 24-27, 2024, Proceedings, Part II, volume 14682 of Lecture Notes Computer Sci- ence, pp. 249-264. Springer, 2024.
gence, volume 38(18). pp.
Lasserre, J. B. Simple formula for integration of polyno- mials on simplex. BIT Numerical Mathematics, 61)) 523-533, 2021.
Raghu, M., Poole, B., Kleinberg, J., Ganguli, S., and Sohl- Dickstein, J. On the expressive power of deep neural networks. In Precup, D. and Teh, Y. W. (eds.), Proceed- ings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learn- ing Research, pp. 2847-2854. PMLR, 06-11 Aug 2017.
Xu, K., Shi. Z., Zhang, H., Huang, M., Chang, K., Kailkhura, B., Lin, X., and Hsich, C. Automatic perturba- tion analysis on general computational graphs. Lawrence Livermore National Lab. (LLNL), Livermore, CA (United States), 02 2020.
Rao, R. R. Relations between weak and uniform conver- gence of measures with applications. The Annals of Math- ematical Statistics, 33(2):659- 6800 1962.
Yurtsever, E., Lambert, J., Carballo, A., and Takeda, K. A survey of autonomous driving: Common practices and emerging technologies. IEEE Access, 8:58443-58469, 2020.
Rudin, W. Principles of mathematical analysis. McGraw- Hill New York, 3rd edition, 1976.
Shi, Z., Jin, Q., Kolter, Z., Jana S., Hsieh, C.-J,, and Zhang, H. Neural network verification with Branch-and- -ouudd for general nonlinearities. In Gurfinkel, A. and Heule, M. (eds.), Tools and Algorithms for the Construction and Analysis of Systems, pp. 315-335, Cham, 2025. Springer Nature Switzerland.
Zhang, B. and Shin, Y. C. An adaptive gaussian mixture method for nonlinear uncertainty propagation in neural networks. Neurocomputing, 458:170- 170-183, 2021.
Zhang, H., Weng, T.-W., Chen, P.-Y,, Hsieh, C.-J., and Daniel, L. Efficient neural network robustness certifica- tion with general activation functions. In Proceedings of the 32nd International Conference on Neural Information Processing Systems, NIPS' 18, pp. 4944 4953, Red Hook, NY, USA, 2018. Curran Associates Inc.
Sicking, J., Akila, M., Schneider, J. D., Hüger, F., Schlicht, P., Wirtz, T., and Wrobel, S. Tailored uncertainty estima- tion for deep learning systems. CoRR, abs/2204.13963 2022.
Zhang. H., Chen, H., Xiao, C., Gowal, S., Stanforth, R., Li, B., Boning, D., and Hsich, C.-J. Towards stable and efficient training of verifiably robust neural networks. In International Conference on Learning Representations, 2020.
Sudjianto, A., Knauth, W., Singh, R., Yang, Z., and Zhang, A. Unwrapping the black box of deep relu networks: Interpretability, diagnostics, and simplification. ArXiv, abs/2011 1. 1,40411, 2020.

---

Exact Upper and Lower Bounds for the Output Distribution of NNs with Random Inputs
Ziegler, G. M. Lectures on polytopes. Springer-Verlle Verlag, New York, 1995.
12

---

Exact Upper and Lower Bounds for the Output Distribution of NNs with Random Inputs
A. Proof of Theorems
A.1. Theorem 3.8
Suppose the activation function in the prediction NN 2 is ReLU and no and nz are the number of input and output neurons, respectively. The integral of a function over a given domain can be expressed as the sum of integrals over a partition of the domain (disjoint subdomains whose union constitutes the original domain). To compute the edf of Y - (x) at y. (v) Pr[f2. (x) S y], compute the sum of (X(NNN(x) S y E Pyl] each of which is the integral of the pdf of the input over the given polytopes subject to NN (x) S These sets of polytopes {P3} and the corresponding local affine transformations [NN' always exist, as shown in (Raghu et al. (017]).
Pj is a convex polytope and can be represented as the intersection of halfspaces (see (Ziegler 1995 ). The set (x (x) = ed + V S y is defined as the intersection of halfspaces
no
NN) c + Σ
Vt
which when intersected with P, and k; defines the reduced complex polytope Pj. The desired local probability, P(NN'(x) y x E Pyl, is the integral I (x): of the pdf (x(x) over the reduced polytope Pj...
Using the Delaunay triangulation (Delaunay 1934) one can decompose any convex polytope Pj into disjoint set of simplices Te. This triangulation allows us to compute the integral over the polytope as a sum of integrals over each simplex. Assuming that the pdf of the input is piecewise polynomial allows us to use the algorithm from (Lasserre 2021) to compute exact integrals over simplices. The sum of all these localized integrals (probabilities) is the exact edf value at pointy. y.
A.2. Theorem 3.11
Since Y is feedforward NN with continuous activation functions on compact support C R"", for any > O, let {en). > 6 > 0 be decreasing sequence. By the UAT (Hornik al. 1989), there exists sequence of ReLU networks (Ye.). such that
sup fx() Ye, (R Cz. K
Setting (x(()) = YY,,
(x) (x) 6n and (,(x) = Ye, (x) + en. we have
(x, - (x) S (x).
Now, we let Y, (x) = maxigiga Y", (x), which is still not greater than Y(x), and Y,,))) minigisn Y, (x), which is still not smaller than (x)). One can see that ( and (Y,) are monotonically decreasing and increasing, respectively. should be noted that the min and max operators can be represented as ReLU networks (see Fig. 4 and the composition of ReLU networks is itself ReLU network.
A.3. Theorem 3.16
According to the UAT (Hornik et al. 1989 for any 6 > O there exist one-layer networks Y, T with ReLU activation function, such that
sup || W(x) Y(x) || x
sup || d(x) - x(x) |<E.
xek
Define (x) = (xx) - e, which also a NN. Similarly, for Y of G.. Then,
W(x) - 2e S = (x) es W(x) < (xx) + W(x) +e o(x) - 2 < (x) σ(x) - o(x)= < x(x) + (x) ((x) + 2e,
Letting e + 0 results in Y,, Y, + W and Φ, 2 + の uniformly on a compact domain while guaranteeing that they be lower/upper bounds.
13

---

Exact Upper and Lower Bounds for the Output Distribution of NNs with Random Inputs
Let
/
F.( min 1,
(x)dx
F,(()) = max 0,
2, (x)dx
(x0)))
(xxEKnY (x))
The limit edf is
F(y)= =
o(x)dx
(xXEX)W(X)5)
Since 2, (x) & 6(x) S (x) and Y, (x) z W(x) & (x) for any X K. xx: KnY, S 5 fx:
Knw(() S y) and fx: xx K Kny, v) C (x: S v). Since 0 S F(y) < for all y.
E.(0) S F()) s Fv
for y E (W(x) xx K}.
Now let us fix an arbitrary y W(x) for x E K, such that y is continuity point of F.
/
(r) FF)) s
(x)x)d
((x)dx
(xxE
(xXEKNWO)SSU)
/
(Φ,(x) o(x))dxx+
o(x)dr
((x)dr
(xxEKnY,,xx55)
(xxEKnY,
(xxEKnW(x)s))
Fly) - E.( -
o(x)da
L (x)xx
(xxx
(xxxKnY (1))sw))
/
/
(6(x) - ②, (x))
((x)dx
Φ(x)dx
=
(xxX
The left integrals in both equations (A and C) converge to zero due to the uniform convergence to zero the integrands over the whole set K. The second differences and D) converge to zero, since the superset {x x E Kny, (x) S y) and subset (x K n Y (x) s y) of the limits of integrals, respectively, converge the true limit set {x due to continuity. We have proven pointwise convergence for every point of continuity of the limiting edf.
K n W(x) S y)
Requiring & means that the limiting distribution has no point mass; i.e., it is continuous. The support of Y is compact because is the continuous image of the compact set K. We can then apply Polya's 's theorem (Rao) 1962) that the convergence of both bounds is uniform as sequences of monotonically increasing functions converging pointwise to a continuous function on a compact set.
B. Convergence of the ReLU bounds
B.1. Preliminary results
Lemma B.I. Let ake dx+1 R with 0x41 dn and assume f (ax, ax+1] R is twice continuously differentiable, monotone increasing and strictly convex (concave) on [ak, ax+1]. Then the values of the linear interpolation and piecewise tangent approximation at boundary points coincide with the original function. That is, Jlin int (ax) fiie. ton (ax) f(x)) and film (ax+1) = foie.tan (ax+1) (ax+1))
14

---

Exact Upper and Lower Bounds for the Output Distribution of NNs with Random Inputs
Proof.
Pln inn.itt (ax))= ((ax) + (ax -xx)m ((s))
plinine (ax+1)) = f(a + (ax+1 - Okk'
1intte )152
lin int) (gin f(ax+1) - f(a iin_inty
[(akiinng int) + (ax+1 ok
ak+1 0k
=f(ax+1),
frie (ax) = f(ax)) + ft (ax))(a -ax) f(ak);
+ fniosan (akt1)= = ((ax+1) + f' (ax+2) - ax+1) f(ax11).
Lemma B.2. Let (ak, ax+1] E R be a closed interval with ax+1 > ak, and let f [ak, ox+1] + R be twice continuously differentialle, monotonically increasing, and strictly convex (or concave) on (ag, ak+1I. Then, the intermediate points lie strictly within the interval [ay,au ak+1]. ag 0k lin.int < ak+1 and as dg
pie ak+1, and the functions of linear
interpolation and piecewise tangent approximation are continuous [ak, ak+1)
Proof. For linear interpolation, the statement follows immediately from the definition. Specifically, for the midpoint interpolation, we have
lin.int
0kk + ax+1
d'g
2
ak + ax+1
< ax+1-
2
For piecewise tangent approximation, we define d'
fan as
f(ax) - f(ax+1) (f4(aa)aa - 8- )ax+1)a)))1
&k pie
(ax+11) f4((x))
Consider first the case where f is convex on (ak, -1- By the convexity of f. we have the inequality
(as (x+1)f_((H++ < f(ax) f(ax+1) < (ax ax+1)f4
Adding the terms f' (ax+1)ax41 - (ax)ax to each part of the inequality, we get
an (ax+1) an.f (ax) < f(an) - f(a++)) f4 + f' (ax+1)ax41
< axt1ff' (ax+1) (0++1f((ax))
i.e., l- (ak+1)) f4 (ax), is strictly positive by convexity, dividing the entire inequality Since the denominator of pie tan d' by this denominator yields
ag < o
pie.tan
< <a++1-
In the case where f is concave on [ag, ds+1], we have similar inequality:
(a++1 - ((x4)())))))) < f(x+11) - f(ax) < a4)
Adding the terms F4 (as )a - ff (ax+))))aa to each part the inequality, we get
axf% (ax) a+f" (a++1) < f(ax+1) f(xx) + (ac)ax f- (as+2)0441
< ax+1 (ax) ax+1f- (f
15

---

Exact Upper and Lower Bounds for the Output Distribution of NNs with Random Inputs
Since the denominator I'_ (= by this negative denominator yields
- -f- (ax) of pie strictly negative for concave functions, dividing the entire inequality dg
pie.tan Ok'
ag
<
Next, we show that the interpolation functions are continuous at the intermediate point. For linear interpolation, we check the continuity by verifying that the parts meet at ag afimsint We have
Mak Nin int
f(ak fin int)
lin snt + (ak
Iên int )ml right plän int akk
(ake lin int).
fla (int) f((k))
left f(ax) + (aim-int aK)K1 f(ax) + (ak (aim-sme ag)
Ok
aj
Thus, the two parts of the linear interpolation meet at ain and finsint is continuous at liriintt
For piecewise tangent approximation, we check that the left and right parts meet at picston From the left, we have
juic tan (an(T)= = [[(u)) + (a(())(T ax), E [ax, ok'
pie tan Jan].
Substituting T = ag', we get
fosee fere.tun ( s(ax) + (4 (a<)))
- ox)-
From the right, we have
frieson (T) f(ak+1) + (ax+1)(T ax+1),
-a++1), TE lak pie ak+1].
Substituting ak's we get
jrinsan (ax)= S - ax+1).
Setting these equal, we have
f(xx) + S (ax) - (a axt1).
Solving for aw, we obtain
f(ax) - S(a+1)) (a= ax+15- (ax+1))
pin_tan
ag
(_((x+1) (ux)
dg
This confirms the continuity of fpie.tar(T) at a
tan
Lemma B.3, Let 0k+1 > a% and (ak: ax+1] C R. and let f [ax, ax+1] + IR be a twice continuously differentiable, monotonically increasing, and strictly convex (or strictly concave) function on [ake 4k-1). Then, the estimating functions defined by linear interpolation fliminnt and piecewise tangent approximation fricton are non-decreasing [ak, d++1].
Proof. We prove that both local estimators non-decreasing functions.
Case 1: Linear Interpolation. The slopes of the linear interpolation segments are given by
f(atininty
lin f(as) O, 162 f(ax+1) 1) f( (aat) (kk
K1
int
0.
0K+1
0
lin ak+1, as established by Lemma B.2 both slopes are Since the original function f is non-decreasing and aj dg non-negative, ensuring that the interpolated function is non-decreasing.
Case 2: Piecewise Tangeni Approximation. The derivatives of both the left and right segments of the piecewise tangent approximation are non-negative due to the increasing nature of the approximated function. Furthermore, by Lemma B.2 we have < dg pie_tan < ax+1, and the approximation remains continuous everywhere. Since both segments are non-decreasing linear functions, their combination also results in a non-decrend non-decreasing estimator over [ak; ax+1].
Thus, both estimation methods preserve the monotonicity of f-
16

---

Exact Upper and Lower Bounds for the Output Distribution of NNs with Random Inputs
Lemma B.4. Let a a IR be closed interval, and let f c IR be a continuous function satisfying C..u. w. (le, a)). Assume that there exist points @ =1 02 <
0n+1 for some E N such that Twice continuously differentiable, monotonically increasing, and either strictly conves, strictly concave, linear each subinterval [ax, ax+il for 1 S k S n.. Then, the approximation method described in Section 3.2 constructs a continuous, non-decreasing estimating function (x) over the entire interval (a, a.
Proof. The claim follows from the following observations:
a) Each subinterval [as, 4x+1] is suitable for approximating f using either linear interpolation, piecewise tangent approxi- mation, or local linear approximation.
b) All the methods mentioned in (a) ensure that the estimator matches the original function at the boundary points of each subinterval, by Lemma B.
c) The approximations described in (a) are continuous within their respective subintervals (see Lemma B.2).
d) The methods in (a) produce non-decreasing functions within each subinterval (see Lemma B.3
By sequentially linking the estimators across all segments {[ax, (k+1]], we obtain a continuous, non non-decreasing, piecewise linear function fxx) over (a, a].
Lemma B.5 (Image of the local estimator). Let f la, a) TR. be continuous, e Cp. (la, al) with a1 02 a a for n E N, such that I'I [1 is twice continuously differentiable, monotonic increasing and either strictly convex (concave), or linear on [ak, 0x-1 S k S n. Then the image of the estimating function f(x) defined by the approximation method in Section 3.2 coincides with the image of the target function that is f((a.a)) f(la, a)).
compact [g, a], and monotonic, then it maps [g, a] into the closed interval (compact) [f()), [()] E IR. But the estimator f also continuous monotonic function [a, a). and by Lemmas B. B.2 f(a) (a) and f(a) f(a). That is why, the ranges of values of f and f on [g. ब coincide and equal to f(aa..
Proof. Since function is continuous
R.
Theorem B.6 (Convergence of piecewise linear interpolation (Driscoll & Braun 2018). Ch.5), Suppose that f(x) has continuous second derivative in [ag, as+1), that is C2 (aaR,, ax+1)) Let Prina (x) be the piecewise linear interpolant of (ok, flay,) )) for i O...., Dint, where
ak agk + ih, h 0h+1 Ok
nint
Then, the error bound satisfies
Il Primeel loc
max x xE(ap,a2+1)
1f(x) Prine (x)) Mh²
where
M max [ar.0241]
f*)
Theorem B.7 (Convergence of piecewise tangent approximation). Suppose that f E C2 ([ax, ax+1)) and is strictly convex (or concave) in [ak, ax+1]. Let fini tan (x) be the piecewise tangent approximator over subsegments (ak, akmn for i == ....., Ntans where
0 0k
ak; = =0g + ihi h
Usan
Then, the error bound satisfies
||
Fric tan
Rean
||
max x xE(ap.aat11)
[f(x) Jiie tan (x) Mh?
where
M = max (x)
(6441]
17

---

Exact Upper and Lower Bounds for the Output Distribution of NNs with Random Inputs
Proof. Each element of the piecewise tangent approximation the Taylor series expansion of the first order around the boundary point of the subsegment. We consider the double Taylor series approximation the refinement (ak, Ok4+1) of the segment [ak, dx+11 for O,... Rint, where
h 0k++ 0k
ak; = + th,
Rian
Since for the twice continuously differentiable in [ake aka+1 1976), which represents eπrоr term, be bounded with
the Lagrange Remainder of the Taylor series expansion (Rudin
pic_tan akk )2
f(xx) (o))) s Mi,
(ak)
(ak; (ka)) )2
max
S M,
Abe a aiictton pictton
2
2
(aka+) pic_tan dK
max
f(x) - -fiietan(a) S Mir
Mi
[ak;
Okk. )2
piestan bf
2
2
where
Ms,- = max
(²) max
fz(x) = M
M2
max
(z)
f"(x) M
plccton
The maximum M exists and is attainable due to the continuity the second derivative on compact [ak. 01-1]. That is, the maximum error on the whole segment of approximation can be bounded as
M
ak+1 ak
Mh²
ntan
Mean -0000
The following lemma demonstrates that the approximation procedure presented in Section 3.2 generates sequences of estimators that: i) serve as valid bounds for the target function, and ii) converge monotonically to the target function. This implies that each new estimator can only improve upon the previous onc.
Lemma B.8. Suppose that f C² at+1l) is monotonic increasing and strictly convex (or concave) [ak, ak++)). Let iliniint (x:) be the piecewise linear interpolant and [gie_tan (x:) be the piecewise tangent approximator over subsegments [kk, ak.o for O..... 2n, where
0k+ Ok
ak; = ag + ilt h
2
with n. E N. Then the estimating functions defined by the linear interpolation Flanant and piecewise langent approximation prie.ton define upper (lower) and lower (upper), respectively, bounds the target function f. Moreover, flin (x) hm and {f ton (z) are non-increasing (non-decraasing) and non-decreasing non-increassing sequences. respectively.
Proof. By the definition of the convex (concave) function,
f(am (1 a))2)) < (2)af(m)) + (x)
for any a E [0, 1] and 1 2y from the region of convexity, and plot of the linear interpolant between any 3 31, lies above (below) the plot of the function. That is, linear interpolation is always an upper (lower) approximation of the convex (concave) function.
18

---

Exact Upper and Lower Bounds for the Output Distribution of NNs with Random Inputs
On the other hand, any tangent line lies below (above) the plot of the convex (concave) function. Indeed, since on a convex segment the derivative of the function increases, that is S4 (ak,) (x) I f''''''aaaa (ak for all [ak; 0k,-1 then
f(x) = f(aa)) + / ((t)dt 2 f(ak)) + / f4 fa(k))dt
42;
= =f(ax, + f+ fa(m - and), x xe [ax, axl
/ f(t)dt f(aker)) /
f(z) = f(an.1))
f" (atls))dt
aRL+1
=f(akut) ('(oort))) aketl); (arrakues
Similarly, we can show that the piecewise tangent upper bounds the true concave function. Without loss of generality, we consider the case of a convex segment. We fix n. The current element of the sequences of piecewise linear interpolants (l)n.in((:)) includes the local linear item based on the interval [2]. The case of a concave segment is analogous. We define the current local linear approximation of the element of sequence of upper approximation as
(12) - ft1)
Piiy.in (x) f(t)) + (x t1) (1)
to -1
Let be such that 0 S < 1 and define new point é as a convex combination t ati + (1 a2)2. Let show that (lninn (x) S Ain.int( (z) for [ti, t2). where
f(ti) -x-tt) f(e)) 6-
x [f1,t]
f(t) + -x-x-x00
- 1) 1(2)-111))
x x tt,
Taking into account the convex segment,
f(h) + (a 1n) 1(1) ((1) 1())) =
f(t) + (x- -t) 11) f(at (1 (1-a))) s(4)) at + (1 (1- attz t1
t- t1
f((h)) + (s - 11) af(t)) + (1- (t) f(t)
(1 (1-a) - (1)
= -((h) + (s -1) f(2)] - ((4)
t2 - t
f(1) + (x [(12) s(1)
s(at1 a1)a))
t22
+ (z at + (1 - a)t2) f(22) f(at1 + (1 a)t2)
l2 a4i + (1 a1
$ aa((4) + - (1-a),(1) (b2)
+ + 41) af(()) + (1 - ()) - f(t2)
4h + +(1 (1-0)12 12
+ +( - th + (1 - aa)2) al(()) + (1 - -a) s(12)
at + (1 a)tt - 12
= f((h) + (ar - (1) f(12) - f(h))
4 - 4
19

---

Exact Upper and Lower Bounds for the Output Distribution of NNs with Random Inputs
Since the refinement on each subsegment leads to the reduced next element of the sequence, the sequence is decreasing the whole convex segment.
We next consider the piecewise tangent approximation on [t,. t2l.
Jpie.tan (z) f(tt)) f4 t1), ((0z) (((())((x= 12).
x [1,
x (2r,
where Ep is the point of intersection the tangents. If we choose some parameter a, where & a $ 1, and define the corresponding intermediate point as t, atı + (1 a2) then the refined approximation is given by:
f(t,) + (4)((x +(t))+ -otottt,)) to), x E [t b2+l, ((12) (())))(((()))) (2),
x E [11, 41=,
tan (x) Aereer
[ 12),
where i is the intersection point of the left and middle lines, and 12+ the intersection point of the middle and right lines. We aim to show that ppictton
[1)m+1 (x) 2 [lpp. (x) for x (2).
Thus, the plot of the linear function corresponding to the middle curve lies above that of the left curve for x > t1s and above that of the right curve for x 1z, due to the monotonicity of the estimator, by Lemma B.3 Consequently, the refined estimator f, coincides with the previous estimator f, on the left and right segments, i.c., for [t,, and [t2. while it takes higher values for [[1,, /2.]. To confirm this, remains show that Lt. I by 12..
First, we note that t S tr I t2. by Lemma B.2 This automatically Ieads to 1 S t1 I t. S 22 S tz.
Consider the function
fua(z) xf-(()) - f(x) -C
f=(x)) K
defined on (41,22) where C= (1)((()) 5(4) and K (11).. Its derivative is given by
fea(x)) (x) + C- Kx)
(f((z)- - K)2
The numerator simplifies to
Kx) -V (f(ti) f4(t)((-
which is positive due the convexity of f. This implies that shifting the right boundary 12 to the left, reaching position t.,, also shifts the intersection point t the left, reaching f1.. A similar argument holds for the right boundary.
The same reasoning applies to the concave segment.
Lemma B.9. Let Jn As Ag be continuous functions, uniformly convergent to continuous function f A, - Ag on a compact interval A₁ C R. and let 9n Ag + A be continuous functions, uniformly convergent to continuous function funnciionnnnn gn (f(())) converges uniformly to g: Ag + A on a compact interval Ag C R. Then the sequence of composition functions g(f(x)) on Ay with n + S0.
functioonn
Proof. An outer limit function, g. is uniformly continuous by the Heine-Cantor theorem (Rudin, 1976), since it is continuous and defined on the compact Ag. That is, for any é1 > there exists e2 >> O such that
1 g(12)1 1l
whenever ø1, V/2 Ag and 1y1 y2l < E2.
Since the sequence (f)( converges uniformly to fxx on Ay. and (gn (v) converges uniformly to g on Ags we can conclude that for any 6 and é2 > 0, there exists e N such that for all n 2 N, we simultaneously have
(v) 9(y)) < C1 2
for all Ag.
[fx(z) f(x) E2,
for all x E Aç.
20

---

Exact Upper and Lower Bounds for the Output Distribution of NNs with Random Inputs
Since the range of possible values of Jn (x)) coincides with the range of values of (xx which equals [1(A,) Ag. the
domain of g. obtain
9n (fn -
uniform convergence of the outer
uniform continuity of the outer
61 E1
=1
Thus, the uniform convergence of the composition follows.
Lemma B.10. Let fn : R be continuous functions that converge uniformly to continuous functions f : AA + on compacts A; C R. and let a; E I be constants for i = 1...., ny. Then linear combination of sequences, fa, E. onJn defined on the direct product A A x x Ani C _ converges uniformly to fa Σi ogS A + R. where a = (001,-..,00 ); that is,
sup =EA
|a (x) - fa(()) n-000
The images of Sa (A) and (A) compact.
Proof. Since (A) and (AA linear combinations of continuous functions, defined compact sets, they continuous functions. Also, since direct product of compact sets is compact, by the continuous mapping theorem (Rudin) 1976), the images of (AA and (A4 are also compact.
If every sequence of functions {fx (x')) converges uniformly to the corresponding limit function f((x), then for any >> 0,
there exists an integer N E N such that for all n 2 N, we have
((('))(((x))) fax
sup z'EA:
for all I.... ny. Consequently, the supremum of the differences in the linear combination can be bounded
sup #EA
(x)) - [a(()) = sup =EAA
Σ (
72
Lal sup #E
[Pa(x*) (1(at) <c
This establishes the uniform convergence of the linear combination of functions.
Lemma B.11. Let F. S fr A. + R. be continuous functions compact intervals A, R. satisfying
[()) ( (()))
for all x² E A, and for every i 1....,
Suppose that the index sets I and J are disjoint and their union forms the full sequence:
IUJ [1....,mu}
Then, for any coefficients a4,, By 2 0.for i J, the following inequality holds:
Eaf((z) Σ EBF (x') < Σ aufi(z²) - E8,((+))
S Eaj
- (,f((')))
Σ
RET
JEJ
jEJ
GET
jEJ
for all xt E As and at E Aj.
21

---

Exact Upper and Lower Bounds for the Output Distribution of NNs with Random Inputs
Proof. For any non-negative coefficients 04,, By and given that [(xt) [i())) f,(x9), we derive the following inequalities:
af = a.f ax)), (x'), -B,7 - -3,5( (x))
Summing these inequalities over all indices completes the proof.
Lemma B.12. Let J. A be continuous funicttionn on compuct interval A C R, satisfying
f(z) S f(x) s f(z) for all x E A.
Furthermore, let 9.9.9 + R be continuous and monotonically increasing functions on a compact interval B R. satisfying
g(v) glu) S gly) for all y B. Then, for all x E A, the following inequality holds:
(((=)) (zf()))
Proof. By assumption, for any y E B,
g(v) g()) S g(y)-
Since g(v), I(v), and g(v) are monotonically increasing, it follows that for any 31: V2 E B such that y1 S S y2.
g(v)) $ g(v) S g(v) - g(v) S g(22).
Setting f(),, in [(x), and f(x), and using that [(a) S f(x) - f(x) for all x E A, we obtain the desired result:
g(f(z)) s <0 (x))..
B.2. Proof of Theorem 3.14
We establish the uniform monotonic convergence of the bounds by considering several key properties.
First, analyze the approximation of neuron domains. Using the concept of over-approximmting the input domain of each neuron (as in IBP (Gowal et al. 2019)), we define Q2 superset of the true input domain Ω for each neuron i. It is well known that if sequence of approximations converges uniformly on L2 it must also converge uniformly on any subset of Ω4 including Ω4 Therefore, we perform all neuron-wise approximations over the supersets of their original domains.
By Lemma B.2, both the upper and lower bounds for each neuron every layer are continuous functions over a bounded domain. Additionally, by Lemma [.5] the images of these bounds coincide with the image of the activation function. Since the activation function and its estimators continuous, and the domain is compact, the output range remains compact throughout the approximation process. This compactness follows from the continuity of the mapping (Rudin) 1976).
Furthermore, the eтror between the linear interpolation (or piecewise tangent approximation) and the original activation function converges to zero as the approximation grid is refined, as established in Theorems B.6 and B.7 Consequently, the sequence of bounds on the activation function converges uniformly to the target activation function for neurons in the network.
By Lemma B.5 the image of each activation function is preserved by both the upper and lower bounds, ensuring that the domain of uniform convergence for the bounds of the outer functions is also preserved. Moreover, by Lemma B. the linear combination of these estimators converges uniformly to the corresponding linear combination of the true activation functions, thereby preserving the image of the original linear combination. As result, Lemma B.9 guarantees the uniform convergence of the composed bounds on the outer functions and the linear combinations of the inner functions.
Next, we establish monotonicity. By Lemma B.8 the sequences of upper and lower bounds on the activation functions are monotonie: the upper bounds (f, (x)) are monotonically decreasing, while the lower bounds { (z)) are monotonically
22

---

Exact Upper and Lower Bounds for the Output Distribution of NNs with Random Inputs
increasing. Additionally, by construction, these sequences bounded by the target function itself, ensuring they remain within the correct range. By Lemma B.11 the linear combination of neurons' estimators forms the overall upper and lower bounds for input arguments in subsequent layers.
Moreover, the compositions of the upper and lower bounds for the inner and outer functions provide valid upper and lower bounds for the composition of the target inner and outer activation functions. Since the outer activation function (and consequently its bounds, by Lemma B.3) is non-decreasing, Lemma B.12 ensures these bounds are valid. Finally, these bounds converge monotonically by Lemmas B.8 and B.12
Combining all the results above, conclude that the sequences of upper and lower ReLU bounds for the network converge uniformly and monotonically to the target network. The monotonicity of the sequences, the uniform convergence of individual approximations, and the preservation of continuity and compactness through compositions collectively ensure the uniform convergence of the entire network. This completes the proof.
C. Subnetwork equivalent for specific functions for upper/lower approximation
C.1. Square
The quadratic function 2 not monotone on an arbitrary interval. But x² is monotonic on Rzo. We can modify the form of the function by representing it as subnetwork to be valid set of sequential monotonic operations. Since x² | E I. and the output range of |a:] is exactly R2o. we represent |uc] as combination of monotonic ReLU functions, |x| = (x) + ReLUU-=). The resulting subnetwork is drawn in Figure 3a
ACLV/D
y
-1
Rell((y)
X1
Y2 = x2
ya
Y
x1
X2
ACLVIT
ReLULY
X2
AcLopy
40.5
(a) Subnetwork equivalent to operation of taking a square.
RAIRTT
(b) Subnetwork equivalent to product operation.
Figure 3: Subnetworks for square and multiplication operations.
C.2. Product of two values
To find a product of two values 31 and 22 one can use the formula x1 22 = 0.5 ((21 + x222 x7 -x3). This leads us to the feedforward network structure in Figure 3b
C.3. Maximum of two values
The maximum operation can be expressed via a subnetwork with ReLU activation functions only, as follows. Observing that max( [1, x2) 0.5 (z) + 22 + [ - .22)) results in the corresponding network structure in Figure
C.4. Softmax
The function softmax transforms a vector of real numbers to a probability distribution. That is, if (x1,..., In) E IR.". then there is multivariate function Sf ax _ _ so that
SfMar, softmax(x.)
7
Σ
6
j=1
23

---

Exact Upper and Lower Bounds for the Output Distribution of NNs with Random Inputs
0.5
0.5
X1
.0.5
Re10(P)
-0.5
= max(x
0.5
X2
.0.5
0.5
0.5
Figure 4: Subnetwork equivalent to maximum of two values.
Then, log(SfMaz.) = log Σ e which is composition of monotonic functions. This leads to the feedforward
network structure in Figure 5
x
Xn
softmax(ss)
vnftmmex(sg)
voftmen(s
Figure 5: Subnetwork equivalent to one softmax node.
D. Description of Iris Experiments
We trained fully connected [3 x 12) ReLU NN with a final [1 x 3] linear layer, as well as fully connected [3 x 12] tanh NN with the same final [1 3] linear layer, on the Iris dataset. The networks classify objects into three classes: Setosa, Versicolor, and Virginica, using two input features: Sepal Length and Sepal Width. The allocation of the data for these two variables in the three classes is shown in Figure 6al The input data were rescaled be within the interval [0. 1].
The ReLU network was pre-trained. We next introduced randomness to the input variables by modeling them as Beta-distributed with parameters (2, 2) and (3, 2), respectively. The pdfs of these input distributions are shown in Figure 6b The first is symmetric about 0.5 and the second is left-skewed.
Experiment 1: ReLU-Based Network with Random Inputs.
In our first experiment, Example 3.10 we computed the exact cdf of the first output neuron (out of three) in the ReLU network before applying the softmax function. Due to the presence of final linear layer, the output may contain negative values. To validate computation, simulations, where the empirical cdf was estimated using 10 samples. As shown in Figure I both edf plots coincide. The cdf values were computed at 100 grid points across the estimated support of the output, determined via the IBP procedure.
compared against conditional ground truth obtained via extensive Monte Carlo
For further comparison, Figure 6c presents an approximation of the output pdf based on the previously computed edf values. This is compared to a histogram constructed from Monte Carlo samples. Additionally, we include a Gaussian kernel density estimation (KDE) plot obtained from the sampled data using a smoothing parameter of 0.005. The results indicate that
24

---

Exact Upper and Lower Bounds for the Output Distribution of NNs with Random Inputs
our pdf approximation better represents the underlying distribution compared to KDE and tracks the histogram closely.
Experiment 2: Bounding Tanh-Based Network with ReLU Approximations. In our second experiment, used pre-trained NN with similar structure but replaced the ReLU activation functions in the first three layers with tanh. To approximate this network, we constructed two bounding fully connected NNs- one upper and one lower-using only ReLU activations.
We conducted computations two regimes: one using 5 segments and another using 10 segments, into which both convex and concave regions of the tanh activation function at each neuron the first three layers were divided. Over each segment, we performed piecewise linear approximations according to the procedure described Section 3.2 and combined these approximations into three-layer ReLU networks with an additional final linear layer for both upper and lower bounds. The results of these approximations are shown in Figure 2 In the 10-segment regime, both the upper and lower approximations closely align with the original NN's 's output.
146
10-
Class
PDF Value
Class 2: versicolor
Class virginica
Peeature (Sepal Width)
83-
Bexd2 sloest 6
Beng 2): go 12 -xy)
supmort
26-
(b) Plots of the Bela(2, 2) and Beta 2) pdfs, resp.
Semppted hittggun Inttraned pdl
PDF Value
02-
14
60-
0.
62
0.4
06
Feature (Sepal Length)
first element NN
(c) Estimated output pdf compared with KDE with smoothing parameter h = 0.005 and histogram of MC simulations.
(a) Plot of Sepal Width vs Sepal Length with class indication.
Figure 6: Iris Dataset
25


Reviews:
[{"review_id": "vqbDPSQztZ", "reviewer": "Reviewer_3TKr", "summary": "This paper proposes a novel method to compute exact upper and lower bounds for the cdf of a neural network’s output, assuming stochastic inputs. Key contributions include:\n- A method to compute the exact cdf of the output of ReLU networks under inputs with piecewise polynomial pdf over compact hyperrectangle.\n- A method to compute tight upper and lower bounds of the cdf of the output of neural networks with piecewise differentiable activations under inputs with a pdf with compact support.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "*Theoretical Claims:*\n\nYou are right that Theorem 3.11 follows from the universal approximation theorem. But it also **constructs** the sequences whereas the UAT simply shows their existence.  \nIndeed the significance of Theorem 3.13 is that the bounds are built using ReLU functions, which is our main innovation as it also renders our approach applicable to general feedforward NNs. If we have the opportunity to revise, we will restate the theorem as you suggest. \nTheorem 3.8 obtains a more general result even though we apply it in the proof of Theorem 3.13 only with piecewise constant pdfs.\n\n*Claims And Evidence:*\n\nIn Table 1 we see that our bounds differ from the competing PLT and MC methods. You are right that it provides no evidence for the correctness of ours. We show formally that our bounds are correct. Table 1 serves to illustrate that the other two methods, PLT and MC, are less accurate. In effect, if we *significantly reduce* the accuracy of our bounds then they include the PLT and MC approximations. In other words, we can recover the PLT and MC estimates by making our method less precise.\nIn the Iris data, the inputs are perturbed by Beta noise since in this case (the Beta pdf is a polynomial) we can compute exact bounds. Of course, we could have had used other random perturbations as well.", "claims_and_evidence": "Correct me if I misunderstood but it seems that Table 1 merely shows that MC and PLT disagree with your bounds, but it does not provide evidence to support the correctness of your bounds?\n\nIn Figure 1, it says that you are assuming beta-distributed inputs, but how can you assume that if the input is the Iris dataset? Are you talking about random perturbations to the Iris dataset?", "methods_and_evaluation_criteria": "see claims and evidence", "theoretical_claims": "I did not check the correctness of the proofs, only looked at the statements.\n\nTheorem 3.13 seems to be a strong result but it is not well presented. For any F (y), of course there always exists upper and lower bounds F_n such that they converge to F (y): the interesting part of the result is that the authors are talking about particular F_n that can be constructed, so the last sentence \"The sequence {F n} can be constructed by bounding the distributions of\nsequences of ReLU NNs.\" should be the core of the theorem. The subsequent two paragrpahs (l.289-304 right block) where the authors describe how to construct F_n should be part of the statement, and not using words but mathematical notations... Also, from the text in the second paragraph (l.295-304 right block), it seems that polynomial pdfs are never used, so it would be perhaps be less confusing if Theorem 3.8 was stated for piecewise constant pdfs.\n\nI'm wondering if Theorem 3.11 does not follow directly from universal approximation theorems?", "experimental_designs_or_analyses": "see claims and evidence", "relation_to_broader_scientific_literature": "This paper contributes to an important area of the literature on estimating uncertainty in neural networks.", "essential_references_not_discussed": "Not that I know of", "other_strengths_and_weaknesses": "Strengths: strong theoretical contribution\nWeaknesses: see other comments", "other_comments_or_suggestions": "see other comments", "questions_for_authors": "if you could address my concerns about the experiments, and answer my comments on the theoretical contributions, I would be willing to increase my score.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "q5x2hr2U9k", "reviewer": "Reviewer_xfm3", "summary": "The paper addresses the challenge of uncertainty quantification in neural network (NN) outputs by deriving exact upper and lower bounds for the cumulative distribution function (CDF) of NN outputs under noisy (stochastic) inputs. The method is designed to apply to feedforward NNs and convolutional NNs (CNNs) using continuous monotonic piecewise differentiable activation functions, such as ReLU, tanh, and softmax. A novel aspect is the use of ReLU NNs to bound general NNs, enabling the computation of these bounds, which converge to the true CDF as resolution increases.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "*Claims And Evidence and Theoretical Claims:*\n\nThe proofs are not detailed as most of the steps are based on well known calculus facts and are easy to obtain. The challenge is that they are numerous and this is why we provided only the sketch of the proof. We appreciate that for the sake of clarity and completeness all the details should be made available and we will include them in the appendix if we are given the opportunity to revise. We will also place them in the arXiv version of the paper (available after the final decision).\n\n*Weaknesses:*\n\nTheorem 3.8 derives the exact cdf of a ReLU NN when the input's pdf is a piecewise polynomial. Theorem 3.11 follows from the universal approximation theorem, but the proof we provide **constructs** the sequences whereas the UAT simply shows their existence. Theorem 3.13 is our main result as it constructs bounds for the cumulative distribution function of the output of any feed forward NN whose inputs are continuous random variables with compact support. The proof of Theorem 3.13 is provided in detail. We fail to see how it is a direct corollary of the UAT and we would be interested to see the justification of this claim.", "claims_and_evidence": "I find a lot of the proof details missing. For example, the proof of Theorem 3.8 is hand-wavy, and I don't see all the steps. Even more problematic is the proof of Theorem 3.11, where the authors only provide proof sketches. I don't know precisely how all these fall into place without knowing the details.", "methods_and_evaluation_criteria": "yes", "theoretical_claims": "yes. See my comment above about the proofs. The proofs in these papers are just not explicitly written out, and the quality of the proofs is something that I would not even tolerate in an undergraduate proof-based class, let alone a research paper. \n\nWhile I understand the proofs are mostly analytical and require calculation, the authors must do this to the point where it is clear how the proofs follow. If they use any well-known theorems as black boxes, then they must explicitly cite a formal statement used, or even better include it as a preliminary subsection of the appendix, and invoke it from here.", "experimental_designs_or_analyses": "No issues.", "relation_to_broader_scientific_literature": "Not sure. See weakenesses.", "essential_references_not_discussed": "NA.", "other_strengths_and_weaknesses": "See summary.\n\nWeaknesses:\n\nSee my complain about the paper. Overall, I don't find any details to verify the proofs. Even the statements themselves look trivial to me. For example, 3.13 just seems to be a direct corollary of the classical universal theory of approximation of NN. Theorem 3.11 is perhaps the most interesting of the three theorems, but I don't find details to understand this for myself.", "other_comments_or_suggestions": "Please provide the proof details to the point where one can verify the correctness of the claims.", "questions_for_authors": "I would want to understand each theorem-by-theorem, why the statement is non-trivial, why previous works could not solve this, and what the key difficulty addressed by this work is.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "H9aFCcLXgg", "reviewer": "Reviewer_WtAa", "summary": "This paper provides upper and lower bounds on the cdf of a ReLU neural network which converge to the exact cdf as the granularity increases. Other monotonic piecewise differentiable activation functions can in turn be approximated using ReLU activation functions, extending the results. Experiments validate the theoretical results.", "questions": "", "limitations": "", "rating": 1, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "*Relation To Broader Scientific Literature and Questions For Authors:*\n\nOur approach aims to approximate as accurately as possible the **cdf** of the output in a probabilistic NN. Krapf et al (2024) estimate the **pdf** of a NN output. Our Theorem 3.8 derives the exact cdf (not pdf) of a ReLU NN when the input's pdf is a piecewise polynomial. That is, given a ReLU NN and an input with piecewise polynomial pdf, our tool (implementing Theorem 3.8) computes the cdf values **exactly** (not approximately). Theorem 2 in Krapf requires the input distribution be defined as what they call a piecewise pdf (different pdfs on the partition sets of the input support). The two theorems are *not* directly analogous. \nFurthermore, Theorem 2 in Krapf et al (2024) gives the formula of the probability density function (pdf) of the output of a feedforward NN subject to random perturbation of the inputs. The formula (Eqn (2) in Krapf et al (2024)) contains an integral that is never explained how to compute. In effect, Krapf et al (2024) do not apply Theorem 2 in their paper but rather its corollary (Theorem 3 in Krapf et al) which is a discrete approximation of the integral in Theorem 2. The latter is the reason we compare our approach to Theorem 3 in Krapf et al, which we could actually apply using the code in the Krapf et al GitHub page.  \nTo make things clear, in our Example 3.10, Krapf et al's method cannot compute the exact output pdf in finite time. Our approach computes the exact cdf of the output in less than a second. \n\nRegarding question 2, Theorem 3 and our approach yield similar results in the sense that the estimated cdf lies within our bounds (once the pdf is transformed to a cdf) provided the input pdf is very smooth (not wiggly).  \n \n*Methods And Evaluation Criteria:*\n\nFor MC, 100,000,000 sample points were used, which far exceeds the sample size in Krapf et al. The triangulation of the domain was upper bounded by 50,000 points, an ad hoc selection. The grid size we used was 20 raised to the number of classes of the output values. In the case of continuous outcome, the grid size was 20. For the Iris and Wine datasets, it was 20^3. We applied the provided PLT code from Github. For example, for the Iris dataset, we follow the settings in Krapf et al. in order to compare the performance of our method with theirs and MC. We report here computing details when we use 1, 2, 3 input variables. The emerging pattern is that our method is far more accurate and typically faster than PLT. We would also like to thank the referee for giving our method a concise name, namely U/L. \n\n*1 input variable: 25% of the data are imputed*\n\n**Computed MC in 1.873 sec**\n\n**Computed ours (U/L) in 5.700 sec**\n\n**Computed TLP in 0.348 sec**\n\nPLT out of bounds count = 2179 / 8000\n\nMC out of bounds count = 378 / 8000\n\n*2 Input variables: 25% of the data are imputed*\n\n**Computed MC in 2.114 sec**\n\n**Computed ours (U/L) in 19.663 sec**\n\n**Computed TLP in 46.043 sec**\n\nPLT out of bounds count = 0 / 8000\n\nMC out of bounds count = 53 / 8000\n\n*3 Input variables: 50% of the data are imputed*\n\n**Computed MC in 2.309 sec**\n\n**Computed ours (U/L) in 97.876 sec**\n\n**Computed TLP in 240.542 sec**\n\nPLT out of bounds count = 7 / 8000\n \nMC out of bounds count = 9 / 8000", "claims_and_evidence": "The theoretical claims are adequately proved.", "methods_and_evaluation_criteria": "There is a lack of detail provided about the experiments, making it difficult to evaluate the methodology. For instance, how much compute/granularity was provided for the MC, PLT, and U/L methods? Why should numerical estimation inaccuracy affect PLT but not the proposed U/L method? How many points were tested for the OOB calculations (is it the grid size)?", "theoretical_claims": "I did not carefully check the correctness, but the results are unsurprising in light of prior work.", "experimental_designs_or_analyses": "See \"Methods And Evaluation Criteria\".", "relation_to_broader_scientific_literature": "This paper is very much related to [Krapf 24]. One of the main theorems, Theorem 3.8, appears to be subsumed be Theorem 2 in Krapf.\n\n\nKrapf, T., Hagn, M., Miethaner, P., Schiller, A., Luttner, L., & Heinrich, B. (2024). Piecewise Linear Transformation – Propagating Aleatoric Uncertainty in Neural Networks.", "essential_references_not_discussed": "Even though [Krapf 24] is cited, it appears that the relationship between Theorem 2 in [Krapf 24] and Theorem 3.8 in this work is not discussed.", "other_strengths_and_weaknesses": "This paper extends prior work on computing cdfs to providing upper and lower bounds, instead of just an approximation, which better quantifies the uncertainty of the estimate.", "other_comments_or_suggestions": "N/A", "questions_for_authors": "1. Why use Theorem 3 in [Krapf 24] as a comparison for theoretical and experimental results, when Theorem 2 also applies?\n2. Is there a major difference between Theorem 3 in [Krapf] and the upper/lower bounds provided? Looking at equation 6, it appears that the upper and lower bounds follow from just upper and lower bounding the integrand.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: The paper looks at estimation of CDF for neural networks with
stochastic inputs. This is applied to ReLU and some other
activations. The topic is of interest in uncertainty quantification.
Most of reviewers agree that the contribution is novel and has
theoretical potential.  At the same time, some reviewers are concerned
with lack of experimental detail and insufficient discussion regarding
relation to [Krapf 24].  Some reviewers also point out lack of detail,
rigor, and clarity in the proofs.  Authors are strongly encouraged to
address these limitations.
Decision: Accept (poster)