According to the perturbation range \(x\in [x_0-\epsilon ,x_0+\epsilon ]\)  of input layer, we compute the linear bounds of the first hidden layer. Then, by backsubstitution [1]}from current layer to the input layer and leveraging eq(1),(2) we compute the bounds \([l^1,u^1]\)  of \(F^1(x)\) . Next, based on \([l^1,u^1]\) , we can compute linear bounds of the second hidden layer and then compute the bounds \([l^2,u^2]\)  by backsubstitution. This process is repeated until the global bounds \([l^K,u^K]\) (a.k.a \([\gamma _{L},\gamma _{U}]\) ) of output layer is produced. To better illustrate the workflow of computing the global bounds, we cite a toy example from DeepPoly [2]} of how it works on a simple FNN with one hidden-layer, two nodes of every layer and the bias set as zero(Figure REF ). In Figure REF , \(x^k_i\)  represents the \(i^{th}\)  Neuron at the \(k^{th}\)  layer. We use \(u^k_i,l^k_i\)  to denote the upper and lower bounds of \(x^k_i\) , and \(u^k_i(\cdot ),l^k_i(\cdot )\)  to denote the upper and lower linear bounds, for example, \(u^2_2(x^1_1,x^1_2)=0.5x^1_2+1,l^2_2(x^1_1,x^1_2)=0\) .
