Shifting-based Optimizable Linear Relaxations for General Activation Functions

Unknown

Conceptio › Archive › arXiv CS

arXiv CSopen access

Shifting-based Optimizable Linear Relaxations for General Activation Functions

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

machine learning, deep learning, neural networks

arXiv:2606.20292v1 [cs.LG] 18 Jun 2026

Shifting-based Optimizable Linear Relaxations for General Activation Functions Philipp Kern

László Antal and Erika Ábrahám

Carsten Sinz

Karlsruhe Institute of Technology [email protected]

RWTH Aachen University {antal,abraham}@cs.rwth-aachen.de

Karlsruhe University of Applied Sciences [email protected]

Abstract—The use of neural networks (NNs) is rapidly increasing, including in safety- and security-critical domains. To provide formal guarantees about NN behavior, many verification methods rely on optimizable linear relaxations of activation functions. However, existing techniques depend on hand-crafted relaxations for each activation function. Extension to state-of-the-art activation functions therefore requires substantial manual effort. In contrast, our approach SLiR (Shifting-based Linear Relaxations) is broadly applicable, requiring only a Lipschitz constant or a set of critical points. SLiR parameterizes relaxations by their slope and computes the corresponding offset via a shifting procedure that ensures sound upper and lower bounds over the input domain, enabling efficient optimization while maintaining correctness. Our experiments show that SLiR produces tight relaxations across a wide range of practical activation functions and enables verification of up to 7.8× more properties compared to state-of-the-art methods.

I. I NTRODUCTION Bound-propagation-based techniques are among the most successful approaches for scalable neural network verification. These methods overapproximate the nonlinear behavior of a network by propagating symbolic linear bounds through its layers, enabling efficient reasoning about safety or robustness. A key factor behind the effectiveness of recent approaches is the use of optimizable linear relaxations of activation functions, whose parameters (e.g., slope or a tangent point) are tuned to minimize concrete upper bounds on the output. Stateof-the-art verifiers such as α-CROWN [1] and its successors rely heavily on gradient-based optimization. Although this paradigm is effective for networks with ReLUs (and some other well-studied activation functions), extending optimizable linear relaxations to general activation functions remains a challenge. Modern NNs increasingly rely on complex activation functions (see Fig. 2) such as GELU (used in GPT [2], BERT [3], and Vision Transformers [4]), S WISH and its variants (appearing in models such as LLaMA [5], YOLOv4 [6], and non-edge versions of YOLOv7 [7]), M ISH, or adaptively learned functions (used in physics-informed neural networks (PINNs) [8]), for which no simple closed-form relaxations are available. For such functions, existing approaches typically construct static linear relaxations using sampling, SMT solving, or Lipschitz-based global optimization. Although these techniques can yield tight relaxations, they are fundamentally ill-suited for extending to slope optimization: soundness must be re-established from

scratch for each candidate slope, rendering gradient-based optimization prohibitively expensive. Our key observation is that, for any activation function f (x), a sound upper linear relaxation can be obtained by fixing a slope m and computing an offset b(m) = max f (x) − mx x∈[l,u]

(1)

that shifts the line mx to upper-bound f for its domain [l, u]. This can be computed efficiently, either given a closed form for the critical points1 of Eq. (1), or by optimizing over a sound piecewise linear (PWL) surrogate that we need to construct only once before gradient-based optimization. We show that closed-form solutions are available for a broad class of activation functions, and for the remaining cases PWL surrogates can be computed using Lipschitz-based optimization. Since any choice of the slope m yields a valid upper relaxation (or a lower relaxation when taking the minimum), we can optimize over the slope parameter m to obtain relaxations that lead to tighter bounds of the network output. Our Shifting-based Linear Relaxations (SLiR) approach, illustrated in Fig. 1, integrates directly with modern bound-propagation frameworks and supports slope-optimized verification for activation functions unsupported by existing methods. Our contributions in this paper are as follows: (C1) We introduce a novel framework for optimizable linear relaxations, based on parameterization and a shifting formalism, which is applicable to general activation functions (Section III). (C2) If no closed form for the critical points is available, we replace the original activation function by a piecewise linear envelope once using Lipschitz optimization and then optimize over that (Section III-B) for efficiency. (C3) Improvements for the Lipschitz optimization method tailored to our problem setting (Section III-C), resulting in tighter overapproximation and faster convergence. (C4) We derive an admissible range of slope values from the convex hull of the graph of the activation function and apply parameter transformation to significantly improve (or enable in the first place) convergence of the output bound optimization presented in this work (Section IV). 1 A critical point of Eq. (1) is a point x ∈ [l, u] where the derivative of f (x)−mx is zero, i.e. where the slope of f equals the candidate slope m.

Inputs

Input bounds

Neural network

Output constraints

Layer-wise initialization: Layer i

Exact optimization possible (Sec. III-A) Piecewise linear approx. of activation function (Sec. III-B)

Network-level optimization Neuron-level optimization: Neuron j

Sound lower and upper piecewise linear envelope (Sec. III-B)

Linear output relaxation

Minimize volume (local) (Sec. IV-A)

Compute volume (local) (Sec. IV-A)

Compute shift

Choose slopes m1 , ..., mk

Choose slope mj

Compute shift

Compute output bounds (Sec. II)

Optimize output bounds (global) (Sec. IV-B)

Linear relaxation

Verification result

Fig. 1: Overview of our Shifting-based Linear Relaxation (SLiR) approach for neural network verification. Name

Definition

S WISH GELU

x · σ(x), where σ(x) = (1 + e−x )−1 x·Φ(x), where Φ(x) is the CDF of the Gaussian distribution x · tanh ln(1 + ex ) x · tanh(x) arctan(x)2 − x x 1 − e−e

M ISH L I SHT ATAN S Q L OG L OG

S WISH GELU

L I SHT ATAN S Q

M ISH L OG L OG

Fig. 2: Activation function definitions (left) and their visual illustration (right).

Proofs for lemmas and propositions stated in our paper can be found in Appendix D. An artifact with our code and experimental evaluation is available online at: https://zenodo.org/ records/20133926. Our approach significantly reduces the effort required for implementing new activation functions — requiring only a formula for computing the critical points (in the exact case) or an overapproximation of the Lipschitz constant: while relaxations for the set of shared univariate activation functions require 1,340 lines of code in α-CROWN, only 154 lines of code are sufficient using our approach. Moreover, handcrafting custom relaxations can be complex and difficult to handle, even for domain experts. For the GELU and tanh functions, our approach yields better initializations than the hand-crafted implementation in α-CROWN (see Appendix F3).

Paulsen et al. [12] also employ a sampling-based approach combined with SMT-based verification to ensure soundness. Although applicable to a broad class of activation functions, this procedure is computationally expensive, as soundness must be re-established for each relaxation. In [13], the authors propose a template-based approach combined with machine learning. However, this method requires predefined ranges for the bounds and may fall back to less precise techniques outside these ranges. Biktairov et al. [14] synthesize linear overapproximations that are close to optimal, but compute only a single relaxation rather than a family of parameterized bounds. Similarly, Ma et al. [15] construct convex hulls for multi-neuron constraints, but do not support slope optimization. Related but more distant are synthesis- and SMT-based approaches such as [16], which construct linear or mixed-integer overapproximations but are not designed for efficient optimization. Overall, existing methods for constructing optimizable linear relaxations remain difficult to apply in practice and often require substantial expert knowledge (see, e.g., [17]). Decomposition-based approaches, such as [1], enable optimization by breaking functions into simpler components, but can introduce loose bounds and increased computational cost. Methods based on sampling, SMT solving, or repeated (expensive) soundness checks are likewise poorly suited for gradientbased optimization, since soundness must be re-established at

a) Related Work: Constructing tight linear relaxations has received significant attention in recent years, as the effectiveness of many NN verification algorithms depends on the quality of these bounds. Various works focus on synthesizing non-optimizable relaxations for specific activation functions. In particular, [9], [10], [11] follow a two-stage approach: they first approximate the function over the input domain using sampling-based techniques, and then restore soundness using dedicated correction procedures. While these methods can yield tight bounds, they are typically tailored to individual activation functions and do not directly support optimization.

every iteration. In contrast, our approach decouples soundness from parameter optimization by construction. By parameterizing relaxations via their slope and computing sound offsets through a shifting procedure, we obtain a family of valid relaxations that can be optimized efficiently for general activation functions. In the next section we provide background on linear relaxations for NN verification. SLiR extends this method by changing the construction of linear upper and lower relaxations for activation functions in Eq. (3).

0.0

0.1

0.3 0.4

II. BACKGROUND a) Neural Networks: A feed-forward NN is a function N : Rd0 → RdL with an input-, an output- and L − 1 hidden layers. The k-th layer consists of dk neurons and its activation values nk and preactivation values n̂k can be computed as n̂k = Wk nk−1 + bk

x∈X

Fig. 3: Illustrating our approach for M ISH: (i) Fit a PWL approximation of the function, (ii) extend it to lower/upper envelopes, and (iii) use them to construct linear relaxations.

nk = f (n̂k ) ,

where Wk and bk are the weight matrix and bias of the kth layer, and f : R → R is a non-linear activation function applied element-wise. To compute y = N (x), we identify the input x ∈ Rd0 with n0 and the output y ∈ RdL with n̂L . This representation also covers convolutional networks, since convolutions can be represented as matrix multiplications [18]. b) NN Verification: The NN verification problem that we consider consists of proving properties of input-output relations of the network. These properties are expressed as linear inequalities aT y ≤ d over the output y = N (x), which have to hold for all inputs x within a hyper-rectangle X ⊆ Rd0 . Instead of directly proving these universally quantified properties, we consider the equivalent optimization problem v ∗ = max aT y − d

mish pwl approximation lower envelope upper envelope lower relaxation upper relaxation

0.2

such that

y = N (x) ,

same procedure recursively to bound preactivation values of earlier layers — once from above and once from below. We denote this procedure to compute an upper bound v ∗ of the maximum violation v ∗ as v ∗ = BACKSUBSTITUTE(N , X ) .

no conclusions can be drawn. Due to the recursive bounds computation, the runtime of this approach is in O(L2 ). Concerning the linear upper and lower relaxations in Eq. (3), computed for each neuron, most methods ([12], [13], [14], [19]) choose α, β, α and β such that the area between αx + β and αx + β over a given bounded range of x ∈ [l, u] is minimized. However, while these relaxations are locally optimal, they may not lead to the tightest possible concrete upper bound on v ∗ achievable by linear relaxations [1]. Some approaches therefore parameterize the relaxations as

(2)

that finds the maximum violation v ∗ of the property aT y ≤ d (properties can be checked independently). If v ∗ is nonpositive, the property is proven safe; otherwise, it is violated. c) Linear Relaxations for NN Verification: This optimization problem, due to the activation functions, is highly non-convex in general and therefore difficult to solve exactly. Thus, many approaches focus on solving a relaxation of Eq. (2), where the activation functions are overapproximated using linear constraints, yielding a sound but incomplete linear programming (LP) formulation. However, modern NN verifiers ([19], [20], [21]) do not rely on standard LP solvers: Restricting the overapproximation of each activation function to exactly two constraints f (x) := αx + β ≤ f (x) ≤ αx + β =: f (x)

(4)

If v ∗ ≤ 0, then the property is proven successfully; otherwise,

f (x | θ) = α(θ) x + β(θ)

f (x | θ) = α(θ) x + β(θ) ,

where α(θ), β(θ), α(θ) and β(θ) yield slope and bias coefficients for valid lower and upper relaxations for each parameter θ. The search for the smallest possible upper bound can then be expressed as another optimization problem Pback over the stacked parameters θ for all neurons in N : Pback :

min BACKSUBSTITUTE(N , X | θ) . θ

(5)

This optimization problem can then be solved using gradientbased methods. If multiple linear properties need to be verified, Pback minimizes the sum of their violations.

(3)

III. G ENERAL A PPROACH Solving problem Pback requires a parameterization of relaxations f (x | θ) = α(θ)x + β(θ) and f (x | θ) = α(θ)x + β(θ) in the first place. In our approach, we only optimize over the slope values α(θ) and α(θ) ensuring soundness of the relaxations via shifting up or down by a sufficient amount

