SOL: Sampling-based Optimal Linear bounding of arbitrary scalar functions
Authors: Yuriy Biktairov
Jyotirmoy Deshmukh
Venue: NeurIPS 2023 poster (2023)
Keywords: Neural network verification
Robustness
Linear bounding

Abstract: Finding tight linear bounds for activation functions in neural networks
is an essential part of several state of the art neural network robustness 
certification tools. An activation function is an arbitrary, nonlinear,
scalar function $f: \mathbb{R}^d \rightarrow \mathbb{R}$. In the existing work 
on robustness certification, such bounds have been computed using human 
ingenuity for a handful of the most popular activation functions. While a 
number of heuristics have been proposed for bounding arbitrary functions,
no analysis of the tightness optimality for general scalar functions has been 
offered yet, to the best of our knowledge. We fill this gap by formulating a concise 
optimality criterion for tightness of the approximation which allows us to
build optimal bounds for any function convex in the region of interest $R$. For
a more general class of functions Lipshitz-continuous in $R$ we propose a 
sampling-based approach (SOL) which, given an instance of the bounding problem, 
efficiently computes the tightest linear bounds within a given $\varepsilon > 0$ 
threshold. We leverage an adaptive sampling technique to iteratively build a set
of sample points suitable for representing the target activation function. While 
the theoretical worst case time complexity of our approach is
$O(\varepsilon^{-2d})$,
it typically only takes $O(\log^{\beta} \frac{1}{\varepsilon})$ time for some $\beta \ge 1$ and is
thus 
sufficiently fast in practice. We provide empirical evidence of SOL's practicality
by incorporating it into a robustness certifier and observing that it
produces similar or higher certification rates while taking as low as quarter of the time compared to the other methods.

SOL: Sampling-based Optimal Linear bounding of arbitrary scalar functions
Yuriy Biktairov
Jyotirmoy Deshmukh
University of Southern California (biktairo, jaee.ddesssoe jdeshmuk) Gusc. edu
Abstract
Finding tight linear bounds for activation functions in neural networks an essen- tial part of several state of the art neural network robustness certification tools. An activation function an arbitrary, nonlinear, scalar function f IR + R. the existing work on robustness certification, such bounds have been computed using human ingenuity for a handful of the most popular activation functions. While number of heuristies have been proposed for bounding arbitrary functions, no analysis of the tightness optimality for general scalar functions has been offered yet, to the best of our knowledge. We fill this gap by formulating concise optimality criterion for tightness of the approximation which allows us to build optimal bounds for any function convex in the region of interest R. For more general class of functions Lipschitz-continuour in R we propose sampling-base bbaedd approach (SOL) which, given an instance of the bounding problem, efficiently computes the tightest linear bounds within given ε > O threshold. We leverage an adaptive sampling technique to iteratively build set of sample points suitable for representing the target activation function. While the theoretical worst case time complexity of our approach is ( it typically only takes Ologg 1) time for some 3 > I and is thus sufficiently fast in practice. We provide empirical evidence of SOL's practicality by incorporating it into a robustness certifier and observing that it produces similar higher certification rates while taking as low as quarter of the time compared to the other methods.
Introduction
Neural networks (NNs) have become the pervasive machine learning model for various tasks including classification, function approximation, and generative modeling. In typical classification task using supervised learning, the NN given as input a nigh-dimensional data object (such as an image or a large piece of text) and its corresponding class label (from a finite set of labels). The NN classifier deemed robust if for reasonable perturbations to a given input, the predicted label for the input does not change. It has been previously shown that NNs susceptible adversarial perturbations in inputs. These approaches take advantage of the non-robustness of a trained NN; for example, in image classification, by injecting a small amount of noise imperceptible human, the predicted label for an input image can be changed. In safety-critical applications where object detection and classification may be used input for decision-making (e.g., in autonomous mobile robot), such misclassifications can have safety implications. Similarly, when used in a function approximation task (e.g., for regression), is important that small changes to the NN input do not lead large changes in the predicted output. Checking the robustness of neural networks has thus emerged as an important problem in recent years 28
23
number of different approaches have been developed for certifying the robustmess of the neural networks with popular architectures. The majority of these approaches fall into one of three overar-
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

