Recent years have seen a resurgence in popularity of Automatic Differentiation (AD) and Differentiable Programming ś not only from machine learning applications requiring derivatives (e.g., for training neural networks [Abadi et al. 2016]), but also in areas as diverse as Computer Graphics [Li et al. 2018], Physical Simulation [Hu et al. 2020] and Scientific Computing [Walther and Griewank 2012]. However, despite the widespread adoption of AD, as pointed out in Hückelheim et al. [2018], there is surprisingly little work on formally verifying that these programs compute correctly. While there have been several works [Di Gianantonio and Edalat 2013;Krawiec et al. 2022;Sherman et al. 2021] on proving correct the concrete semantics of AD, there is almost no work on defining different abstract AD semantics for the purpose of verifying program properties expressed over derivatives.
Abstractly interpreting AD programs is further complicated by the fact that these programs do not have the typical concrete semantics of standard numerical programs. One either overloads all operators to have a non-standard interpretation over a tuple of numbers corresponding to the derivatives, such as with dual numbers [Griewank and Walther 2008], or one defines the semantics of the differentiable program using reverse-mode program execution. Additionally, most semantics for computing higher-order derivatives are defined over complicated structures, such as derivative towers [Karczmarczuk 2001] or lazily evaluated data structures as in Pearlmutter and Siskind [2007]. Hence, given the complexities in existing AD program semantics, adapting existing static analysis techniques, specifically abstract interpretation [Cousot and Cousot 1977], to soundly reason about differentiable programs is a notable challenge.
While there has been limited work for soundly reasoning about AD, the few techniques that have been developed in this area are significantly restricted: they are either limited to the interval domain [Bendtsen and Stauning 1996;Deussen 2021;Laurel et al. 2022a;Sherman et al. 2021] or do not support higher-order derivatives [Jordan and Dimakis 2021;Laurel et al. 2022a]. Thus, precisely tracking relational information across higher derivatives is out of the realm of existing abstract interpreters, despite the fact that this precision is needed in tasks such as improving the performance of verified ODE solvers [Immler 2018]. Furthermore, fundamental questions that arise when constructing an abstract semantics, such as how to systematically construct sound transformers for broad classes of functions or understanding the differences in the precision of different domains, have not been studied ś especially for AD involving higher-order derivatives.
Challenges. We focus on developing a general framework for systematically constructing abstract semantics for AD with arbitrary order derivatives. The first key challenge is that we need to first define a concrete semantics of higher-order AD that lends itself to precise and scalable abstract interpretation ś particularly with higher derivatives. This is challenging since the same abstract interpreter can have significantly different cost and precision for different syntactic representations of the same computation.
The second key challenge is ensuring that this abstract reasoning can then be done generally with expressive relational abstract domains and is not just constrained to the simpler interval domain, as in Bendtsen and Stauning [1996]; Deussen [2021]; Laurel et al. [2022a]. This problem is challenging since unlike in conventional programs, reasoning about higher-order AD requires precise abstract transformers for an exponential number of non-linear assignment statements which result from all the different partial derivative expressions. Designing these abstract transformers manually would require considerable expertise as well as substantial AD domain knowledge to know how the functional forms of different derivatives can be leveraged for precision.
Therefore, we need a generic construction that can systematically construct sound transformers for all higher-order derivatives by composing the transformers for a small number of primitive functions in our language. We also want to systematically leverage analytical properties of these derivatives to further refine the precision of the abstraction in a way that goes beyond naively composing existing numerical domains' abstract transformers together. Simultaneously addressing all of these challenges is beyond the scope of any existing work.
Our Work. We present a novel construction for abstractly interpreting higher-order AD to address these challenges. Our construction involves defining a concrete semantics for arbitrary order AD that semantically desugars the computation of derivatives into primitive operations on single variables, and that also leverages data-dependence across derivatives in the program's syntactic representation. Thus, given only a small set of sound transformers for the elementary arithmetic primitives, we can immediately abstractly interpret these transformed AD programs, and still attain precision due to the dependency-aware syntactic representation of the AD program. Hence, we reduce the problem of formally reasoning about complicated AD semantics and having to manually design custom transformers for all possible partial derivatives to reasoning about (straight-line) programs using standard abstract domains. This technique allows us generality in the abstraction ś we can easily instantiate the construction with relational domains. We also show how to combine the abstract transformers with domain-specific knowledge about derivatives to further improve the precision of the analysis.
Contributions. The paper makes the following contributions:
(1) Concrete AD Semantics: We present a novel concrete semantics for higher-order AD. By designing our concrete semantics to be forward-mode and imperative, and by exposing each variable in a way that captures data dependence across derivatives, we enable one to intuitively define precise abstractions of these concrete semantics.
(2) Generic Abstract Interpretation of AD: We provide the first generic construction to allow one to abstractly interpret forward-mode AD with arbitrary order derivatives. This construction is fully general and can be instantiated with both relational and non-relational abstract domains. (3) Implementation: We implement our formal construction into a practical tool and instantiate it with the interval and zonotope domains for first and higher derivatives, thus providing the first implementation of higher-order AD with relational abstract domains. Our implementation is publicly available at https://github.com/uiuc-arc/AbstractAD [Laurel et al. 2022c]. (4) Evaluation: We empirically show the advantages of using a relational abstract domain and the benefits of abstractly interpreting higher-order derivatives through two case studies:
(1) computing robust interpretations of regression neural networks and (2) bounding the Lipschitz constant of classification neural networks with respect to a semantic perturbation.
For the robust interpretation task, we are the first to bound both first and second derivatives of a neural network with respect to an input region. First and second derivatives computed via zonotope AD are up to 4.76× and 6.98× more precise, respectively, compared to interval AD. For Lipschitz certification, we obtain bounds that are up to 11,850× more precise with zonotopes, compared to the state-of-the-art interval-based tool of Laurel et al. [2022a].
We start with a simple illustrative example that highlights the generality of our construction.
Running Example. Fig. 1 presents an example neural network (the reader can find a simpler example in the Appendix [Laurel et al. 2022b]). The network takes two inputs 𝑥 1 and 𝑥 2 , propagates them through an affine layer to get the values of hidden neurons 𝑥 3 and 𝑥 4 , applies a SoftPlus activation to get the values of 𝑥 5 and 𝑥 6 , then passes those results through a final affine layer to get the outputs 𝑥 7 and 𝑥 8 . Our goal is to compute second derivatives of the network's outputs (𝑥 7 and 𝑥 8 ) with respect to the inputs (𝑥 1 and 𝑥 2 ). These derivatives are useful for a variety of applications, such as if we want to explain the network by computing interactions between inputs (as in Janizek et al. [2021]) or if we want to certify input regions where the function is locally concave or convex [Deussen 2021]. For this example, we study the impact of the interaction between 𝑥 1 and 𝑥 2 on the
𝑤 57 = 0.7074 𝑤 5 8 = -0 .7 0 6 8
Fig. 1. Concrete 2 𝑛𝑑 order forward-mode AD evaluation network output over a local region. The computation of the neural network output can be described by a straight-line imperative program, where each of the eight neurons is a program variable. However, when interpreting this program with forward-mode AD semantics, our interpreter will add additional augmented variables to keep track of all the derivatives (described in Section 4).
Second Derivative Analysis. While our prior work [Laurel et al. 2022a] can propagate 1 𝑠𝑡 derivatives through a neural network using the interval domain, for computing higher-order interactions of multiple variables, we need to compute higher derivatives with respect to the inputs, in this case 2 𝑛𝑑 derivatives. To compute 2 𝑛𝑑 derivatives, we instantiate a 2 𝑛𝑑 -order instance of our construction.
A key contribution of our approach is that it provides a general construction for giving both a concrete and abstract semantics for AD of arbitrary order. For this example, each variable 𝑥 𝑖 in the original program will have four associated components in the forward-mode AD interpretation. The first of the four components is the łreal partž, denoted as 𝑥 𝑖 [∅], which intuitively corresponds to a ł0 𝑡ℎ ž derivative. Equivalently, 𝑥 𝑖 [∅] is just the value of 𝑥 𝑖 , if the program were run under standard semantics instead of forward-mode AD semantics.
The next two components are 𝑥 𝑖 's first derivative terms, 𝑥 𝑖 [{1}] and 𝑥 𝑖 [{2}], which respectively index the first derivatives with respect to variable 1 and variable 2. Just as dual numbers (the canonical approach for first-order AD) only support differentiation of one variable at a time, our second-order analysis only supports differentiation of two variables at a time, which we explain formally in Theorem 4.11. In this example, variable 1 is just 𝑥 1 and variable 2 is just 𝑥 2 , however this need not always be the case. Hence, for each 𝑥 𝑖 , 𝑥 𝑖 [{1}] = 𝜕𝑥 𝑖 𝜕𝑥 1 and 𝑥 𝑖 [{2}] = 𝜕𝑥 𝑖 𝜕𝑥 2 . The last of the four components is the 2 𝑛𝑑 derivative term 𝑥 𝑖 [{1, 2}], which tracks the derivative with respect to the first variable, then with respect to the second; equivalently, 𝑥 𝑖 [{1, 2}] = 𝜕 2 𝑥 𝑖 𝜕𝑥 1 𝜕𝑥 2 . Furthermore, because derivatives are symmetric (e.g., 𝜕 2 𝑥 𝑖 𝜕𝑥 1 𝜕𝑥 2 = 𝜕 2 𝑥 𝑖 𝜕𝑥 2 𝜕𝑥 1 ), we do not need to worry about terms such as 𝑥 𝑖 [{2, 1}]. In our example, we can think of the forward-mode AD semantics as propagating coefficients of a 2 𝑛𝑑 degree Taylor polynomial that has been truncated to only include these terms. A key point is that our forward-mode AD semantics will automatically transform the original program to include additional variables to track these additional derivatives, which we will also call augmented variables. These derivative will be evaluated not symbolically, but rather at specific points ś in this example, all derivatives will be concretely evaluated at the scalar point (𝑥 𝑖𝑛 1 , 𝑥 𝑖𝑛 2 ) ∈ Rfoot_1 shown in Fig. 1. Lastly, to properly initialize the input state so that the concrete evaluation produces correct derivatives with respect to 𝑥 1 and 𝑥 2 , we initialize the first derivative terms of the inputs as
This initialization is the secondorder analog of how one initializes states for 1 𝑠𝑡 -order forward-mode AD with dual numbers. Forward-Mode AD. The question then becomes how to correctly propagate these additional terms corresponding to first and second derivatives of each 𝑥 𝑖 through the differentiable program (in this example a neural network). For first derivatives, there are the canonical dual numbers which give a standard way to do this. However, for higher derivatives, we must lift all standard calculus rules to higher-order derivatives. For instance, the chain rule must be generalized using the Faa di Bruno formula [Di Bruno 1857], and the product rule must be generalized using the general Leibniz rule.
We next compute the values for 𝑥 3 and 𝑥 4 . Because of the linearity of derivatives, both the first and second derivative terms for 𝑥 3 and 𝑥 4 are also linear combinations of the respective first and second derivative terms of 𝑥 1 and 𝑥 2 . Hence 𝑥 3
Upon computing all real and derivative coefficients for 𝑥 3 and 𝑥 4 , we then apply a SoftPlus activation to obtain 𝑥 5 and 𝑥 6 , respectively. Computing the first derivative terms ś 𝑥 5 [{1}], 𝑥 5 [{2}], 𝑥 6 [{1}], and 𝑥 6 [{2}] ś is done virtually identically as with dual numbers. We compute the first derivative of the SoftPlus function which is the sigmoid function, 𝜎 (𝑥), and evaluate it at both 𝑥 3 [∅] and 𝑥 4 [∅], saving these results in unique intermediate variables, 𝑑𝑥 5 [1] and 𝑑𝑥 6 [1].
However, for first derivatives, the chain rule tells us that we need not only the derivatives of the sigmoidal function (𝑑𝑥 5 [1] and 𝑑𝑥 6 [1]) we are composing with, but also the first derivative of the inputs (𝑥 3 [{1}] and 𝑥 4 [{1}]), hence why we must also evaluate 𝑥 5
Computing the higher-derivative terms for composition with the SoftPlus function is more involved. While computing the higher-derivative terms for affine combinations of variables is easy due to linearity, higher-derivatives of functional composition require Faa di Bruno's formula. In this example, we need to compute the value of the second derivative of the SoftPlus function, which is 𝜎 (𝑥) • (1 -𝜎 (𝑥)). Since we have already computed 𝜎 (𝑥 3 [∅]) and 𝜎 (𝑥 4 [∅]), we can reuse these results (respectively stored in 𝑑𝑥 5 [1] and 𝑑𝑥 6 [1]), hence why the second derivatives are computed as
It is important to note that there is a data-dependence across derivatives: for instance, 𝑑𝑥 6 [2] has data-dependence on 𝑑𝑥 6 [1]. This data-dependence is a direct result of how we have set up our concrete AD semantics, as one could have naively had the transformed program recompute terms instead of producing a program that exposes data dependence between lower and higher-order derivatives as we do. We will later see in Section 7 that transforming the program using this technique improves the precision when performing abstract interpretation with a relational domain.
Lastly, we again apply an affine layer to get the final output of the network. As before, we take linear combinations of not just the łrealž parts (𝑥 5 [∅] and 𝑥 6 [∅]) but also their derivative parts, since derivatives still follow linearity. The final outputs are now the derivatives with respect to the original inputs, 𝑥 1 and 𝑥 2 . Thus, 𝑥 7
= 0 .7 0 6 8 𝑤 24 = 0.7074
, where the notation 𝜕𝑥 𝑖 𝜕𝑥 𝑗 ,...,𝜕𝑥 𝑘 | (𝑥 1 ,...,𝑥 𝑚 ) ∈R 𝑚 denotes the partial derivative of some variable (𝑥 𝑖 ) with respect to other variables (𝑥 𝑗 , ..., 𝑥 𝑘 ) evaluated at some concrete point (𝑥 1 , ..., 𝑥 𝑚 ) ∈ R 𝑚 . Likewise, we also know that 𝑥 8
2 ) ∈R 2 . Abstract Interpretation. Having shown how the concrete semantics are interpreted and expose each derivative explicitly as a separate variable in memory (indexed by a set), we can think about how we might abstractly interpret these semantics to get formal bounds on all of these derivatives inside a region. Getting bounds on these derivatives is useful as (1) derivatives are canonically used for explaining the behavior (e.g. 1 𝑠𝑡 -order feature attribution or 2 𝑛𝑑 -order interactions) of a neural network and (2) prior work has shown that such explanations evaluated at scalar points are not robust [Alvarez-Melis and Jaakkola 2018], hence why recent work [Fel et al. 2022] has focused on computing provably robust explanations using interval abstractions. Abstractly interpreting the AD used to compute the derivatives for these explanations gives us these provable bounds.
For our example, we instantiate our construction with the zonotope abstract domain [Ghorbal et al. 2009], which associates to each variable an affine form written as 𝑐 0 + 𝐾 𝑖=1 𝑐 𝑖 𝜖 𝑖 , where each noise term satisfies 𝜖 𝑖 ∈ [-1, 1]. Since multiple variables can share noise terms, dependencies across variables are preserved to a high-degree; therefore, the zonotope domain is a relational domain.
For the purposes of this example, we want to compute provable bounds on the derivatives of all variables for the local region
, its affine form is .5 + .5𝜖 2 , noting that it has an independent noise term (𝜖 2 instead of 𝜖 1 ) since the variables are not correlated. Even for the abstract analysis, we still set 𝑥 1 [{1}] = 1 and 𝑥 2 [{2}] = 1 if we wish to differentiate with respect to 𝑥 1 and 𝑥 2 , hence why there are still constants (constants are represented exactly in the zonotope domain).
For AD, affine transformations on variables correspond to affine transformations on all of their derivatives. Thus, because we will be computing not just one, but 4 affine transformations per variable (one for each component), it is advantageous to have an abstract domain that is as precise as possible for affine transformations. The zonotope domain is exact for affine transformations, thus making it attractive for abstract interpretation of AD. For the łrealž parts of 𝑥 3 and 𝑥 4 , we Applying any non-linear function, such as a 𝜎 (𝑥) or SoftPlus(x), to a zonotope or affine form is more challenging, as we must come up with an abstract transformer that employs a sound linearization of that function. Furthermore, for AD, we need a sound linearization for not just the activation function, but also for all of its derivatives. Hence, we need abstract transformers for SoftPlus(x), 𝜎 (𝑥), and multiplication. In Section 6.2, we describe how to automatically construct zonotope transformers for functions like SoftPlus(x) that other works [Jordan and Dimakis 2021;Singh et al. 2018a] cannot support, using the Chebyshev construction [Stolfi and de Figueiredo 2003]. The core idea is to add new noise symbols for each application of a non-linear function. The resulting affine forms after applying the SoftPlus(x) are shown above variable 𝑥 5 and below variable 𝑥 6 in Fig. 2. Since the 𝑑𝑥 𝑖 [1] and 𝑑𝑥 𝑖 [2] are merely intermediate variables for the computation of each 𝑥 𝑖 [{1}], 𝑥 𝑖 [{2}], and 𝑥 𝑖 [{1, 2}], we omit them for simplicity. However, it is important to recall that because of how we set up the semantics earlier, there is a data dependence between 𝑑𝑥 𝑖 [1] and 𝑑𝑥 𝑖 [2]. The zonotope domain will be able to leverage this dependence for improved precision.
We can now finally propagate the zonotope through the last affine layer. As mentioned, for the zonotope domain, this step is exact. The final affine forms (which collectively form the zonotope) for all derivatives are shown in Fig. 2. Furthermore, we can ultimately convert this output zonotope to interval bounds for each variable. Fig. 3a presents in red the bounds for all first derivatives computed from the zonotopes and likewise Fig. 3b shows in red the bounds on the second derivatives. We also show the result of performing the abstraction instantiated with the standard interval domain, denoted by the black bars in the respective plots. Lastly, we show the derivatives at randomly sampled points (𝑥 1 , 𝑥 2 ) ∈ [0, 1] × [0, 1] using only the concrete semantics (colored points). Both plots show that zonotopes are more precise for bounding both the first and second derivatives, and both abstractions soundly enclose the derivatives computed at scalar points. There is slightly more over-approximation for the bounds in Fig. 3b, as second derivatives require more computations, compounding the over-approximation inherent to the abstraction. For illustration simplicity, the example of Fig. 2 corresponds to a single pass of (abstract) forward-mode AD. However, to obtain all entries plotted in Figs. 3a and3b, we need to rerun forward-mode AD for multiple passes. Theorem 4.11 describes how many separate passes of forward-mode AD are required.
By understanding bounds on the derivatives, we can interpret and explain the magnitude of the contribution of each variable (1 𝑠𝑡 derivatives) or the interactive combinations of two variables (2 𝑛𝑑 derivatives) to the final output, and do this provably for an entire input region. We empirically evaluate this approach further in Section 7.2.
We now detail the necessary mathematical preliminaries to describe forward-mode automatic differentiation, particularly for higher derivatives.
As many of our operations involve working with sets we first detail our set-based notation. We let 𝑛 𝑘 represent the scalar binomial coefficient 𝑛! (𝑛-𝑘)!𝑘! . We will let P (𝑆) denote the powerset of a set 𝑆. Furthermore, we will write |𝑆 | to denote the cardinality of 𝑆.
We will also write P 𝑘 (𝑆) to denote the collection of all subsets of some set 𝑆 that have a specific cardinality 𝑘, equivalently {𝐴 ∈ P (𝑆) : |𝐴| = 𝑘 }. For instance, P 1 ({1, 2}) = {{1}, {2}}. For finite sets of integers in a given range, we may write {1, ..., 𝐷 }, e.g. {1, ..., 4} = {1, 2, 3, 4}. To denote the set of all possible partitions of a finite set of integers, we will write Part (𝑆).
The most basic implementation of first-order forward-mode AD is operator overloading with dual numbers [Griewank et al. 2000]. Intuitively, dual numbers are two-dimensional numbers that correspond to the value of a function at a point, and the value of the function's derivative at that point. Dual numbers are canonically defined as numbers of the form 𝑎 + 𝑏𝜖 with 𝜖 being an infinitesimal part, similar to the imaginary part of complex numbers. However, because our work later uses 𝜖 for denoting zonotope noise symbols, and because it is more intuitive when generalizing to higher derivatives, we present dual numbers simply as two-element tuples.
Definition 3.1. Dual numbers, are numbers of the form (𝑎, 𝑏) where 𝑎, 𝑏 ∈ R and for any realvalued differentiable function 𝑓 : R → R, we can lift the interpretation of 𝑓 to dual numbers by defining 𝑓 (𝑎, 𝑏) = 𝑓 (𝑎), 𝑓 ′ (𝑎) • 𝑏 .
One may also define all the arithmetic operations over dual numbers. For example for two dual numbers 𝑥 = (𝑎 1 , 𝑏 1 ) and 𝑦 = (𝑎 2 , 𝑏 2 ) one defines 𝑥 +𝑦 = (𝑎 1 +𝑎 2 , 𝑏 1 +𝑏 2 ), which intuitively encodes the notion of linearity. Rules for multiplication and division are similar, implicitly encoding the product and quotient rules. Hence, given a set of primitive functions 𝑓 𝑖 whose derivatives 𝑓 ′ 𝑖 we know analytically, we can propagate dual numbers through arbitrary compositions of these functions including ones involving arithmetic operations to compute first derivatives compositionally.
Higher-Order Derivatives. Just as forward-mode automatic differentiation can compute first derivatives using an alternate interpretation of functions and arithmetic, it can be extended to compute higher-order derivatives, as noted in Griewank et al. [2000]. However, higher-order derivatives are more challenging, because the formulas such as the product or chain rules must be lifted to their higher-order versions. For first derivatives there are 2 elements in the dual number tuple, hence one might expect that for 𝐷 th derivatives the generalization of a dual number will be a R 2 𝐷 tuple, which we will indeed show is the case. Intuitively, for a 𝐷 th -order derivative, each of the 2 𝐷 tuple elements represents a partial derivative with respect to a subset of the program variables. However, there is no loss of generality; indeed, we will show that this idea still allows one to differentiate with respect to the same variable multiple times.
This stems from the fact that in forward-mode AD, one is essentially propagating coefficients of a truncated 𝐷 degree multi-variate Taylor polynomial, as the coefficients (without the integer factorial division) are precisely the values of the derivatives. We define this intuition formally. Definition 3.2. For a variable 𝑥 𝑜 corresponding to some arithmetic function of variables 𝑥 1 through 𝑥 𝐷 , a truncated 𝐷-degree multi-variate Taylor polynomial is a 2 𝐷 length tuple of the form
..𝜕𝑥 𝐷 ) containing the 0 th through 𝐷 th -order partial derivatives with respect to 𝑥 1 ...𝑥 𝐷 . Furthermore, each element will not be symbolic but rather the result of evaluating that derivative at a given input point, hence the tuple is an element of R 2 𝐷 .
For a fixed 𝐷, there will be 𝐷 0 zeroth derivatives (the single łrealž part 𝑥 𝑜 ), 𝐷 1 first derivatives ( 𝜕 𝜕𝑥 1 through 𝜕 𝜕𝑥 𝐷 ) all the way up to 𝐷 𝐷 derivatives of order 𝐷 ś hence we can only compute one 𝐷 th -order derivative at a time, which is why the polynomial tuple is considered truncated. Thus, for any 𝑑 ∈ {0, 1, ..., 𝐷 }, there will be exactly 𝐷 𝑑 derivatives of order 𝑑 stored in the tuple. Furthermore, these derivatives will only be with respect to a predefined set of input variables.
Encoding Truncated Taylor Polynomials. Since each entry of a truncated 𝐷-degree Taylor polynomial is a derivative evaluated at a concrete real value, we will ultimately reduce reasoning about complicated polynomials to reasoning about their individual coefficients. However, to łdepackagež a tuple of length 2 𝐷 where each element has a specific semantic meaning (some particular partial derivative), we need an indexing scheme to access the individual elements. This will form a core development of our concrete semantics and will later allow us to capture data dependencies between program variables (corresponding to derivatives) when abstracting the semantics.
Variable Indexing. As pointed out by Paszke et al. [2021], maintaining clarity in how one indexes program variables corresponding to the various derivatives in AD is difficult. Therefore, a part of our contribution is to develop a novel set-based indexing scheme, generalizing existing higherorder AD implementations [Fike and Alonso 2011]. For a 𝐷 th -order truncated Taylor polynomial tuple, each of the 2 𝐷 elements corresponds to a particular element of the power-set P ({1, ..., 𝐷 }). Furthermore, that element of the power set is exactly the set of variables the partial derivative is with respect to. For instance, the empty set ∅ corresponds to the real part (e.g. 𝑥 𝑜 ) since it is the partial derivative with respect to no variables. Likewise, the singleton set {1} corresponds to 𝜕 𝜕𝑥 1 and the full set {1, ..., 𝐷 } corresponds to the 𝐷 th derivative 𝜕 𝐷 𝑥 𝑜 𝜕𝑥 1 ...𝜕𝑥 𝐷 . Hence, this indexing scheme respects the semantic meaning of the entries. Lastly, while the natural idea is to associate a unique program variable to each number in the indexing set, we can also associate the same program variable to multiple numbers of the set. Hence, we could associate a variable 𝑥 𝑖 to both 1 and 2, meaning the set {1, 2} corresponds to 𝜕 2 𝜕𝑥 𝑖 𝜕𝑥 𝑖 . In these cases, some of the tuple elements become redundant copies of the same derivative. More details on input state initialization are in Section 4.4.
Example 3.3. For the following truncated degree 2 Taylor polynomial, where 𝑥 1 is associated to 1 and 𝑥 2 is associated to 2, given by: 𝑇
Figure 4 presents the syntax of the language. Our language is imperative and syntactically supports arithmetic operations and differentiable function primitives. We will denote the variables in the
One can translate a pure mathematical function that uses only the differentiable primitives in our syntax into a well-formed program using an ANF conversion [Flanagan et al. 1993].
We now define a meta-semantics, which is a construction that takes as input a maximum derivative order 𝐷 to be computed and gives a concrete forward-mode AD semantics for computing derivatives up to that order using truncated Taylor polynomials. While we could map each syntactic variable 𝑥 𝑖 in the original program to a 2 𝐷 length real-valued array storing its truncated Taylor polynomial, we instead syntactically desugar this mapping so that each individual coefficient in the Taylor polynomial is identified by its own unique variable. Our novel set-based indexing scheme (Section 3.2) is used to distinguish these individual variables. We call this variable set the augmented variables, or AugVars, to distinguish it from SynVars, and we will access these augmented variables using the indexing scheme. The added benefit is that by semantically exposing each entry as an intuitively indexed, augmented variable, we can (a) reduce reasoning about complicated mathematical objects (Taylor polynomials) to reasoning about standard real-valued programs and (b) we can more easily track correlations between individual variables when later using a relational abstract domain.
Hence, for forward-mode AD, our concrete states 𝜎 will map 𝜎 : AugVars → R. However this means that if we have 𝑚 syntactic variables in the program source code (e.g. 𝑥 1 , ..., 𝑥 𝑚 ), then AugVars
Having now defined a concrete AD program state 𝜎 (which is what the semantics will ultimately return), we can formally define the concrete domain upon which the meta-semantics will be built.
Definition 4.3. The concrete domain D is P (AugVars → R) but this set is isomorphic to P (R |AugVars | ), hence we will also say that D = P (R |AugVars | ).
Intuitively, for a concrete state 𝜎 we can just enumerate all O (𝑚 • 2 𝐷 ) augmented variables into a single state vector with their corresponding Taylor coefficient value, thus this state vector 𝜎 is equivalently an element of R |AugVars | . Because our AD states are now standard mappings of (augmented) variables to real numbers, we can construct an interpreter for • 𝐷 using only a standard interpreter for real-valued imperative programs • : (𝑃, D) → D. We now formally define • 𝐷 in terms of the base interpreter • .
Definition 4.4. The base interpreter • : (𝑃, D) → D for assignment statements, and 𝐸𝑥𝑝𝑟 : (𝐸𝑥𝑝𝑟, D) → R for arithmetic sub-expressions, is given by the rules of Fig. 5.
), the evaluation of • returns ⊥. Likewise, any arithmetic operation with ⊥ returns ⊥. We can now use • as a subroutine to build an interpreter for performing concrete 𝐷 th degree forwardmode AD. Lastly, while • is defined in Fig. 5 only for elementary expressions (binary arithmetic and unary functions) we will still notationally write 𝐸𝑥𝑝𝑟 when 𝐸𝑥𝑝𝑟 is a non-elementary expression, as these are desugared to sequences of elementary operations by merely adding more augmented variables for storing all intermediate sub-expressions.
Definition 4.5. The concrete meta-semantics for performing forward-mode AD are a parametric semantics given by • 𝐷 : (𝑃, D) → D, where the parameter 𝐷 ∈ N >0 is the order of the highest derivative the AD semantics can compute. Each statement's rule is shown in the following section.
We now give the formal rules for evaluating • 𝐷 . We immediately see why the concrete metasemantics are parametric in 𝐷 ś choosing a 𝐷 governs how many times • will be called. This is why we refer to these as a meta-semantics ś choosing a different 𝐷 instantiates a different concrete semantics for each statement and expression. We also mention that in all cases, despite the complex form of the right-hand side expressions used in assignments, the entire expression can always be unpacked into the language primitives of Fig. 4, even if it requires introducing additional intermediate variables to store intermediate results. Lastly, this construction is defined imperatively, as these are a state-transformer semantics which update 𝜎 after each application of • .
Addition. The rule for addition is simple, as higher-order derivatives also follow linearity. The meta-semantics are given, where for all coefficients 𝑥 𝑖 [𝑆] in the truncated Taylor polynomial, we compute that coefficient by adding the corresponding coefficients of syntactic variables 𝑥 𝑗 and 𝑥 𝑘 .
Subtraction. Subtraction follows nearly identically to addition, with the only difference being the state gets updated via
Multiplication. The rule for multiplication follows directly from the Generalized Leibniz formula (𝑆 \ 𝑃 is set subtraction of 𝑃 from 𝑆).
While the summation term,
is given in a general form, the only primitive operations involved are addition and multiplication, hence when instantiating a concrete instance of these semantics for a fixed derivative order 𝐷 ∈ N >0 we can unroll the entire summation into primitive additions and multiplications of existing variables, and thus we can still evaluate each assignment in the for loop with the base interpreter of Fig. 5.
Constants. Constants are simple, as all their higher-order derivatives are zero. Hence, for any non-empty set 𝑆, the coefficient in 𝑥 𝑖 's Taylor polynomial indexed by 𝑆 will necessarily be zero.
Unary Functions. For compositions with unary functions 𝑓 : R → R, we will need to use Faa di Bruno's formula [Di Bruno 1857] as a generalization of the Chain rule to compute higherorder derivatives. The multivariate Faa di Bruno formula for a 𝐷 th -order derivative of a composite function is given below in Def. 4.6: We will see that for the subsequent unary function semantic rules, there is a direct correspondence between each term in Def. 4.6 and the respective rule; in particular, each possible partial derivative needed by the Faa di Bruno formula will be stored in a separate variable. Furthermore, as with the multiplication rule, though the notation is condensed, in the subsequent semantic rules, all primitive arithmetic operations are just additions (from 𝑃 ∈𝑃𝑎𝑟𝑡 (𝑆) ), multiplications (from 𝐵 ∈𝑃 ), variable lookups (e.g. due to 𝑥 𝑗 [𝐵]) or elementary function evaluations of 𝑓 and its analytically known derivatives. Further, when fixing a particular order of derivative, 𝐷, all possible partitions of all possible sets (e.g. 𝑃 ∈ 𝑃𝑎𝑟𝑡 (𝑆)) can be precomputed and enumerated (as they are finite), as can |𝑃 |. Hence, these seemingly complicated expressions still reduce to elementary operations.
The core requirement is that we are able to analytically compute the first through 𝐷 th derivatives of 𝑓 evaluated at the real part of 𝑥 𝑗 (itself stored in 𝑥 𝑗 [∅]). Thus, we restrict 𝑓 to only elementary functions whose derivatives of all orders can be known analytically as combinations of the other elementary functions in the language, hence why the language is restricted as shown in Fig. 4. Each of the 𝐷 derivatives of 𝑓 will be stored in a unique 𝐷 element array: 𝑑𝑥 𝑖 [0] through 𝑑𝑥 𝑖 [𝐷] within 𝜎, but as with the Taylor coefficients, each of these elements will be exposed as its own variable. Lastly, this is the only part that changes when defining the meta-semantics of different functions.
Division. Since Division is composition with the function 𝑓 (𝑥) = 1 𝑥 , we can encode it using the above technique. Furthermore, the derivatives of 1 𝑥 have a known elementary form that can be encoded using only primitive arithmetic operations (multiplications and divisions), which specifically is
which allows us to capture data dependencies across derivative terms.
Exp and Log . As before, for log and exp the only difference will be the computation of the 𝐷-derivatives array, which for the exponential function is given as
For the log function, the first derivative is just 1 𝑥 , hence we can reuse the functional form from the division rule to compute all further derivatives.
Sigmoid and SoftPlus. The sigmoid function 𝜎 𝑎 (𝑥) = 1 1+𝑒 -𝑎𝑥 is such that all of its derivatives can be given in terms of an elementary combination of other sigmoid functions. The 𝑑 th derivative is given as 𝑑 𝑑 𝑑𝑥 𝑑 𝜎 𝑎 (𝑥) = 𝑑 𝑘=0 (-1) 𝑑+𝑘 (𝑘!)(𝑆 𝑑,𝑘 )(𝑎 𝑑 )𝜎 𝑎 (𝑥)(1 -𝜎 𝑎 (𝑥)) 𝑘 , where 𝑆 𝑑,𝑘 are Stirling numbers of the second kind (which are just constants). Hence for any fixed 𝑑 ∈ N, we can combinatorially enumerate all the terms in order to have an analytical closed-form expression for the 𝑑 th derivative that uses only elementary operations (sums, products and constants) and the sigmoid function itself. This means that we need only compute the sigmoid once, which we do by having the base interpreter evaluate 𝑥 𝑖 [∅] = 𝜎 𝑎 (𝑥 𝑗 [∅]) (𝜎), after which we can fetch the variable 𝑥 𝑖 [∅] when computing each 𝑑𝑥 𝑖 [𝑑] augmented variable. This helps us to capture the dependency across derivatives.
Furthermore, because the sigmoid function is just the first derivative of the SoftPlus function, the 𝑑 th derivative of a SoftPlus is just the 𝑑 -1 th derivative of a sigmoid, hence we can leverage the exact same functional form when computing 𝑑𝑥 𝑖 [𝑑] for the SoftPlus function.
Sequencing. Having defined the construction to produce a semantics for computing the 𝐷 th derivative for individual assignment statements, we now note how sequencing of multiple statements works. We also note that the simplicity in the sequencing rule stems from the fact that programs are assumed to be in SSA form, hence no program variable will be overwritten.
Remark 1. When 𝐷 = 1, • 𝐷 computes exactly the dual numbers [Griewank and Walther 2008] and when 𝐷 = 2, • 𝐷 computes exactly the hyper-dual numbers [Fike and Alonso 2011].
A key contribution of our concrete semantics is to ensure that the transformed program which computes all derivatives is generated in such a way that the later abstract interpretation will be precise. For AD, the computation of derivatives involves heavy reuse of common sub-expressions, hence by sharing these sub-expressions across multiple variables, the data dependence between different derivatives can be made explicit. For instance, in the case of the exp function, all derivatives, 𝑑𝑥 𝑖 [𝑑], will be computed as 𝑑𝑥 𝑖 [𝑑] = 𝑥 𝑖 [∅] (𝜎) where 𝑥 𝑖 [∅] already stores the result of computing exp(𝑥 𝑗 ). Likewise, one can see similar data dependencies made explicit in the rules for 1 𝑥 and the 𝜎 𝑎 function. While this will not change the execution result of the concrete semantics, it does substantially improve the precision of the abstract interpretation, as will be seen in Section 7.
We now formally state the correctness as well as how one initializes the state to compute derivatives with respect to specific variables.
Input Variables. We must first decide which input variables we wish to differentiate with respect to. However, we first note the following important remark that provides a limitation on how many input variables we can differentiate with respect to in a single forward execution 𝑃 𝐷 (𝜎).
Remark 2. For any set of syntactic input variables 𝑥 𝑖𝑛 𝑖 ,..., 𝑥 𝑖𝑛 𝑗 when evaluating 𝑃 𝐷 (𝜎), we can only differentiate with respect to up to 𝐷 distinct syntactic variables.
This means that if the program 𝑃 has more syntactic input variables than the maximum derivative order 𝐷, some variables will necessarily be treated as constants (meaning not differentiated). Additionally, if we differentiate with respect to the same syntactic variable more than once (e.g. to compute 𝜕 2 𝜕𝑥 𝑖 𝜕𝑥 𝑖 ) the number of distinct variables we will be differentiating with respect to will be strictly less than 𝐷, again meaning some input variables will be treated as constants. However, as we will see, if one wishes to compute more derivatives, one may always rerun • 𝐷 with different input variables initialized to obtain their derivatives. Input State. After deciding which of the (up to 𝐷) syntactic input variables we wish to differentiate with respect to, we must properly initialize each of those syntactic input variables' corresponding augmented variables in the input state, as AD only produces correct derivatives if input states are initialized correctly [Griewank and Walther 2008]. We generalize this idea to arbitrary order derivatives, as described below.
Definition 4.7. For • 𝐷 , an initial state 𝜎 is valid if for each singleton set 𝑆 ∈ P 1 ({1, ..., 𝐷 }), exactly one augmented variable satisfies 𝑥 𝑖𝑛 𝑖 [𝑆] = 1 while all others satisfy 𝑥 𝑖𝑛 𝑗 [𝑆] = 0. Furthermore, for each 𝑆 ∈ P 𝑘 ({1, ..., 𝐷 }) where 𝑘 ≥ 2, all augmented variables satisfy 𝑥 𝑖𝑛 𝑖 [𝑆] = 0. For different singleton sets 𝑆, 𝑆 ′ ∈ P 1 ({1, ..., 𝐷 }), where 𝑆 ≠ 𝑆 ′ , the same input variable 𝑥 𝑖𝑛 𝑖 can satisfy both 𝑥 𝑖𝑛 𝑖 [𝑆] = 1 and 𝑥 𝑖𝑛 𝑖 [𝑆 ′ ] = 1. Setting two different augmented variables that are both associated with the same syntactic variable to 1 allows us to differentiate with respect to the same variable twice (e.g.
) instead of only being able to differentiate with respect to different variables. We now illustrate an example of a valid initial state.
Example 4.8. For 𝐷 = 2 and 𝑃 ≜ 𝑥 3 = 𝑥 𝑖𝑛 1 + 𝑥 𝑖𝑛 2 , and state 𝜎 given as
]] = 0 , then we have that 𝜎 is a valid input state to 𝑃 𝐷 that can be used to compute
]] ≠ 0 then 𝜎 would no longer be a valid initial state.
We next detail a lemma that says a correctly initialized input state computes valid derivatives of the input variables. This lemma also relates the syntactic variables of the program to the augmented variables computed by the interpreter and serves as the base case of our correctness theorem. Lemma 4.9. Let 𝜎 be a valid input state. Then for any syntactic input variable 𝑥 𝑖𝑛 𝑖 and any non-empty set 𝑆 ∈ P (𝐷) \ ∅, we have that
where 𝐴𝑐𝑡𝑖𝑣𝑒 (𝑆) is the list of all input variables 𝑥 𝑖𝑛 𝑗 satisfying 𝑥 𝑖𝑛 𝑗 [{𝑠 𝑗 }] = 1 for some 𝑠 𝑗 ∈ 𝑆. Proof. See Appendix. □
Computing Derivatives. Having established that a valid input state stores derivatives correctly, the idea is to now show that after executing each statement, the derivatives are still correct.
Theorem 4.10. Let 𝜎 be a valid input state and 𝑃 be a well-formed program, and let 𝜎 ′ = 𝑃 𝐷 (𝜎). Then for any variable 𝑥 𝑖 in 𝜎 ′ and any non-empty set 𝑆 ∈ P (𝐷) \ ∅, we have that
Proof. (Sketch) We show select cases with all other cases shown in the Appendix.
• Constants. For non-empty set 𝑆, and constant 𝑐 ∈ R,
) ∈R 𝑚 = 0. However, for any non-empty 𝑆, the interpreter always assigns 0 to 𝜎 ′ [𝑥 𝑖 [𝑆]].
• Addition. As the current state 𝜎 has been correctly initialized, before executing 𝑥 𝑖 = 𝑥 𝑗 +𝑥 𝑘 𝐷 (𝜎)
we inductively assume (by Lemma 4.9) that for each
After executing the meta-semantic rule for addition [𝑆]] by the rules of the base interpreter in Fig. 5. Thus, by substitution
Complexity. While Theorem 4.10 tells us that for the predetermined set of input variables, the derivatives will be correctly computed, as mentioned these derivatives are only with respect to a subset of the input variables. To compute all possible derivatives with respect to all possible input variables (for the full Jacobian or Hessian), we must rerun 𝑃 𝐷 multiple times for different (valid) input states 𝜎. In the simplest case, for 𝐷 = 1 (dual numbers), one needs to rerun 𝑃 1 exactly 𝑚 times for functions of the form 𝑓 : R 𝑚 → R 𝑛 to obtain the entire Jacobian. We can generalize this result for arbitrary 𝐷.
Theorem 4.11. For a program 𝑃 corresponding to a function 𝑓 : R 𝑚 → R 𝑛 where 𝑚 > 𝐷, to compute all possible derivatives up to the 𝐷 𝑡ℎ order, requires evaluating 𝑃 𝐷 exactly 𝑚+𝐷-1 𝐷 times.
Proof. (Sketch) We provide the basic intuition here. First, 𝑚+𝐷-1 𝐷 is the number of ways to select 𝐷 items from a set of size 𝑚, with replacement when order does not matter. The 𝐷 items we select are precisely the 𝐷 variable we wish to differentiate with respect to. The reason why replacement is allowed is because we can differentiate with respect to the same variable multiple times. Likewise, the reason why order does not matter is because 𝜕 𝜕𝑥 𝑖 𝜕𝑥 𝑗 = 𝜕 𝜕𝑥 𝑗 𝜕𝑥 𝑖 , hence we do not need to recompute them separately. □
As we can see, for dual numbers (𝐷 = 1) this reduces to 𝑚 1 = 𝑚 independent forward passes, as is well-known. While the complexity for higher-order AD seems expensive, it is the same complexity given in Griewank and Walther [2008], and more recent work [He 2019] has noted that specifically for higher-order derivatives, forward-mode AD is better suited than reverse mode.
Our construction is the first to provide a generic approach to abstractly interpret the semantics of higher-order AD. The key benefit of our construction is that by exposing each derivative term explicitly in a memory state as an augmented variable, we reduce the problem of reasoning about complex mathematical objects used in AD to reasoning about standard program states. Thus, we can readily apply existing numerical abstract domains. Furthermore, because our formulation is imperative, and our transformed AD program captures data-dependence across derivative terms, we can leverage the state-based semantics during analysis to conveniently use relational abstract domains that also track correlations between variables. This allows us to precisely track correlations across derivatives (including different orders). We first define the preliminaries of abstract interpretation suitable for our setting.
Definition 5.1. [Abstract Interpretation] The abstract interpretation primitives that our construction requires are given by the following symbols: (D # , 𝛾, ⊥ # , ⊤ # , • # ) where D # represents the set of abstract program states, 𝛾 : D # → D is a concretization function mapping abstract states to sets of concrete states. Further, ⊥ # and ⊤ # are the abstract domain's respective least and greatest elements. Lastly, • # are sound abstract transformers for statements in the language of Fig. 4.
Given an abstract domain, we can interpret the same program syntax over the abstract domain in a way that over-approximates the concrete program semantics, which are computed by • . Definition 5.2. [Soundness] An abstract transformer • # for statements is sound if the following holds: For any 𝜎 # ∈ D # and any valid syntactic expression 𝐸𝑥𝑝𝑟 in the language of Fig. 4:
where
Hence, we say that the abstract interpreter • # soundly over-approximates the semantics of the base interpreter • of Fig. 5.
To construct an abstract interpreter • # 𝐷 to soundly over-approximate • 𝐷 , we use an existing abstract interpreter • # defined over a numerical abstract domain that soundly over-approximates the semantics of the base interpreter • . Having defined the primitives needed for the construction in Defs. 5.1 and 5. (4) For any arithmetic expression 𝐸𝑥𝑝𝑟 in Fig. 4, we require a sound transformer for assignments with that expression, 𝑥 𝑖 = 𝐸𝑥𝑝𝑟 # : D # → D # . (5) (Optionally) we need an abstract test transformer 𝑥 𝑖 > 𝑐 # : D # → D # for refining abstract states with linear constraints. We also require that 𝛾 ( 𝑥 𝑖 > 𝑐 # (𝜎 # )) ⊆ 𝛾 (𝜎 # ).
We provide detailed explanation of Items 4 and 5 below, but we first note a (non-exhaustive) list of abstract domains satisfying these requirements.
Remark 3. The Interval domain, Octagon domain, Zonotope domain, DeepPoly domain, and Polyhedra domain can all be endowed with the necessary transformers to satisfy the requirements.
Assignment of Arithmetic Expressions. The abstract domain must have sound transformers for assignment with arithmetic expressions of Fig. 4, which includes arithmetic primitives (e.g., addition and multiplication) and constants but also differentiable unary functions (e.g., exp, log).
Arithmetic Primitives. The abstract domain must provide sound transformers for all primitive arithmetic operations: +, -, •, and /, as well as assignments with constants. For the interval domain, there are sound transformers for all of these. For domains like zonotopes, octagons, or polyhedra, one can always employ linearization of the non-linear operations of • and / as in Stolfi and de Figueiredo [2003] and Miné [2006] to construct the necessary transformers. Shi et al. [2019] has also shown how to construct DeepPoly abstract transformers for • and /.
Differentiable Unary Functions. As our language (Fig. 4) includes unary functions, the chosen abstract domain must provide sound transformers for the following functions: 𝜎 𝑎 , SoftPlus, exp, and log. As mentioned in Section 4.5, the derivatives of each of these functions are given in terms of elementary arithmetic combinations of functions already in the language (e.g., 𝜎 ′ 𝑎 (𝑥) = 𝜎 𝑎 (𝑥) • (1 -𝜎 𝑎 (𝑥)), hence this set of functions is łclosedž under n th derivatives. As with the arithmetic primitives, for the interval domain, there are sound transformers for all of these, while for zonotopes and polyhedra, one would need to soundly linearize these non-linear functions, such as the Chebyshev method for zonotopes [Stolfi and de Figueiredo 2003], the polyhedral linearization of Miné [2006], or using the methods of Ryou et al. [2021] for the DeepPoly domain.
Automatically Improving Precision. Beyond naively composing existing abstract transformers from a numerical domain, we also want to devise a method for automatically improving the precision of the analysis. Having an optimal transformer for one primitive function, does not necessarily imply one has an optimal transformer for all derivatives of that function, particularly if the derivative requires the composition of multiple primitive function transformers (e.g., 𝜎 ′ 𝑎 = 𝜎 𝑎 (𝑥) • (1 -𝜎 𝑎 (𝑥)) requires composing transformers for 𝜎 𝑎 and multiplication), since naive, off-the-shelf composition of numerical abstract transformers can be non-optimal.
To avoid this issue, we leverage analytical properties of the derivatives of the functions to systematically improve the precision of • # 𝐷 . We now describe an optional requirement on the domain that can improve the abstraction's precision. By using an abstract test operator 𝑥 𝑖 > 𝑐 # (𝜎 # ), we can refine the abstract state 𝜎 # using analytical information about the range of the function's derivatives. This allows us to enforce constraints, such as how every odd derivative of 1 𝑥 is necessarily negative. However, one can always define 𝑥 𝑖 > 𝑐 # (𝜎 # ) = 𝜎 # (the identity), hence why this step is optional. By requiring 𝛾 ( 𝑥 𝑖 > 𝑐 # (𝜎 # )) ⊆ 𝛾 (𝜎 # ), this step never loses precision.
Remark 4. We do not need a widening operator, since all the programs expressible in the syntax of Fig. 4 are loop-free. We also do not require a join operator ⊔ or meet operator ⊓.
We now provide an abstract meta-semantics which shows how to leverage an existing numerical abstract domain satisfying our criteria to construct a sound abstract interpreter for 𝐷 th -order AD. We highlight that the ease in defining our abstract semantics stems from design choices made in constructing our concrete semantics (Definition 4.5). However, merely using the abstract version of each arithmetic operator naively is imprecise, hence we will also see how the abstract meta-semantics leverages mathematical properties of differentiable functions to improve precision.
Definition 5.3. The abstract meta-semantics for performing 𝐷-degree Taylor polynomial forwardmode abstract AD are a parametric (abstract) semantics given by • # 𝐷 : (𝑃, D # ) → D # , where 𝑃 is a program in the syntax of Fig. 4. The abstract meta-semantics are parametric in both the order of derivative 𝐷, as well as the underlying numeric abstract domain D # and its associated abstract transformers • # . The abstract meta-semantics ultimately return an abstract state, 𝜎 # at the program exit. This is in contrast to the concrete meta-semantics of Definition 4.5 which are parametric only in 𝐷, (they have a fixed base interpreter, • ). This is because • # 𝐷 can use any set of abstract transformers • # satisfying the criteria of Section 5.1, thus the construction is configurable. This also allows the correctness to be easily established ś proving that the abstract meta-semantics soundly over-approximate the concrete meta-semantics reduces to showing that the instantiated numerical abstract domain over-approximates the base interpreter at each step.
Addition. The meta-semantics of abstract addition follow very similarly to the concrete case, by adding respective derivative terms (via linearity) albeit using the abstract interpreter • # of the chosen numerical abstract domain instead of the base interpreter • , noting that by the assumptions of Section 5.1, • # has sound transformers for (abstract) addition of two variables.
Multiplication. The rule for multiplication follows an abstracted form of the generalized Leibniz formula, again noting by assumption that • # has sound transformers for assignments involving the sum and product of variables and that 𝑆 \ 𝑃 represents set subtraction.
Constants. For this case, we need to assign each augmented variable to a constant (either 𝑐 or 0), hence all this requires is that the abstract domain has sound transformers constant assignment.
Unary Functions. As in the concrete meta-semantics, with unary functions, we must use Faa di Bruno's formula (Def. 4.6). In the computation of the first through 𝐷 th derivatives of 𝑓 (stored in 𝑑𝑥 𝑖 [1] through 𝑑𝑥 𝑖 [𝐷]), since we know the derivatives' analytical forms, we also know their valid ranges, hence we can refine the abstract state 𝜎 # to use this knowledge via the 𝑥 𝑖 > 𝑐 # abstract test transformer. This is helpful for enforcing domain-specific knowledge. Since the refinement depends on the particular function, we only show cases for specific unary functions.
As noted, division is composition with the function 𝑓 (𝑥) = 1 𝑥 , hence we again use Faa di Bruno's formula. Additionally, every odd derivative of 1 𝑥 is necessarily negative, hence we can apply the abstract test operator to refine 𝜎 # . Lastly because the derivatives (𝑑𝑥 𝑖 ) that are used to compute 𝑥 𝑖 [𝑆] for 𝑆 ≠ ∅ are given in terms of the real part 𝑥 𝑖 [∅], when instantiated with a relational domain like zonotopes, we are tracking correlations across derivative orders.
Exp and Log. Since all 𝑛 𝑡ℎ derivatives of exp are also exp, and the exp function is strictly positive, we can always refine the abstract state with this constraint. Additionally, we can exploit the fact that all derivatives of exp are also exp. Therefore, we only need to compute this once via the
), and then we can merely copy the result into the variables corresponding to each derivative 𝑑𝑥 𝑖 . However, the benefit of this is not only computational savings, but also the potential for improved precision if instantiating the construction with a relational abstract domain. Since all the 𝑑𝑥 𝑖 variables are directly used to compute the 𝑥 𝑖 [𝑆] for 𝑆 ≠ ∅, if • # is instantiated with a relational domain, one is able to track correlations across those derivatives.
Abstract Sequencing. Abstractly interpreting the sequencing of multiple statements follows similarly to the concrete meta-semantics.
Proc. ACM Program. Lang., Vol. 6, No. OOPSLA2, Article 161. Publication date: October 2022.
Remark 5. Since all variables and their derivatives are included in the abstract state 𝜎 # , when instantiated with a relational domain (e.g., zonotopes), this construction produces an abstract interpreter that relationally tracks dependencies across derivatives and derivative orders.
Intuitively the soundness follows from construction. By composing abstract transformers at each step that are sound for each primitive, the end construction soundly over-approximates the original concrete semantics. However, we now state this formally.
Theorem 5.4. (Soundness of Abstraction) For any program 𝑃 expressible in the syntax of Fig. 4, 𝐷 ∈ N, 𝜎 # ∈ D # and any 𝜎 ∈ 𝛾 (𝜎 # ), then we have that 𝑃 𝐷 (𝜎) ∈ 𝛾 ( 𝑃 # 𝐷 (𝜎 # )). Thus, we can compute a sound over-approximation of the set of values that the (possibly higherorder) derivatives take, provided all the requirements of Section 5.1 are satisfied.
Proof. (Sketch) We provide sketches of representative cases (others detailed in the Appendix):
• Addition: We need to show that for any 𝐷, 𝜎 # ∈ D # with 𝜎 ∈ 𝛾 (𝜎 # ) that the following holds:
in the for loop of the meta-semantic rule for 𝑥 𝑖 = 𝑥 𝑗 +𝑥 𝑘 𝐷 (𝜎), there is the corresponding application
), thus by the initial assumption that 𝜎 ∈ 𝛾 (𝜎 # ), and that fact that • # soundly over-approximates • (also by assumption) for every application (regardless of 𝐷) we know 𝑥 𝑖
• Division: We need to show that for any 𝐷, 𝜎 # ∈ D # with 𝜎 ∈ 𝛾 (𝜎 # ) that the following holds:
The idea is the same, we note that for each application 𝜎 = 𝑑𝑥 𝑖 [𝑑] = -𝑑 • 𝑥 𝑖 [∅] • 𝑑𝑥 𝑖 [𝑑 -1] (𝜎) in the concrete meta-semantics, that the corresponding application
) in the abstract metasemantics soundly over-approximates it. To address the application of 𝑑𝑥 𝑖 [𝑑] < 0 # (𝜎 # ), we note that in the concrete semantics, 𝑑𝑥 𝑖 [𝑑] will never be greater than 0 if 𝑑 is odd. Thus, even after applying this refinement, we still have that 𝜎 ∈ 𝛾 (𝜎 # ). Likewise, for each application
We formally describe instantiations of our framework, which we will later evaluate in order to study the effects of first vs. higher derivatives as well as relational vs. non-relational abstract domains.
While interval arithmetic has been applied to AD, our construction is more general, thus using the interval domain is but one instantiation of our framework.
Interval Domain. The interval domain, denoted D # = P (IR |𝐴𝑢𝑔𝑉 𝑎𝑟𝑠 | ), is among the simplest numeric abstract domains where in an abstract state 𝜎 # a variable is mapped to an interval [𝑎, 𝑏] ∈ IR, hence 𝜎 # : AugVar → IR. The interval domain is non-relational.
First Derivatives. Our construction can be used to produce the Dual Interval domain of Laurel et al. [2022a]. This is an instantiation of our framework where 𝐷 = 1 with the interval domain.
Higher Derivatives. We may also encode an interval abstraction of the hyperdual numbers of Fike and Alonso [2011] using our framework, which allows us to obtain interval bounds on second derivatives. This is done by instantiating our construction with 𝐷 = 2 and the interval domain.
To the best of our knowledge, prior to our work, there has been no general zonotope abstract interpretation for forward-mode AD, especially for higher-order derivatives.
Zonotope Abstract Domain. We first describe the zonotope abstract domain, which we denote with D # = Zono. The abstract state 𝜎 # ∈ Zono maps each variable to an affine form, where an affine form is a tuple (𝑐, 𝑔) with center 𝑐 ∈ R and 𝑔 ∈ R |generators | where |generators| is the number of noise symbols (also called generators). For a variable 𝑥, to denote its affine form, we will write
𝑔 𝑖 𝜖 𝑖 , where each noise symbol 𝜖 𝑖 ∈ [-1, 1]. To index a variable's affine form in the abstract state 𝜎 # , we may write
to access the coefficients of the center and noise symbol terms. Because the 𝜖 𝑖 are shared across variables, this is a relational domain; furthermore, the set of states encoded by 𝜎 # is exactly a zonotope. This can be seen from the concretization function 𝛾 : Zono → P (R |AugVars | ) defined as:
As mentioned in Albarghouthi [2021], |generators| grows as the program executes since new noise symbols (the 𝜖 𝑖 ∈ [-1, 1]) are dynamically added with each non-linear operation (e.g., multiplication, 𝜎 𝑎 (𝑥), etc.). However, the zonotope domain loses no precision when encoding affine transformations, since these can be done precisely in the domain.
One can also convert an affine form 𝑥 = 𝑐 + |𝑔 𝑖 |. This will prove useful when constructing sound transformers for arithmetic primitives.
Differentiable Zonotope Transformers. While the construction of Section 5 shows us how to produce a 𝐷 th -order forward-mode AD abstract interpreter • # 𝐷 using the chosen abstract domain D # , it assumes the existence of sound transformers, • # , for the arithmetic functions and gives no way to produce them. However, for zonotopes, prior works [Jordan and Dimakis 2021;Singh et al. 2018a] give sound transformers for some (e.g., tanh), but not all of the functions we support. Therefore, when instantiating with zonotopes, we need an automatic construction for sound differentiable function transformers. To do so, we use the Chebyshev construction [Fryazinov et al. 2010;Stolfi and de Figueiredo 2003 . We refer to this formula as the Chebyshev construction for 𝑓 . The Chebyshev construction gives us sound zonotope transformers for the following functions: SoftPlus a (𝑥), log(𝑥), exp(𝑥), and 1 𝑥 . For the other functions, we can use existing zonotope transformers (e.g., those in Singh et al. [2018a]) or just default to the interval domain transformers. We can then collectively use all of these abstract transformers in the general construction of • # 𝐷 . Theorem 6.2. When 𝐷 = 2 and when the abstract domain is the zonotope domain, • 𝐷 computes a zonotope abstraction of the hyper-dual numbers of Fike and Alonso [2011].
Precision. For many functions, the zonotope transformers are strictly more precise than the interval ones. As mentioned in [Fryazinov et al. 2010;Stolfi and de Figueiredo 2003], the Chebyshev construction is optimal in the input-output plane, hence it gives the same guarantees as Singh et [2018a]. The zonotope transformers in this case are incomparable to the interval ones. We also found that existing abstract test transformers for zonotopes [Jeannet and Miné 2009;Singh et al. 2017] do not boost precision, but increase cost. Hence, when instantiating our construction with zonotopes, we use the identity function 𝑥 𝑖 > 𝑐 (𝜎 # ) = 𝜎 # for the tests.
We now present a set of case studies that highlight the benefits of an abstraction supporting both higher-order derivatives and more precise abstract domains.
We briefly describe the experimental setup and the network architectures used in our experiments. We ran our experiments on a 3.70 GHz Intel Xeon W-2135 CPU with 32 GB of main memory. For the first case study on the robust interpretation of first-and second-order effects, we trained a 3-layer neural network with 5 inputs, 64 neurons in the first hidden layer, 10 neurons in the second hidden layer, and 1 output neuron to approximate a 5-input, 1-output synthetic function polynomial of degree two with five interactions; the network uses the SoftPlus activation function after every layer. For the Lipschitz certification case study, we trained four fully connected networks on the MNIST [LeCun 1998] image classification dataset (which contains 70,000 28 × 28 images of handwritten digits) ś three are smaller networks with 3, 4, and 5 layers having 100 neurons in each hidden layer, and one is the FFNNBig architecture from Singh et al. [2019] consisting of 4,106 neurons (4 hidden layers of 1,024 neurons each). Each hidden layer is followed by a SoftPlus activation function. All networks attain over 98% accuracy on the test set. For the perturbation function, we consider the Haze perturbation studied in Paterson et al. [2021]. For all experiments, we also present runtimes for the analyses.
Implementation. The instantations of our framework for both first and second derivatives on the interval and zonotope domains were all implemented in PyTorch.
As neural networks are notoriously hard to understand, significant work has focused on constructing more interpretable explanations, in order to understand which features are most relevant to the network's outputs. A standard way of interpreting and explaining neural networks involves computing derivatives of the network's outputs with respect to their inputs [Ancona et al. 2018;Janizek et al. 2021]. However, prior work has shown that it is not enough to generate explanations for scalar points; rather, explanations should be robust, hence why Fel et al. [2022] used robust interval explanations. Our case study serves to show how our construction allows a user to compute provably robust interpretations via derivative-based first-order feature attributions and secondorder feature interactions. For our experiments, we uniformly generate 5 random inputs in the range [0, 1) and enclose each input in an interval of ±0.01. We then soundly bound the first and second derivatives of the neural network (trained to approximate a polynomial as described in Section 7.1) with respect to these input ranges. We repeat this experiment for 5 different seeds.
Table 1 and Table demonstrate the precision improvement of the zonotope domain on bounding the Jacobian and the Hessian, respectively. In both tables, each cell presents the geometric mean over 5 trials of the ratio of the interval-bounded derivative's width over its zonotope-bounded counterpart. For zonotopes and intervals, the average runtimes across all trials are 0.014 and 0.009 seconds for the first-order information and 0.23 and 0.11 seconds for the second-order information, respectively. This demonstrates that the substantial increase in precision when using the zonotope domain does not incur a large runtime overhead. Furthermore, performing the separate passes to compute each Jacobian or Hessian term takes virtually the same time, meaning that the input one differentiates with respect to does not affect the runtime.
We next show detailed results for the trial that yields the median improvement in precision when using zonotopes over intervals (the other trials show the same trend). Fig. 6 presents the bounds on the first derivatives (the Jacobian). The entries on the x-axis correspond to the first derivative of the network's output node with respect to each of the 5 input variables. The value of the y-axis denotes the value of the derivative. We compare derivatives computed at sampled scalar points with both an interval and zonotope instantiation of our framework for first derivatives. The zonotope bounds on the Jacobian entries are much tighter than the bounds computed via the interval domain, both of which enclose the sampled points (shown in multiple colors).
We also compute the second derivatives to study the magnitude of interactions between input variables over an entire region of the input space. Fig. 7 presents the results (for the same trial as in Fig. 6). The entries along the x-axis represent Hessian terms instead of just individual input derivatives (as in Fig. 6). As before, the y-axis represents the magnitude of the (second) derivative. Zonotope bounds on the second derivatives are also much more precise at enclosing the sampled points than standard interval domain bounds, which stems from the fact that relational domains allow us to preserve correlation across derivative orders, thus improving precision.
Our second case study involves bounding the local Lipschitz constant of neural networks with respect to semantic perturbations, as in [Laurel et al. 2022a;Yang et al. 2022]. We study networks trained on MNIST data and consider the Haze perturbation, defined as the pixelwise transformation 𝑝 𝛼 (𝑥 𝑖 ; 𝛼) = (1 -𝛼)𝑥 𝑖 + 𝛼, where 𝑥 𝑖 is the value of each pixel and 𝛼 represents the haze amount. Let 𝑓 denote a neural network. Then, given an image 𝑥, we compute a sound upper bound on the (local) Lipschitz constant of 𝑓 (𝑝 𝛼 (𝑥; 𝛼)) with respect to 𝛼 over a range of haze values, specifically where 𝛼 lies in [0, 𝛼 𝑚𝑎𝑥 ]. The parameter 𝛼 𝑚𝑎𝑥 represents the maximum haze amount, and we consider values of 𝛼 𝑚𝑎𝑥 ∈ {10 -𝑘/4 • 2 : 𝑘 ∈ [2, 18]}. We evaluate on the zonotope domain (enabled by this work) and compare with Laurel et al. [2022a] (as no other work can handle this setting).
Fig. 8 shows the Lipschitz certification results for the 3-layer (left), 4-layer (center), and 5-layer (right) MNIST networks. Fig. 8a presents the computed Lipschitz constant bounds. The x-axis shows the value of 𝛼 𝑚𝑎𝑥 , and the y-axis shows the upper bound on the Lipschitz constant computed with the zonotope and interval domains (smaller is better). Fig. 8b presents the increase in precision of the zonotope domain. The x-axis again shows the value of 𝛼 𝑚𝑎𝑥 , and the y-axis shows the ratio of the interval-bounded Lipschitz constants over the zonotope-bounded Lipschitz constants (larger is better). For each plot, the x-axis uses a logarithmic scale, while the y-axis uses a linear scale. Each data point is the average over the first 1,000 correctly classified test set images. For the 3-, 4-, and 5-layer networks, the average zonotope runtimes are 2.8, 4.1, and 5.9 milliseconds per image and the average interval runtimes are 2.9, 3.2, and 4.0 milliseconds per image, respectively. These runtimes showcase how our analysis, even when using a precise relational domain, is fast and scalable.
AD with zonotopes is always more precise than AD with intervals ś the zonotope-bounded Lipschitz constants are always smaller. Compared to intervals, zonotopes are up to 6×, 31×, and 216× more precise and the computed Lipschitz constants are up to 3374, 39105, and 382755 smaller for the 3-, 4-, and 5-layer networks, respectively. This showcases how AD with intervals is much more over-approximate for deeper networks, while zonotope AD retains much more precision.
For each network, we see a similar trend across values of 𝛼 𝑚𝑎𝑥 . Fig. 8a shows that as 𝛼 𝑚𝑎𝑥 (i.e., the range of perturbation) increases, the absolute difference between the interval-and zonotopebounded Lipschitz constant increases significantly. For small perturbation ranges (when 𝛼 𝑚𝑎𝑥 < 10 -3 ), the Lipschitz constant between the two domains are similar in magnitude. Fig. 8b shows that for small values of 𝛼 𝑚𝑎𝑥 , both domains are quite precise, since the input range is very small. For large values of 𝛼 𝑚𝑎𝑥 , both domains incur larger degrees of over-approximation (though intervals are much more over-approximate). This is why the relative increase in precision of zonotopes over intervals is not as great for both ends of the input range. Zonotopes are relatively most precise for an input range that is around 10 -1.5 . The maxima in Fig. 8b exactly correspond to the points in Fig. 8a where the zonotope curve already asymptotes (i.e., induces little over-approximation), while the interval curve is still very over-approximate. Zonotopes allow us to obtain much tighter bounds on neural networks' Lipschitz constants with respect to semantic perturbations than the previous interval-domain work of Laurel et al. [2022a]. Our approach can also scale to larger networks: Fig. 9 shows the Lipschitz certification results on the FFNNBig architecture. All axes are in log scale for clarity; as before, each data point is the average over the first 1,000 correctly classified test images. The average zonotope and interval runtimes for this experiment are 0.58 and 0.12 seconds per image, respectively. Compared to the interval domain, our zonotope domain construction obtains Lipschitz constants that are up to 11,850× more precise, and the computed Lipschitz constants are up to 4.67 × 10 8 smaller. Again, we can observe that using the more precise zonotope domain (versus just intervals) increases analysis precision significantly for larger networks, highlighting the importance of our analysis.
To the best of our knowledge, there is not any work that provides a general construction for abstract interpretation of higher-order AD.
Higher-Order AD. Higher-order forward-mode AD has been explored going back to Griewank et al. [2000] and Karczmarczuk [2001]. Pearlmutter and Siskind [2007] also proposed a forward-mode scheme; then, later implementations in JAX [Bettencourt et al. 2019], 𝜆 𝑆 [Sherman et al. 2021], and the work of Huot et al. [2021] also followed in this spirit. In particular, 𝜆 𝑆 also supports Clarke Jacobians. However, all of these works only give concrete (and not abstract) semantics for AD. Additionally, unlike our concrete semantics, Huot et al. [2021]; Karczmarczuk [2001]; Krawiec et al. [2022]; Pearlmutter and Siskind [2007]; Sherman et al. [2021] all define their concrete semantics in functional styles, whereas we use an imperative semantics (that generalizes Fike and Alonso [2011]) for ease in building the formalism. Lastly, while we refer to our semantics as propagating tuples of Taylor coefficients, these are technically closer to Tensor coefficients of Griewank and Walther [2008], as our tuple coefficients are unscaled by factorial terms.
Abstract Interpretation of AD. As pointed out in Laurel et al. [2022a] and Vassiliadis et al. [2016], there is very little work on Abstract Interpretation for AD. Most recently, Laurel et al. [2022a] developed an interval abstract interpretation for bounding Clarke Jacobians; however, they cannot support relational domains or higher derivatives. Jordan and Dimakis [2021] provide a way to propagate zonotopes through vector-Jacobian products, but their formalism, soundness guarantees, and implementation are only valid for first derivatives. The benefit of their work is that beyond improving precision, by capturing linear dependencies between derivative terms, they can solve linear programming optimization problems defined over the derivative terms in closed form.
Other works like Vassiliadis et al. [2016] and Mangal et al. [2020] use the interval domain to bound AD, but their works cannot use any relational domain (e.g., zonotopes or polyhedra). Deussen [2021] can abstract higher-order derivatives and in fact use derivatives up to third order for abstract sensitivity analysis (applied to approximate computing); however, like the other works, they are also restricted to only the interval domain. Immler [2018] develops a solver to bound the solutions of differential equations using interval and zonotope bounds on first and second derivatives. However, all their derivatives are evaluated purely symbolically, and so they make no use of automatic differentiation; thus, their work suffers scalability problems. Despite their derivative computations being purely symbolic, evaluating those symbolic expressions using zonotopes to capture linear dependencies between first and second-order derivative terms does greatly improve the precision of their verified second-order Runge Kutta solver. Also motivated by verified ODE solutions, Bendtsen and Stauning [1996] and Stauning [1997] develop a validated AD framework for computing interval arithmetic liftings of the higher-order derivatives needed to compute higher-order Taylor series expansions of the solution to an ODE, but their work is restricted to only the interval domain.
Meta-Abstract Interpretation. There are a few works which provide meta-abstract interpretations [Cousot et al. 2019;Reps and Thakur 2016;Singh et al. 2018b;Sotoudeh and Thakur 2020], meaning they do not propose a new abstract domain, but rather a new construction that can be instantiated with different abstract domains. However, none of these works target AD. Additionally, there are techniques to automate the construction of transformers for a variety of numerical abstract domains [Miné 2006], but the precision and scalability of these transformers are usually suboptimal.
We developed the first general construction for abstract interpretation of automatic differentiation that supports both higher-order derivatives and virtually any numerical abstract domain. We instantiate our method with both intervals and zonotopes, and in the latter case, we show how to relationally leverage dependencies across derivatives to further improve precision. Our evaluation demonstrates the applicability and scalability of our technique.
Proc. ACM Program. Lang., Vol. 6, No. OOPSLA2, Article 161. Publication date: October 2022.
) ∈R 2 and Proc. ACM Program. Lang., Vol. 6, No. OOPSLA2, Article 161. Publication date: October 2022.
Received 2022-04-15; accepted 2022-09-01 Proc. ACM Program. Lang., Vol. 6, No. OOPSLA2, Article 161. Publication date: October 2022.