valid for all preactivation values x ∈ [l, u], enables them to instead use a specialized bound propagation procedure. This technique, referred to as backsubstitution ([19], [20], [22]), computes the same bounds on v ∗ as a standard LP solver on the two-constraint LP ([23], [24]) under mild conditions, but is significantly more efficient (see Appendix A for more details). The preactivation bounds necessary to construct the linear relaxations given in Eq. (3) can be computed by applying the

β(θ) = max f (x) − α(θ) · x , x∈[l,u]

(6)

where [l, u] are the preactivation bounds for the neuron that needs to be relaxed. Similarly, we set β(θ) to the corresponding minimum value for the lower relaxation. In cases — as discussed in Sec. III-A — where the critical points of f (x) − α(θ) · x for x ∈ [l, u] can be computed exactly, we can directly solve Eq. (6) to obtain sound relaxations. Otherwise, the approximate solution to Eq. (6) must guarantee soundness. Prior work has already followed a similar approach using interval ([12], [13]) or Lipschitz optimization [14] as overapproximating methods. However, calling these procedures for every neuron during every iteration of the gradient-descent process for Pback is prohibitively expensive. Our solution to this problem is illustrated in Fig. 3: Instead of calling these solvers in every iteration, we compute a PWL approximation and extend it to sound PWL lower and upper envelopes before the first pass of the gradient-descent iteration (Section III-B). Then, we use the lower or upper envelopes as surrogates for the original activation function in Eq. (6). This approach is efficient, since optimization over PWL functions is cheap and guarantees soundness, as the lower and upper envelopes enclose the original activation function. In both cases, (1) exact and (2) overapproximate optimization, our approach is easy to implement and to extend to new functions. In the latter case, we only require a Lipschitz constant. In the first case, we only need a formula (see Section III-A for further details) for computing the critical points. In the remainder of this section, we use a simplified form of Eq. (6), disregarding the parametrization: b = maxx∈[l,u] f (x) − mx, where f is the activation function, m a fixed candidate slope, and [l, u] the input domain.

B. Verified Piecewise Linear Envelopes If there is no closed form for f ′−1 (m), then we replace the function f in Eq. (6) with a sound upper PWL envelope fpwl , and maximize fpwl (x) − m · x. As a first step, we construct a PWL approximation f˜pwl (that is not necessarily sound for verification) and then apply a sufficiently large shift to each linear segment to achieve a valid upper envelope fpwl . Our approach works for any f˜pwl , however, its tightness will affect the tightness of the upper envelope and thus the linear relaxations. Additionally, Eq. (6) needs to be solved many times during gradient-based search for the best slopes. Here, fewer linear segments lead to higher efficiency. Therefore, we use an approach that tries to construct f˜pwl with linear segments that are maximally long while approximately maintaining a target approximation error (see Appendix C). Construction of Sound Envelopes: Given any PWL approximation f˜pwl (x) of an activation function f : R → R, specified by the set of segment boundaries Ξ = {ξi | i = 1, . . . , n + 1}, we compute an overapproximation ŷi∗ ≥

max x∈[ξi ,ξi+1 ]

f (x) − γi · x,

∀i = 1, . . . , n

(7)

of the shifting distance for γi · x to be a valid overapproximation of f (x) over the interval [ξi , ξi+1 ], where γi · x + δi is the i-th linear segment of f˜pwl over the interval [ξi , ξi+1 ]. The extension of the approximation f˜pwl to a verified upper envelope is then  ∗  γi · x + ŷi , x ∈ [ξi , ξi+1 ), i ≤ n (8) fpwl (x) = γn · x + ŷn∗ , x = ξn+1   ⊥, otherwise .

A. Exact Shifting

Similarly, we use an underapproximation of the minimum to construct a valid lower envelope. Note that, due to different ŷi∗ in neighboring segments, fpwl is not necessarily continuous (see Fig. 3). This does not affect correctness, since verification only requires pointwise overapproximation of the activation function. To compute the overapproximate shift ŷi∗ , we use an extension to the Lipschitz-based Piyavskii algorithm [26], which we introduce in Section III-C. This kind of overapproximate optimization methods can often compute tighter bounds on smaller intervals. When the segment bounds [ξi , ξi+1 ] of f˜pwl are smaller than the domains defined by the pre-activation bounds [lj , uj ] of each neuron j, our approach can benefit from tighter bounds compared to approaches that require verified optimization of f (x) − m · x for each neuron in a NN ([12], [14]). Additionally, the number of optimization problems we have to solve only depends on the width [l, u] = [minj lj , maxj uj ] of the pre-activation bounds of all neurons j in a layer of the NN and the corresponding number of linear segments of f˜pwl . This can be advantageous when the number of neurons is large compared to the number of linear segments of f˜pwl . Optimization over PWL Envelopes: Because fpwl is not necessarily continuous, it is not sufficient to consider just {fpwl (x) − m · x | x ∈ {l, u, ξ1 , ..., ξn+1 } ∩ [l, u]}, which only evaluates each linear segment for its left boundary point, as

For some functions, we can exactly solve Problem (6). This is the case for functions f where the pre-image f ′−1 (y) = {x | f ′ (x) = y} of the derivative of f is computable: Proposition 1. Let f : R → R be differentiable over the interval [l, u] and let f ′−1 (y) = {x | f ′ (x) = y} be the preimage of the derivative of f . Then b = maxx∈[l,u] f (x) − m · x = max f (x∗)−m·x∗ | x∗ ∈ {l, u}∪ f ′−1 (m)∩[l, u] . If there is a closed form expression for f ′−1 (m) and the set f (m) ∩ [l, u] is finite, then we can compute the maximum by just evaluating a finite number of points. Many functions fit this category. Examples include sigmoid, tanh and ELU [25]: ( x, x>0 1 , tanh(x), ELU(x) = x σ(x) = −x 1+e e −1, otherwise ′−1

where the critical points can be computed using q σ ′−1 (m) = −log( s1i −1) s1,2 = 21 ± 14 − m, m ∈ (0, 41 ] √ tanh′−1 (m) = arctan(si ) s1,2 = ± 1 − m, m ∈ (0, 1] ELU′−1 (m) = {0} ∪ {log(m) | m > 0} . Further functions of this category are listed in Appendix B.

Algorithm 1 Pseudocode of the Piyavskii method. Require: Function g(x), domain [l, u], Lipschitz constant Lg Ensure: Sound underestimation ŷi∗ 1: xl , zl , xr , zr , xm , zm ← INITIALIZE(g, Lg , l, u, nini ) 2: for k ← 1, . . . , max_iter do i 3: i ← argmini zm ▷ index of minimal segment i i 4: if zm − g(xm ) ≤ ε then i 5: return min g(l), zm , g(u)

candidates for the optimal value of Eq. (6). Instead, the set containing evaluations of both boundary points for each segment {γ1 · l + ŷ1∗ − m · l, γn · u + ŷn∗ − m · u}

(9)

∪ {γi · ξi + ŷi∗ − m · ξi | i ∈ 1, ..., n}

(10)

∪ {γi · ξi+1 + ŷi∗ − m · ξi+1 | i ∈ 1, ..., n}

(11)

is guaranteed to contain the optimal value. Here, the indices 1, ..., n correspond to the segments intersecting with [l, u]. Using that the linear terms (γi −m)·x+ ŷi∗ are monotonic, we can halve the number of boundary points that need to be evaluated. The maximum value over each linear segment is then ( (γi −m)· ˆl+ ŷi∗ , γi −m ≤ 0 ∗ (12) max (γi −m)·x + ŷi = (γi −m)· û+ ŷi∗ , otherwise x∈[l̂,û]

6: IMPROVE(i, g, Lg , xl , zl , xr , zr , xm , zm ) 7: end for i 8: i ← argmin i zm i 9: return min g(l), zm , g(u)

where [ˆl, û] = [l, u]∩[ξi , ξi+1 ]. Minimization is analogous.

The IMPROVE step (Line 6) refines the current approximation by splitting segment i into two subsegments via inserting a new trial point at the minimizer of the current minorant, denoted by xim . The objective function is evaluated at this point, i = g(xim ), and the original segment [xil , xir ] is yielding zm replaced by the two segments [xil , xim ] and [xim , xir ]. For each of the newly created segments, the corresponding minorant k ) is computed according to Eq. (15). intersection point (xkm , zm The algorithm iteratively applies this refinement by selecting, at each step, the segment i whose current minorant minimum i is smallest. The algorithm terminates either when a predezm fined maximum number of iterations is reached or when the i and difference between the minimum value of the minorant zm the objective function g is at most ε. Upon termination, the algorithm returns a sound lower bound i obtained as the mini(Lines 5 and 9) for mum of the minorant, that is min zm i , g(u) if the normal Piyavskii optimization, but min g(l), zm we consider our improved envelope described below. a) Improved Piyavskii: We add two improvements to the Piyavskii method that increase its efficiency. First, we allow the practitioner to provide a piecewise constant overestimation of the Lipschitz constant of the activation function f

C. Piyavskii Optimization For computing ŷi∗ of Eq. (7), we resort to Lipschitz optimization methods. These methods assume that the objective function g is Lipschitz-continuous over the domain [l, u], i.e., |g(x) − g(y)| ≤ Lg |x − y|

∀x, y ∈ [l, u] ,

(13)

for some known constant Lg ∈ R≥0 . The optimization method implemented in our framework is a tensor-friendly and efficient implementation of the Piyavskii–Schubert method (sawtooth method) [26]. Using this method, by underapproximating the minimum (overapproximating the maximum), we compute a sound lower (or upper) bound ŷi∗ for g(x) = f (x) − γi x over [ξi , ξi+1 ]. In the remainder of this section we focus on the minimization case — maximization can be achieved by minimizing g(x) = γi x − f (x). The intuition behind the Piyavskii method is that, given the Lipschitz constant Lg (or a sound overapproximation thereof), one can construct a straightforward lower Lipschitz envelope between two sample points xil < xir via linear lower bounds g(x) ≥ −Lg (x−xil )+g(xil )

g(x) ≥ Lg (x−xir )+g(xir ) (14)

(and similarly for the upper envelope) that forms a sound underpproximation of the objective function (see left sub-figure of Fig. 4) without requiring an in-depth analytical treatment of g. Based on this envelope, a lower piecewise linear approximation, referred to as a minorant (see right sub-figure of Fig. 4), can be constructed over a finite decomposition of the domain [l, u]. The Piyavskii algorithm iteratively refines this minorant until its minimum is ε-close to the function g. We provide pseudocode in Alg. 1. The INITIALIZE method (Line 1) builds an initial minorant with nini − 1 segments by sampling nini uniformly spaced points in the interval [l, u]. Each minorant segment i is characterized by its left and right endpoints xil and xir , their corresponding trial values zli = g(xil ) and zri = g(xir ), and the i minorant intersection point (xim , zm ), computed as the intersection of the linear functions in Eq. (14) xim =

z i − zli xil + xir − r , 2 2Lg

( Lif Lf (x) = ⊥

, x ∈ [ai , ai+1 ) , otherwise ,

(16)

for some interval bounds ai , rather than a single global constant over the entire domain. That is, for each segment [xil , xir ], a tighter local Lipschitz constant maxx∈[xil ,xir ] Lf (x) can be used, yielding a tighter minorant (right sub-figure of Fig. 4). Secondly, we exploit that we apply Lipschitz optimization to the specific case of optimizing the difference of a general function and a linear function: g(x) = f (x) − γi x. Traditional Lipschitz optimization would infer the global Lipschitz constant for g as Lg = Lf +|γi |, where Lf is a Lipschitz constant of f (see Proposition 2.3.3 of [27]). However, in our case, we can derive tighter lower bounds for x ∈ [xl , xr ] directly from the Lipschitz envelope of f :

i zm = Lg (xim − xir ) + zri . (15)

f (x) ≥ −Lf (x − xl ) + f (xl )

f (x) ≥ Lf (x − xr ) + f (xr ) ,

(ℓ, g(ℓ))f (ℓ) f(

ℓ)

−

ℓ+

(

+ f (ℓ) −

)( γ|

− L f (x

x−

ℓ)

g(x) f( u)

ℓ) −

−

x2 L

γx

−L f (u) ℓ)

g(x) +

γx )−