ching categories: the algorithms using approximate calculations in order to speed the certification time up
25
33 23 [3]. the algorithms using exact calculations via some sort of constraint solver
10
11
29
and the algorithms taking advantage of both techniques [26] 24
30 32)
The most
practical approaches which can handle the largest neural networks are usually built on top of the approximate techniques and rely heavily on being able to find tight linear bounds of non-linear scalar functions.
A lot of research has been dedicated to increasing the quality of linear bounds of the most popular
activation functions such as ReLU sigmoid. hyperbolic tangent and others 32 25 34 But not much has been done in the direction of efficiently bounding arbitrary functions. This led to the current situation where a number of perfforin novel activation functions like GeLU [1] and Swish[22] and the neural networks containing them are effectively uncertifiable due to their linear bounding not being supported by the robustness certifiers.
36 25
To this end we introduce a novel approach for efficient optimal linear bounding of general scalar functions. The considered optimality with respect to the conventional measure of tightness- the volume of the discrepancy region between the function and the linear bound. At the core of our approach is a simple idea of leveraging the Lipschitz continuity of the activation function in order to build a near-optimal bound using appropriately constructed sample of points from the region of interest. The major contributions of the presented paper are the following
we analyze the problem of optimal linear bounding of scalar functions and derive simple optimality criterion which allows us to optimally bound arbitrary functions convex or concave in the region of interest;
we propose a novel optimal linear bounding approach called SOL which unlike any other known method efficiently produces linear bounds which are arbitrarily close optimum for the general class of functions Lipschit in the region of interest;
we show that using SOL for the robustness certification is viable by benchmarking its performance on synthetic dataset of bounding problems and by comparing it to the known alternatives on real neural network certification tasks.
2
Related work
To the best of our knowledge, there are only two approaches presented in the literature capable of producing linear bounds for arbitrary scalar functions. The first one is implemented inside the AutoLiRPA (35 framework for the robustness certification. The second is the LinSyn 20 approach.
AutoLiRPA supports bounding several general binary operations: addition, subtraction, multiplication and division. It also capable of bounding compositions these operations, which thus allows it to bound an arbitrary function decomposable into sequence of some elementary operations. While this class of functions is indeed quite general and contains all the novel activation functions, the decompositional nature of the approach inevitably leads to producing bounds which are significantly suboptimal.
The authors of LinSyn, on the other hand, use an idea reminiscent of what we propose in this paper: they take uniform sample of points from the region of interest, find the tightest bound for this finite set of points and then repeatedly adjust the obtained bound until its soundness is successfully certified by an SMT solver 4 Unlike SOL, however, their approach does not provide any optimality guarantees for the generated bounds and heavily relies on the use of the SMT solver, which might limit the method's s running time performance.
We provide an empirical robustness certification performance comparison between these two ap- proaches and SOL in section 5.3 of this paper.
3 Optimal linear bounding
In this section we give formal definition of the optimal linear bounding problem and derive the optimality criterion. We then use this criterion to give a recipe for the optimal linear bounding of functions which are convex or concave in the region of interest.

---

((x)
convaf(a))
f())
[[s())]
f(a))
Figure I: Linear bounding example.
Figure 2: Optimal bounding criterion.
Without the loss of generality from now on we are only going to discuss the optimal upper bounding. The lower bounding problem can easily be reduced the former by negating the function and then negating the values of the bound.
Definition (Optimal linear bounding problem). Given a function f R - R defined in a convex bounded region R C RR, aim to find such linear function g(æ²)
b, that
1. g(x) upper-bounds the original function: g(x) 2 f(x) for every R:
2. the volume Y [9(x)) = JZER
dx the discrepancy region between the function graphs is as small possible: 9opt (x) arg g()) TTa-aaa 2+b (l)())).
The value of the discrepancy volume defined above is the conventionally accepted measure of tightness of the linear bound and is known to correlate well with the performance of robustness The notion is illustrated in the figure I We are going call any linear function satisfying the two conditions an "optimal bound". Upon detailed examination it can be noted that the discrepancy volume may be expressed in a more convenient form.
certification [36 25 13
Proposition For any function g(z) upper-bounding the target function f(x) the equality [[g[=)] = g(xx) - V + C' holds, where V = SaER dæ is the volume of R, Ie = V Jae xdæ is the region's 's center of mass and C is constant independent of the function g(x)).
1
Proof. Indeed, for any upper bound g(æ) the integrand in the definition of [[((x)] can simplified
[s(x) - g(x))((()))) dx =
(g(#)) f(x)) dæ
g(x)dxx + C,
JEE
ØER
ZER
where
Smen (x)dx does not depend on g(ar). The straightforward integration
J g(x)dx
(a + b)dx a x.V + bV =g(x)) V
MER
JER
then yields the desired expression.
Since R is convex, H. always lies in R. Then, g(x) being an upper bound for the target function gives us lower bound on the discrepancy volume [[g(=)] = -x.) +C' 2 V-f(x,) (=f(x)++C +C. For a concave target function f(x:) this inequality concludes the problem of optimal linear bounding. Indeed, for concave target function there always exists such an upper bound g" (xx)) that g" (xx) = f(x)). This
upper bound is then optimal due to the inequality.
Now, conv I' ((x) be the upper convex hull the function ((x) in the region R a function of representing the upper boundary of the convex hull of f's graph as is illustrated in figure 2 by the dashed line. The following then holds for arbitrary, possibly non-concave, functions / (ax).
Proposition linear function g(xr) hoper-bounds the target function f (a) in the region iff upper-bounds the convex hull cony R fxx of the target function in the region R.

---

Proof. If ((x) does not bound the target function f(x), then it does not upper-bound its convex hull either since conv R f(x) 2 f(z) for every æ E Conversely, if ( does not upper-bound the convex hull there is at least one point mo such that / (To) >> g(zo). Otherwise the whole convex hull would lie below the linear function.
This proposition allows us to reduce the problem of bounding arbitrary function ((x) to bounding a concave function conv ʃ (x). This reduction combined with the reasoning suggested for the concave target functions above gives us the following concise optimality criterion for the linear upper-bounding problem.
Theorem 1. A linear upper bound g(z) of a function f(x) is optimal iff g(xc) conv R fxx))
Figure 2 Jillustrates the criterion showing a portion couvR R f(xee) distinct from f()) with a dashed line. Theorem iimmediately provides us with two efficient optimal bounding algorithms for concave and convex functions.
For concave and differentiable f(x) the convex hull conv f fx) coincides with the function itself. Finding an optimal upper bound is then simply a matter of finding a hyperplane tangent to f((x) at the center of mass Ze. This only takes (dd operations if we have the access to the function's 's gradient.
For convex (ax) and R being polytope one can show that only the values f(x;) of the function in the vertices 2; of It matter. The problem of upper-bounding then essentially becomes discrete linear program with d variables and [x:3] constraints. In the typical case of a cubic this means having (22) constraints which can still be viable for reasonably high- dimenssinni functions. The process of solving this kind of linear programs efficiently is discussed in details and analyzed experimentally in sections 4.1 and 5.1 of this paper.
The two algorithms combined form a complete optimal upper- and lower-bounding suite to be used with arbitrary convex concave differentiable activation functions. Optimal bounding of both sides of the function can be done exactly in this case with no need for approximations. Among the suitable popular activation functions are: softplus 377, ELU and SELU[12 with appropriate choices of parameters. The approach can easily be generalized to piecewise-difiee differentiable convex or concave functions as well. This extends the applicability to such functions as [eLU[[9], leaky ReLU]16). SELU with a wider set of suitable parameters and the max pooling function (although not an activation function, also frequently used in neural networks).
4 SOL
In order to solve the problem of optimal linear bounding for the functions which are neither convex nor concave in R we introduce the SOL approach: Sampling-basee based Optimal Linear bounding algorithm. In contrast to the approaches discussed above where an exactly optimal solution could be found in finite time, here only aim synthesizing target of the optimum in terms of the discrepancy volume. In return, we get an algorithm which efficiently upper-bounds a very general class of activation functions functions Lipschitz-continuous in R. Some of the popular non-convex and non-concave Lipschitz-continuous activation functions are: Sigmoid, Tanh, GeLU [ Swish 22
upper bound within certain required tolerance
Log Log(6) Mish [8] ELU and SeLU with a >>
The key idea of SOL is to perform multiple iterations alternating between constructing an appropriate sample of points S C R and solving the discrete version of the optimal bounding problem for the constructed sample until properly adjusted solution of the discrete problem meets the required accuracy target. During these iterations the generated sample progressively adapts the target function increasing its density near the points where the optimal bound "touches" the function which allows for higher accuracy. Figure illustrates the adaptive nature of these iterations by showing sample S in the first 3 iterations of SOL run.
4.1 Discrete problem
We define the discrete linear upper-bounding problem mentioned above in the following way.
Definition (Discrete optimal linear bounding problem). Given a set of points S = {(z,, y:)} where ar E d and ya E IR, and a center point Be E Rd. we aim to find such linear function g(æ)) = a x. b.b. that

---

1. g(z) upper-bounds the sample points: g(x,) 2 Ys for every (z; V.) S;
2. the value of the bound at De is small as possible: gopt arg g(=)= T a+ (Ze).
Solving an instance of this problem with appropriately set ae can be seen as relaxing an instance of the continuous optimal bounding problem to only contain a finite subset of the constraints while maintaining the original objective.
Clearly, the discrete problem is an instance the linear programming (LP) optimization problem. It can be solved quite efficiently with one of many known approaches. We benchmark several popular general-purpose LP solvers as well as our implementation of two specific algorithms in the section 5.1 of this paper. It is worth noting that when the number of variables is bounded one solve an LP in time linear in the number of constraints [7]. For our particular application means that an instance of the discrete optimal linear bounding can be solved in O([S)) for fixed d, which makes it as easy asymptotically as simply drawing the sample.
4.2 Soundness conditions
Relaxing the constraints of an optimization problem generally leads to solution which is not feasible with respect to the original constraints. In case of SOL we hope reason about the soundness of the bound in a continuous region based on the values of the target function in finite set of points. One way enable such reasoning limit the analysis functions conforming to a certain smoothness criterion.
If the target function is "smooth enough" then given subregion U C R, a point x within it and
slope a of the solution to the discrete problem (s(x) there exists a certain threshold AUU f) > 0 such that big enough gap 9s(x*) f(xx) > A implies gs (ar) > f(x) for every
U. Having
such a threshold allows us to obtain a sound upper bound in U by appropriately shifting gs (æ) up if
needed. We analyze two different classes of target functions and show how the value of the threshold
A depends on the parameters of the region and the smoothness of the target function.
4.2.1 f is Lipschitz continuous in R
Let Li > 0 be an upper bound on the Lipschitz constant of the function f, meaning that
(x()) - fxx2) s L1 x1 x2l
for 21.22 E R. Let g(x) = a. x: b be the solution of the discrete problem associated with particular sample of points C containing x=" can be shown that |a| S L1 - G, where G is a geometric factor determined by the shape of the region R and the structure of the initial sample S. Importantly, in one die-dimensional case (e.g.. all traditional activation functions) the factor exactly I for any region and any reasonable sample S. For an arbitrary point æ E U the following inequalities then hold
g(x) - f(x)= 2 g(x") f(x)) |a|| x01 Li x x= 2 g(x") f(x") - (1 +1+ G)
where 8 supzEU xx" - xe| is the "width" of the region. Hence, as long as
g(x²) f(x A (1 (1
the bound is sound for the whole subregion U.
The connection between |a| and L1 originates from the following inequality which is based on g(æ=) being a discrete upper bound for the L1 Lo-ooniinuuus function f(x)
x
24
lal < L1
inf sup V (2.),.eeS.)
G(R,
x
xx
where the infimum is taken over all finite subsets {a;} C R such that the center o. lies in the subset's convex hull and all unitary vectors The factor G(R, S) depends on the geometry of the target region R and on the structure and the density of the sample S. Note that this factor decreases monotonically when add new points to the sample, reaching certain limit G(R) when S = R.
5

---

Algorithm Adaptive SOL
- 50
[U+] uniform partition into cubes of size 8
Ze + center of mass of R
K <
do
{U:)
(PLITCUBES((U)).. K)
a. <- SOLVEDISCRETE( (#; ((x;)),,
{A;} <- GETGAPTHRESHODL.S (UU)): ( a.
K < {i]A: (95(z))
- f(x))) > }
while K o 0
= b+ V
return a.
4.2.2 V f is Lipschitz continuous in R
be upper bound on the Lipschitz constant of the gradient V of the target function. Given a bounded R this immediately implies the Lipschitzness of the target function f itself. Let L be upper bound the Lipschitz constant of f-
Let La
The functions fx) of this type behave "approximately linearly" near the sample point x=. denote the linearized version of fx) as fi(æ) f(x") + ff(x") (æ xx), the inequality [f())
- [a)) S L2 x x²12 holds for an arbitrary point æ.
As consequence, the following inequality holds for the discrete bound g(x)
g(x) f(x) 2 g(xx) - f(x') f((*)) Vg(l")]6 L20².
Combined with the previous condition which is still applicable in this setting the inequality gives us the following expression for the gap threshold
min((1 + G)L,, Vf(=*) Vg(2= L207)
Intuitively, the introduction of this additional condition matters the most near the "touching points" between the function and the optimal linear bound. Around these points the gradients the expression are going to be close to each other, so the second term in the minimization will allow for significantly smaller gaps than those of the first term.
4.3 The algorithm
Now we have all the parts necessary to define the SOL algorithm as well as its simplified version.
At all times during the runtime of the algorithm
going to maintain the partition of the original region into collection of cubes U, = {x x Rd \ - 4 x < x + +4, with side lengths of 27, such that R S U. U,.
U. The sample 5 = { we are going to use for the formulation of the discrete problem at particular iteration is going to consist of the center points of these cubes.
Given the solution gs gs(x)) of the discrete problem we are going to calculate the gap thresholds A; associated with all (U;, æx) pairs according to one of the soundness conditions. The shift of value 17 = max; ( (gs (2.) - fxx)))] can then be applied to the discrete bound to obtain a bound sound in the entire target region R: (z(x)) gs(x) 17.
Importantly, since (x) is solution of the relaxed problem its discrepancy is a lower bound for the optimal discrepancy. On the other hand, (x:) is a feasible bound for the continuous optimization problem, hence its discrepancy an upper bound for the optimal discrepancy. Combining these
two observations with the proposition we can conclude that (x) is within 17 V of the optimum discrepancy volume-wise. The value of The shift can directly controlled by varying the density of the
cubes the partition since 7) S max; A, & (1+G)T, max; 6g where = supzeU, |æ:-œr a| = Lvd is the size of cube U..

---

This observation inspires the uniform version of SOL. It only performs single iteration with uniform partition of R into equally-sized cubes of size 6 I [I+G)L,V where is the desired accuracy target. While being attractively simple this version turns out to be too slow in practice. Indeed, in typical case would only need such high density of sample near the touching points of the optimal upper bound. Having sampled the whole region with a uniform density we end up with an intractably big discrete problem.
The solution to this issue to start with a sparse regular lattice of cubes as a partition and gradually change it according to the values of A, obtained after solving discrete problem. Specifically, we would want increase the density of cubes such subregions U/ that
A. (gs( f(x))) V
(I)
as these regions require too big of shift in order to guarantee the soundness. In order locally increase the density we split the target cube U; into 2 evenly spread cubes of half the size of the original. An example of the split procedure can be seen figure This modification concludes the adaptive version of the SOL algorithm algorithm
4.4 On the complexity of SOL
Note that the cubes of size 8 < 5, [I+G)L cannot violate the inequality ① and, consequently, never split. The maximum total number of points in S achievable until the accuracy target is met is then defined by the critical cell size d. as Ne = 0 gg O(E where we only track the complexity's dependence on e. Given that each iteration increases the sample size by at least 1 and solving the discrete problem as well as updating S can be done in ((S)), we can bound the worst case complexity of SOL by ( O(==((22)) 22).
While this runtime bound may seem excessively conservative, there exist linear bounding problem instances which SOL is going take similar O(=->) amount of time for certain values of ry. One particular example is bounding a linear function, where SOL. would have to split all cubes down to the size 8, or 8% V L depending on which soundness criterion used. This example would
give us the lower complexity bounds of 2(= and A(s ś) respectively just for the final iteration of SOL alone.
The situation is quite different, however, for many practical instances of the linear bounding problem. In typical scenario the optimal bound Sopt (z)) only "touches" the target function f(x) separate points and f(x has nondegenerate Hessian in each of them. The high sample density of 8- or (6")=dis is then only necessary in the vicinity of these touching points. Taking the continuous limit one can show that verifying an s-optimal linear bound in this case only requires Ollog _) points in the sample when using the L2 soundness condition. This logarithmic estimate of the number of points required for the verification is the primary reason we introduce the L2 condition. Provided that the adaptive sampling procedure "finds" the optimal sample "fast enough" by splitting the cells in appropriate regions, one can hope achieve the total running time of (log { for some 8 2 1 in arbitrary dimensionality d. We show that this is indeed the case in practice in section 5.2 where we measure the running time of SOL in a synthetic dataset of bounding problems involving popular activation functions.
5 Evaluation
In this section we study the performance of SOL empiricaly by measuring its runtime and the quality of the linear bounds on synthetic dataset containing randomly generated instances the bounding problem. We also show the practical applicability of SOL by incorporating into AutoLiRPA and measuring its robustness certification performance on several image classification neural network architectures.
5.1 Benchmarking discrete problem solvers
Each run of SOL may involve dozens of invocations of the discrete problem solver. There are numerous different LP solvers that can be used for handling the discrete problem. Their performance
7

---

Ma(())
Gurobi
Scipy
2.0
1D bisect
Megiddo
1.5
t. t.ms
2=(#)
0.5
ss(e)
92
0.0
1000
10
100 10-7
10
1001 10
1000
10
10
E
Figure 3: Adaptive sampling.
Figure Runtime of LP solvers.
Algorithm 2 ID bisect algorithm
max + (in
<-0 0
while - 1>> do
1+u m 2
Imax max
m-
24 Fe
Rc 1o
Tmax f max
U--m m
34 > 1e
F1-Be
if Tmaa Imax then
u x m
Lmao trmax
else
+ m
end if
end while return
known to vary drastically between the different algorithms used at their and the distributions of problems targeted. We benchmark 4 approaches on synthetic dataset of the discrete problem instances to find the one most suitable for our applications.
The dataset was generated by sampling 2000 discrete problem instances for each of the following activation functions: GeLU, Log Log, Swish. We use these three activation functions they have become quite popular in recent years(8) 21 22 27 but don't have optimal linear bounds handcrafted for them specifically. Each problem instance has the region's 's R [1,r] boundaries sampled uniformly from the [-2, 2) interval and contains 500 points sampled uniformly from the region. This number of points was chosen for the experiments corresponding to the typical final sample size of SOL in application to the robustness certification. Although the uniform distribution points might not give an accurate representation of the samples encountered by SOL, the difference between the running times the solvers is pronounced enough to ignore this inaccuracy.
The approaches we investigate include two well-known optimization libraries capable of solving LP: Gurobi and SciPy [3] The other two approaches are our implementation Megiddo's linear-time algorithm[17 17 and a specific bisection procedure which takes advantage of the linear ordering of the sample points in the one-dimensional case. The idea behind the bisection procedure is that particular g(x.) h guess on the objective value can be validated in linear time by checking whether

---

300
LinSyn
regression
15
-10
10
(a) Running time comparison.
(b) Optimality gap.
(c) Running time in log-scale.
Figure 5: SOL benchmarking results.
(Zey h)) lies inside the convex hull of the sample or not. A pseudocode implementing the procedure is shown in the algorithm 2 Note that this procedure is only applicable for d = 1.
The benchmarking results shown in the figure The plot shows the dependence of solver's 's average runtime across all the problems in the dataset on the targeted optimality tolerance E. Interestingly, both general LP solvers show slower performance than the alternatives, which highlights the specific nature linear programs encountered by SOL. Both their running times seem almost independent of the optimality tolerance required. One possible reason would be the presence of significant overhead of setting sophisticated optimization routines.
The results of the Megiddo's 's algorithm show expectedly saturating runtime dependence supposed find the exact solution in O([S)) time. The only dependence on e is due to the presence of an early stopping upon reaching the target accuracy. The bisection approach shows the best results with a significant margin along the entire bandwidth asymptotically worse than that of the Megiddo's algorithm O !) O(1) the constant factor of the latter appears to too high for it to outperform the bisection procedure at practical values of optimality tolerance. Taking into consideration this result, the bisect algorithm the discrete problem solver in all of the remaining evaluations.
since it's
values. Although its dependence on is
5.2 Benchmarking SOL
In order to study the average runtime and the bounding performance of SOL on the linear bounding problem instances encountered in practice we generate another synthetic dataset. This time sample 1000 instances of the problem for each of the same three functions. The region bounds for the problems are sampled uniformly from the [-3,3] interval.
We compare the SOL's performance on this benchmark to the performance of LinSyn n 20 as their approach also aims build bounds close the optimal ones. does not, however, give any guaranties on the optimality gap and does not have parameters controlling the accuracy target. Its metrics are thus independent of the optimality target
Figures [egures[5ajadae and [5]]shww the average runtime and the average optimality gap respectively as functions of the optimality target e. The optimality gap is estimated by taking the most accurate linear bound produced by SOL and shifting in down appropriately to obtain a tight underestimate of the optimal value of objective. Notably, there appears to be wide range of accuracy targets & [10-7 10 where SOL gives better average bounds while taking less time then LinSyn. For example, at & = 10 the SOL bounds are 10 times tighter than these of LinSyn while the runtime is 3 time shorter.
While figure 5a already hints on the form the asymptotic dependency of the running time on the accuracy target e, we make the analysis clearer by investigating the relation between log and log log instead. This dependence shown in figure 5c turns out to be quite close to linear with the slope of a z= 2.0 estimated by the linear regression. The linear relation with this slope corresponds to the complexity of Olog" !) and agrees with our hypothesis on the typical running time of SOL.

---

Table 1: Certification performance comparison. The upper number is the fraction of properties certified, the lower is the estimate of the running time (with std < 3%). Neural network architectures are denoted by the activation function used and their respective depths.
Dataset
Approach
MNIST [13]
CIFAR 14
Ge
Log Log 41
Swish 41
GeLU 51
Log Log 51
Swish 51
0.01
0.0
0.34
0.0
0.59
0.03
AitA
3
230s
23s
9s
700s
50s
35s
0.72
0.23
0.76
0.31
0.69
0.35
LinSyn 20
430s
450s
410s
580s
360s
550s
0.73
0.23
0.76
0.31
0.69
0.37
SOL (5 10-3)
100s
96s
100s
150s
87s
150s
0.73
0.24
0.76
0.31
0.69
0.37
SOL (5 10-5)
140s
140s
130s
190s
95s
170s
SOL (5 = 10-")
0.73
0.24
0.76
0.31
0.69
0.37
310s
300s
290s
400s
160s
340s
5.3 Robustness certification performance
Finally, we study the performance of SOL as part of a robustness certification pipeline by using SOL as an activation function bounding subroutine in the AutoLiRPA (35] framework. We follow the protocol introduced in LinSyn and compare the robustness certification performance of our approach against that of LinSyn and that of AutoLiRPA with their default decompositional linear bounding attempt the robustness certification of 6 image classification
approach. In these experiments neural networks with convolutional architectures. Three of them were trained on MNIST while the other three CIFAR[14].
Each certification task consists of proving that any -bounded perturbation a particular correctly- classified input image retains the same label prediction as the image itself. For each neural network we randomly sample 100 images from the test part the dataset, filter out the misclassified ones and attempt to certify the remaining images. We use perturbation magnitude of 8/255 for networks trained on MNIST and 1/255 for CIFAR networks as is common in the literature. We have also conducted experiments with larger perturbations of 2/255 and 4/255 for the CIFAR models, but found the performance too weak for a proper analysis. Qualitatively the results for these larger perturbations were similar to the presented results.
The results are shown I Both LinSyn and SOL show much higher certification rates compared to Auto since they aim at producing tight linear bounds while the default bounding procedure of AutoLiRPA systematically ends up with loose bounds.
The SOL variants show similar or slightly higher certification rates than LinSyn while taking as low a quarter of the time. Lowering the accuracy target down from 10-5 noticeably increases the runtime of SOL without any improvement to the success rate of the certification. Notably, although not depicted in the table, the average sizes of the final bounds on the logit outputs of the networks only improve by a fraction of 1% upon changing E from 10 to 10 which is negligible compared the improvement of x2-x5 from switching from AutoLiRPA to either LinSyn or SOL.
Conclusion
In this paper we have introduced SOL sampling- based approach for finding linear bounds arbitrarily close to optimum in terms of tightness for Lipschitz-continuouss functions. We have shown experimentally that despite having rather pessimistic theoretical running time guarantees, it is quite fast in practice. As part of neural network robustness certification framework increases the certification rates and takes significantly smaller time compared to the alternative approaches.
10

---

Acknowledgments
The authors would like to thank the anonymous reviewers for their feedback. This work was supported by the National Science Foundation through the following grants: CAREER award (SHF-2048094), CNS-1932620, CNS-2039087, FMitF- 13, CCF-SHF-19326200, funding by Toyota R&D and Siemens Corporate Research through the USC Center for Autonomy and AI, an Amazon Faculty Research Award, and the Airbus Institute for Engineering Research.
The authors would also like to acknowledge early discussions with Prof. Chao Wang and Dr. Brandon Paulsen that influenced this work.
References
[1] M. Alzantot, Y. Sharma, Elgohary, B.-J. Ho, M. Srivastava, and K.-W. Chang. Generating natural language adversarial examples. arXiv preprint arXiv:1804.077 079988 2018.
[2] D.-A. Clevert, T. Unterthiner, and Hochreiter. Fast and accurate deep network learning by exponential linear units (elus). arXiv preprint
2015.
[3] L. Engstrom, B. Tran, D. Tsipras. L. Schmidt, and A. Madry. Exploring the landscape of spatial robustness. In International conference on machine learning. pages 1802-1811. 1. PMLR, 2019.
[4] S. Gao, S. Kong, and E. M. Clarke. dreal: An smt solver for nonlinear theories the reals. In Automated Deduction-CADE-24 24th International Conference on Automated Deduction, Lake Placid, NY, USA, June 9-14, 2013. Proceedings 24, pages 208-214. Springer, 2013.
[5] T. Gehr, M. Mirman, D. Drachsler- P. Tsankov, S. Chaudhuri, and M. Vechev. Ai2: Safety and robustness certification of neural networks with abstract interpretation. In 2018 IEEE symposium on security and privacy (SP), pages 3-18. IERE, 2018.
[6] G. S. d. S. Gomes and T. B. Ludermir. Complementary log-log and probit: activation functions implemented in artificial neural networks. 2008 Eighth International Conference on Hybrid Intelligent Systems, pages 939-9499999999 992.9.222222 IEEE, 2008
[7] Gurobi Optimization, LLC. Gurobi Optimizer Reference Manual, 2023. gurobi com
[8] D. Hendrycks and K. Gimpel. Gaussian error linear units (gelus). arXiv preprint arXiv: (6066084155 2016.
[9] C. Kanbak, S.-M. Moosavi- Moosavi-Dezfooli, and Frossard. Geometric robustness of deep networks: analysis and improvement. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 4441-4449, 2018.
[10] G. Katz, C. Barrett, D. L. Dill, K. Julian, and M. J. Kochenderfer. Reluplex: An efficient smt solver for verifying deep neural networks. In Computer Aided Verification: 29th International Conference, CAV 2017, Heidelberg, Germany, July 24-28, 2017, Proceedings, Part I 30, pages 97-117. Springer, 2017.
[11] G. Katz, D. A. Huang, D. Theling, K. Julian, C. Lazarus, R. Lim, P. Shah, S. Thakoor, H. Wu, A. Zeljié, al. The marabou framework for verification and analysis of deep neural networks. In Computer Aided Verification: 31st International Conference, CAV 2019, New York City, NY, USA, July 15-18, 2019, Proceedings, Part 31, pages 443-452. Springer, 2019.
[12] G. Klambauer, T. Unterthiner, A. Mayr, and S. Hochreiter. Self-normalizing neural networks. Advances in neural information processing systems, 30, 2017.
[13] C.-Y. Z. Lyu, L. Weng, L. Danicl, N. Wong, and D. Lin. Popqorn: Quantifying robustness of recurrent neural networks. In International Conference on Machine Learning, pages 3468-3477. PMLR, 2019.
[14] A. Krizhevsky, G. Hinton, et al. Learning multiple layers of features from tiny images. 2009.

---

[15] Y. LeCun, L. Bottou. Y. Bengio, and P. Haffner. Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86(11):2278-333) 2324, 1998.
[16] A. L. Maas, A. Y. Hannun, A. Y. Ng. et al. Rectifier nonlinearities improve neural network acoustic models. In Proc. icml, volume 30, page 3. Atlanta, Georgia, USA, 2013.
[17] N. Megiddo. Linear programming in linear time when the dimension is fixed. Journal of the ACM (JACM), 31(1):114-127, 1984.
[18] D. Misra. Mish: A self regularized non-man oo-moootonic activation function. arXiv preprint arviv: 08681, 2019.
[19] V. Nair and G. E. Hinton. Rectified linear units improve restricted boltzmann machines. In Proceedings of the 27th international conference machine learning (ICML-I0), pages 807- 814, 2010.
[20] B. Paulsen and C. Wang. Linsyn: Synthesizing tight linear bounds for arbitrary neural network activation functions. In Tools and Algorithms for the Construction and Analysis of Systems: 28th International Conference, TACAS 2022, Held as Part of the European Joint Conferences on Theory and Practice of Software, ETAPS 2022, Munich, Germany, April 2-7, 2022, Proceedings, Part I, pages 357-376. Springer, 2022.
[21] A. Radford, K. Narasimhan, T. Salimans, I. Sutskever, et al. Improving language understanding by generative pre- pre-training. 2018.
[22] P. Ramachandran, B. Zoph, and Q. V. Le. Searching for activation functions. arXiv preprint arXiv: 2017.
[23] Z. Shi, Zhang, K.-W. Chang, M. Huang. and C.-J. Hsieh. Robustness verification for transformers. arXiv preprint arXiv:2002.06622 2020.
[24] G. Singh, R. Ganvir, M. Püschel, and M. Vechev. Beyond the single neuron convex barrier for neural network certification. Advances in Neural Information Processing Systems, 32, 2019.
[25] G. Singh, T. Gehr, M. Püschel, and M. Vechev. An abstract domain for certifying neural networks. Proceedings of the ACM on Programming Languages, 3(POPL):1330 2019.
[26] G. Singh, T. Gehr, M. Püschel, and M. Vechev. Boosting robustness certification of neural networks. In International conference on learning representations, 2019.
[27] V. Singla, S. Singla, S. Feizi, and D. Jacobs. Low curvature activations reduce overfitting in adversarial training. In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 16423-16433, 2021.
[28] C. Szegedy, W. Zaremba, I. Sutskever, J. Bruna, D. Erhan, I. Goodfellow, and R. Fergus. Intriguing properties of neural networks. arXiv preprint arXiv.1312.6199,22 2013.
[29] V. Tjeng, K. Xiao, and R. Tedrake. Evaluating robustness of neural networks with mixed integer programming. arXiv preprint arXiv:
2017.
[30] H.-D. Tran, D. Manzanas Lopez, P. Musau, X. Yang, I. V. Nguyen, W. Xiang, and T. T. Johnson. Star-based reachability analysis of deep neural networks. In Formal Methods-The Next 30 Years: Third World Congress, FM 2019, Porto, Portugal, October 7-11, 2019, Proceedings 3, pages 670-686. Springer, 2019.
[31] Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. Nelson, F. Jones, R. Kern, E. Larson, J. Carey, i. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribciro. F. Pedregosa, van Mulbregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientitic Computing in Python. Nature Methods, 17:261-272, 2020. doi: 10.
[32] Wang, K. Pei, J. Whitehouse, J. Yang, and Jana. Efficient formal safety analysis of neural networks. Advances in neural information processing systems, 31, 2018.
12

---

[33] L. Weng, H. Zhang, H. Chen, Z. Song, C.-J. Hsieh, L. Daniel, D. Boning, and I. Dhillon. Towards fast computation of certified robustness for relu networks. In International Conference on Machine Learning, pages 5276-5285. PMLR, 2018.
[34] Y. Wu and M. Zhang. Tightening robustness verification of convolutional neural networks with fine men-grained linear approximation. Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, pages 11674-1 1667-11681, 2021.
[35] K. Xu, Z. Shi, H. Zhang, Y. Wang, K.-W. Chang, M. Huang, B. Kailkhura, X. Lin, and C.-J. Hsieh. Automatic perturbation analysis for scalable certified robustness and beyond. Advances in Neural Information Processing Systems, 33:1 121-119-144 141, 2020.
[36] H. Zhang, T.W. Weng, P.-Y. Chen, C.-J Hsieh, and L. Daniel. Efficient neural network robustness certification with general activation functions. Advances in neural information processing systems, 31, 2018.
[37] H. Zheng. Z. Yang. W. Liu, J. Liang, and Li. Improving deep neural networks using softplus units. In 2015 International joint conference neural networks (IJCNN), pages 1-4. IEEE, 2015.
13


Reviews:
[{"review_id": "vi9Li04JS2", "reviewer": "Reviewer_Uxf6", "summary": "This paper introduces a new method for finding tight (in the l1 error sense) linear bounds for a general non-linear activation functions. Under their definition of tightness, the authors propose a method for finding optimal bounds for any convex function. For non-convex functions with Lipschitz continuity, they propose a sampling based method that efficiently computes the linear bounds within some error threshold. The performance of the proposed methods are benchmarked in robust certification tasks.", "questions": "Please see Weakness section above.", "limitations": "I do not see any limitations.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "This is an interesting paper tackling an important problem (linear bounding of arbitrary activation functions) in the field. It formulates and solves the problem in a rigorous fashion, and a systematic benchmarking proves the utility of the new method.", "weaknesses": "I do not see particular weaknesses in this paper, but I do have a question out of theoretical interest. The optimality criterion is based on the L1 error of between the bounds g(x) and the activation function f(x). And this L1 error is important in the derivation of Proposition 1. I wonder how the landscape would change if we change the loss function to, e.g., L2 or other meaningful functions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review and for the interesting question.\n\n[significance of optimizing $L_1$ error] Right now we don't see any straightforward way to generalize our approaches to alternative tightness measures.\nMoving from $L_1$ measure to, for example, $L_2$ poses several immediate problems:\n1) As you have already mentioned, Proposition 1 no longer works. What is worse, now the loss function depends on $\\int_R x g(x)$ and $\\int_R g(x)$ in a non trivial way.\nGenerally, there is no analytical expression for these two.\n2) The analogue of our discrete problem is no longer an LP, but a QP, which is harder.\n3) Shifting the bound up to adjust for possible unsoundness of the discrete problem's solution no longer corresponds directly to incrementing the loss function.\n\nOn the other hand, we expect that, at least in 1d, in many situations (when the optimal $L_1$ bound touches $f(x)$ in two distinct points) the optimal bound stays the same for\na wide spectrum of discrepancy measures. Therefore, the exact choice of measure may not influence the performance of the robustness certification too much.\nThen we might as well choose a measure which facilitates the optimization."}, {"review_id": "e9neFYo1DH", "reviewer": "Reviewer_173q", "summary": "This paper focuses on finding tight linear bounds for the activation functions in neural networks, which are used for certifying the robustness of a neural network.\n\nIt aims to provide optimality guarantees for the tightness of the bounds of a scalar function represented by the neural net.\nTwo settings are considered:\n1. for functions that are convex in some region $R$, optimal bounds are obtained through an optimality criterion for the tightness of the approximation in $R$;\n2. for functions that are Lipschitz continuous in some region $R'$, a sampling-based approach is proposed called SOL. Given an instance of the bounding problem and a positive scalar $\\epsilon$, SOL efficiently computes the tightest linear bounds within the $\\epsilon$ threshold.  \n\nThe empirical simulations show that the proposed method SOL typically takes a quarter of the time other methods take.\n\n\n", "questions": "1. What is the average time complexity of SOL, and does it relate better with the observed time? If not, could you discuss what are the limitations / worst-case examples that yield the worst-case time complexity?\n", "limitations": "Please see above.", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "**Novelty to more systematically analyze the linear bound tightness.**\nAlthough I am not working in this area, it seems this paper is among the first to study the tightness of the linear bounds for the activation functions used for neural nets.\n\n**Elegant/Simple approach for the convex setting.**\nThe results for convex functions are very simple, but to me, that makes the approach very nice: it shows that it suffices to estimate the center of the mass, which is a point in the given region from the domain (in practice estimated after sampling), and compute the gradient at that point to obtain a linear bound. \n\n\n", "weaknesses": "**Incomplete results for $L$-Lipschitz setting.**\nWhile I personally enjoyed the overall idea, I think the results for the $L$-Lipschitz setting are incomplete.\nIn particular:\n- these should be more rigorous (see also comment below), and moreover\n- I expected to see a theorem that relates the tightness of the bounds estimated with sampling with the $L$-constant, the sample size, the dimensionality of the problem, etc. (even if this is done for some specific non-convex function). \n\n\n**Writing.**\nThe writing is, in general, easy to follow.\nHowever, parts of the text are relatively ambiguous or less formal. \nFor example, the following should be improved:\n- *relating to motivation*: since robustness certificates seem to be the main motivation mentioned in the paper for linear bounds, explain how the latter gives the former (and if it is the most consuming step)\n- *Prop. 1*. The way it is written, the proposition is not a full statement -- write it more formally to include the previous definition of the volume or point to it when making the statement.\n- *Theorem 1*. Similarly, theorem 1, as stated, is incomplete -- the setting and the assumptions should be repeated here. \n- *Missing proof of Thm. 1*. \n- *the discussions for the $L$-Lipschitz setting should be made much more rigorous.* (and stated as lemmas/theorems).\n- Fig. $3$ should be explained better in the caption.\n- etc. \n\n\n\n\n## Minor\n\n- Abstract: *tightness optimality* is hard to understand; elaborate better as it seems central\n- line 157: *the* sample $S$", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for such a detailed feedback. We will definitely address clumsiness in formulations and unintentional ambiguities.\n\n[average time complexity of SOL] Unfortunately, we cannot provide a rigorous analysis of the average time complexity, since the dynamics of the estimated linear bound\nbetween the iterations seems too complex. However, the proposed intuitive reasoning (lines 252-261) aligns pretty well with what we get in practice (fig. 5a, 5c).\nThe empirical complexity of $O(\\log^2 \\frac{1}{\\varepsilon})$ can be attributed to:\n1) 1d bisect taking $O(n \\log \\frac{1}{\\varepsilon})$ time to solve each discrete problem with $n$ points;\n2) the final discrete problem typically having $O(\\log \\frac{1}{\\varepsilon})$ samples concentrated near the touching point as per the reasoning;\n3) the number of points growing fast enough for the last iteration's runtime to be dominant.\n\n[worst-case examples] Bounding linear functions gives complexities close to our worst-case bound of $O(\\varepsilon^{-2d})$ (lines 246-251). Currently, we don't know whether closer\ncomplexities are achievable on not.\n\n[theorem that relates the tightness of the bounds ...] Such a relation is given implicitly for the uniform (not adaptive) SOL at line 227 by defining the cell size\nin terms of the desired optimality and the Lipschitz constant. We should definitely make this relation more explicit for the sake of clarity.\nFor the adaptive version, again, no rigorous connection of this sort is known to us."}, {"review_id": "dmCSxgpPYb", "reviewer": "Reviewer_hjHu", "summary": "[Context]\nNeural network verification algorithms are based on bounding the activations and outputs of neural networks.\nThis is achieved by propagating linear bounds through the network. \nFor each non-linear operation in the network, it is required to provide linear bounds of the operation: given the range of inputs that the function accepts, give a linear mapping that upper (or lower) bound the function. \nUsing these linear bounds, propagation algorithm like the variants of Crown / Lirpa are able to bound the whole network.\nThese linear bounds of operations are usually manually derived.\n\n[Contribution]\nThe authors define an optimality criterion in the context of convex or concave activation function [Theorem 1].\nThe main contribution of the paper is an algorithm to computationally derive linear bounds for functions that are not convex or concave, as long as they are lipschitz continuous. The algorithm consists in sampling a finite number of points, solving a LP to obtain the optimal linear bound that upper bound the sampled points [4.1] , and using smoothness criterion to adjust these bounds and guarantee that it will be a valid upper bound [4.2]. This procedure can be re-iterated until the desired accuracy gap (to the optimal linear bound) is achieved (4.3)\n\n\n[Opinion of the paper]\nI think that the paper is very interesting and the proposed algorithm can be very useful. The general idea of the algorithm is clearly presented, however some important details are missing which makes the current description of the paper insufficient to reproduce / re-implement it. If my main weaknesses comments ([Missing Description], [Mising Experiment])  and questions ([Computational aspects] [Details about the experiments]) are addressed, I would be happy to increase my Soundness to good/excellent and my General Rating to Weak Accept / Accept.\n", "questions": "[Results of Theorem 1]\nThis might be a result of the definition of \"optimal\", but in the case where conv_R(f) is not differentiable, but only subdifferentiable, then it is possible that we have several functions which are all \"optimal\" at the same time, is that correct?\n\n[Claim in l.136 to l.142] \nUnless I'm mistaken, both Softmax and SELU are neither convex nor concave, so should not be in the list of functions that are directly supported before introducing SOL.\n\n[Claim in l.190 to l.192]\n\"It can always be chosen in such a way that $|a| \\leq L_1$\"\nThis claim doesn't seem justified by either a proof sketch or a reference?\n\n[Computational aspects]\nIt's unclear to me on what plaftorm is every operation run? I assume that Gurobi and Scipy are running on CPU, but what about the other algorithms? Is the implementation of \"1D bisect\" amenable to vectorization?\n\n[Details about the experiment]\nThe paper is missing some information about what type of networks are used for the \"end to end\" test including the method in AutoLIRPA. What activations were those network using? What size are they? How were they trained? Are they from a standard \"NN robustness verification\" benchmark such as the one from VNN-comp?\n\n[Suggestion for future work - Ignore this for the review process]\nIt seems like this method could potentially be improved quite significantly by introducing some sort of smart caching. Most of a time, a network will use the same activation function through the network and it seems non optimal to restart from a uniform splitting of the domain every time, even if certain points (inflection points, local extrema) are much more likely to be important to include in the sample. \n\n[Notes]\n* Throughout the whole paper, Lipschitz is spelled Lipshitz. I did not find any examples of that spelling anywhere so I'm wondering if it's a typo or a different accepted spelling.", "limitations": "The paper would be improved by adding an experiment showing the cost it introduced vs. deriving hand crafted optimal bounds. See the [Missing experiment] comment in Weaknesses.", "rating": 7, "confidence": 5, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "* The problem solved is important, and the proposed solution is interesting. Having had to derive manually convex relaxations of activation functions for some project, I definitely appreciate the value of an automated algorithm to handle this problem.\n\n* The validation of the proposed measure is broad. \nIt is both performed on small problems where internal design choices (picking the LP solver algorithm) can be validated [5.1]  and compared to the most relevant baseline [5.2], but also on more (I assume, see Questions) large scale problems.\n\n* The paper is very well structured, framing first the problem and their work in the context of the existing literature, introducing the framework in which they operate, and then building the description of their algorithm part by part in a logical fashion. While some elements are not fully described (see Weaknesses + Questions), I believe that this can be remedied for the final version.", "weaknesses": "[Definition of the Optimal Linear Bounding Problem]\nThe paper motivates the choice of the \"Minimum volume criterion\" by citing [36, 25, 13]. \nI could not find anything relevant to that point in the Popqorn paper [13]. In the DeepPoly paper [25] and Crown paper [36], the choice of lower bound for ReLU is indeed done by minimizing the area of the relaxation created, but without actually showing that it \"correlates well with the performance of robustness certification\" (except by showing that it performs by the old \"parallel bounds\"), and we know that those choices can be improved (as evidenced by paper like alpha crown that give better robustness results that those choices). In addition, no point is made about this for general functions beyond ReLU.\nI completely agree with the authors that this is probably a good criterion, but if this is to be a basis for the method, I think the paper would be much stronger if this point was more strongly defended. \nSome possible way this could be done:\n- Can we link final robustness results to the volume criterion in some way? Maybe in expectation over some random weights, or at least clarifying that on a restricted type of activation function (convex? Monotonously increasing?), the volume optimal bound can be shown to be dominating other bounds.\n\n[Missing description]\nThe \"1D bisect\" algorithm seems critical for the performance of the method (around twice as fast as the best alternative option), but it isn't described anywhere. The LP to solve need to return a slope vector $\\mathbf{a}$ and an intercept $b$, but the only explanation given is in line 288. and 289 which does not describe how to implement those values.\n\n[Missing experiment]\nOne experiment that is missing is the comparison to a network where we *know* what the correct bounds should be. If you take a network trained with sigmoid or tanh activation, the handcrafted convex hull / optimal linear bounds have been derived and should be tested (in the setting of Table 1). It is expected that the proposed SOL method would lose to that baseline, but knowing the gap in performance would be quite valuable, so that the \"cost\" of not deriving manually optimal bounds is known.\n ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review and the detailed feedback.\n\n[Missing description] We should definitely describe the 1d algorithm in more details. The idea of the bisection is as follows\n1) guess the value $g(x_c) = y$ which we want to optimize;\n2) check its feasibility by finding $a = \\min_{x_i < x_c} \\frac{y - f(x_i)}{x_c - x_i}$ and $b = \\max_{x_i > x_c} \\frac{f(x_i) - y}{x_i - x_c}$ in $O(|S|)$ time;\n3) if $a < b$, the guess is infeasible since to line passing $(x_c, y)$ upper-bounds both argmaxes from the previous step;\n4) otherwise, the guess is feasible and a valid solution would be $g(x) = y + (x - x_c) (a + b) / 2$.\n\nSo, starting with the initial feasible guess of $\\max f(x_i)$ and the initial infeasible guess of $\\min f(x_i)$ we can find the solution within $\\varepsilon$ of optimum in $O(|S| \\log \\frac{1}{\\varepsilon})$ time.\n\n[Missing experiment] Such an experiment would, indeed, be a great addition to the paper and we will gladly include it.\nRight now similar quantities can be estimated by comparing:\n1) Certification rates of SOL with higher eps and the rates of SOL with low enough eps for rates to be saturated (1e-7).\nThe saturated rates should be close enough to the hypothetical rates of returning the exact optimal bounds.\n2) Running times of SOL and the running times of the default AutoLIRPA. For simpler functions decompositional bounding of AutoLIRPA\nshould take time similar to how much calculating optimal hand-crafted bounds takes.\nHigher times for gelu networks are supposedly caused by massive decompositional structure of the function.\n\n[Definition of the Optimal Linear Bounding Problem] Connecting the choice of the tightness measure to the robustness properties would be very appropriate for the paper.\nWe will look into it. Right now the motivation is more practical:\n1) other researchers have tried it and didn't find it pathological;\n2) this exact tightness measure is convenient to optimize.\n\nAlso, we believe that for many practical scenarios the exact choice of a measure (withing a certain spectrum) might no affect the optimal bound much.\n\n[Results of Theorem 1] Correct, any subderivative passing through $(x_c, conv_R(f)(x_c))$ would be optimal.\nE.g., optimal bounds for $-|x|$ in [-1, 1].\n\n[Claim in l.136 to l.142] For SELU we meant to refer to the case when the scale parameter of the linear portion of the function is set high enough for the function to be convex. Softmax, indeed, should not be mentioned there. What should have been mentioned instead is max_pool – which is a convex function of several variables.\n\n[Claim in l.190 to l.192] Thank you for noticing this one! A few weeks ago we gave this statement a little more thought and realized that it's not as straightforward as it had seemed.\n1) Given an arbitrary solution to the discrete problem we can always \"perturb\" it in such a way that the corresponding hyperplane passes through $d + 1$ points $[(x_i, f(x_i))]$ of the sample and the center of mass $x_c$ lies in the closed simplex $[x_i]$.\n2) For the one-dimensional case this is enough to satisfy the statement: $|a| = \\frac{|f(x_1) - f(x_2)|}{|x_1 - x_2|} \\le L_1$.\n3) In higher dimensions it does not necessarily hold as is. Note that for an arbitrary point $x \\in S$ the upper-bounding condition requires\n       $$f(x_i) + a(x - x_i) \\ge f(x) \\ge f(x_i) - L_1 |x - x_i|$$\ntherefore\n        $$a \\frac{x - x_i}{|x - x_i|} \\ge - L_1 \\Rightarrow (-a) \\frac{x - x_i}{|x - x_i|} \\le L_1 \\Rightarrow |a| \\le L_1 / (\\frac{-a}{|a|} \\cdot \\frac{x - x_i}{|x - x_i|}).$$\nThis gives us the general bound of\n        $$|a| \\le L_1 / \\inf_{[x_i], v} \\sup_{i, x \\in S} (v \\cdot \\frac{x - x_i}{|x - x_i|}),$$\nwhere infimum is taken over all simplexes $[x_i]$ such that $x_c$ lies in the simplex and $v$ are all possible unitary vectors corresponding to the direction of the gradient.\nIt's easy to see that that this min-max optimization depends on two things:\n1) the density of points in the sample: higher density gives lower bound;\n2) the geometry of R which determines the limit value of the bound when density approaches infinity\n            $$q =  \\inf_{[x_i]\\subset R, v} \\sup_{i, x \\in R} (v \\cdot \\frac{x - x_i}{|x - x_i|}).$$\nThe geometric factor $q$ is quite important in practice. Intuitively, $R$ having \"acute angles\" on the boundary may make $q$ smaller than 1.\n\nOne example of such situation is in 2d: let $R$ be the triangle {(0, 0), ($w$, 1), (-$w$, 1)}, for which $q \\rightarrow 0$ as $w \\rightarrow \\infty$. If we chose the target function as $f(x,y) = |x|$, the Lipschitz constant is always $L_1 = 1$, but the optimal upper bound is $g(x,y) = wy$ with $|a| = w \\rightarrow \\infty$.\nAnd the estimated discrete optimal bounds might have similarly big |a|.\n\nOn the other hand, $R$ with \"no angles\" – a closed sphere seems to have $q = 1$ in any dimension, so we can guarantee $|a| <= L_1 (1 + \\varepsilon)$ for any positive $\\varepsilon$ if the initial sample size is big enough.\n\nAll in all, these factors only introduce an additional constant factor into the complexity.\n\n\n[Computational aspects] All methods only use CPU. We run the experiments from the paper in a VM with 3 cpus allocated.\nAutoLIRPA runs several instances of the bounding problem solver concurrently.\nThe calculation of two $\\frac{f(x_i) - y}{x_i - x_c}$ arrays in 1d bisect algorithm should be vectorizable. We will look into it.\n\n\n[Details about the experiment]\nThese details should definitely be included in the paper:\n1) We run AutoLIRPA with the simple \"CROWN\" method for bounding.\n2) The networks with \"4l\" in the name have 2 convolutional layers followed by 2 fully-connected layers, the \"5l\" networks – 3 conv + 2 fc.\n3) The activation layer corresponding to the name of the model is present after each main layer except the last one.\n4) All of the networks are trained by optimizing traditional cross-entropy. Similarly to how it was done in LinSyn paper.\n\n[Lipschitz vs Lipshitz] It seems to be a typo in the spell-checker on our side."}, {"review_id": "bvuY0DfWyV", "reviewer": "Reviewer_GjZS", "summary": "This paper describes an approach for finding a linear upper bound of scalar functions that are Lipschitz. To find a bound that approximately optimizes discrepancy with the target function, it samples points to construct LP instances whose solutions bound the corresponding \"discrete\" bounding problem, and continues until finding a bound that matches an optimality target. This technique can be applied as a primitive in neural network verification routines, and the evaluation shows that when it is used in place of prior work, it leads to faster verification times with no penalty on accuracy for MNIST and CIFAR10 models.", "questions": "1. Please provide information on the missing details discussed in (3) from above.\n2. If you have done experiments with different robustness radii or architectures, but did not report them, please describe the results.", "limitations": "Aside from the questions about experimental details discussed above, limitations were adequately addressed.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "This work addresses a specific problem with well-known applications, introducing a new technique that improves measurably on prior work. Notably, while the performance of this approach when used in robustness verification is quite a bit better than the recent best-in-class techniques, it does not seem to impose any additional restrictions that would limit its applicability to networks with different architectures or activation functions, and it does not seem to degrade the precision of the verifier.\n\nUnlike the most closely related work, this approach provides a parameter that can be tuned (i.e., the approximation threshold) to trade performance for precision. While this isn't explored too much in the paper, this flexibility could potentially be used by verifiers to adaptively scale to larger problems by offering progressively weaker guarantees.\n\nThe writing is clear and understandable, although the main contribution takes some time to get through. The writing might be improved with a  more focused exposition of the algorithm. Currently, some of the intuition is given in section 3, some at the beginning of section 4, and the full algorithm is finally presented in 4.3. Soundness conditions and a brief diversion into LP are discussed in between these sections, which was distracting.", "weaknesses": "While the performance gains are apparent (modulo some experimental concerns discussed below), the significance of this approach may be somewhat limited, as it targets a specific application and offers incremental improvements, not new capabilities. Nonetheless, this will be of interest to those working on neural network verification.\n\nThe experimental analysis is limited, and focuses more on microbenchmarking LP solvers and SOL's isolated performance than showing gains in its most useful application, robustness verification.\n1. Verification results are shown over two models, configured with three different activation functions. It's interesting to see the differences across activations, but this doesn't show consistent improvements as architectures grow deeper or wider. Additionally, MNIST results don't offer too many insights at this point, as prior work on certified training yields results than are difficult to improve on.\n2. The experiments don't vary the robustness radius, so we can't determine how this impacts runtime or certification rates. In particular, it would be good to see that this approach's reliance on sampling doesn't break down for stronger guarantees. The 1/255 radius for CIFAR is not really standard anymore; although not a decisive argument, note that the [certified robustness community leaderboard](https://sokcertifiedrobustness.github.io/leaderboard/) doesn't have papers in this category, as most evaluate 2/255 or 8/255.\n3. Important details of these models are not given: how many parameters, which layers, and how were they trained? Importantly, were these trained using techniques that would make them efficiently certifiable, or adversarial training, or something else? This will have a significant effect on the results. If these models were not trained for certification, then the experiments should include some models that were, as AutoLiRPA's default bound and LinSyn are likely to show improvement, and we would hope to see a comparable improvement in this approach.\n\nThe limited experiments make it difficult to judge the significance of this work, which is positioned as a drop-in replacement for existing bound approximation methods within verifiers.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review and the feedback.\n\n[exposition of the algorithm] Thank you for the suggestion. We will try to rearrange some of the paragraphs to make the narrative more coherent.\n\n[information on the missing details discussed in (3)]\nHere are the details. We will make them more clearly specified in the paper:\n1) All of the models were trained by optimizing a simple cross-entropy. No robust training techniques were used. This is the setup\nLinSyn paper uses for their experiments, so we chose a similar approach.\n2) The networks with \"4l\" in the name have 2 convolutional layers followed by 2 fully-connected layers, the \"5l\" networks – 3 conv + 2 fc.\n3) The activation layer corresponding to the name of the model is present after each main layer except the last one.\n\n[experiments with different robustness radii] We conducted additional experiments running SOL on CIFAR with perturbation radii 2/255, 4/255 and 8/255.\nSimilarly to the results for 1/255 presented in the paper the certification rates for any other radius do not depend on whether\nthe optimality target is set to 1e-7, 1e-5 or 1e-3. This means that for bigger radii the saturation still happens somewhere above 1e-3.\nThe certification rates drop significantly with increased radius as can be seen here\n```\n      gelu, loglog, swish\n1/225 0.31  0.69    0.37\n2/225 0.    0.38    0.02\n4/225 0.    0.08    0.\n8/225 0.    0.03    0.\n```\nSupposedly, the drop is so steep due to the models themselves not being very robust to begin with.\nThe running times also increase with increased perturbation scale. The change is more pronounced for the small eps = 1e-7\n```\n      gelu, loglog, swish\n1/225 400s  160s    340s\n8/225 700s  312s    630s\n```\nthen for the larger eps=1e-3\n```\n      gelu, loglog, swish\n1/225 150s  87s     150s\n8/225 179s  90s     165s\n```\n\nWe hope to get the LinSyn results for the increased radii soon, so that we can make the comparison complete.\nHowever, so far nothing indicates that the comparison might be qualitatively different from the 1/255 comparison presented in the paper."}, {"review_id": "U6ebWEXVfX", "reviewer": "Reviewer_8jE6", "summary": "This paper proposed a method to upper and lower bound a scalar function within a convex set using linear functions. The proposed method works by solving a sequence of discrete linear bounding problems with an increasing number of sample points.", "questions": "1. What kind of training method is used to train the models in Table 1? Are those models robustly trained to against $\\ell_\\infty$ attacks?\n2. Table 1 suggests that lowering the $\\epsilon$ down from $10^{-5}$ would only increase the runtime of SOL without gaining any improvement to the \"fraction of properties certified\". Is this empirical finding consistent among larger $\\ell_\\infty$ bound perturbations? Because the $\\ell_\\infty$ bound perturbation used in Table 1 is quite small, 8/255 for MNIST and 1/255 for CIFAR-10. Would SOL require smaller $\\epsilon$ in order to achieve a similar \"fraction of properties certified\" as LinSyn?\n3. It would be helpful to also include the upper bound on the \"fraction of properties certified\" in Table 1 in order to gauge how robust the models are. Such upper bound can be easily computed using projected gradient descent or algorithm for finding the attacks for the neural network.\n", "limitations": "The authors have adequately addressed the limitations.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The proposed method works for neural networks with general activation functions. Most robust verifier in the existing liturature can only handle popular activation functions such as ReLU, Sigmoid and Softmax.\n2. This method achieve similar performance as the current state-of-the-art, LinSyn, while requires significantly less computational time. ", "weaknesses": "1. The experimental section is a little bit weak. The proposed method does not seem to improve the \"fraction of properties certified\" compared to LinSyn in Table 1. From Table 1, it seems like the proposed method is simply a more efficient implementation of LinSyn.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review and the feedback.\n\n[proposed method does not seem to improve the \"fraction of properties certified\" ...] We would like to point out that this can be easily explained by how LinSyn was tuned.\nWhile their method does not provide any tuning parameters for the user, it has a number of hardcoded parameters determining the trade-off between the running time and\nthe accuracy. These parameters might have been chosen in a such a way as to make the bound tight enough for the certification rates to be close to saturation. This\nhypothesis seems quite realistic considering that the experimental setup we use is exactly the same.\n\n[training method] We train the models by optimizing a simple cross-entropy without implementing any robust training techniques. Comparison on the models trained\nspecifically to be efficient at robustness certification would definitely be a useful extension.\n\n[scale of perturbation] We conducted additional experiments running SOL on CIFAR with perturbation radii 2/255, 4/255 and 8/255.\nSimilarly to the results for 1/255 presented in the paper the certification rates for any other radius do not depend on whether\nthe optimality target is set to 1e-7, 1e-5 or 1e-3. This means that for bigger radii the saturation still happens somewhere above 1e-3.\nThe certification rates drop significantly with increased radius as can be seen here\n```\n      gelu, loglog, swish\n1/225 0.31  0.69    0.37\n2/225 0.    0.38    0.02\n4/225 0.    0.08    0.\n8/225 0.    0.03    0.\n```\n\nSupposedly, the drop is so steep due to the models themselves not being very robust to begin with.\nThe running times also increase with increased perturbation scale. The change is more pronounced for the small eps = 1e-7\n```\n      gelu, loglog, swish\n1/225 400s  160s    340s\n8/225 700s  312s    630s\n```\nthen for the larger eps=1e-3\n```\n      gelu, loglog, swish\n1/225 150s  87s     150s\n8/225 179s  90s     165s\n```\n\nWe hope to get the LinSyn results for the increased radii soon, so that we can make the comparison complete.\nHowever, so far nothing indicates that the comparison might be qualitatively different from the 1/255 comparison presented in the paper.\n\n[optimization-based bound on the certification rates] Thank you for the suggestion! It might, indeed, give us new insights. We'll try to implement such bounds in a couple\nof days."}]

Author rebuttal:

Meta-review: This paper introduces a technique for estimating linear bounds on neural networks for use in certification settings. The approach is able to deal with arbitrary lipschitz activation functions without requiring manual user construction of bounds, thus suggesting a path towards automation of the procedure. The reviewers recognize the utility of this result and recommend acceptance.
Decision: Accept (poster)