Recent years have seen growing adoption of machine learning (ML) models in several safetycritical domains, including autonomous driving [Bojarski et al. 2016] and healthcare [Esteva et al. 2019]. The first derivatives specified via the Jacobian matrix are the backbone of many prominent learning paradigms and are used in all facets of the machine learning pipeline, from training to testing. Automatic Differentiation (AD), and more broadly, Differentiable Programming, have been developed to offer a principled, language-based method of compositionally computing derivatives. However, the needs of ML researchers have rapidly outpaced formal development on the programming languages side. For instance, ML techniques regularly must differentiate through functions that may have points of non-differentiability, such as the commonly used ReLU activation in neural networks, which is not differentiable at 0. In practice, for such a situation, existing AD frameworks just return any arbitrary number in the interval [0, 1]. While such an ad-hoc approach may suffice in many practical scenarios, for critical domains requiring formal certification of properties defined over the first derivatives, this lack of rigor is troubling. To resolve this limitation, programming language researchers have begun using generalizations of the Jacobian for describing AD semantics [Sherman et al. 2021], such as the Clarke Generalized Jacobian [Clarke 1990] for non-smooth, but Lipschitz continuous functions.
However, these developments still cannot provide the desired formal guarantees for practical verification tasks, such as obtaining formal bounds on the Lipschitz constant (which measures how rapidly a function's output changes) of the composition of a non-differentiable perturbation and a differentiable network within an input region. Lipschitz constants can be used as a metric to compare the stability and smoothness of the output of neural networks prior to deployment, as a network with a smaller constant is often preferable [Lin et al. 2019]. Formal bounds on the Lipschitz constant can also be used during training to learn classifiers that are certifiably robust to adversarial perturbations [Tsuzuku et al. 2018], robust to quantizations [Lin et al. 2019], or to improve interpretability by making network explanations themselves more robust [Alvarez-Melis and Jaakkola 2018]. Further, analyzing the Lipschitz constant has direct applications in algorithmic fairness [Dwork et al. 2012] and differential privacy [Dwork et al. 2006], where fairness and privacy are established by certifying a small Lipschitz constant. Beyond using the Jacobian for formally bounding Lipschitz constants, a Jacobian analysis can also be used to formally reason about the local geometry of ML models [Zhang et al. 2019], which has applications in explainability, e.g., why a network works well for certain inputs but not for others.
Challenges. We focus on designing a static analysis that can formally reason about not only Jacobians of functions that are differentiable (e.g., tanh), but can also handle non-differentiable behavior due to functions like ReLU and the branching that can arise in differentiable programs. Furthermore, we must also handle high-dimensional computational graphs with arbitrary arithmetic (instead of, e.g., solely neural networks that avoid non-scalar multiplication and division). Simultaneously handling all of these requires a theoretical formalism that is beyond the scope of the existing works [Edalat and Maleki 2017;Jordan and Dimakis 2020;Mangal et al. 2020;Zhang et al. 2019]. One of the main difficulties arises from the fact that generalized notions of Jacobians do not always obey the same rules as classical Jacobians. For example, the Clarke Jacobian cannot simply be computed by concatenating partial Clarke derivatives into a matrix [Clarke 1990;Khan and Barton 2013]. In this same light, we want our desired generality to also come with compositionality: we want to be able to combine different functions together instead of restricting ourselves to a non-compositional analysis limited to a specific, fixed type of function (e.g., a DNN) or input perturbation. More practically, we want this analysis to be scalable and fully end-to-end, specifically for real-world problems such as analyzing Lipschitz robustness with respect to multiple input perturbations. This has proven challenging, as most formalisms that aim for broad theoretical generality cannot scale beyond toy examples [Di Gianantonio and Edalat 2013], do not have implementations [Edalat and Lieutier 2004;Edalat et al. 2013;Edalat andMaleki 2017, 2018], or lack end-to-end integration with specific analyses for practical problems. Moreover, the Lipschitz robustness analyses that go beyond toy examples, such as RecurJac [Zhang et al. 2019] and ProLip [Mangal et al. 2020], are either tailored for handling fully-connected architectures (like RecurJac) and therefore cannot immediately analyze the state-of-the-art convolutional architectures, or are heavily restricted in the activation functions supported (like ProLip). These issues substantially limit the practical applicability of these tools.
Results. We implement DeepJ and evaluate it on two tasks:
• Lipschitz Robustness: Certifies bounds on the local Lipschitz constant of a given input region.
• Local Optimization Landscape: Uses the Clarke Jacobian to analyze a function's local geometry in a specified input region to determine the absence of stationary points. We apply DeepJ to neural networks with both fully-connected and convolutional layers that are trained on the CIFAR10 [Krizhevsky et al. 2009] and MNIST [LeCun et al. 1998] datasets. For each of our analyses, we compose the neural networks with three perturbation functions: haze (which models images as foggy), contrast variation (which accentuates the differences between bright and dark pixels), and rotation (which rotates the image by a specified angle 𝜃 with bilinear interpolation). Analyzing the Jacobian of a network with respect to these perturbations is out of reach of the existing techniques.
For each perturbation, DeepJ is able to leverage its fully localized analysis, which computes Clarke Jacobians solely for the specified input region, to obtain local Lipschitz constants that can be up to several orders of magnitude smaller than a baseline analysis based on Gouk et al. [2021] (which multiplies a network's global Lipschitz constant by the perturbation's local Lipschitz constant instead of being fully localized). DeepJ can also extend the analysis to compositions of multiple perturbations. Furthermore, the localized Jacobian analysis can certify the absence of stationary points in a network's optimization landscape. DeepJ's parallel CPU-based implementation is efficient: it can precisely analyze 100 CIFAR10 images on our largest convolutional network containing > 62K neurons within a median time of 15 seconds each for haze and contrast variations, and under 1.4 minutes for rotation. It can also compute precise results for 10 CIFAR10 images on the same network within a median time of 9 seconds for contrast variation composed with haze, and under 51 seconds for contrast variation or haze composed with rotation. Our largest network has the same architecture as the one commonly handled by state-of-the-art CPU-based robustness verifiers [Singh et al. 2018[Singh et al. , 2019;;Urban and Miné 2021]. The differences in the computed constants between the DeepJ versions with and without sound floating-point rounding are negligible, with 2.8-4.1x execution time overhead for the sound version.
Contributions. This paper makes the following main contributions:
(1) A new dual number-based interval abstraction for analyzing Clarke Jacobians. Our domain soundly over-approximates the Clarke Generalized Jacobian of locally Lipschitz and piecewise differentiable functions, allowing it to soundly handle functions like max, ReLU, and limited branches. Furthermore, as our abstraction is defined for sets of points, we can analyze local properties of locally Lipschitz functions beyond the scope of prior work. (2) A novel Clarke Jacobian analysis and Lipschitz certification of neural networks composed with non-differentiable perturbations. (3) A scalable, optionally floating-point sound, implementation of our method which supports both convolutional and fully-connected neural networks. (4) An extensive evaluation against multiple perturbation types on several deep neural networks, showing that DeepJ can (a) achieve orders of magnitude tighter bounds on local Lipschitz constants compared to a baseline analysis and (b) certify the absence of stationary points within a given input region.
DeepJ is available at https://github.com/uiuc-arc/DeepJ. Our implementation exploits CPU-level parallelism. We believe that DeepJ can be easily parallelized over GPUs to boost the scalability to even larger architectures, such as those considered by GPU-based verifiers [Müller et al. 2021]. The GPU extension of DeepJ may help train networks to be robust to semantic non-differentiable perturbations like rotation and contrast variation, which is beyond the reach of existing robust training methods [Balunović et al. 2019;Mirman et al. 2018Mirman et al. , 2019;;Zhang et al. 2021Zhang et al. , 2020]].
In this section, we start with a small illustrative example that showcases a real-world use of our abstraction for a scenario not handled by any prior work.
Running Example: Contrast Variation Perturbation. We consider the simple fully-connected network shown in Fig. 1. For simplicity, the network takes two inputs in the input layer, and we assume the network has been fully trained. The network contains a single hidden layer and all activation functions in both the hidden and output layers are tanh (with no biases). Most importantly, we compose the network with a perturbation function modeling contrast variation. We are interested in knowing precisely how sensitive the network's outputs are to inputs that are perturbed by contrast variation, which often arises when the image passed to a network was obtained with a fixed aperture lens [Paterson et al. 2021]. Hence, instead of passing input image pixels 𝑥 1 and 𝑥 2 directly into the input layer, we pass their perturbed values, 𝑥 ′ 1 and 𝑥 ′ 2 . The contrast variation
Hidden Layer Output Layer tanh tanh Fig. 1. Composition of a Neural Network with the Contrast Variation Perturbation perturbation for a pixel 𝑥 𝑖 is given in Paterson et al. [2021] as
where 𝛼 specifies the amount of contrast variation. We apply this perturbation to every input pixel (there are only two in this example). We can thus think of the perturbation as a function of the perturbation parameter 𝛼.
Jacobian Analysis. We will measure the sensitivity of the network with respect to the perturbation by computing the Jacobian of the composition of the network with the perturbation. This allows us to analyze how sensitive and robust the network's outputs are with respect to a change in 𝛼 in some local region. For this example, the local region we are interested in is when 𝛼 ∈[0, 0.1]. This allows us to analyze what happens when the amount of perturbation ranges from none up to a modest amount. Because of the combination of both non-differentiable (max and min) and differentiable (tanh) functions, as well as the division in the perturbation function, computing a Jacobian for the composition of the network and the perturbation is beyond the capabilities of existing frameworks.
Abstract Domain. To perform the analysis, we need to compute bounds on the Jacobian of the composition of the network and the perturbation. However, this is complicated by the fact that (a) we cannot settle for the derivatives at a single point (as standard automatic differentiation gives) and instead need a bound on the derivatives for an entire input region, (b) max and min are Lipschitz, but not differentiable, thus we need a more general notion of differentiation that works for such functions, and (c) the analysis cannot be restricted to only neural networks, since it needs to be able to compositionally handle arbitrary combinations of functions. As mentioned, prior work cannot address these challenges, thus our solution necessitates a novel abstract domain. The full formalism is described in Sections 5 and 6. Our abstract domain associates to each variable an interval bounding the variable itself and an interval bounding that variable's Clarke derivative via a dual interval of the form [𝑎, 𝑏] + [𝑐, 𝑑]𝜖, where 𝑎, 𝑏, 𝑐, 𝑑 ∈ R ∪ {±∞}, 𝑎 ≤ 𝑏 and 𝑐 ≤ 𝑑. We will call [𝑎, 𝑏] the real part and [𝑐, 𝑑] the dual part. Dual intervals are an adapted interval-domain abstraction of the canonical dual numbers of forward-mode automatic differentiation. The key benefit of this approach is that we can leverage an existing numerical system to track derivatives instead of relying on non-extensible, ad-hoc approaches as in Edalat et al. [2013]; Mangal et al. [2020]; Zhang et al. [2019]. However, a naive adaptation is not sufficient, as one still has to contend with non-differentiable functions like max. Furthermore, all primitive operations must be reinterpreted for dual interval arithmetic.
Abstract Interpretation of the Perturbation and Network. We now step through the abstract interpretation of the composition of the neural network and the contrast variation perturbation, which is detailed in Fig. 2. As every term in the abstract domain must be a dual interval, to perform the analysis, we first must lift the constant edge weights 𝑤 𝑖 to dual intervals of the form [𝑤 𝑖 , 𝑤 𝑖 ] + [0, 0]𝜖, as shown in Fig. 2. Upon lifting all constants and edge weights to the abstract domain, we start from the inputs 𝑥 1 , 𝑥 2 and the perturbation parameter 𝛼. As our goal is to compute the sensitivity solely with respect to 𝛼, we set its dual part to [1, 1]𝜖. This means that we treat 𝑥 1 and 𝑥 2 as constant; since we do not need to compute derivatives with respect to 𝑥 1 and 𝑥 2 , the analysis sets their dual part to [0, 0]𝜖. Furthermore, as 𝑥 1 and 𝑥 2 are pixel values of some fixed image, their real parts are just the degenerate interval of their pixel intensities: [0.95, 0.95] and [0.2, 0.2], respectively.
We begin by propagating the dual interval inputs through the contrast variation perturbation function in Eq. 1. This requires that all arithmetic operations (+, -, •, /) and function primitives (min and max) be redefined for dual intervals. Section 6 presents a full formalization of dual intervals.
For each pixel 𝑥 𝑖 , we first compute 𝑥 𝑖 -0.5•𝛼 1-𝛼 inside the min function. Though not shown in the function, the constants 0.5 and 1 are actually interpreted as the dual intervals [0.5, 0.5] + [0, 0]𝜖 and [1, 1] + [0, 0]𝜖, respectively. The numerator 𝑥 1 -0.5•𝛼 abstractly evaluates to [0.9, 0.95] + [-0.5, -0.5]𝜖 and likewise 𝑥 2 -0.5 • 𝛼 evaluates to [0.15, 0.2] + [-0.5, -0.5]𝜖. These follow from the rules of dual interval arithmetic: scaling by a constant scales a term's real and dual part, and dual interval addition between terms adds their respective real and dual components. Therefore, we have implicitly encoded the notion of linearity of the derivative.
To compute the entire quotient for each variable, we next perform dual interval division (described in Section 6.1). The terms 𝑥 1 -0.5•𝛼 1-𝛼 and 𝑥 2 -0.5•𝛼 1-𝛼 ultimately evaluate to [0.9, 1.05]+[0.4, 0.62]𝜖 and [0.15, 0.22] + [-0.43, -0.25]𝜖, respectively. Dual interval division implicitly encodes the quotient rule of differentiation.
Non-Differentiable Functions. Upon computing 𝑥 1 -0.5•𝛼 1-𝛼 and 𝑥 2 -0.5•𝛼 1-𝛼 , we now must take the min of each with [1, 1] + [0, 0]𝜖. This is highly challenging as min is not differentiable. To resolve this, we must use a more general notion of differentiation, specifically the Clarke Jacobian [Clarke 1990], which has also recently emerged in the programming languages literature [Di Gianantonio and Edalat 2013;Sherman et al. 2021]. The Clarke Jacobian can compute a generalized derivative for non-differentiable, but Lipschitz functions like min, max, ReLU, abs, and even functions defined by conditional branching statements (provided one proves such functions are continuous and piecewise differentiable). The Clarke Jacobian does so by returning a convex set of points as the "derivative" instead of just a single point. To perform the Clarke differentiation through the min function, one takes the standard derivative of whichever of the min function's two arguments attains the minimum. The caveat is that if both arguments attain the minimum, one must take the convex hull of both arguments' derivatives. However, this becomes even more difficult for the interval domain, as due to the inherent uncertainty when two intervals overlap, it is possible that either could attain the minimum. Thus, if the real parts of two dual intervals overlap, our analysis must take the convex hull of their respective dual parts.
Recall that for 𝑥 1 we must take the min of [1, 1] + [0, 0]𝜖 and 𝑥 1 -0.5•𝛼 1-𝛼 =[0.9, 1.05]+[0.4, 0.62]𝜖. In this case, the real parts overlap, so we must take the convex hull (denoted by co) of their respective dual parts: [0, 0] and [0.4, 0.62], as shown in Eq. 2.
min 1,
For 𝑥 2 , when computing the min of [1, 1] + [0, 0]𝜖 and 𝑥 2 -0.5 Finally, as the contrast variation perturbation is also composed with the max function, we must repeat the same procedure, albeit with max instead of min. For 𝑥 1 and 𝑥 2 , we compute: The dual part of the result is also multiplied by the dual part of the input ([-0.43, 0.37]), implicitly encoding the chain rule. In this example, 𝑓 3 evaluates to [0.78, 0.84] + [-0.17, 0.14]𝜖. Likewise for 𝑓 4 , the input to the tanh function is the sum of 𝑥 ′ 1 and 𝑥 ′ 2 scaled by the edge weights, which is just [0.67, 0.85] + [0.25, 1.05]𝜖. Hence, 𝑓 4 evaluates to [0.59, 0.69] + [0.13, 0.68]𝜖. As our analysis is fully compositional, we easily repeat this procedure for the subsequent (final) layer. Multiplying 𝑓 3 and 𝑓 4 by their respective edge weights, then passing these values to the tanh activations, allows us to determine that 𝑓 5 =[0.88, 0.91]+[-0.01, 0.19] 𝜖 and 𝑓 6 =[0.09, 0.24]+[-0.84, 0.012]𝜖.
Interval Clarke Jacobian. Upon computing the outputs 𝑓 5 and 𝑓 6 , we take their dual parts as the Interval Clarke Jacobian. This is because we show in Section 6.3 that abstractly evaluating functions with dual intervals is equivalent to computing the Interval Clarke Jacobian, which in turn soundly over-approximates the true Clarke Jacobian. Since we are modeling the composition of the network and the perturbation as a function of only 𝛼 (while holding all other inputs fixed), the Interval Clarke Jacobian is a 2 × 1 interval matrix, as the network has 2 outputs. In this case, it is [[-0.01, 0.19],[-0.84, 0.012]] 𝑇 . Practical Applications. Upon computing this over-approximation of the Clarke Jacobian, we can use it for several practical applications. For instance, we can compute the local Lipschitz constant in the region 𝛼 ∈ [0, 0.1] by taking the maximum norm of the Interval Clarke Jacobian, which intuitively gives us a point summary of the network's sensitivity to perturbations by 𝛼 in this local region. For this example, the Lipschitz constant (with respect to the ℓ ∞ -norm) evaluates to 0.84.
We can also use the over-approximation of the Clarke Jacobian to analyze the local landscape of the composition of the network and the perturbation for stationary points. If a point is a local extremum, then the Clarke Jacobian at that point contains 0. Therefore, if any entry of the Interval Clarke Jacobian does not contain 0, it certifies that no point in the input range is a local extremum. In this example, as both entries contain 0, the analysis determines that the input region 𝛼 ∈ [0, 0.1] could still contain a local extremum.
We define all the mathematical preliminaries needed to describe automatic differentiation, as well as the Clarke Jacobian. We start with the definition of the standard Jacobian.
Definition 3.1. The Jacobian of a function 𝑓 :
When 𝑚 = 𝑛 = 1, the Jacobian is merely the classical derivative: J(𝑓 , x 0 ) = 𝑑 𝑓 𝑑𝑥 | 𝑥=x 0 . Dual Numbers. The question then arises, how do we automatically compute this Jacobian? The most popular method is via Automatic Differentiation, or AD. AD has two modes: reverse and forward. The former recursively evaluates derivatives of sub-expressions, and can be thought of as a generalization of backpropagation. We focus on the latter, as forward-mode automatic differentiation is much easier to implement and excels when the function's input dimension is small (as will be in our use cases). To implement forward-mode automatic differentiation, one may overload all primitive arithmetic operations to work on dual numbers, which we now describe. Definition 3.2. Dual numbers are numbers of the form 𝑎 + 𝑏𝜖, where 𝑎, 𝑏 ∈ R and 𝜖 is a symbolic variable (akin to 𝑖 for imaginary numbers). We denote the set of all dual numbers as D. Dual number arithmetic is given by the following rules:
The above arithmetic rules for dual numbers implicitly encode linearity, the product rule, and the quotient rule in the computation of their dual part. To access the real part of a dual number, we write fst(𝑎 +𝑏𝜖) = 𝑎, and likewise the dual part is accessed by snd(𝑎 +𝑏𝜖) = 𝑏. For any differentiable function 𝑓 : R → R, we use the standard lifting of 𝑓 to dual numbers 𝑓 : D → D given as
Therefore, the dual part of a dual number corresponds to the value of the function's derivative evaluated at the real part, 𝑎. Further, multiplying the derivative 𝑓 ′ (𝑎) by the existing dual part 𝑏 implicitly encodes the chain rule of calculus.
A Dual Number Abstraction for Static Analysis of Clarke Jacobians 56:9
We subsequently show how to extend the previous concepts to non-differentiable but locally Lipschitz functions. Hence, we first define the local Lipschitz property.
where ∥ • ∥ 𝛼 and ∥ • ∥ 𝛽 are arbitrary 𝑝-norms over R 𝑚 and R 𝑛 , respectively. Furthermore, if for a given point 𝑥 0 ∈ R 𝑚 , there exists a positive real 𝛿 > 0 such that 𝑓 is locally Lipschitz within a ball of radius 𝛿 centered at 𝑥 0 , we say 𝑓 is Lipschitz near 𝑥 0 .
Lipschitz Constant. The constant 𝐾 𝛼,𝛽 is called the (local) Lipschitz constant, which provides a formal bound on how much a function's output can change (measured by the ∥ • ∥ 𝛽 norm) given a change in input (measured by the ∥ • ∥ 𝛼 norm). The constant can easily be obtained once one has the Jacobian J. For a Lipschitz function 𝑓 : R 𝑚 → R 𝑛 that is differentiable, one can compute the local Lipschitz constant on a region 𝑋 by taking the maximum dual norm of the Jacobian over 𝑋 :
where the dual norm of any matrix 𝑀 ∈ R 𝑛×𝑚 is given as ∥𝑀 ∥ 𝛼,𝛽 = sup ∥𝑣 ∥ 𝛼 ≤1 ∥𝑀𝑣 ∥ 𝛽 . For common values of 𝛼 and 𝛽, there is a closed-form expression for ∥𝑀 ∥ 𝛼,𝛽 . For instance, ∥𝑀 ∥ 1,1 = max 1≤ 𝑗 ≤𝑚 ( 𝑛 𝑖=1 |𝑀 𝑖,𝑗 |), thus we simply take the norm of the Jacobian (over all points in 𝑋 ) that has the maximum absolute column sum.
However, the following question arises: what if we need to compute the derivative at a point where it is not defined, such as ReLU(𝑥) at 𝑥 = 0? Can one extend Jacobians (and methods to compute them) to non-differentiable functions? We follow recent work by Sherman et al. [2021] and employ the notion of the Clarke Generalized Jacobian [Clarke 1990] for this extension. Intuitively, this generalizes the notion of a Jacobian to non-differentiable, but still Lipschitz continuous functions.
Convexity. The Clarke Jacobian of a function 𝑓 : R 𝑚 → R 𝑛 evaluates to a convex set of 𝑛 × 𝑚 matrices, hence we define the following operators. Let Co(R 𝑛×𝑚 ) denote all convex sets of 𝑛 × 𝑚 real matrices. Further, let co : P (R 𝑛×𝑚 ) → Co(R 𝑛×𝑚 ) be the convex hull operator, which given a set of matrices, takes their convex hull. We can now define the Clarke Generalized Jacobian. Definition 3.4. (Clarke Thm. 2.5.1 [Clarke 1990]) Let 𝑓 : R 𝑚 → R 𝑛 be a function that is locally Lipschitz at x 0 ∈ R 𝑚 , where the set of non-differentiable points of 𝑓 has Lebesgue measure 0 (we denote that set as 𝑆 ⊂ R 𝑚 ). The Clarke Generalized Jacobian of 𝑓 at the point x 0 is denoted with the following signature
, and is given as:
We first detail this definition when x 0 is a point of non-differentiability of 𝑓 . Intuitively, since we cannot compute the actual Jacobian J at a point of non-differentiability, we instead take any sequence of points that converges to that point of non-differentiability, such that the Jacobian is defined for all points in that sequence; we then compute what the Jacobians evaluated at those points in the sequence converge to. The Clarke Jacobian is then just the convex hull over all such sequences' respective limiting Jacobians. Thus, we obtain a convex set of matrices. When the function is differentiable at the point x 0 , this limit reduces to exactly the standard Jacobian at x 0 (and the convex hull becomes superfluous), hence 𝜕 𝑐 (𝑓 , x 0 ) will be a singleton set containing exactly J(𝑓 , x 0 ).
We describe our differentiable programming language, which is based upon 𝜆 𝑆 [Sherman et al. 2021] but with additional branching primitives. Our language is first-order and purely functional (no side-effects), yet expressive enough to encode neural networks and locally Lipschitz perturbations.
Figure 3 presents the syntax of the language. DeepJ syntactically supports standard arithmetic operations, differentiable function primitives, and limited branching for encoding Lipschitz but non-differentiable functions like min and max. One can encode neural networks in DeepJ; however, instead of encoding them as a series of edge-weight matrix multiplications, our syntax constructs them inductively via compositions: •, Cartesian products: ×, and branching (e.g., for ReLU networks).
Arithmetic Operations. We support all of the key arithmetic primitives: addition, multiplication, and division. While syntactically speaking, these are only defined for functions of a single output variable, one can easily encode multi-variable versions (e.g., vector addition) by also making use of the Cartesian product, ×, where
Differentiable Functions. Syntactically, we support as primitives all the standard primitive differentiable functions (e.g., 𝑒 𝑥 , sin(𝑥), tanh(𝑥), etc.). Hence, we denote these as 𝐶 1 because each function is 𝐶 1 -smooth, meaning the function is continuous and differentiable everywhere on its domain, and the first derivative is also continuous everywhere on its domain.
Non-differentiable, Lipschitz Functions. We support a branching primitive, 𝑓 0 > 0 ? 𝑓 1 : 𝑓 2 , to implement Lipschitz, but not necessarily differentiable functions such as abs, ReLU, min, and max (and by extension max-pooling). For example, abs can be implemented as 𝑥 > 0 ? 𝑥 : -𝑥. However, one could easily define a discontinuous function such as 𝑥 > 0 ? 1 : 0 using our syntax. Thus, for the computed Jacobian (and by extension Lipschitz constant) to be semantically meaningful and valid, we restrict the type of branching we support. We will later show that checking for these restrictions is undecidable in Section 4.2, and thus the responsibility of ensuring that the restrictions are satisfied rests upon the programmer or an (incomplete) program analysis. Lastly, even though we restrict the type of branches we support, more complex branches with arbitrary Boolean predicates can be systematically desugared into simpler ones. For example, the branch 𝑐 1 < 𝑥 ∧ 𝑥 < 𝑐 2 ? 𝑓 1 : 𝑓 2 can be desugared into 𝑐 1 < 𝑥 ? 𝑥 < 𝑐 2 ? 𝑓 1 : 𝑓 2 : 𝑓 2 . Disjunctions and negations of Booleans can be encoded similarly.
Type System. As our language only employs real-valued functions, the typing rules for a function are simple and based on standard real-valued arithmetic. For instance, when adding two functions, their dimensions must agree. Likewise, when dividing by a function 𝑓 1 (as in 1/𝑓 1 ), the output dimension of the function 𝑓 1 must be 1. The full typing rules can be seen in Fig. 3, where Γ corresponds to the typing context that maps all arithmetic expressions (including intermediate ones) to their types.
As our language is an augmented differentiable programming language, the semantic interpretation of a function 𝑓 is its derivative, which in our case corresponds to its Clarke Generalized Jacobian 𝜕 𝑐 (𝑓 , •). Our language is inspired by Sherman et al. [2021], hence we follow their convention of lifting the Clarke Jacobian to become a total function by defining the result to be ⊤ whenever 𝜕 𝑐 would be undefined, such as trying to evaluate the Clarke Jacobian of log(𝑥) at 𝑥 = 0. Because of this, ⊤ corresponds to the entire space of all real 𝑛 × 𝑚 matrices (R 𝑛×𝑚 ).
The semantic rules for recursively defining the Clarke Generalized Jacobian of each language primitive are shown in Figure 4. Unlike the regular Jacobian, these rules use ⊆ instead of equality. This means that the exact Clarke Jacobian is not computable; however, our end goal is to compute a sound over-approximation.
Convex Arithmetic Operations. We denote + 𝑐 : Co(R 𝑛×𝑚 ) × Co(R 𝑛×𝑚 ) → Co(R 𝑛×𝑚 ) to be the addition of two convex sets (Minkowksi addition), where 𝐴 + 𝑐 𝐵 = {𝑎 + 𝑏 : 𝑎 ∈ 𝐴, 𝑏 ∈ 𝐵}. Likewise, we denote
𝑣 • 𝑎 : 𝑎 ∈ 𝐴} to be convex scalar multiplication and division, respectively. Lastly, as mentioned, we let co : P (R 𝑛×𝑚 ) → Co(R 𝑛×𝑚 ) be the convex hull operator that takes the convex hull of a set of matrices. We now detail the rules of Fig. 4.
Variables and Constants. The Clarke Jacobian of a single variable 𝑥 𝑖 is {1}, but if we compute it with respect to x 0 ∈ R 𝑚 , the Clarke Jacobian will be 𝑚-dimensional, hence denoted as {e 𝑖 }, where the 𝑖 𝑡ℎ entry of e 𝑖 is 1 and all other 𝑚 -1 entries are 0. Equivalently, 𝑓 ≜ 𝑥 𝑖 can be thought of as a projection function going from R 𝑚 → R that takes the 𝑖 𝑡ℎ component of x 0 . Similarly, the Clarke Jacobian of any constant is just the vector of 𝑚 zeroes, denoted as {0}.
Cartesian Product. The Clarke Jacobian of the Cartesian product of two functions (defined over the same input) is a subset of the matrix concatenation of each function's respective Clarke Jacobians. This follows directly from proposition 2.6.2 (e) of Clarke [1990].
Addition. The Clarke Jacobian does not obey exact linearity; however, the Clarke Jacobian of the sum of two functions is contained in the Minkowski sum of each function's respective Clarke Jacobian. This rule follows directly from Proposition 2.3.3 of Clarke [1990].
Multiplication and Division. The Clarke Jacobian follows both a product and quotient rule, but as with the other rules, the relationship is of containment instead of strict equality. These follow directly from Propositions 2.3.13 and 2.3.14 of Clarke [1990].
Composition. To ensure our language is fully compositional, we can exploit the fact that the Clarke Jacobian follows a chain rule, when the outermost function is 𝐶 1 -smooth. This result follows directly from Theorems 2.3.9 and 2.6.6 of Clarke [1990].
𝐶 1 Functions. For the 𝐶 1 primitive functions, the Clarke Jacobian reduces to the standard Jacobian J(𝐶 1 , •). Furthermore, one can catch erroneous behavior by leveraging knowledge about the known domain of each 𝐶 1 primitive function. For instance, if one tries to evaluate 𝜕 𝑐 (sqrt(-𝑅𝑒𝐿𝑈 (𝑥 𝑖 )), 1), this would require evaluating J(sqrt, -1). However, since -1 lies outside the domain of sqrt, this result is undefined; thus, the expression will evaluate to ⊤, and errors will be caught at this step. One can then propagate ⊤ up the remainder of a function's expression tree, since all other rules in Fig. 4 first check if any sub-expression evaluates to ⊤ (in which case they will also evaluate to ⊤).
Branching. Branching is absent in Clarke's original formulation. Furthermore, 𝜆 𝑆 [Sherman et al. 2021] and Di Gianantonio and Edalat [2013] also do not support a branching primitive. This is because a branch can introduce the possibility of encoding a discontinuous function, such as 𝑓 (𝑥) = 𝑥 > 0 ? 1 : 0. As a branch can be thought of as "splitting" a function into two separate pieces, we need to ensure that the entire function is still locally Lipschitz continuous (on the region of interest containing x 0 ) for it to have a well-defined Clarke Jacobian at x 0 . Thus, to formally establish the necessary conditions for the well-definedness of the Clarke Jacobian of a branching function, we use results from the theory of piecewise differentiable functions [Khan and Barton 2013; Scholtes 2012]. We first state a useful lemma. Lemma 4.1. Let 𝑓 : R 𝑚 → R 𝑛 be a function expressible using the syntax of Figure 3 that does not contain branches. For any x 0 where J(𝑓 , x 0 ) is defined, we have that 𝜕 𝑐 (𝑓 , x 0 ) = {J(𝑓 , x 0 )}.
Proof. (Sketch) Since a function that does not have branches is a constant or only uses addition, composition with a 𝐶 1 function, multiplication, division, or the Cartesian product, one can directly compute J(𝑓 , x 0 ) using linearity, chain rule, product or quotient rule, or the direct computation of the derivative (for 𝐶 1 functions and constants), provided J(𝑓 , x 0 ) is well-defined (e.g., there is no division by 0). Furthermore, by Proposition 2.2.4 of Clarke [Clarke 1990], if a function has a well-defined Jacobian at a point, the Clarke Jacobian reduces to that value.
□
We now formally define piecewise differentiability.
Definition 4.2. (Piecewise Differentiability [Scholtes 2012]) A function 𝑓 : R 𝑚 → R 𝑛 is piecewise differentiable on an open set X ⊆ R 𝑚 if 𝑓 is continuous on X and for every 𝑥 ∈ X there exists an open neighborhood 𝑂 ⊆ X and a finite number of differentiable functions, denoted {𝑔 1 , ..., 𝑔 𝑘 }, such that for any 𝑥 𝑜 ∈ 𝑂, 𝑓 (𝑥 𝑜 ) ∈ {𝑔 1 (𝑥 𝑜 ), ..., 𝑔 𝑘 (𝑥 𝑜 )}. We will refer to the set of differentiable functions {𝑔 1 , ..., 𝑔 𝑘 } as the selection set.
Any function 𝑓 that is differentiable on a set X ⊆ R 𝑚 (meaning J(𝑓 , •) exists) is trivially piecewise differentiable on X, albeit with a single piece. A piecewise differentiable function has a well-defined Clarke Jacobian that can be given in terms of the convex hull of the standard Jacobian of constituent pieces, provided an active set is known a priori [Scholtes 2012]. If 𝑓 has the selection set of differentiable functions {𝑔 1 , ..., 𝑔 𝑘 }, then (by Proposition 4.3.1 of Scholtes [2012]):
where 𝐴(𝑓 , x 0 ) ⊆ {1, ..., 𝑘 } is the active set, which denotes at x 0 which of the 𝑘 selection functions satisfy 𝑔 𝑖 (x 0 ) = 𝑓 (x 0 ). To adapt this to our setting, we start with the simplest branching function 𝑓 ≜ 𝑓 0 > 𝑐 ? 𝑓 1 : 𝑓 2 where 𝑓 1 , 𝑓 2 are branch-free (meaning they are compositions, sums, or products of 𝐶 1 functions). Since 𝑓 1 and 𝑓 2 do not contain branches, one can compute J(𝑓 𝑖 , x 0 ) directly in accordance with Lemma 4.1. Hence, 𝑓 1 and 𝑓 2 are the selection set, as they are differentiable (by assumption) and the value of 𝑓 0 > 𝑐 ? 𝑓 1 : 𝑓 2 for some 𝑥 0 will necessarily be in {𝑓 1 (𝑥 0 ), 𝑓 2 (𝑥 0 )}. Therefore, when writing a branch, the selection set is known by construction. So all that remains for 𝑓 to be piecewise differentiable (and have a well-defined 𝜕 𝑐 ) is to ensure that it is continuous, which is true provided 𝑓 1 (𝑥) = 𝑓 2 (𝑥) for {𝑥 : 𝑓 0 (𝑥) = 𝑐}.
We do not need to restrict to cases where 𝑓 1 and 𝑓 2 are branch-free. We can nest branches arbitrarily deeply, provided they agree on the decision boundary {𝑥 : 𝑓 0 (𝑥) = 𝑐}, as all this does is increase the number of selection functions (or pieces) by finitely many. This is because we do not allow infinite recursion or while loops that would permit expressing countably infinite possible branches. Intuitively, we recursively "unroll" the nested branches into all possible innermost functions (which themselves will no longer contain branches). This ultimately yields a finite set of branch-free functions 𝑓 𝑖 for which we can compute J(𝑓 𝑖 , •) directly (as in Lemma 4.1). Using these notions, we can now formally describe necessary conditions for the Clarke Jacobian of the branching primitive to be well-defined. Theorem 4.3. (Well-Definedness of the Clarke Jacobian of a Branch) The function given by the branch 𝑓 0 > 𝑐 ? 𝑓 1 : 𝑓 2 is piecewise differentiable on an open set X ⊆ R 𝑚 if 𝑓 1 is piecewise differentiable on {𝑥 ∈ X : 𝑓 0 (𝑥) ≥ 𝑐}, 𝑓 2 is piecewise differentiable on {𝑥 ∈ X : 𝑓 0 (𝑥) ≤ 𝑐}, and 𝑓 1 (𝑥 0 ) = 𝑓 2 (𝑥 0 ) for all 𝑥 0 ∈ {𝑥 ∈ X : 𝑓 0 (𝑥) = 𝑐}.
Proof. (Sketch) Since 𝑓 1 is piecewise differentiable on {𝑥 ∈ X : 𝑓 0 (𝑥) ≥ 𝑐}, it has some selection set {𝑔 1 , ..., 𝑔 𝑘 }; hence, 𝑓 0 > 𝑐 ? 𝑓 1 : 𝑓 2 is also piecewise differentiable on {𝑥 ∈ X : 𝑓 0 (𝑥) ≥ 𝑐}, with the same selection set on that region. Likewise, since 𝑓 2 is piecewise differentiable on {𝑥 ∈ X : 𝑓 0 (𝑥) ≤ 𝑐}, it has some selection set {ℎ 1 , ..., ℎ 𝑙 }; thus, 𝑓 0 > 𝑐 ? 𝑓 1 : 𝑓 2 is also piecewise differentiable on {𝑥 ∈ X : 𝑓 0 (𝑥) ≤ 𝑐}, with the same selection set. Furthermore, since {𝑥 ∈ X :
□
We now present the rule for computing the Clarke Jacobian for a branching function, using the notions of piecewise differentiability. Formally, for 𝑓 0 : R 𝑚 → R, 𝑓 1 , 𝑓 2 : R 𝑚 → R 𝑛 , and x 0 ∈ R 𝑚 :
When 𝑓 0 (x 0 ) > 𝑐, the active set is 𝐴(𝑓 0 , x 0 ) = {𝑓 1 }; similarly, when 𝑓 0 (x 0 ) < 𝑐, the active set is 𝐴(𝑓 0 , x 0 ) = {𝑓 2 }. Along the decision boundary, both 𝑓 1 and 𝑓 2 are in the active set, hence the Clarke Jacobian is the convex hull, co, of both of their respective Clarke Jacobians. In the rule shown in Eq. 5, we evaluate 𝜕 𝑐 (𝑓 1 , x 0 ) and 𝜕 𝑐 (𝑓 2 , x 0 ) instead of J(𝑓 1 , x 0 ) and J(𝑓 2 , x 0 ) as in Eq. 4. This recursive definition allows us to capture the notion of unrolling a nested branch, as once 𝑓 1 and 𝑓 2 are themselves branch free, 𝜕 𝑐 (𝑓 𝑖 , x 0 ) and J(𝑓 𝑖 , x 0 ) become equivalent (by Lemma 4.1); thus, this equation would coincide with Eq. 4.
Example 4.4. We can encode ReLU(𝑥) ≜ 𝑥 > 0 ? 𝑥 : 0, and hence 𝜕 𝑐 (𝑥 > 0 ? 𝑥 : 0
Despite the elegant theory of Scholtes [2012], which allows us to characterize the conditions and well-definedness of the Clarke Jacobian of a piecewise differentiable function (or branch in our language), the problem of statically checking that these conditions are satisfied is undecidable. This is because even the smaller problem of ensuring that a branch is continuous along the decision boundary is known to be undecidable [Chaudhuri et al. 2010;Di Gianantonio and Edalat 2013;Griewank 2013]. As those works do not target the exact framework we focus on, we offer a short self-contained result: Lemma 4.5. (Undecidability of Ensuring Piecewise Differentiability) Let 𝑓 1 : R 𝑚 → R 𝑛 and 𝑓 2 : R 𝑚 → R 𝑛 be arbitrary functions in the language, and let X ⊆ R 𝑚 be an arbitrary set on which 𝑓 1 and 𝑓 2 are both defined. Checking if 𝑓 1 (𝑥) = 𝑓 2 (𝑥) for all 𝑥 ∈ X where 𝑓 0 (𝑥) = 𝑐 is undecidable.
Proof. (Sketch) Reduction using Richardson's Theorem [Richardson 1969]. □
Existing works can heuristically check for this condition using SMT solvers [Chaudhuri et al. 2010], or restrict the language to only pre-specified operators (e.g., max, min, ReLU and abs) [Beck and Fischer 1994;Di Gianantonio and Edalat 2013;Sherman et al. 2021]. Despite the difficulties introduced by a branching primitive, it will prove immensely useful in making the static analysis more precise, particularly for functions such as bilinear interpolation.
Having now defined the syntax and the standard Clarke Jacobian, we now formalize a computable, sound over-approximation: the Interval Clarke Jacobian, 𝜕 𝐼𝑛𝑡 . We will then show how to scalably implement our abstraction with an equivalent formulation that leverages a sound abstraction of forward-mode automatic differentiation.
To soundly approximate the Clarke Jacobian, we use the interval domain as it is fully computable and can soundly abstract the convex set of matrices corresponding to the original Clarke Generalized Jacobian, since interval matrices are necessarily convex sets.
Preliminaries. Denote the set of real-valued intervals of the form [𝑎, 𝑏] where 𝑎, 𝑏 ∈ R ∪ {±∞} and 𝑎 ≤ 𝑏 as IR. The set of 𝑛 × 𝑚 matrices of intervals is denoted as IR 𝑛×𝑚 . We will use the notation x 0 instead of x 0 to distinguish matrices and vectors whose entries are intervals instead of scalars. Similarly, to denote the evaluation of a function 𝑓 where all of its operations are lifted to interval arithmetic, we will write 𝑓 ( x 0 ). To denote the lower and upper bounds of an interval x 0 = [𝑎, 𝑏], we will write 𝑙𝑏 ( x 0 ) and 𝑢𝑏 ( x 0 ), respectively. We denote + IR , • IR , and / IR as the interval arithmetic versions of addition, multiplication, and division, respectively. Likewise, we denote ⊔ : IR×IR → IR to be the standard interval join, which returns the smallest interval enclosing both arguments. We may also apply ⊔ element-wise to matrices in IR 𝑛×𝑚 . Lastly, we denote ⊤ = IR 𝑛×𝑚 , which is the entire space of 𝑛 × 𝑚 interval matrices. We now define the Interval Clarke Jacobian, 𝜕 𝐼𝑛𝑡 .
Definition 5.1. The Interval Clarke Jacobian 𝜕 𝐼𝑛𝑡 : R 𝑚 → R 𝑛 × IR 𝑚 → IR 𝑛×𝑚 for a function 𝑓 : R 𝑚 → R 𝑛 and interval vector x ∈ IR 𝑚 is denoted 𝜕 𝐼𝑛𝑡 (𝑓 , x 0 ) and is given by the rules of Fig. 5.
The interpretation of the Jacobian of a function 𝑓 is now an interval matrix that over-approximates the convex set of matrices corresponding to the original Clarke Jacobian. This will be necessary to be able to define the local neighborhood of points for which we want to abstractly analyze or compute a Lipschitz constant. We now detail each abstract transformer.
Constants and Variables. The Interval Clarke Jacobian of a constant function is nearly identical to the respective Clarke Jacobian (in that 𝜕 𝐼𝑛𝑡 of a constant is 0); it is lifted to become the constant interval [0, 0] of the same dimension as x 0 , which we denote as 0 ∈ IR 𝑚 . Likewise, the Interval Clarke Jacobian of a single variable 𝑥 𝑖 is also the constant vector where the 𝑖 𝑡ℎ component is 1 and all other 𝑚 -1 components are 0, again lifted to become a constant interval, which we denote as e 𝑖 ∈ IR 𝑚 .
Cartesian Product. The Interval Clarke Jacobian of a Cartesian product of functions is just the concatenation of their respective Interval Clarke Jacobians.
Addition. The Interval Clarke Jacobian obeys linearity exactly.
Composition. For compositions with 𝐶 1 -smooth functions, the Interval Clarke Jacobian satisfies a chain rule as well, where J 𝐼𝑛𝑡 is the classical Jacobian of the 𝐶 1 -smooth function, just interpreted using interval arithmetic instead of ordinary arithmetic.
Example 5.2.
Just as with 𝜕 𝑐 , primitive 𝐶 1 functions are where we can check if the Jacobian is undefined (in which case 𝜕 𝐼𝑛𝑡 evaluates to ⊤), which we likewise propagate up the function's expression tree. Branching. If the lower bound of 𝑓 0 ( x 0 ) is larger than 𝑐, we definitively know that 𝑓 1 is the only possible function of the active set. Similarly, if the upper bound of 𝑓 0 ( x 0 ) is smaller than 𝑐, we definitively know that 𝑓 2 is the only possible function of the active set. If 𝑐 ∈ 𝑓 0 ( x 0 ), then it is possible that both 𝑓 1 and 𝑓 2 are in the active set, thus we have to consider both possibilities.
In either case, we can also refine our information about x 0 . For example, when evaluating 𝜕 𝐼𝑛𝑡 of the function 𝑓 (𝑥) ≜ 𝑥 > 0 ? 𝑓 1 : 𝑓 2 on the interval x 0 = [-1, 1], conditional upon entering the true branch, we should only evaluate 𝜕 𝐼𝑛𝑡 (𝑓 1 , [0, 1]). Likewise, conditional upon entering the false branch, we should only evaluate 𝜕 𝐼𝑛𝑡 (𝑓 2 , [-1, 0]). We denote the refinement of x 0 for a Boolean guard of the form 𝑓 0 > 𝑐 as 𝑓 0 > 𝑐 ( x 0 ). However, the general problem of determining how to optimally refine the input x 0 , conditioned on the information of the branch 𝑓 0 > 𝑐, is undecidable. Therefore, we follow the approach of Miné [Miné 2017] and refine x 0 when the function 𝑓 0 in the Boolean guard has a simple form: a single variable 𝑥 𝑖 , and then use the fallback transformer (the identity function) for all other cases (which can then be simplified to get the rule from Fig. 5):
We now state the soundness of 𝜕 𝐼𝑛𝑡 . We first define the concretization 𝛾.
Definition 5.3. Define the concretization 𝛾 : IR 𝑛×𝑚 → Co R 𝑛×𝑚 for an interval matrix 𝑆 ∈ IR 𝑛×𝑚 as follows:
This is because an interval matrix is already by definition a convex set of matrices. The soundness theorem is now given as follows.
Theorem 5.4. (Soundness) Let 𝑓 : R 𝑚 → R 𝑛 be any function constructible according to Fig. 3 andx
Proof. See Appendix A.2 (Laurel et al. [2022]). □
While we could compute 𝜕 𝐼𝑛𝑡 recursively in a reverse-mode pass, a key contribution of our work is to provide an equivalent, forward-mode version of the static analysis based on dual numbers. We can lift the ordinary arithmetic operators to dual intervals. We define dual interval addition, + ID : ID × ID → ID as follows:
It will also be useful to define a join for dual intervals,
We now describe how to perform a forward-mode abstract evaluation where all operations are lifted to operate on dual intervals using the abstract interpreter 𝐸𝑣𝑎𝑙 ID .
Definition 6.2. The abstract interpreter 𝐸𝑣𝑎𝑙 ID : R 𝑚 → R 𝑛 × ID 𝑚 → ID 𝑛 takes a real-valued function 𝑓 : R 𝑚 → R 𝑛 and evaluates it abstractly by lifting the interpretation of all operations to dual interval arithmetic, as shown in Fig. 6.
Because our language is functional, there is no "state. " However, the second argument to 𝐸𝑣𝑎𝑙 ID is the input, which serves the same purpose. We now detail the rules shown in Fig. 6.
Constants and Variables. The abstract evaluation of a constant is that constant, albeit abstracted to a (degenerate) dual interval. Likewise, the evaluation of a single variable 𝑥 𝑖 is the 𝑖 𝑡ℎ element of the input dual interval (denoted by the [𝑖] accessor).
Cartesian Product. To abstractly evaluate the Cartesian product of two functions 𝑓 1 and 𝑓 2 on a given input x 0 + y 0 𝜖, we abstractly evaluate each one, then take the Cartesian product of the results.
Addition. To abstractly evaluate the sum of two functions, we abstractly evaluate each on the given input x 0 + y 0 𝜖, then take the dual interval sum + ID of the respective results.
Multiplication and Division. As with addition, to abstractly evaluate the product of two functions or their quotient, we first abstractly evaluate the individual functions on the input, then use the dual interval forms of multiplication • ID and division / ID .
Composition. The abstract evaluation of the composition of any function 𝑓 1 with a 𝐶 1 function for an input dual interval x 0 + y 0 𝜖, is the composition of the successive abstract evaluations: 𝑓 1 is abstractly evaluated on the input, then the 𝐶 1 function is abstractly evaluated on that result.
Branching. When evaluating a branch abstractly, we only use the real part of the dual interval, denoted 𝑓 𝑠𝑡 (𝐸𝑣𝑎𝑙 ID (𝑓 0 , x 0 + y 0 𝜖)), to select which branch to take (or whether to abstractly evaluate both). While this may seem strange, one will note that in Eq. 5, the Clarke Jacobian is only computed for the piece that is ultimately chosen -𝑓 1 or 𝑓 2 (possibly both) -and not the threshold function 𝑓 0 . In fact, 𝑓 0 is only used to select which branch to take, hence its derivative information (which corresponds to its dual part) is unnecessary. If the lower bound 𝑙𝑏 (𝑓 𝑠𝑡 (𝐸𝑣𝑎𝑙 ID (𝑓 0 , x 0 + y 0 𝜖))) is larger than 𝑐, we definitively know to take only the true branch. Conversely, if the upper bound of the real part 𝑢𝑏 (𝑓 𝑠𝑡 (𝐸𝑣𝑎𝑙 ID (𝑓 0 , x 0 + y 0 𝜖))) is less than 𝑐, we only take the false branch. Otherwise, we abstractly evaluate both branches then take their join ⊔ ID . In all cases, we can refine the information of the real part x 0 but not the dual part y 0 , as there is no way of knowing which sub-regions of x 0 correspond to which sub-regions in y 0 , since the interval domain is non-relational.
We now prove that 𝐸𝑣𝑎𝑙 ID can be used to compute exactly the same result as 𝜕 𝐼𝑛𝑡 .
Theorem 6.3. (Equivalence of 𝜕 𝐼𝑛𝑡 and 𝐸𝑣𝑎𝑙 ID ) Let 𝑓 : R 𝑚 → R 𝑛 and let x 0 ∈ IR 𝑚 . Then
Proof. We show the sketch for select cases; the full proof is in Appendix A.4. Base Cases. Constants and single variables are straightforward.
𝑠𝑛𝑑 𝐸𝑣𝑎𝑙 ID (𝑐,
Arithmetic Operations. Equivalence for arithmetic primitives is straightforward, though it requires the inductive hypothesis and using the definitions of + ID , • ID , / ID . We now detail addition:
= 𝑠𝑛𝑑 𝐸𝑣𝑎𝑙 ID (𝑓 1 , x 0 + e 1 𝜖) 𝑇 × ... × 𝐸𝑣𝑎𝑙 ID (𝑓 1 , x 0 + e 𝑚 𝜖) 𝑇 (Ind. Hyp.)
The case for the Cartesian product proceeds analogously. Composition also uses the inductive hypothesis and relies upon the fact that
). Branching. As stated in Lemma 4.1, up to this point, without branching the language is just a standard differentiable programming language, thus the result is not surprising. However, with branching, establishing this equivalence is more challenging. We split the proof into three cases:
(1) 𝑙𝑏 (𝑓 ( x 0 )) > 𝑐, (2) 𝑢𝑏 (𝑓 ( x 0 )) < 𝑐, and (3) 𝑐 ∈ 𝑓 ( x 0 ). Starting with 𝑙𝑏 (𝑓 0 ( x 0 )) > 𝑐:
And for any 𝑖 ∈ {1, .., 𝑚}, 𝑙𝑏 (𝑓 0 ( x 0 )) > 𝑐 iff 𝑙𝑏 (𝑓 𝑠𝑡 (𝐸𝑣𝑎𝑙 ID (𝑓 0 , x 0 + e 𝑖 𝜖))) > 𝑐, hence:
By induction and the definitions of 𝐸𝑣𝑎𝑙 ID (𝑓 0 > 𝑐 ? 𝑓 1 : 𝑓 2 , x 0 + e 𝑖 𝜖) and 𝜕 𝐼𝑛𝑡 (𝑓 1 , 𝑓 0 > 𝑐 ( x 0 )):
The case where 𝑢𝑏 (𝑓 0 ( x 0 )) < 𝑐 proceeds exactly the same by symmetry. When we assume 𝑐 ∈ 𝑓 0 ( x 0 ), we must take the join. Starting from the desired right-hand side:
of 𝜕 𝐼𝑛𝑡 and Ind. Hyp.) □
The key insight of Theorem 6.3 is that, as with standard forward-mode AD, the dual interval abstraction requires as many forward passes as there are inputs. Therefore, for 𝑓 : R 𝑚 → R 𝑛 , we require 𝑚 independent passes of 𝐸𝑣𝑎𝑙 ID in order to compute the full interval over-approximation of the Clarke Jacobian. However, they can be computed in parallel (as we do). Hence, as with standard forward-mode AD, our method is better suited for functions where the input dimension is small. Additionally, the real part for each pass will be the same, and need only be computed once. Furthermore, for the basic arithmetic functions, dual interval arithmetic requires more primitive operations than the same operation defined over the reals. For example, multiplication of two dual intervals as shown in Section 6.1 requires 12 primitive multiplications to perform all the interval arithmetic multiplications for both the real and dual parts, as well as 2 primitive additions. However, for all arithmetic operations {+, -, •, /}, the amount of additional primitive operations required for the dual interval version (vs. the real-valued version) is still just a constant factor more.
While it has been established that for standard AD, the amount of primitive operations of a single pass of forward-mode evaluation with dual numbers is at most a constant amount (typically 2 -3x) more than the number of operations to evaluate the function itself [Griewank and Walther 2008], this is not the case for our analysis. A key complexity issue arises due to our support of branching. As the abstract evaluator 𝐸𝑣𝑎𝑙 ID could potentially evaluate both branches of a conditional, we may have exponentially many more evaluations compared to evaluating the function with ordinary real numbers. Thus, if 𝑂𝑃𝑆 (𝑓 ) denotes the number of primitive arithmetic operations to evaluate 𝑓 (not its Jacobian) over real numbers, then O (2 𝑂𝑃𝑆 (𝑓 ) ) is the upper bound on the number of primitive arithmetic operations required to evaluate a single pass of the dual interval lifted version of 𝑓 . For the evaluation, the ReLU function as well as the contrast variation and rotation perturbations employ branching. Yet, despite the theoretically worst-case exponential operation complexity, in practice, DeepJ's implementation of 𝐸𝑣𝑎𝑙 ID was still quite fast. Lastly, each term will always have only two intervals associated to it (the real and dual parts). Hence, unlike affine interval arithmetic [De Figueiredo and Stolfi 2004] or tools like Rosa [Darulova and Kuncak 2017], the number of intervals we must track per variable does not grow as a function of the expression size.
We next describe the experimental setup, including the perturbation functions and networks used in our experiments. We also detail the training procedures and accuracies of these networks. We ran our experiments on a 2.20 GHz 14 core Intel Xeon Gold 5120 CPU with 256 GB of main memory.
Implementation. We implemented our analysis in a tool called DeepJ, written in C++. DeepJ employs operator and function overloading to allow dual interval inputs to be propagated through arbitrary perturbation functions and neural networks. Even though our formalization is from a purely functional view, one can easily use an ANF conversion [Flanagan et al. 1993] to produce imperative code, since our programs do not have recursions or while loops. All code, neural networks, and results are available at https://github.com/uiuc-arc/DeepJ. Perturbation Functions. We consider three previously studied image perturbations affecting pixel intensity and image geometry, as well as compositions of these perturbations. We assume that the images' pixel values are in the range [0, 1].
(1) Haze [Paterson et al. 2021 4) Composition: We look at three ways of composing the above functions -haze followed by rotation, contrast followed by rotation, and contrast followed by haze.
The contrast, rotation, and composite perturbations described above are non-differentiable but Lipschitz continuous, and cannot be handled by prior work [Edalat and Maleki 2017;Jordan and Dimakis 2020;Mangal et al. 2020;Zhang et al. 2019]. In our experiments, we consider 𝛼 𝑚𝑎𝑥 ≤ 0.63 for the haze and contrast perturbations, and 𝜃 𝑚𝑎𝑥 ≤ 0.32 radians (which corresponds to approximately ±18 • ) for rotation.
Network Architectures. We trained three ReLU networks each for the CIFAR10 and MNIST datasets. The first network (FFNN) is a 7-layer fully-connected architecture from RecurJac [Zhang et al. 2019]. The other networks are convolutional networks from Mirman et al. [2018] -ConvMed features two convolutional and two fully-connected layers, while ConvBig features four convolutional and three fully-connected layers. Details on these networks are in Appendix A.5. Our largest network is the CIFAR ConvBig network containing > 62 K neurons. These network sizes are comparable to those of other state-of-the-art verification techniques [Singh et al. 2019]. Furthermore, for the local optimization landscape experiment, we trained 7 fully-connected networks on the MNIST dataset, varying the total number of layers from 3 to 9. Each hidden layer contains 30 neurons, and a ReLU activation is applied after every layer (including the final layer).
Data Transformation. For all MNIST networks, we transformed the training set so that for each image, we padded it by 4 and took a random 28x28 crop of the resulting 32x32 image. For the CIFAR10 networks, we introduced random cropping with padding 4 and randomly flipping each image horizontally with probability 0.5. Afterwards, we normalized the images using 𝜇 = 0.1307, 𝜎 = 0.3081 for MNIST and 𝜇 = (0.4914, 0.4822, 0.4465), 𝜎 = (0.2023, 0.1994, 0.2010) for CIFAR10.
Training Hyperparameters. For the 7 networks trained for the local optimization experiment, we used the Adam optimizer [Kingma and Ba 2014] with a learning rate of 10 -4 and L2-regularization with 𝜆 = 0.001; we trained with a batch size of 500 for 60 epochs, using 6,000 images from the training set for validation. For all other networks, we used the Adam optimizer with a learning rate of 10 -3 on MNIST and 10 -4 on CIFAR10. We trained the MNIST networks with a batch size of 500 for 30 epochs, using 6,000 images from the training set for validation; we trained the CIFAR10 networks with a batch size of 64 for 60 epochs, using 5,000 images from the training set for validation. In all cases, we then picked the model that attained the highest validation accuracy out of all epochs.
Network Accuracies. The fully-connected networks trained for the local optimization landscape experiment all attain test accuracy of at least 92%. The classification accuracies for the networks in the Lipschitz experiments are shown in Table 1.
We evaluate the effectiveness of our approach on two tasks: (i) Lipschitz robustness certification and (ii) optimization landscape analysis. Both tasks are defined across a variety of datasets, perturbation functions, and neural network architectures, demonstrating our language's flexibility and scalability.
We compute an upper bound on the Lipschitz constant of functions of the form 𝑓 • 𝑛 • 𝑝, where 𝑝 denotes one or more perturbation functions, 𝑛 is an input normalization function (which simply rescales images by a constant), and 𝑓 is a neural network. The correctness of bounding the Lipschitz constant via the Interval Clarke Jacobian is proved in Appendix A.3. We use the ℓ ∞ -norm for the calculation of all Lipschitz constants. For both individual and composite perturbations, we consider five versions of our tool: DeepJ with the given range for the input perturbation parameter(s) and four versions of DeepJ where the input intervals are subdivided uniformly, denoted DeepJ 𝑘x, where 𝑘 represents how many subintervals are used per input parameter. When using splitting, we compute the upper bound on the Lipschitz constant separately for each split and take the maximum constant across the splits. Each pass of our analysis per image and split is completely independent of one another, and is thus parallelized in our implementation. Since no existing work can handle functions of the form 𝑓 • 𝑛 • 𝑝 considered in our evaluation, we employ a baseline combining global and local Lipschitz analysis. The baseline computes 𝐿 𝑙𝑜𝑐 (𝑝) • 𝐿 𝑔𝑙𝑜𝑏𝑎𝑙 (𝑛) • 𝐿 𝑔𝑙𝑜𝑏𝑎𝑙 (𝑓 ), where 𝐿 𝑙𝑜𝑐 and 𝐿 𝑔𝑙𝑜𝑏𝑎𝑙 are the local and global Lipschitz constants of their corresponding functions. As 𝑝 is not globally Lipschitz, we obtain 𝐿 𝑙𝑜𝑐 (𝑝) with DeepJ. We compute 𝐿 𝑔𝑙𝑜𝑏𝑎𝑙 (𝑛) as the reciprocal of the standard deviation (MNIST) or the reciprocal of the smallest standard deviation of the three RGB channels (CIFAR10) used for normalizing the input. We calculate 𝐿 𝑔𝑙𝑜𝑏𝑎𝑙 (𝑓 ) by multiplying the norm of each layer's weights using the tool from Gouk et al. [2021].
Results for Individual Perturbations. input to our function is in radians and covers the same range of interval sizes as the other two perturbations, but we convert the units to degrees for clearer presentation. In all cases, DeepJ is more precise than the baseline, often obtaining Lipschitz bounds orders of magnitude smaller than the baseline. The computed bounds become larger as the perturbation size is increased. It can also be seen that increased splitting leads to more precise results, with DeepJ 25x achieving much lower bounds than vanilla DeepJ (which does not do splitting).
Table 2 shows statistics on the runtime of the different methods for the same networks as in Figure 7. We consider all values of the parameters on the x-axis shown in Fig. 7. The runtimes for the remaining networks can be found in Appendix A.5. We report the minimum, median, and maximum runtimes across the 100 images for each function. DeepJ 25x takes the longest to run due to more splits, while the baseline usually runs the fastest, except for rotation on large angles. In some cases, splitting may be faster than the baseline or vanilla DeepJ, since without splitting, intervals can become too over-approximate; evaluating conditionals on over-approximate intervals often requires evaluating both branches, which may lead to the exponential blowup phenomenon discussed in Section 6.4. For rotation with bilinear interpolation, we empirically observed this beyond ±0.1 radians. The computation of 𝐿 𝑙𝑜𝑐 (𝑝) via our method contributes the most to the runtime of the baseline. For all versions of DeepJ, the rotation perturbation has the highest runtime. For a given perturbation type, the runtime increases with the size of the network. On the most expensive rotation perturbation with the CIFAR ConvBig network, the median time for DeepJ to finish is under 1.4 min per 100 images. Finally, the precision of our analysis can be improved by increasing the number of splits (as seen in Figure 7) at the cost of additional runtime (as seen in Table 2).
Results for Composite Perturbations. Figure 8 shows the upper bounds on the Lipschitz constant computed for compositions of perturbations on the same networks as in Figure 7. The results for the remaining networks can be found in Appendix A.5. Haze-Rotation denotes a haze perturbation composed with a rotation, in that order; the terminology is similar for the other composite perturbations. We use the same interval width when perturbing each parameter independently; the x-axis shows this width. For compositions that involve rotation, if the interval size on the x-axis is denoted 𝑠, we utilize the real interval [-𝑠/2, 𝑠/2] for rotation and the interval [0, 𝑠] for the other perturbation. The y-axis shows the upper bound on the Lipschitz constant computed by each method. Again, both axes use logarithmic scales. Each data point is the average over 10 correctly classified images, taking the first correctly classified test-set image from each output category. We use 𝑠 ∈ {10 -𝑘/4 • 2 | 𝑘 ∈ {4, 7, 10, 13, 16, 19}} for all experiments.
As with individual perturbations, DeepJ is better than the baseline in all cases, obtaining upper bounds on the Lipschitz constants orders of magnitude smaller than the baseline, as seen in Figure 8. DeepJ 9x with the largest number of splits is the most precise. The Lipschitz constants for both the Haze-Rotation and Contrast-Rotation perturbations are nearly identical, since the entries in the Clarke Jacobian corresponding to the rotation variable have much higher magnitudes. Thus, these entries dominate the ℓ ∞ -norm, overshadowing the effect of the other perturbation. Table 3 shows the runtime statistics for the different methods in the same way as Table 2. The results for the remaining networks can be found in Appendix A.5. We observe similar trends in the relative runtimes of the different methods as with individual perturbations. For all versions of DeepJ and for the baseline, the compositions that take the longest time to analyze are those containing rotations. On the most expensive composite perturbation (Haze-Rotation) with the CIFAR ConvBig network, the median time for DeepJ to finish is under 51 seconds per 10 images. As with individual perturbations, the precision of DeepJ can be improved by considering more splits at the cost of additional runtime.
Our tool, though implemented via floating-point, assumes real arithmetic. To ensure this assumption does not lead to substantially different results, we also implemented a floating-point sound version using the techniques in [Miné 2004]. We ran floating-point sound experiments for vanilla DeepJ. Table 4 shows that the difference in the computed Lipschitz constants is negligible (< 10 -10 ) in all cases. However, ensuring floating-point soundness adds up to 4.1x runtime overhead.
Obtaining an Interval Clarke Jacobian allows us to study the local geometry of 𝑓 • 𝑛 • 𝑝. We focus on finding the largest input regions where no stationary point exists. To prove that a given region does not contain a stationary point, we check if there exists an interval entry [𝑙 𝑖,𝑗 , 𝑢 𝑖,𝑗 ] in the Interval Clarke Jacobian such that 𝑙 𝑖,𝑗 > 0 or 𝑢 𝑖,𝑗 < 0. If this holds, then we can guarantee that the Clarke Jacobian matrix will not become zero, and therefore no stationary point exists within the given input region (which follows from Theorem 2.3.2 of [Clarke 1990]). Hence, similar to RecurJac [Zhang et al. 2019], we study the relationship between network depth and the maximal interval size for which the absence of stationary points can still be guaranteed. However, RecurJac cannot handle functions of the form 𝑓 • 𝑛 • 𝑝 considered in our work. We define the set of input intervals to analyze as {[0, 10 -𝑘/4 • 2] | 𝑘 ∈ [2, 21]} for haze and contrast and {[-10 -𝑘/4 , 10 -𝑘/4 ] | 𝑘 ∈ [2, 21]} for rotation. We consider 100 test images from MNIST that are correctly classified, using the first 10 correctly classified test images from each category. For each combination of perturbation type and network, we compute the largest input region where a stationary point does not exist. For each image, we uniformly split every input interval into 25 subintervals, computing the Interval Clarke Jacobians separately for each subinterval. Next, we identify the subintervals where the entry in the Interval Clarke Jacobian corresponding to the correctly classified class is either strictly positive or negative. We sum the widths of all such subintervals, then compute the maximum of these sums across all candidate widths. Figure 9 shows that as the depth of the fully-connected network increases, the maximal interval size that guarantees a non-zero Jacobian decreases for all perturbations. Further, we use the same procedure to obtain maximal interval sizes for the networks used in Section 8.1, as shown in Table 5. We observe that rotation requires noticeably smaller interval sizes as it is a more complex perturbation involving interpolation, whereas the other two are simpler pixel-wise operations. Furthermore, we can certify the absence of stationary points for larger regions on MNIST compared to CIFAR10, due to the former's smaller input dimension.
While there is considerable work on the static analysis of input-output properties of ML models specified via constraints on the network inputs and outputs [Balunović et al. 2019;Ehlers 2017;Huang et al. 2017;Katz et al. 2017Katz et al. , 2019;;Singh et al. 2018Singh et al. , 2019;;Sotoudeh and Thakur 2020;Urban and Miné 2021], such as for robustness and safety, much less work has been done on formally analyzing the Jacobian matrix.
From the programming languages (PL) and automatic differentiation (AD) literature, prior works [Di Gianantonio and Edalat 2013;Edalat et al. 2013;Khan and Barton 2013;Sherman et al. 2021] have examined AD through the Clarke Jacobian [Clarke 1990]. However, many define their semantics for computing the Clarke Jacobian at a single point [Khan and Barton 2012Barton , 2013]], instead of an abstraction of points. Additionally, the works that can formally analyze Jacobians for sets of input points are insufficient for our tasks. Di Gianantonio and Edalat [2013] restrict functions' domains to [-1,1] and require all Lipschitz constants be less than 1, thus their analysis cannot be used for local Lipschitz certification. The work by Edalat and Lieutier [2004] is restricted to functions of a single variable with input domain on [0, 1]. Follow-up works [Edalat et al. 2013;Edalat andMaleki 2017, 2018] all suffer other restrictions, namely requiring the output dimension be one and only supporting limited arithmetic operations (e.g., no division), which render them unable to analyze both neural networks and our perturbations. Furthermore, these techniques provide only a theoretical discussion, with no implementation of their approaches. 𝜆 𝑆 [Sherman et al. 2021] presents semantics for the Clarke Jacobian for concrete input points, extended partially to intervals (described in their appendix) as compared to our approach, which presents a sound abstract interpretation of the Clarke Jacobian that we compute with a set of fully compositional dual interval domain transformers, whose behavior is specified exactly for all language primitives. Another closely related PL work by Chaudhuri et al. [2011] aims to establish Lipschitz robustness of programs, but they use the classical Jacobian and do not support AD. Follow-up works showed how to analyze Lipschitz robustness of non-differentiable programs by analyzing their smooth approximation [Chaudhuri and Solar-Lezama 2011], stemming from the fact that several PL works focus on differentiable approximations of non-differentiable programs and ML models [Chaudhuri and Solar-Lezama 2010;Laurel and Misailovic 2020]. In contrast, DeepJ can analyze Lipschitz properties of programs with points of non-differentiability directly without using approximations, by employing the more general Clarke Jacobian.
In addition, while the practical problem of Lipschitz certification of neural networks has been studied [Jordan and Dimakis 2020;Scaman and Virmaux 2018;Weng et al. 2018a,b;Zhang et al. 2019], to the best of our knowledge, none of these works can bound local Lipschitz constants for composite, non-smooth perturbations. Both Mangal et al. [2020] and Jordan and Dimakis [2020] use proof techniques that require the classifier function analyzed be piecewise linear, hence they cannot support arbitrary activations. Additionally, RecurJac [Zhang et al. 2019] cannot support arbitrary arithmetic primitives like non-scalar multiplication or division. Hence, none of these works can reason about non-differentiable input perturbations, such as those generated by rotation [Balunović et al. 2019] or contrast variation [Paterson et al. 2021].
Our work also bears similarity with Rosa [Darulova and Kuncak 2017], as they also bound Jacobians to compute Lipschitz constants. However, unlike us, they do not support bounding Clarke Jacobians. Further, while their abstract domain tracks "𝜖-terms," the semantics of these terms do not correspond to first derivatives, but rather numerical round-off error effects.
We developed a novel abstraction for bounding the Clarke Jacobian of a Lipschitz, but not necessarily differentiable function for local input regions. Our domain, based upon dual numbers, soundly over-approximates all first derivatives needed to compute the Clarke Jacobian. We implemented our analysis in tool named DeepJ and showed that it can efficiently compute Lipschitz bounds and analyze the local geometry of multiple deep neural networks with respect to multiple nondifferentiable input perturbations. Our work is the first to address the problem of local Lipschitz certification of non-smooth perturbations, such as haze, contrast variation, rotation with bilinear interpolation, and their compositions.
Proc. ACM Program. Lang., Vol. 6, No. POPL, Article 56. Publication date: January 2022.
𝐶 1 ::= sin(𝑥) | cos(𝑥) | 𝑒 𝑥 | log(𝑥) | sqrt(𝑥) | 𝜎 (𝑥) | tanh(𝑥) Γ ⊢ 𝑐 : R 𝑚 → R Γ ⊢ 𝑥 𝑖 : R 𝑚 → R Γ ⊢ 𝑓 1 : R 𝑚 → R 𝑛 Γ ⊢ 𝑓 2 : R 𝑚 → R 𝑝 Γ ⊢ 𝑓 1 × 𝑓 2 : R 𝑚 → R 𝑛+𝑝 Γ ⊢ 𝑓 1 : R 𝑚 → R Γ ⊢ 𝑓 2 : R 𝑚 → R Γ ⊢ 𝑓 1 + 𝑓 2 : R 𝑚 → R Γ ⊢ 𝑓 1 : R 𝑚 → R Γ ⊢ 𝑓 2 : R 𝑚 → R Γ ⊢ 𝑓 1 • 𝑓 2 : R 𝑚 → R Γ ⊢ 𝑓 1 : R 𝑚 → R Γ ⊢ 1/𝑓 1 : R 𝑚 → R Γ ⊢ 𝑓 1 : R 𝑚 → R Γ ⊢ 𝐶 1 • 𝑓 1 : R 𝑚 → R Γ ⊢ 𝑓 0 : R 𝑚 → R Γ ⊢ 𝑓 1 : R 𝑚 → R 𝑛 Γ ⊢ 𝑓 2 : R 𝑚 → R 𝑛 Γ ⊢ 𝑓 0 > 𝑐 ? 𝑓 1 : 𝑓 2 : R 𝑚 → R 𝑛