−

−γ − x)

γu

x) −

γ x(u, g(u)) u) − x )(x

γ|

Fig. 4: Lipschitz envelope and minorant of an objective function g(x) = f (x) − γx. Left: Construction of a straightforward Lipschitz envelope (purple/solid) and of the improved envelope (orange/dashed). Right: Lipschitz minorant constructed using a global Lipschitz constant (purple/dashed) and local Lipschitz constants (orange/dashed). The next trial points x4 obtained from both minorants are shown. Dotted lines (purple and orange) show the minorant segments of the next iteration.

and subtracting γi x from both sides we obtain the following valid lower bounds for g(x): g(x) = f (x) − γi x ≥ −Lf (x − xl ) + f (xl ) − γi x

(17)

g(x) = f (x) − γi x ≥ Lf (x − xr ) + f (xr ) − γi x .

(18)

IV. PARAMETER O PTIMIZATION Following existing work [1], we divide the process of finding parameters that lead to a good solution of the optimization problem Pback into two phases: First, we heuristically set initial slopes for the lower and upper linear relaxations. Then we use these slopes as starting point for gradient-based optimization of Pback .

These linear bounds define an envelope (see left sub-figure of Fig. 4) that provably improves tightness:

A. Parameter Initialization

Proposition 2. Let f : [xl , xr ] → R be Lipschitz continuous with constant Lf and g(x) = f (x) − γx. Let the straightforward Lipschitz envelope use Lg = Lf + |γ|, and let the improved envelope be obtained by applying the Lipschitz bounds to f and subtracting γx. Then, for all x ∈ [xl , xr ], the improved lower envelope g(x) ≥ f (xl ) − Lf (x − xl ) − γx =: ℓimp (x) l

(19)

g(x) ≥ f (xr ) + Lf (x − xr ) − γx =: ℓimp r (x)

(20)

When optimizing Pback , it is beneficial to start with initial parameters that already yield close bounds. Many NN verifiers ([12], [13], [14], [19], [28]) choose parameters that minimize the area between the lower and upper linear relaxation over the pre-activation bounds [l, u] for each neuron, by solving Z u αx + β − (αx + β) dx . (23) argmin α,α,β,β

While other verifiers rely on manually crafted closed-form solutions for specific activation functions ([19], [28]), sampling and linear programming ([12], [14]), or precomputation [13], we use gradient descent to find slopes α, α that result in small areas. Note that via Eq. (6), β and β are functions of the slopes in our case, reducing Eq. (23) to a bivariate optimization problem. Examples for initial relaxations found by SLiR are shown in Fig. 5a.

is pointwise at least as tight as the straightforward one: g(x) ≥ g(xl ) − Lg (x − xl ) =: ℓstd l (x) std g(x) ≥ g(xr ) + Lg (x − xr ) =: ℓr (x) ,

(21) (22)

imp std i.e., ℓimp (x) ≥ ℓstd l (x) and ℓr (x) ≥ ℓr (x), ∀x ∈ [xl , xr ]. l Moreover, if γ > 0, the bottom-right segment ℓimp r (x) is strictly tighter while the bottom-left coincide; if γ < 0, the bottom-left segment ℓimp (x) is strictly tighter while the l bottom-right coincide.

B. Parameter Optimization After initial slopes have been found for all neurons, we use them as starting point for solving Pback . Thereby, the choice of the parametrization functions α(θ), α(θ), β(θ) and β(θ) is highly relevant. However, a naive choice of

A similar statement holds for the upper improved Lipschitz envelope. In this case, the top-left segment is strictly tighter for γ > 0, while for γ < 0, the top-right segment is strictly tighter than the standard Lipschitz envelope. Note that, in the standard Piyavskii method (Eq. (14)) the bounds have different slope signs, which guarantees that the minimum of the minorant is at the intersection point xm of the bounds. However, in Eq. (17) and (18) we might have matching slope signs, so the true minimum of the minorant might be also at the endpoints g(xl ) or g(xr ).

α(θ) = θ

β(θ) = min f (x) − α(θ) · x x∈[l,u]

α(θ) = θ

β(θ) = max f (x) − α(θ) · x x∈[l,u]

(24) (25)

for parameter θ ∈ R leads to unstable optimization performance. As shown in Fig. 5c (topmost two lines for learning

1.75

loglog convex hull SLiR lower relaxation SLiR upper relaxation

1.50 1.25

loglog convex hull loose relaxation tight relaxation

2.0 1.5

65 Loss

1.00

1.0

0.75

slope slope (lr = 0.05) slope_bounds slope_bounds (lr = 0.05) SLiR SLiR (lr = 0.05)

0.50

0.5

0.25 0.00 3

1 x

(a) SLiR relaxation after area minimization.

0.0

1 x

(b) Tight relaxation with α = amax and loose relaxation with α = amax + 0.05.

10 12 Iteration

Fig. 5: Effect of using the convex hull’s lower and upper envelopes to bound the admissible slopes of the linear relaxations.

rates 0.1 and 0.05), the loss curve spikes sharply after the first gradient descent iteration. We identify two reasons for this behaviour: High sensitivity to the learning rate and suboptimal relaxations achievable via these parameters. While SLiR generates valid linear relaxations for any slope α, α ∈ R, many of these relaxations are not tight: For some slopes a the corresponding upper relaxations U (x) = a · x + b are dominated by linear upper relaxations Û (x) that are strictly better over the interior of the whole domain x ∈ [l, u] (see Fig. 5b). To address this issue, we restrict the admissible slopes to the range defined by the upper (or lower) convex hull of f (x) over [l, u]. This approach is justified by the following result:

In contrast to parameterized lower relaxations for the ReLU function, where the interval [0, 1] always produces tight lower relaxations for any unstable ReLU, admissible ranges for other activation functions can be narrower (for the loglog activation function shown in Fig. 5b the admissible range is just [0.0571, 0.1248] for the depicted domain) and vary between each neuron (a different domain would result in different admissible ranges for the slopes). To ensure a smooth mapping from parameter θ ∈ R to slope values within [amin , amax ], we parameterize slopes via the sigmoid function σ(·): α(θ) = almin + (almax − almin ) · σ(θ)

(26)

Lemma 1. Let l < u, f : R → R, amax the maximum slope of its upper convex hull H(x) over x ∈ [l, u], ϵ > 0 and U (x) = (amax + ϵ)x + b ≥ f (x) be a linear upper relaxation over [l, u]. Then there is another linear overapproximation Û (x) with slope amax , s.t. f (x) ≤ Û (x) < U (x) over x ∈ (l, u].

α(θ) = aumin + (aumax − aumin ) · σ(θ) ,

(27)

almin , almax , aumin and aumax are the minimum and maximum slope of the lower and upper convex envelopes. Initial slopes α are mapped to θ by application of the inverse sigmoid function α − amin θ = log , (28) amax − amin

A similar statement holds for the minimum slope amin of the upper convex hull (see Lemma 2 in Appendix D). Hence, we do not compromise quality of the achievable linear relaxations by clamping the admissible slope values for the upper relaxations to α ∈ [amin , amax ] (and similarly for lower relaxations and the lower convex hull). We obtain amin and amax , by computing the convex hull of sampled points {(xi , f (xi )) | xi ∈ [l, u]} (if we computed a PWL approximation fˆpwl for f , we use the function values at the boundary points {(l, f (l)), (u, f (u))} ∪ {(ξi , f (ξi )) | ξi ∈ [l, u]}) using the monotone chain algorithm [29]. This algorithm directly returns the lower and upper convex envelope and runs in O(n) time when the input points are sorted first by their x-values and then by their y-values. However, our vectorized implementation may not take full advantage of the linear runtime. Since the local input bounds [l, u] can get tighter with each gradient-descent step for Pback , we have to recompute the convex hull for each neuron in each iteration. As illustrated in Fig. 5c (center two lines), restricting the range of the slope values drastically reduces the deterioration of the loss. However, the optimization process is still highly sensitive to the learning rate.

where we clamp the argument of the log to [ϵσ , 1−ϵσ ] to avoid numerical issues and amin and amax correspond to almin and almax or aumin and aumax respectively. Illustrated by the lines labeled SLiR in Fig. 5c, this transformation allows consistent improvement in the loss curves, even for higher learning rate. V. E VALUATION We evaluate (i) how our approach compares to prior work and (ii) the impact of our optimizations on performance. We report the number of verified properties (# certified) and total runtime (in seconds) at initialization and after 20 gradientbased optimization steps for each lower and upper output bounds. All experiments were run on a cluster with multiple nodes of Intel Xeon Platinum 8358 processors (64 cores, 256 GiB RAM, 2.6 GHz), with each run limited to 4 CPU cores. Further evaluation results can be found in Appendix F. Tools. We implemented our approach in PyTorch [30] and integrated it into the backsubstitution-based framework AU TO L I RPA [20]. We compare against α-CROWN [1], the leading bound propagation method with optimizable relaxations

and the basis of α-β-CROWN, winner of VNN-COMP 20212025 ([31], [32], [33], [34], [35]). While α-CROWN relies on custom relaxations or decomposition into elementary operations, we also compare to methods for general activation functions: L IN S YN [12] and SOL [14], which synthesize relaxations at runtime, and the example-guided approach [13], which precomputes relaxations and falls back to interval overapproximation outside known bounds.

depends strongly on pre-activation bounds and approximation quality: more segments lead to more critical points to evaluate. Comparison to Prior Work. Tab. II compares SLiR to SOL [14], L IN S YN [12], the example guided (E. Guided) approach [13] and α-CROWN [1]. We used adversarial attacks (the PGD method [45]) to generate counterexamples for the verification problems. The top row indicates instances with no found counterexample, serving as an upper bound on the number of provable instances. For all approaches, we report the number of verified instances and runtime (if source code was available). If parametric relaxations are used, we show both metrics for the initial relaxations and after optimization. Our experiments show that optimizing the relaxations significantly boosts verification rate — with the effects being more pronounced for the MNIST than for the CIFAR10 benchmark. A possible reason could be that the pre-activation ranges, and thus the range of admissible slope values, of the neurons were larger in the MNIST than in the CIFAR benchmark. Notably, SLiR verifies more instances — both after initialization and after optimization — than α-CROWN even for GELU, where α-CROWN uses a custom relaxation. Similarly, we can verify significantly more properties than E. Guided for the L I SHT and M ISH MNIST NNs. While times for L IN S YN have to be taken with a grain of salt2 , repeating their approach for all 40 gradient descent iterations will clearly exceed SLiR’s total time after optimization. As expected, SLiR’s initialization times are slower than for α-CROWN and E. Guided, as these approaches use composition of custom relaxations or precomputation, whereas SLiR has to generate a PWL approximation and extend it to an overapproximation for each instance. However, SLiR is only 1.4 − 1.9× slower than α-CROWN on MNIST (except for M ISH, where it is significantly faster). On CIFAR10, performance varies: Notably, α-CROWN struggles with M ISH because decomposing it into a large number of elementary functions increases the NN’s effective depth (see complexity of Eq. (4)).

Benchmarks. We use benchmarks from the literature with general activation functions ([12], [13], [14]), which are convolutional NNs trained on MNIST [36] and CIFAR10 [37] by Paulsen et al. [13]. The MNIST NNs have input dimension 28 × 28, hidden layers of size 1568, 784, 256, and output size 10, while the CIFAR10 NNs have input dimension 3×32×32, hidden layers of size 2048, 2048, 1024, 256, and output size 10. For both datasets, we consider ATAN S Q [38], GELU [39], L I SHT [40], L OG L OG [41], M ISH [42], and S WISH [39] activations. For cases where the critical points of Eq. (6) can be computed exactly, we also evaluate two NNs with 6 layers of 500 neurons using sigmoid and tanh activations from the ERAN benchmark suite ([43], [19]). We additionally train a NN of the same architecture with ELU activations [25]. We verify adversarial robustness [44] in the L∞ -sense with ϵ = 8/255 for MNIST and ϵ = 1/255 for CIFAR10 on the 100 images used in prior work ([12], [13], [14]), skipping inputs that are already misclassified. For the ERAN MNIST NNs, we follow the same procedure as [43], using the same 100 images and ϵ = 8/255. Tab. I: Results for functions with computable critical points. Total time shows runtime; #certified counts verified instances. ERAN MNIST 6 × 500 †

α-CROWN

init opt

SLiR

init opt

SLiR (exact)

init opt

Metric

sigmoid

tanh

ELU

#certified total time #certified total time

27 16.3s 33 3 185s

10 18.5s 20 3 579s

−† −† −† −†

#certified total time #certified total time

29 257s 31 5 942s

16 541s 21 11 364s

62 352s 64 7 137s

#certified total time #certified total time

29 205s 33 4 966s

16 201s 21 4 953s

62 200s 64 4 730s

α-CROWN cannot decompose ELU.

Method

Ablations. To evaluate our optimizations from Sec. III-C, we conducted an ablation study on the MNIST benchmark. The top row in Tab. III shows the number of verified instances and runtimes with all optimizations enabled. Using only the traditional Piyavskii method for Lipschitz optimization (SLiR (base)), adding just the improved envelopes (SLiR (envelope)) or just the local Lipschitz constants (SLiR (local)) decreases the number of verified properties, showing the effectiveness of each optimization. However, SLiR (envelope) verifies just one instance less than our final approach, showing good performance even if only a global Lipschitz constant is available. Allowing wider bounds for the slope of the linear lower and upper relaxation within the range [minx∈[l,u] f ′ (x), maxx∈[l,u] f ′ (x)] (SLiR (loose-bounds)) of a neuron significantly decreases the number of verified instances after optimization. The difference is especially evident for the NNs with L I SHT and L OG L OG activation functions.

Functions with Computable Critical Points. Tab. I shows that both SLiR configurations verify more instances than αCROWN after initialization, indicating tighter initial relaxations. However, α-CROWN is significantly faster at initialization, as it avoids PWL approximation and gradient-based parameter search. Optimization times for SLiR (exact) are relatively consistent across networks, whereas the PWL variant

2 We had to replace the python module multiprocessing by pathos.multiprocessing to be able to run the code.

Tab. II: Verification results for adversarial robustness. Time (in seconds) shows the total time spent for analyzing all instances while #cert. shows the number of instances verified as safe. UB (upper bound) indicates the number of instances for which the PGD attack was unable to find a counterexample; it therefore is an upper bound on the number of certifiable instances.

Method

Metric

atansq

gelu

lisht

loglog

CIFAR10 CNN 5-Layer (100 instances each)

mish

swish

atansq

gelu

lisht

loglog

mish

swish

PGD

SOL§ L IN S YN¶

#cert. #cert. time #cert. time

−† 73 −† 72 −† 1 345 16 70 138 183

−† −† −† 11 114

24 23 1 393 −† −†

−† −† −† 28 140

76 76 1 338 74 115

−† −† −† 14 214

20 20 1 803 19 263

−† −† −† 0 185

27 27 1 017 −† −†

−† −† −† 19 222

24 23 1 747 23 182

#cert. time #cert. time

1 7.8 3 4 868

64 0 6.7 8.9 70 9 3 488 3 625

−‡ −‡ −‡ −‡

0 45.2 6 16 172

34 9.8 47 3 959

0 13.1 3 29 602

11 0 13.6 15.3 14 0 19 465 17 899

−‡ −‡ −‡ −‡

0 451 1 76 707

2 14.5 5 16 339

#cert. time #cert. time

18 958 25 6 808

72 81 647 950 75 86 5 794 6 940

23 308 32 4 844

63 857 66 6 213

76 567 76 5860

15 785 16 19 300

20 0 466 953 21 0 17 880 21 328

26 173 28 11 888

19 462 20 18 806

23 434 23 18 372

E. Guided α-CROWN

init opt

SLiR (ours)

init opt

Tab. III: Effects of different contributions (ablation study). Top lines show # cert. instances, bottom lines show times (s).

Replacing the sigmoid-based parametrization of Eq. (28) with direct optimization of the slope values (SLiR (directslope)) likewise reduces the number of verified instances, indicating that both the convex-hull restriction and the parameter transformation contribute substantially to overall performance.

Metric

atansq

gelu

lisht

loglog

mish

swish

#cert. time #cert. time

18 666 25 5,763

72 490 75 5 095

81 681 86 5 939

23 270 32 4 433

63 579 66 5 105

76 437 76 5 259

#cert. time #cert. time

18 72 80 829 641 777 25 74 86 6 299 5 439 6 634

23 62 615 676 31 64 4 880 5 650

76 718 76 6 023

SLiR SLiR (envelope) (local) Section III-C Section III-C

init

#cert. time

18 926

72 675

81 799

23 649

63 727

76 729

opt

#cert. time

25 7 064

74 5 822

86 6 759

32 5 286

65 6 084

77 6 166

init

#cert. time

18 1010

72 811

80 953

23 732

62 917

76 962

opt

#cert. time

25 6 556

74 5 617

86 6 593

32 5 333

64 6 230

76 6 481

SLiR (loose bounds) Section IV

init

#cert. time #cert. time

18 72 74 870 544 808 22 74 79 4 005 3 389 3 672

23 63 272 763 25 64 3 543 3 812

76 524 76 3 819

SLiR (direct slope) Section IV

MNIST CNN 4-Layer Method

init

#cert. time #cert. time

18 72 81 878 585 867 22 74 86 6 665 5 727 6 759

23 63 282 811 24 64 5 294 6 414

76 533 76 6 331

SLiR

init opt

SLiR (base)

init opt

opt

† Activation function not supported by the tool. ‡ α-CROWN encountered NaN in the computation. ¶ We had to replace the multiprocessing library, runtimes may be affected. § Results taken from [14].

MNIST CNN 4-Layer (100 instances each)

VI. C ONCLUSION We presented SLiR, a shifting-based framework for constructing parameterized, optimizable linear relaxations for general activation functions. By separating slope optimization from synthesizing valid, tight relaxations through minimal vertical shifting computation, SLiR enables automatic and principled construction of sound relaxations without requiring expert knowledge and manual activation-specific derivations. When the shifting distance can be computed in closed-form, the practitioner only needs to provide the set of critical points. If no closed-form solution exists, it suffices to provide one global or multiple local Lipschitz constants to enable sound overapproximation. While SLiR is not always as fast as some specialized relaxation techniques, our evaluation shows that it can verify a larger number of properties on the overwhelming majority of tested benchmarks. Acknowledgements. The authors acknowledge support by the state of Baden-Württemberg through bwHPC. R EFERENCES [1] K. Xu, H. Zhang, S. Wang, Y. Wang, S. Jana, X. Lin, and C. Hsieh, “Fast and complete: Enabling complete neural network verification with rapid and massively parallel incomplete verifiers,” in Proc. of the 9th International Conference on Learning Representations(ICLR’21), Virtual Event, Austria, May 3-7, 2021. OpenReview.net, 2021. [2] A. Radford, K. Narasimhan, T. Salimans, and I. Sutskever, “Improving language understanding by generative pre-training,” OpenAI, 2018.

Runtimes are faster for SLiR (loose-bounds), because no convex hulls have to be computed. The number of initially verified instances can also differ, as the slopes of the initial linear relaxations are not clipped to the tightest possible range.

[3] J. Devlin, M. Chang, K. Lee, and K. Toutanova, “BERT: pre-training of deep bidirectional transformers for language understanding,” in Proc. of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies (NAACLHLT’19), Minneapolis, MN, USA, June 2-7, 2019. Association for Computational Linguistics, 2019, pp. 4171–4186. [4] A. Dosovitskiy, L. Beyer, A. Kolesnikov, D. Weissenborn, X. Zhai, T. Unterthiner, M. Dehghani, M. Minderer, G. Heigold, S. Gelly, J. Uszkoreit, and N. Houlsby, “An image is worth 16x16 words: Transformers for image recognition at scale,” in Proc. of the 9th International Conference on Learning Representations (ICLR’21), Virtual Event, Austria, May 3-7, 2021. OpenReview.net, 2021. [5] H. Touvron, T. Lavril, G. Izacard, X. Martinet, M. Lachaux, T. Lacroix, B. Rozière, N. Goyal, E. Hambro, F. Azhar, A. Rodriguez, A. Joulin, E. Grave, and G. Lample, “Llama: Open and efficient foundation language models,” CoRR, vol. abs/2302.13971, 2023. [6] A. Bochkovskiy, C. Wang, and H. M. Liao, “YOLOv4: Optimal speed and accuracy of object detection,” CoRR, vol. abs/2004.10934, 2020. [7] C. Wang, A. Bochkovskiy, and H. M. Liao, “Yolov7: Trainable bag-offreebies sets new state-of-the-art for real-time object detectors,” in Proc. of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR’23), Vancouver, BC, Canada, June 17-24, 2023. IEEE, 2023, pp. 7464–7475. [8] H. Wang, L. Lu, S. Song, and G. Huang, “Learning specialized activation functions for physics-informed neural networks,” ArXiv, vol. abs/2308.04073, 2023. [Online]. Available: https://api.semanticscholar. org/CorpusID:260704245 [9] M. Balunovic, M. Baader, G. Singh, T. Gehr, and M. T. Vechev, “Certifying geometric robustness of neural networks,” in Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019 (NeurIPS’19), December 8-14, 2019, Vancouver, BC, Canada, 2019, pp. 15 287–15 297. [10] W. Ryou, J. Chen, M. Balunovic, G. Singh, A. M. Dan, and M. T. Vechev, “Scalable polyhedral verification of recurrent neural networks,” in Proc. of the 33rd International Conference on Computer Aided Verification (CAV’21), Virtual Event, July 20-23, 2021, Proceedings, Part I, ser. Lecture Notes in Computer Science, vol. 12759. Springer, 2021, pp. 225–248. [11] J. Laurel, S. B. Qian, G. Singh, and S. Misailovic, “Synthesizing precise static analyzers for automatic differentiation,” Proc. of the ACM on Programming Languages, vol. 7, no. OOPSLA2, pp. 1964–1992, 2023. [Online]. Available: https://doi.org/10.1145/3622867 [12] B. Paulsen and C. Wang, “LinSyn: Synthesizing tight linear bounds for arbitrary neural network activation functions,” in Proc. of the 28th International Conference on Tools and Algorithms for the Construction and Analysis of Systems (TACAS’22), 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, ser. Lecture Notes in Computer Science, vol. 13243. Springer, 2022, pp. 357–376. [13] ——, “Example guided synthesis of linear approximations for neural network verification,” in Proc. of the 34th International Conference on Computer Aided Verification (CAV’22), Haifa, Israel, August 710, 2022, Proceedings, Part I, ser. Lecture Notes in Computer Science, vol. 13371. Springer, 2022, pp. 149–170. [Online]. Available: https://doi.org/10.1007/978-3-031-13185-1_8 [14] Y. Biktairov and J. Deshmukh, “SOL: sampling-based optimal linear bounding of arbitrary scalar functions,” in Advances in Neural Information Processing Systems 36: Annual Conference on Neural Information Processing Systems 2023 (NeurIPS’23), New Orleans, LA, USA, December 10 - 16, 2023, 2023. [15] Z. Ma, Z. Wang, and G. Bai, “Convex Hull Approximation for Activation Functions,” Proc. of the ACM on Programming Languages, vol. 9, no. OOPSLA2, pp. 1007–1033, 2025. [Online]. Available: https://dl.acm.org/doi/10.1145/3763086 [16] J. Wang, A. Gupta, and C. Wang, “Synthesizing MILP constraints for efficient and robust optimization,” Proc. of the ACM Programming Languages, vol. 7, no. PLDI, pp. 1896–1919, 2023. [17] Z. Shi, Q. Jin, Z. Kolter, S. Jana, C.-J. Hsieh, and H. Zhang, “Neural network verification with branch-and-bound for general nonlinearities,” in Tools and Algorithms for the Construction and Analysis of Systems. Cham: Springer Nature Switzerland, 2025, pp. 315–335. [18] K. Xu, Z. Shi, H. Zhang, M. Huang, K. Chang, B. Kailkhura, X. Lin, and C. Hsieh, “Automatic perturbation analysis on general computational graphs,” CoRR, vol. abs/2002.12920, 2020.

[19] G. Singh, T. Gehr, M. Püschel, and M. T. Vechev, “An abstract domain for certifying neural networks,” Proc. of the ACM Programming Languages, vol. 3, no. POPL, pp. 41:1–41:30, 2019. [20] K. Xu, Z. Shi, H. Zhang, Y. Wang, K. Chang, M. Huang, B. Kailkhura, X. Lin, and C. Hsieh, “Automatic perturbation analysis for scalable certified robustness and beyond,” in Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems (NeurIPS’20), December 6-12, 2020, virtual, 2020. [21] T. Zelazny, H. Wu, C. W. Barrett, and G. Katz, “On optimizing back-substitution methods for neural network verification,” in 22nd Formal Methods in Computer-Aided Design, FMCAD 2022, Trento, Italy, October 17-21, 2022, A. Griggio and N. Rungta, Eds. IEEE, 2022, pp. 17–26. [Online]. Available: https://doi.org/10.34727/2022/ isbn.978-3-85448-053-2_7 [22] E. Wong and J. Z. Kolter, “Provable defenses against adversarial examples via the convex outer adversarial polytope,” in Proceedings of the 35th International Conference on Machine Learning, ICML 2018, Stockholmsmässan, Stockholm, Sweden, July 10-15, 2018, ser. Proceedings of Machine Learning Research, J. G. Dy and A. Krause, Eds. PMLR, 2018, pp. 5283–5292. [Online]. Available: http://proceedings.mlr.press/v80/wong18a.html [23] H. Salman, G. Yang, H. Zhang, C. Hsieh, and P. Zhang, “A convex relaxation barrier to tight robustness verification of neural networks,” in Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada, H. M. Wallach, H. Larochelle, A. Beygelzimer, F. d’Alché-Buc, E. B. Fox, and R. Garnett, Eds., 2019, pp. 9832– 9842. [Online]. Available: https://proceedings.neurips.cc/paper/2019/ hash/246a3c5544feb054f3ea718f61adfa16-Abstract.html [24] Z. Lyu, C. Ko, Z. Kong, N. Wong, D. Lin, and L. Daniel, “Fastened CROWN: tightened neural network robustness certificates,” in The Thirty-Fourth AAAI Conference on Artificial Intelligence, AAAI 2020, The Thirty-Second Innovative Applications of Artificial Intelligence Conference, IAAI 2020, The Tenth AAAI Symposium on Educational Advances in Artificial Intelligence, EAAI 2020, New York, NY, USA, February 7-12, 2020. AAAI Press, 2020, pp. 5037–5044. [Online]. Available: https://doi.org/10.1609/aaai.v34i04.5944 [25] D. Clevert, T. Unterthiner, and S. Hochreiter, “Fast and accurate deep network learning by exponential linear units (ELUs),” in Proc. of the 4th International Conference on Learning Representations (ICLR’16), San Juan, Puerto Rico, May 2-4, 2016, 2016. [26] S. A. Piyavskii, “An algorithm for finding the absolute extremum of a function,” Journal of USSR Computational Mathematics and Mathematical Physics, vol. 12, no. 4, pp. 57–67, 1972. [27] Ş. Cobzaş, R. Miculescu, A. Nicolae et al., Lipschitz functions. Springer, 2019. [28] H. Zhang, T. Weng, P. Chen, C. Hsieh, and L. Daniel, “Efficient neural network robustness certification with general activation functions,” in Advances in Neural Information Processing Systems 31: Annual Conference on Neural Information Processing Systems (NeurIPS’18), December 3-8, 2018, Montréal, Canada, 2018, pp. 4944–4953. [29] A. M. Andrew, “Another efficient algorithm for convex hulls in two dimensions,” Information Processing Letters, vol. 9, no. 5, pp. 216– 219, 1979. [Online]. Available: https://doi.org/10.1016/0020-0190(79) 90072-3 [30] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Köpf, E. Z. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, “PyTorch: An imperative style, highperformance deep learning library,” in Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada, 2019, pp. 8024–8035. [31] S. Bak, C. Liu, and T. T. Johnson, “The second international verification of neural networks competition (VNN-COMP 2021): Summary and results,” CoRR, vol. abs/2109.00498, 2021. [32] M. N. Müller, C. Brix, S. Bak, C. Liu, and T. T. Johnson, “The third international verification of neural networks competition (VNN-COMP 2022): Summary and results,” CoRR, vol. abs/2212.10376, 2022. [33] C. Brix, S. Bak, C. Liu, and T. T. Johnson, “The fourth international verification of neural networks competition (VNN-COMP 2023): Summary and results,” CoRR, vol. abs/2312.16760, 2023.

[34] C. Brix, S. Bak, T. T. Johnson, and H. Wu, “The fifth international verification of neural networks competition (VNN-COMP 2024): Summary and results,” CoRR, vol. abs/2412.19985, 2024. [35] K. Kaulen, T. Ladner, S. Bak, C. Brix, H. Duong, T. Flinkow, T. T. Johnson, L. Koller, E. Manino, T. H. Nguyen, and H. Wu, “The 6th international verification of neural networks competition (VNN-COMP 2025): Summary and results,” CoRR, vol. abs/2512.19007, 2025. [36] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner, “Gradient-based learning applied to document recognition,” Proc. of the IEEE, vol. 86, no. 11, pp. 2278–2324, 1998. [37] A. Krizhevsky and G. Hinton, “Learning multiple layers of features from tiny images,” University of Toronto, Toronto, Ontario, Tech. Rep. 0, 2009. [Online]. Available: https://www.cs.toronto.edu/~kriz/ learning-features-2009-TR.pdf [38] P. Ramachandran, B. Zoph, and Q. V. Le, “Searching for activation functions,” in Proc. of the 6th International Conference on Learning Representations (ICLR’18), Vancouver, BC, Canada, April 30 - May 3, 2018, Workshop Track Proceedings. OpenReview.net, 2018. [39] D. Hendrycks and K. Gimpel, “Bridging nonlinearities and stochastic regularizers with gaussian error linear units,” CoRR, vol. abs/1606.08415, 2016. [40] S. K. Roy, S. Manna, S. R. Dubey, and B. B. Chaudhuri, “LiSHT: Nonparametric linearly scaled hyperbolic tangent activation function for neural networks,” in Proc. of the 7th International Conference on Computer Vision and Image Processing (CVIP’22), Nagpur, India, November 4-6, 2022, Revised Selected Papers, Part I, ser. Communications in Computer and Information Science, vol. 1776. Springer, 2022, pp. 462–476. [41] G. S. d. S. Gomes and T. B. Ludermir, “Complementary Log-Log and Probit: Activation functions implemented in artificial neural networks,” in Proc. of 2008 Eighth International Conference on Hybrid Intelligent Systems, Barcelona, Spain, September 10-12, 2008. IEEE, 2008, pp. 939–942. [42] D. Misra, “Mish: A self regularized non-monotonic activation function,” in Proc. of the 31st British Machine Vision Conference (BMVC’20), Virtual Event, UK, September 7-10, 2020. BMVA Press, 2020. [43] G. Singh, T. Gehr, M. Mirman, M. Püschel, and M. T. Vechev, “Fast and effective robustness certification,” in Advances in Neural Information Processing Systems 31: Annual Conference on Neural Information Processing Systems (NeurIPS’18), December 3-8, 2018, Montréal, Canada, 2018, pp. 10 825–10 836. [44] C. Szegedy, W. Zaremba, I. Sutskever, J. Bruna, D. Erhan, I. J. Goodfellow, and R. Fergus, “Intriguing properties of neural networks,” in Proc. of 2nd International Conference on Learning Representations, ICLR 2014, Banff, AB, Canada, April 14-16, 2014, Conference Track Proceedings, Y. Bengio and Y. LeCun, Eds., 2014. [Online]. Available: http://arxiv.org/abs/1312.6199 [45] A. Madry, A. Makelov, L. Schmidt, D. Tsipras, and A. Vladu, “Towards deep learning models resistant to adversarial attacks,” in Proc. of the 6th International Conference on Learning Representations (ICLR’18), Vancouver, BC, Canada, April 30 - May 3, 2018. OpenReview.net, 2018. [46] Z. Zhou, M. Baratchi, G. Si, H. H. Hoos, and G. Huang, “Adaptive error bounded piecewise linear approximation for time-series representation,” Journal of Engineering Applications of Artificial Intelligence, vol. 126, p. 106892, 2023. [47] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. 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. J. Nelson, E. Jones, R. Kern, E. Larson, C. 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. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python,” Nature Methods, vol. 17, pp. 261–272, 2020. [48] R. P. Brent, Algorithms for Minimization without Derivatives, ser. Dover books on mathematics. Mineola, NY: Dover Publications, 2002.

compute one linear lower and one upper relaxations f (x) and f (x) for that particular neuron such that

A PPENDIX A. Extended Description of the Backsubstitution Procedure

f (x) := αx + β ≤ f (x) ≤ αx + β =: f (x)

In this section, we provide more details on the backsubstitution procedure BACKSUBSTITUTE(N , X ). Fundamentally, backsubstitution is a bound propagation procedure. It computes an overapproximate solution to Eq. (2) by successively computing affine upper and lower bounds aT x + d and aT x + d. These symbolic bounds are computed for both activations and preactivations of each layer, in the order nk → n̂k → nk−1 → n̂k−1 → . . . , going backwards from outputs to inputs. We use the symbols ak and dk in affine forms expressed in terms of preactivation values n̂k , and âk and dˆk expressed in terms of activation values nk ; similarly for lower bounds, replacing the overlines by underlines. In parallel to the symbolic bounds, i.e. affine expressions over earlier-layer variables, the procedure maintains concrete bounds, that is scalar intervals [l, u], for the activation and preactivation values of each neuron, as will be detailed below. The procedure starts with an inital affine function in the NN’s output variables y expressing the violation v ∗ of the output property: aT y − d. This also gives us the first upper bound (in fact, it is not an approximation, but exact) in the form aTL n̂L + dL , with aL = a, n̂L = y, and dL = −d. Note that, for the last layer, we only need the upper bound. For all other layers, upper and lower bounds are needed to compute relaxations of the activation functions. Backsubstitution then computes, for decreasing layers, new values ak−1 and dk−1 as well as ak−1 and dk−1 for affine, overapproximating upper and lower bound functions aTk−1 n̂k−1 + dk−1 and aTk−1 n̂k−1 + dk−1 that only depend on variables n̂k−1 from the previous layer. Analogous computations are interleavingly performed to bound nk−1 . Repeating this, we finally arrive at affine approximations that only contain input variables. The computations, in detail, are as follows: n̂k → nk−1 : Going from preactivation to activation of the T previous layer, defining âk−1 and dˆk−1 by:

for all x ∈ [l, u]. These relaxations, gathered for all neurons on level k, are then employed in Eq. (30). As the backsubstitution process results in a relaxation that only depends on the NN’s input variables, it can, together with the hyperrectangle-constraints of the input region X , be used to compute an upper bound v ∗ of the maximum violation v ∗ . a) Invariant: That v ∗ is indeed an upper bound is justified by the following invariant of the backsubstitution process: aTk n̂k + dk ≥ aT y − d

(32)

aTk−1 n̂k−1 + dk−1 ≥ aTk n̂k + dk , maintaining the invariant for decreasing values of k. We thus finally obtain aT1 n̂1 + d1 ≥ aT y − d, and from this, by Equation (29) and n0 = x, T â0 x + dˆ0 ≥ aT y − d

for all x ∈ X

(33)

Thus, if â0 x + dˆ0 ≤ 0, then so is aT y − d for all x ∈ X . Moreover, a concrete bound, v ∗ , on v ∗ can be computed for X = [x, x] using T v ∗ = [â0 ]T+ x + [â0 ]T− x + dˆ0 ≥ â0 x + dˆ0 .

(34)

If v ∗ ≤ 0, then the overapproximate approach is able to prove the property — otherwise, no conclusions can be drawn. The preactivation bounds necessary for the above steps are computed by recursively applying the steps in Equations (29)(30) and the concretization step (34) for each neuron in the earlier layers — once to compute an upper bound on its preactivation values and once to compute a lower bound. B. Closed-form Solutions for Critical Points In the following, we derive closed-form expressions for the candidate critical points used in our shifting procedure. Given an activation function f and a candidate slope m, we consider g(x) = f (x) − mx and compute the set of critical points X̂ ∗ ⊆ [l, u] such that

(29)

nk → n̂k : From activation to preactivation, defining aTk and dk : T âk nk + dˆk ≤ [âk ]T+ f (n̂k ) + [âk ]T− f (n̂k ) + dˆk = [âk ]T+ α ⊙ n̂k + β + [âk ]T− α ⊙ n̂k + β + dˆk

= aTk n̂k + dk ,

for 1 ≤ k ≤ L,

The invariant obviously holds for k = L by definition. Using Equations (29) and (30) we have

aTk n̂k + dk = aTk (Wk nk−1 + bk ) + dk = aTk Wk nk−1 + aTk bk + dk = âk−1 nk−1 + dˆk−1 .

(31)

max g(x) = x∈[l,u]

max

g(x) .

(35)

x∈(X̂ ∗ ∪{l,u})

For functions with non-constant derivative, X̂ ∗ coincides with the preimage of the derivative in Proposition 1, in other cases, it can also contain the break-points of piecewise functions. Evaluating g(x) at the points in X̂ ∗ determines the minimal shift b such that U (x) = mx + b is a sound upper (or lower) relaxation of f (x) on [l, u]. In some cases, we distinguish maximizers and minimizers by X̂ + and X̂ − , respectively; but in our implementation we simply evaluate all candidates X̂ ∗ = X̂ + ∪ X̂ − ∪ {l, u}.

(30)

where [x]+ = max(0, x) and [x]− = min(0, x), applied element-wise, denote the positive and negative parts of x. The computations for lower bounds are, mutatis mutandis, the same. In this step, we also need symbolic upper and lower bounds f and f of the vectorized activation function f for concrete input ranges of each neuron. So, if a concrete bounding interval [l, u] is available for a neuron’s preactivation value n̂, we can

interval Ii can be either ei (ai ) or ei (bi ), where ai and bi are the endpoints of interval Ii . Then, for indices k1 and k2 , such that l ∈ Ik1 and u ∈ Ik2 , k1 < k2 , the maximum of the error function e(x) over the domain [l, u] is:

1) Piecewise Linear Activation Functions: In this subsection we present a general solution to an arbitrary piecewise linear activation function, consisting of n linear pieces. The function does not need to be continuous, but we expect that it’s domain covers the whole set of real numbers (i.e., f : R → R). We consider a finite decomposition of R into nS∈ N>0 nonn empty disjoint intervals I1 , I2 , . . . , In , such that i=1 Ii = R and Ii ∩ Ij = ∅ for all i, j ∈ {1, 2, . . . , n} with i ̸= j. We use ai and bi to refer to the lower resp. upper bound of the interval Ii . Each of the intervals might be either open or closed in each direction, but we assume w.l.o.g. that the intervals are ordered according to the order on R: I1 is unbounded from below (i.e. I1 is left-open with a1 = −∞), In is unbounded from above (i.e. In is right-open with bn = ∞), and bi = ai+1 for i = 1, . . . , n − 1 (note that the property of decomposition implies that Ii is right-closed if and only if Ii+1 is left-open). A piecewise linear function is composed of n linear function, each of which is applied over one of these intervals.

b = max

max

1 . 1 + e−x

g ′ (x) = σ(x)(1 − σ(x)) − m = 0 2

⇔ − σ(x) + σ(x) − m = 0

(40)

(42) let s = σ(x)

⇔−s +s−m=0

(43) (44)

By the quadratic formula we get the candidate solutions r 1 1 −m, (45) s1,2 = ± 2 4 if m ∈ [0, 1/4]. To ensure there is an x, s.t. si = σ(x), we need to check if si ∈ (0, 1). Note that this condition is fulfilled iff. m ∈ (0, 1/4]. The set of critical points is then X̂ ∗ = {x | si ∈ (0, 1) ∧ x = σ −1 (si )} .

(46)

3) Tanh: The proof for tanh proceeds similarly to the proof for σ. The tanh function is an injection on R and a surjection to (−1, 1). Therefore it also has an inverse function arctanh on (−1, 1). We want to find the extrema of g(x) = tanh(x) − m · x with derivative g ′ (x) = 1 − tanh(x)2 − m. The critical points thus have to satisfy

(38)

ei (x) .

We want to find the extrema of g(x) = σ(x) − m · x with derivative g ′ (x) = σ(x)(1 − σ(x)) − m. The critical points satisfy

g ′ (x) = 1 − tanh(x)2 − m = 0

for some linear functions ei : Ii → R, ei (x) = fi (x) − m · x for i = 1, . . . , n. Our goal is to find the maximum (or minimum for finding a lower relaxation) of the error function over the domain [l, u]: max

∪

It is an injection on R and a surjection to (0, 1). Therefore it has an inverse function on (0, 1): 1 σ −1 (x) = − log −1 . (41) x

In this case, considering the candidate line y = m · x, the error is also piecewise linear function of the following form:

i∈{1,2,...,n} x∈Ii ∩[l,u]

∪

{ek (bk ) | k1 < k < k2 }

σ(x) =

the graph of f over [l, u].

x∈[l,u]

{ek (ak ) | k1 < k < k2 }

2) Sigmoid: The sigmoid function is defined as

for some linear functions fi : Ii → R, fi (x) = mi · x + ni with mi , ni ∈ R for i = 1, . . . , n. For any interval [l, u] ⊆ R, we call F = (x, f (x))T x ∈ [l, u] ⊆ R (37)

b = max e(x) =

∪

{ek2 (ak2 ), ek2 (u)}

Definition 1 (Piecewise linear function). A piecewise linear function is a function of type f : R → R such that there exists a finite decomposition of R into n ∈ N>0 non-empty intervals I1 , I2 , . . . , In with  f1 (x) if x ∈ I1     f2 (x) if x ∈ I2 (36) f (x) = . .. ..   .    fn (x) if x ∈ In

Definition 2 (Piecewise error function).   e1 (x) if x ∈ I1   e2 (x) if x ∈ I2 e(x) = . .. ..   .    en (x) if x ∈ In

{ek1 (l), ek1 (bk1 )}

⇔ −t + 1 − m = 0

(47) (let t = tanh(x))

(48)

with solutions √ t1,2 = ± 1 − m ,

(39)

(49)

which is defined for m ∈ [0, 1]. To ensure that there is an x s.t. ti = tanh(x), we need to restrict our solutions to ti ∈ (−1, 1). Note that this condition is fulfilled iff. m ∈ (0, 1]. The set of critical points is then

Since each piece of e(x) is linear, there are no local extrema inside the domain Ii of a single piece ei (x), thus the maximum of a single piece of the error function can happen only at the endpoints of its domain, i.e. the maximum of ei (x) on the

X̂ ∗ = {x | ti ∈ (−1, 1) ∧ x = arctanh(ti )} .

(50)

4) Abs: Since the absolute value function ( −x if x < 0 abs(x) = |x| = x if x ≥ 0

8) HardSigmoid, ReLU, LeakyReLU, Sign: All these functions are piecewise linear, thus our solution from Section B1 should work here.   ifx ≤ −3 0 (64) HardSigmoid(x) = 1 ifx ≥ +3  x 1 + otherwise 6 2 ( 0 if x < 0 (65) ReLU(x) = x otherwise ( γx if x < 0 LeakyReLU(x) = (66) x otherwise   −1 if x < 0 Sign(x) = +1 if x > 0 (67)   0 otherwise

(51)

is piecewise linear, we can just apply the results of Section B1. 5) Acos: We have g(x) = arccos(x) − m · x and g ′ (x) = −1 √ − m. The critical points can be found by setting the 1−x2 derivative equal to zero: −1 −m=0 1 − x2 −1 ⇔ √ =m 1 − x2 p 1 ⇔ 1 − x2 = − m 1 ⇔ 1 − x2 = 2 m 1 ⇔ x2 = 1 − 2 m r 1 ⇔ x=± 1− 2 m

g ′ (x) = √

(52) (53) , m ̸= 0

(54)

,m < 0

(55)

,m < 0

(56)

,m < 0 .

9) HardSwish: The HardSwish activation is defined as   if x ≤ −3 0 HardSwish(x) = x (68) if x ≥ 3   x(x+3) otherwise. 6

(57)

The function g(x) has no discontinuities or piecewise definitions, so the set of critical points are ( q ± 1 − m12 , m < 0 ∗ X̂ = (58) ∅ , otherwise .

We consider g(x) = HardSwish(x) − mx. On (−3, 3) we have x(x + 3) − mx, 6 2x + 3 g ′ (x) = − m. 6 Thus, any interior critical point must satisfy g(x) =

6) Cos: We want to solve for the critical points of g(x) = cos(x) − m · x, with derivative g ′ (x) = −sin(x) − m. The solution is analogous to the solution for the sine function described in Appendix B10. 7) Exp: We have g(x) = ex −x·m. With derivative g ′ (x) = x e − m. Therefore, critical points are ′

g (x) = e − m = 0

g ′ (x) = 0 ⇔

(59)

⇔ ex = m m>0

(61)

The second order derivative g ′′ (x) = ex , which is purely positive, means that the error function is convex on the whole set of R. Thus, x = log(m) is a global minimum, and it cannot be the x value which maximizes the error function on the bounded domain [l, u]. The function g(x) has no discontinuities or piecewise definitions, so the set of critical points when maximizing the error function is X̂ + = ∅

(62)

Note that for minimizing the error function, the set of critical points is ( X̂

−

log m ∅

,m > 0 , otherwise

(70)

(71)

The point is admissible iff. x∗ (m) ∈ (−3, 3), equivalently m ∈ (− 21 , 32 ). Moreover, g ′′ (x) = 13 > 0 on (−3, 3), hence x∗ (m) is a (strict) local minimum. Therefore, when maximizing g over a bounded interval [l, u], no interior critical point needs to be considered and it suffices to check only the interval endpoints and the breakpoints {−3, 3} (when they lie in [l, u]). Conversely, when minimizing g, the candidate critical point is x∗ (m) whenever m ∈ (− 21 , 32 ) and x∗ (m) ∈ [l, u]. 10) Sin: We have g(x) = sin(x) − m · x with derivative g ′ (x) = cos(x) − m. Assuming that the co-domain of arccos is [0, π] (with [−π/2, π/2] being another popular convention), i.e. arccos : [−1, 1] → [0, π], the critical points inside [0, π] can be computed as

(60)

⇔x = log m,

2x + 3 3 = m ⇔ x∗ (m) = 3m − . 6 2

(69)

g ′ (x) = cos(x) − m = 0

(72)

⇔ cos(x) = m

(73)

⇔ x = arccos(m)

, m ∈ [−1, 1] .

(74)

However, there can be critical points outside of [0, π] as well - to find those, we use symmetry and periodicity of the cos function.

(63)

For every x = arccos(m), also cos(−x) = m since cos(−x) = cos(x). Similarly, if x is a solution, then any x+2π is another valid solution due to periodicity. We can avoid enumerating all infinitely many solutions by taking the bounds [l, u] of the approximation domain into account. Since we are only interested in solutions l ≤ x ≤ u inside the approximation domain, we can restrict our search to the set {x + 2πk | k = kmin , . . . , kmax }

The next three lemmas can be proven analogously to Lemma 2: Lemma 3. Let m < 0, then the sequence of local minima g(xk ) for minimizers xk = − arccos(m)+2kπ,

Lemma 4. Let m > 0, then the sequence of local maxima g(xk ) for maximizers

l−x 2π u−x kmax = . 2π kmin =

(76) xk = arccos(m)+2kπ, (77)

xk = − arccos(m)+2kπ,

1) m ≤ −1: The error function g(x) is monotonically increasing. 2) −1 < m < 0: The error function itself is not monotonous, but the sequence of local maxima (or minima) strictly increases (Lemmas 2,3). 3) m = 0: The error function is the same as the sin(x) function. 4) 0 < m < +1: The error function itself is not monotonous, but the sequence of local maxima (or minima) strictly decreases (Lemmas 4,5). 5) m ≥ +1: The error function is monotonically decreasing.

Lemma 2. Let m < 0, then the sequence of local maxima g(xk ) for maximizers

In cases (1) and (5) due to the error function being monotonic, minimum and maximum are attained at the boundary values l, u and there are no interior critical points, i.e. X̂ + = X̂ − = ∅ In case (3), we all points in M+ = { π2 + kπ | k = kmin , . . . , kmax } attain the maximum value and all points in M− = {− π2 + kπ | k = kmin , . . . , kmax } attain the minimum value. We set X = M+ ∪ {l, u} (for the upper bound) and X = M− ∪ {l, u} (when finding the lower bound). In cases (2) and (4), the error function is not monotonic, but the sequence of maxima and the sequence of minima are monotonic. In case (2), we set X̂ + = {2kmax π + arccos(x)} and X̂ − = {2kmin π − arccos(x)}. In case (4), we use X̂ + = {2kmin π + arccos(x)} and X̂ − = {2kmax π − arccos(x)}. √ 11) Sqrt: q We have g(x)√= x − x · m, with derivative g ′ (x) = 21 x1 − m. Since x needs to be defined, we also

(78)

is (strictly) monotonically increasing in k. Proof. We now show g(xk ) < g(xk+1 ) for all k ∈ Z. For any a ∈ R, we have (79) (80)

⇔ − m · (a + 2kπ) ≤ −m · (a + 2(k + 1)π)

(81)

⇔ a + 2kπ ≤ a + 2(k + 1)π ,

(82)

(85)

Now we analyze five cases based on the value of m:

Before we continue with our analysis, we state four lemmas about sequences of the local minima and maxima.

≤ sin (a + 2(k + 1)π) − m · (a + 2(k + 1)π)

k∈Z

is (strictly) monotonically decreasing in k.

M− = {− arccos(m)+2kπ | k = kmin , . . . , kmax , −1 < m < 1}

sin (a + 2kπ) − m · (a + 2kπ)

(84)

Lemma 5. Let m > 0, then the sequence of local minima g(xk ) for minimizers

M+ = {arccos(m)+2kπ | k = kmin , . . . , kmax , −1 < m < 1}

k∈Z

is (strictly) monotonically decreasing in k.

Indeed, we show that we can restrict our search to evaluation of at most two interior critical points in all cases, where m ̸= 0. As a first step, we can utilize curvature information. The second order derivative is given by g ′′ (x) = −sin(x) and is always positive on the intervals [2kπ, (2k+1)π] and always negative on the domains [(2k+1)π, (2k+2)π], for k ∈ Z. Thus, if arccos(m) ∈ [2kπ, (2k + 1)π], then it is a local maximum, if arccos(m) ∈ [(2k+1)π, (2k+2)π], then it is a local minimum. Since the co-domain of the arccos(·) function is [0, π], we are always in the first case, i.e., x = arccos(m) ∈ [2k ·π, (2k +1)· π] and x+2kπ is therefore a local maximum. Similarly, we are always in the second case for x = − arccos(m) and x+2kπ is therefore a local minimum. Denoting the sets of local maxima and minima as M+ and M− respectively, we obtain

xk = arccos(m)+2kπ,

(83)

is (strictly) monotonically increasing in k.

(75)

for x = ± arccos(m), where

k∈Z

where the last statement is obviously true, Equation (81) holds due to periodicity of the sin function and Equation (82) holds as m < 0. The statement then follows from choosing a = arccos(m).

Algorithm 2 Construction of the PWL approximation. Require: Function f , domain [l, u], error bound ϵtol , N ∈ R+ Ensure: Set Ξ of segment boundaries of f˜pwl such that |f˜pwl (x) − f (x)| ≲ ϵtol (approx. smaller than ϵtol ) 1: function APPROX PWL(f, l, u) 2: lo ← l ▷ lower limit for root finding u−l ▷ upper limit for root finding 3: hi ← N 4: while LINA PPROX E RROR(f, l, hi) < ϵtol do 5: lo ← hi 6: hi ← min (u, hi + (hi − l)) ▷ double step size 7: end while 8: ξ ← ROOT (LINA PPROX E RROR(f, lo, hi) − ϵtol ) 9: ▷ if no root in interior, return hi 10: if ξ < u then return {l, ξ} ∪ APPROX PWL(f, ξ, u) 11: else return {l, u} 12: end function

assume x ≥ 0. The critical points can be derived as r 1 1 ′ g (x) = −m=0 (86) 2 rx 1 1 ⇔ =m (87) 2 r x 1 ⇔ = 2m (88) x 1 ⇔ = 4m2 ,m ≥ 0 (89) x 1 ,m > 0 (90) ⇔x= 4m2 q The second order derivative g ′′ (x) = − 14 x13 is strictly negative, thus the error function is concave on the set of reals 1 R. Thus, the computed extrema x = 4m 2 if m > 0 is the global maximum of the function. Since g(x) has no discontinuities or piecewise definitions, so the set of critical points is just ( 1 if m > 0 2 + X̂ = 4m (91) ∅ otherwise X̂ − = ∅ .

D. Proofs In this section we provide the proofs for Proposition 1, Proposition 2, Lemma 1 and Lemma 2. Proposition 1. Let f : R → R be differentiable over the interval [l, u] and let f ′−1 (y) = {x | f ′ (x) = y} be the preimage of the derivative of f . Then

(92)

C. Piecewise Linear Approximation

y = max f (x) − m · x x∈[l,u]

While the employed widening approach ensures soundness of fpwl for any PWL approximation f˜pwl , we achieve better results, when f˜pwl closely approximates f . However, minimizing the approximation error may result in a large number of linear segments increasing the runtime later on. Therefore, we try to make each linear segment as long as possible while maintaining an approximation error smaller than a given tolerance ϵtol . We are aware that finding good PWL approximations with a low number of segments that maintain an error bound is a known problem [46], however, we did not find the exact algorithm we used in the literature, so we state it here. Our algorithm for constructing such PWL approximations is shown in Alg. 2. Our algorithm maintains bounds [lo, hi], which contain the next segment boundary. It starts at the lower bound l and iteratively evaluates the approximation error over intervals [l, hi], doubling the distance to hi each time until the error exceeds the tolerance ϵtol . At that point, there must exist a value between the previous and current hi where the error equals ϵtol ; this point ξ is selected as a segment boundary of f˜pwl , and the process continues on [ξ, u]. Since the subroutines LINA PPROX E RROR and ROOT can be called many times while constructing f˜pwl , we do not compute them using verified approaches, but instead rely on SciPy’s [47] implementation of Brent’s methods for optimization and root finding [48] which yield good approximations in practice. The linear function in each segment [ξi , ξi+1 ], for the final PWL approximation as well as in LINA PPROX E RROR, is constructed just as the line connecting the end points (ξi , f (ξi )) and (ξi+1 , f (ξi+1 )).

= max {l, u} ∪ f ′−1 (m) ∩ [l, u]

(93) (94)

Proof. The maximum of a function over an interval [l, u] can occur either at the boundary or in the interior of the interval. If the maximum occurs at x∗ ∈ (l, u), then the derivative at x∗ has to vanish: d f (x∗ ) − m · x∗ = 0 (95) dx ′ ∗ ⇔f (x ) − m = 0 (96) ⇔f ′ (x∗ ) = m

(97)

The solution set to the last equation is exactly the pre-image of f ′−1 (m). However, f ′−1 (m) may contain values outside of [l, u]. Since we are only interested in solutions inside [l, u], we need to intersect with the interval. The overall maximum is then max {l, u} ∪ f ′−1 (m) ∩ [l, u]

g(x) ≥ f (xl ) − Lf (x − xl ) − γx =: ℓimp (x) l

(98)

g(x) ≥ f (xr ) + Lf (x − xr ) − γx =: ℓimp r (x)

(99)

Case γ ≤ 0. Then Lg = Lf − γ. Hence

is pointwise at least as tight as the straightforward lower envelope: g(x) ≥ g(xl ) − Lg (x − xl ) =: ℓstd l (x) std g(x) ≥ g(xr ) + Lg (x − xr ) =: ℓr (x) .

std ℓimp r (x) = f (xr ) + (Lf − γ)(x − xr ) = ℓr (x).

(100)

For the left bound,

(101)

imp std Meaning ℓimp (x) ≥ ℓstd l (x) and ℓr (x) ≥ ℓr (x). l Moreover, if γ > 0, the bottom-right segment ℓimp r (x) is strictly tighter while the bottom-left coincide; if γ < 0, the bottom-left segment ℓimp (x) is strictly tighter while the l bottom-right coincide.

(102)

g(x) ≥zr + Lg (x − xr ) =: ℓstd r (x) .

(103)

(106)

ℓimp (x) = f (xl ) − (Lf + γ)(x − xl ). l

(107)

Thus ℓimp (x) − ℓstd l (x) = −2γ(x − xl ) ≥ 0 l

Proof. For an interval [xl , xr ] and g(x) = f (x) − γx with f being Lf -Lipschitz, let zl = g(xl ) and zr = g(xr ). Straightforward envelope. Since g is Lipschitz with Lg = Lf + |γ|, we have for all x ∈ [xl , xr ]: g(x) ≥zl − Lg (x − xl ) =: ℓstd l (x)

ℓstd l (x) = f (xl ) − (Lf − γ)(x − xl ),

for all x ∈ [xl , xr ]. For x > xl and γ < 0, even strict inequality holds. Hence the left-anchored segment is strictly tighter, while the right-anchored one coincides. Showing that the improved upper segments are also as tight as the straightforward ones works analogously. Lemma 1. Let l < u, f : R → R, amax the maximum slope of its upper convex hull H(x) over x ∈ [l, u], ϵ > 0 and U (x) = (amax + ϵ)x + b ≥ f (x) be a linear upper relaxation over [l, u]. Then there is another linear overapproximation Û (x) with slope amax , s.t. f (x) ≤ Û (x) < U (x) over x ∈ (l, u].

Improved envelope. From Lipschitz continuity of f , f (x) ≥ f (xl )−Lf (x−xl ) and f (x) ≥ f (xr )+Lf (x−xr ). Subtracting γx yields two valid lower bounds for g: g(x) ≥ f (xl ) − Lf (x − xl ) − γx =: ℓimp (x) , l

Proof. Recall that the upper convex hull of f (x) over [l, u] is the smallest concave function H(x) ≥ f (x) for x ∈ [l, u]. If U (x) > H(x) for all x ∈ (l, u], we can construct a dominating function Û (x) by using a smaller bias term. Note that U (x) = (amax +ϵ)x+b cannot touch H(x) in two or more points, otherwise it would be a secant and thus tangent to the upper convex hull. It can also not be tangent to H(x) or it would also be part of the upper convex hull. Therefore, U (x) = H(x) is only possible for either x = l or x = u. Note that

and g(x) ≥ f (xr ) + Lf (x − xr ) − γx =: ℓimp r (x) . Now rewrite the improved bounds in terms of f (xl ) and f (xr ). For the left bound, ℓimp (x) = f (xl ) − Lf (x − xl ) − γx l = f (xl ) − γxl − (Lf + γ)(x − xl ) = f (xl ) − (Lf + γ)(x − xl ). For the right bound,

H(x) − H(l) ≤ amax (x − l)

(108)

U (x) − U (l) = (amax + ϵ)(x − l)

(109)

due to amax being the maximum slope of H(x) and U (x) being a linear function. We first assume U (u) = H(u). Let x = u in Equations (108),(109), then Since ϵ > 0 and (u − l) > 0, we know that amax (u − l) > (amax + ϵ)(u − l) and thus

ℓimp r (x) = f (xr ) + Lf (x − xr ) − γx = f (xr ) − γxr + (Lf − γ)(x − xr ) = f (xr ) + (Lf − γ)(x − xr ). We now distinguish two cases.

U (u) − U (l) ≥ H(u) − H(l)

Case γ ≥ 0. Then Lg = Lf + γ. Hence

⇔ U (u) − H(u) > U (l) − H(l) .

ℓimp (x) = f (xl ) − (Lf + γ)(x − xl ) = ℓstd l (x). l

But since U (u) = H(u) this is equivalent to

For the right bound, 0 > U (l) − H(l) ⇔ H(l) > U (l)

ℓstd r (x) = f (xr ) + (Lf + γ)(x − xr ),

(104)

ℓimp r (x) = f (xr ) + (Lf − γ)(x − xr ).

(105)

(110)

which is a contradiction to H(x) being the smallest concave upper bound to f (x). Therefore, U (l) = H(l). Note that, since they touch at l, we have U (x) = (amax + ϵ)(x − l) + H(l). But then, we can construct

Thus std ℓimp r (x) − ℓr (x) = −2γ(x − xr ) ≥ 0

for all x ∈ [xl , xr ]. For x < xr and γ > 0, even strict inequality holds. Hence the right-anchored segment is strictly tighter, while the left-anchored one coincides.

Û (x) = amax (x − l)+H(l) < amax (x − l)+H(l)+ϵ(x − l) = U (x)

for all x ∈ (l, u]. Note that amax (x − l) + H(l) ≥ H(x) ⇔ amax (x − l) ≥ H(x) − H(l) ,

avoids the need for sorting or ordered insertion operations on tensors. During refinement, a selected segment entity is replaced by its left subsegment, while the corresponding right subsegment is appended to the end of the stored lists. This design significantly simplifies the minorant improvement step. Second, since tensor concatenation in PyTorch is a costly operation, we preallocate the tensors xl , zl , xr , zr , xm , zm to their maximum required size, namely n − 1 + max_iter. Although this strategy incurs additional memory usage, it enables substantially faster refinement steps by avoiding repeated reallocation.

(111) (112)

where the right hand side is true by Equation (108). Therefore, Û (x) ≥ H(x) ≥ f (x) and Û (x) is still an upper relaxation. Lemma 2. Let l < u, f : R → R, amin the minimum slope of its upper convex hull H(x) over x ∈ [l, u], ϵ > 0 and U (x) = (amin − ϵ)x + b ≥ f (x) be a linear upper relaxation over [l, u]. Then there is another linear overapproximation Û (x) with slope amin , s.t. f (x) ≤ Û (x) < U (x) over x ∈ [l, u).

F. Supplementary Evaluation Material 1) Parameter Settings: This section summarizes the parameter values used in our experimental evaluation. Unless stated otherwise, all experiments use the same default configuration. a) Initialization: For the gradient-based initialization of the relaxation slopes we perform 100 iterations of gradient descent. b) Optimization of relaxation parameters: During bound optimization we run 20 iterations of gradient-based optimization with learning rate 0.05. c) Piyavskii refinement: For constructing the overapproximation we employ the Piyavskii–Shubert algorithm with a maximum of 1000 iterations or until the tolerance 10−4 is reached. The same tolerance 10−4 is used for the SLiR approximation procedure. d) Additional parameters: Unless explicitly overridden, the following default parameters are used in our implementation: • initialization method: gradient-based • initialization iterations: 100 −4 • overapproximation tolerance: ϵtol = 10 • Piyavksii maximum iterations: max_iter = 1000 −4 • Piyavskii tolerance: ε = 10 • Piyavskii initial points: nini = 100 • local Lipschitz constant estimation enabled • improved Piyavskii envelope enabled −6 • sigmoid transformation epsilon: ϵσ = 10 All remaining parameters use the default values provided by our implementation. 2) Lipschitz Optimization and Bound Tightness: While the improvements to the Piyavskii algorithm described in Section III-C did not lead to a large increase in verified properties for the MNIST benchmark, the effects are noticable, when having a closer look at the verifiable output bounds. Figure 6 shows the mean of the summed width of the output bounds X ui − li , (119)

Proof. By a similar argument to the proof of Lemma 1, we know that U (x) = H(x) is only possible for either x = l or x = u. Now note that H(u) − H(x) ≥ amin (u − x)

(113)

U (u) − U (x) = (amin − ϵ)(u − x)

(114)

We first assume that U (l) = H(l). Let x = l in Equations (113),(114), then as ϵ > 0, we know that amin (u − x) > (amin − ϵ)(u − x) and thus H(u) − H(l) > U (u) − U (l) ⇔ H(u) − U (u) > H(l) − U (l) = 0 which would lead to H(u) − U (u) > 0 ↔ H(u) > U (u) ≥ f (u)

(115)

and thus contradicts the assumption that H(u) is the smallest concave upper bound to f (x). Therefore, U (u) = H(u). Since they touch at u, we can equivalently rite U (x) = (amin − ϵ)(x − u) + H(u). We can then construct a dominating upper bound Û(x) = amin (x − u) + H(u) < amin (x − u) + H(u) + ϵ(x − u) = U (x) , where the last step holds because x − u < 0 for x ∈ [l, u). Note that amin (x − u) + H(u) ≥ H(x)

(116)

⇔ amin (x − u) ≥ H(x) − H(u)

(117)

⇔ amin (u − x) ≤ H(u) − H(x)

(118)

where the right hand side is true by Equation (113). Therefore, Û (x) ≥ H(x) ≥ f (x) and Û (x) is still an upper relaxation. E. Implementation Details of the Piyavskii Method In order to implement the Piyavskii method in a tensorfriendly and computationally efficient manner, we introduce several improvements. First, instead of maintaining an ordered list of trial points, we store a collection of minorant segment entities, each defined by its left and right endpoints and the corresponding minorant intersection point. This representation

where li , ui are the verifiable lower and upper bounds returned by the different configurations of SLiR for the i-th output of a network for both after initialization and after optimization. The mean is taken over all considered images and networks.

sum of widths of output bounds (mean)

initialization — can verify significantly more properties than the example guided approach [13] for the L I SHT and M ISH activation functions for the MNIST benchmark. Similarly, we can verify more properties than α-CROWN for the GELU and tanh activation functions — both after initialization and after optimization — even though, α-CROWN utilizes manually derived custom relaxations in these cases.

width sum (opt) width sum (ini)

260.0 257.5 255.0 252.5 250.0 247.5 245.0

tanh(x) aCROWN lower aCROWN upper lower (ours) upper (ours)

242.5 2

SLiR (local)

SLiR (envelope)

SLiR

SLiR (base)

Fig. 6: Effect of different improvements on the Piyavskii method on the width of the provable output bounds of the MNIST networks. All modifications improve upon the baseline. Combining both local Lipschitz constants (SLiR (local)) and the improved envelope computation (SLiR (envelope)) leads to further improvements (SLiR).

sum of widths of output bounds (mean)

55 54 53 52 SLiR (local)

SLiR (envelope)

In this section, we show specific failure cases, where the (initial) relaxations chosen by the example guided approach and α-CROWN are not as tight with respect to the area enclosed between the lower and the upper relaxation — leading to suboptimal performance of those tools for these benchmarks. While the E-GUIDED approach is generally capable of synthesizing tight linear relaxations, there exist configurations in which the combination of subinterval partitioning and tangent point estimation results in comparatively loose bounds. Figure 10 illustrates representative cases of this behavior. In Figures 10a and 10b, the E-GUIDED method employs the same subinterval partitioning l ∈ [−9.0, −8.0] and u ∈ [−1.0, 0.0], together with the same tangent-line template. For intervals in this partition, the tangent line is first constructed locally based on the predicted tangent point, and is then shifted upward by the maximum violation computed over the entire subinterval in order to ensure soundness. In Figure 10a, the primary source of looseness is the coarse partitioning of the input interval space. Even when the locally constructed tangent is close to tight for the concrete interval, the required global shift—computed with respect to all intervals in the subpartition—introduces a substantial degradation of the upper bound. In contrast, Figure 10b demonstrates a case in which the looseness is mainly attributable to suboptimal tangent point estimation. Although the required soundness shift is relatively small, the predicted tangent point itself does not induce a tight upper relaxation, and the resulting bound remains conservative.

width sum (opt) width sum (ini)

SLiR (base)

Fig. 8: Example where α-CROWN fails to find a good initial overapproximation for the tanh function. Enclosed area: αCROWN: 7.4823, ours: 4.5441

Each of our modifications - using local Lipschitz constants and improved envelopes - leads to SLiR producing tighter output bounds than when just using the baseline Piyavskii method. Moreover, combining our modifications - as used in the final SLiR algorithm - leads to an additional reduction in bounds width. For the CIFAR benchmarks results are similar (see Figure 7).

SLiR

Fig. 7: Effect of different improvements on the Piyavskii method on the width of the provable output bounds of the CIFAR networks. All modifications improve upon the baseline. Combining both local Lipschitz constants (SLiR (local)) and the improved envelope computation (SLiR (envelope)) leads to further improvements (SLiR). 3) Failure Cases of Other Methods: The experimental results in Sections V showed that our approach — even after

1.00 0.75

gelu(x) aCROWN lower aCROWN upper lower (ours) upper (ours)

0.25

0 y

0.50

0.00

0.25

0.50 1.0

0.5

0.0

0.5

1.0

(a) α-CROWN: 2.1758, ours: 1.7261

gelu(x) aCROWN lower aCROWN upper lower (ours) upper (ours)

15 1.5

gelu(x) aCROWN lower aCROWN upper lower (ours) upper (ours)

(b) α-CROWN: 208.0164, ours: 173.3838

Fig. 9: Cases where α-CROWN fails to find good initial overapproximations for the GELU function. The area enclosed between the lower and upper relaxations can be found in the caption of each subfigure.

Moreover, Figure 10c shows that even in situations where the E-GUIDED relaxation appears reasonably tight, our approach can still derive a slightly tighter relaxation, resulting in a smaller enclosed area between lower and upper bounds. Finally, we observe that the E-GUIDED approach may yield particularly loose bounds when the queried interval lies outside the predefined input universe ([−10, 10] in the shared artifact). In this case, the method reverts to its decompositionbased bounding procedure rather than employing the synthesized template-based relaxations. Such decomposition-based bounds are generally more conservative, and therefore tend to be significantly looser than the bounds obtained within the designated input universe.

0.3

0.1 0.2 0.3

(a) Input domain [−8.70, −0.50]. E-GUIDED lower: y = −0.0415x − 0.3624, E-GUIDED upper: y = −0.0279x − 0.0811, E-GUIDED tangent: y = −0.0279x − 0.1729, ours lower: y = −0.0415x−0.3627, ours upper: y = −0.0309x − 0.1876. Relaxation area: E-GUIDED 1.7937, ours 1.0357.

0.05 mish(x) E-GUIDED lower E-GUIDED upper E-GUIDED tangent lower (ours) upper (ours) 8 7 6

mish(x) E-GUIDED lower E-GUIDED upper E-GUIDED tangent lower (ours) upper (ours) 8 6

0.2

0.00

0.0

0.0 0.1

0.05

0.1

0.10

mish(x) E-GUIDED lower E-GUIDED upper E-GUIDED tangent lower (ours) 0.20 upper (ours) 10 9 8 7 0.15

(b) Input domain [−8.01, −0.13]. E-GUIDED lower: y = −0.0454x − 0.3678, E-GUIDED upper: y = −0.0265x − 0.0743, E-GUIDED tangent: y = −0.0265x − 0.1661, ours lower: y = −0.0456x−0.3681, ours upper: y = −0.0100x − 0.0738. Relaxation area: E-GUIDED 1.7066, ours 1.1767.

(c) Input domain [−9.90, −2.29]. E-GUIDED lower: y = −0.0285x − 0.2855, E-GUIDED upper: y = −0.0110x − 0.0810, E-GUIDED tangent: y = −0.0110x − 0.0810, ours lower: y = −0.0289x−0.2863, ours upper: y = −0.0149x − 0.1036. Relaxation area: E-GUIDED 0.7445, ours 0.7409.

Fig. 10: Cases where the E-GUIDED approach fails to obtain a tight overapproximation of the mish(x) activation. The area enclosed between the lower and upper relaxations is reported in the caption of each subfigure. In Figures 10a and 10b, both methods produce identical lower relaxations, while in Figure 10c our method yields a slightly tighter lower bound. The improvement in Figures 10a and 10b stems primarily from significantly tighter upper bounds. In Figure 10c, although the bounds are closer, our approach still achieves a smaller overall relaxation area.

Record · ID 290599 · SHA-256 a141bab12e50f066
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.