\section{Introduction}
   
    \label{sec:intro}
    Nowadays, Deep Neural Networks~(DNNs) are successfully deployed and show remarkable performance in many challenging applications, including facial recognition~\cite{liu2015deep,zhang2016joint}, autonomous driving~\cite{liang2019multi,pang2020clocs}, game playing~\cite{silver2017mastering,vinyals2019grandmaster}, and code completion~\cite{liu2022unified,bui2021infercode}.
    Together with the wide deployment of DNN models, the reliability issue of DNN-based systems has become a serious concern, where malfunctioning DNN-based systems have led to serious consequences such as fatal traffic accidents~\cite{drivingfatal}.
    
    Numerical defects and their resulting numerical failures\footnote{A \emph{numerical failure}~(in short as failure) is a type of failure in DNN-based systems such as a divide-by-zero failure.
    A \emph{numerical defect}~(in short as defect) is a location in a system's code that potentially exhibits a numerical failure such as the location of a division operation. Existing work mainly uses the phrase of  ``numerical bugs'', which indicates either a defect~\cite{zhang2020detecting} or failure~\cite{yan2021exposing}.} are serious threats to the reliability of DNN-based systems,  and deserve attention for two main reasons. 
   
   
   
    First, numerical defects are one of the most frequent defects in DNN-based systems. Empirical studies~\cite{zhang2018empirical,zhang2019empirical} on GitHub and Stack Overflow show that numerical defects, as a representative of DNN defects, widely exist in the code of DNN-based systems.
    Moreover, since numerical defects exist at the architecture level, any model using the architecture extensively inherits these defects. 
    Second, numerical defects can result in serious consequences. Once numerical defects are triggered, the faulty DNN model will output \NaN or \INF instead of producing any meaningful prediction, resulting in numerical failures~\cite{zhang2019empirical}. Thus, numerical defects  hinder the application of DNNs in scenarios with high reliability requirements such as threat monitoring in cybersecurity~\cite{problemaisecurity} and controlling in networks and cloud systems~\cite{sharma2016machine,jay2019deep}.
    
    \begin{figure}[t]
        \includegraphics[width=\linewidth]{figures/workflow_new.pdf}
        \vspace{-2.25em}
        \caption{Workflow for reliability assurance against numerical defects in DNN architectures.
        The left-hand side shows three stages and the right-hand side shows corresponding examples.
        \sysname supports all three stages, and is the first automatic approach for stages \textcircled{2} and \textcircled{3} to our knowledge.}
       
        \label{fig:workflow}
        \vspace{-1.75em}
    \end{figure}
    
   
   
   
    To address numerical defects in DNN architectures, in this paper, we propose a workflow of reliability assurance (as shown in \Cref{fig:workflow}), consisting of three tasks: 
    static defect detection, system test generation, and precondition-fix suggestion, along with our proposed approach to support all the three stages (the approach details are described later in this section).
    
  
   
   
   
   
   
    
    \emph{Static Defect Detection.}
    As the first stage, we need to conduct static analysis to detect all possible numerical defects given a DNN architecture as the input.
    A DNN architecture is usually encoded by program code as the example shown in \Cref{fig:dnn-code-example}.
    First, for training, the user executes the program with the architecture on training/validation data, attains the model weights, and stores them in a weight file.
    The architecture along with the weights is named a model.
    Then, for inference, the user loads the weight file to AI chips or memory, executes the same program with the inference sample and weights as arguments, and gets the model prediction result as program outputs.
    We focus on detecting inference-time numerical defects, which are locations~(operators) that potentially exhibit numerical failures in the inference phase,
    since inference-time numerical defects can be triggered after the model is deployed and thus are more devastating than training-time defects~\cite{odena2019tensorfuzz,zhang2020detecting}.
    
    \begin{discard}
    \textbf{Failure-Exhibiting-Test Generation.}
    Then, we generate failure-exhibiting tests to confirm the feasibility of numerical failures and help fault localization.
    We formulate two types of tests: \emph{training/inference system tests}~(in short as system tests) and \emph{inference unit tests}~(in short as unit tests).
   
    A failure-exhibiting system test is a tuple of training example~($\vxtrain$ in \Cref{fig:workflow}) and inference example~($\vx$ in \Cref{fig:workflow}), such that after using the training example for training the architecture, the resulting model contains a numerical failure which can be triggered by inputting the inference example.
    Therefore, failure-exhibiting system tests can confirm the existence of feasible numerical failures exhibited by real training data and inference data entry.
    A failure-exhibiting unit test is a tuple of model weights ($\vw_1$, $\vw_2$ in \Cref{fig:workflow}) and inference example~($\vx$ in \Cref{fig:workflow}), such that the numerical failure is triggered when the inference example is given to the architecture with these model weights.
    The unit test may not confirm the feasibility of numerical failures since the generated weights may not be feasible via training, but many distinct unit tests can help automatic or manual fault localization.
    \end{discard}
    
    \emph{System Test Generation.}
    After possible numerical defects are detected, in this stage, we need to confirm the feasibility of these numerical defects via generating failure-exhibiting system tests.
    A failure-exhibiting system test is a tuple of a training example\footnote{In real settings, a relatively large number of training examples are used to train the architecture, but generating a single training example to exhibit failures (targeted by our work) is desirable for ease of debugging while being more challenging than generating multiple training examples to exhibit failures.} ($\vxtrain$ in \Cref{fig:workflow}) and an inference example~($\vx$ in \Cref{fig:workflow}), such that after  the training example is used to train the architecture, applying the resulting model on the inference example exhibits a numerical failure.
    Therefore, failure-exhibiting system tests can confirm the existence of feasible numerical failures exhibited by real training data and inference data entry.
   
    
    \begin{discard}
    \textbf{Fault Localization.}
    With failure-exhibiting system and unit tests, we locate the faults.
    A fault corresponds to a program location, which is the output of this stage.
    Both manual inspection and automatic approaches~\cite{wardat2021deeplocalize} can conduct fault localization.
    Specifically, for coverage-based fault localization tools~\cite{abreu2009practical,abreu2009spectrum,abreu2006evaluation,lou2021boosting}, many distinct failing unit tests are required.
    In this paper, we do not focus on this task.
   
   
    \end{discard}
    
    \emph{Precondition-Fix Suggestion.}
    For a feasible numerical defect, in this stage, we need to fix it.
    Given the location where users prefer to impose a fix, we suggest valid fixes for the numerical defect.
    To determine the fix form,
    we have inspected the developers' fixes of the numerical defects collected by \citeauthor{zhang2020detecting}~\cite{zhang2020detecting} by looking at follow-up Stack Overflow posts or GitHub commits.
    Among the 13 numerical defects whose fixes can be located,
    12 fixes can be viewed as explicitly or implicitly imposing interval preconditions on different locations, such as before inputs, before weights, and  before vulnerable operators.
   
    Thus, imposing an interval precondition, e.g., by clipping the input for the vulnerable operator(s), is an effective and common strategy for fixing a numerical defect.
    In this stage, we need to suggest precondition fixes for the given location.
    
    
   
   
    
   
   
   
   
    
   
   
   
   
   
   
   
   
   
   
   
   
   
    
    
    To support all the three stages of the \textbf{R}ealiability \textbf{A}ssurance process against DNN \textbf{NUM}erical defects, in this paper, we propose the \sysname approach, which contains novel technical solutions in two major aspects: 
    (1) a DNN static analysis framework to provide highly efficient and effective support for the first-stage task (static defect detection) than a recent state-of-the-art approach~\cite{zhang2020detecting}, and 
    (2) two novel techniques~(the technique of two-step test generation and the technique of abstraction optimization) to provide automatic support for the second-stage and third-stage tasks.
   
    
   
    
   \textbf{DNN Static Analysis.} 
   The static analysis framework in \sysname provides efficient abstractions for 84 DNN operators, 17 more than those supported by the recent state-of-the-art approach~\cite{zhang2020detecting}.
   Moreover, the framework incorporates two novel abstraction techniques (backward fine-grained node labelling and dynamic partitioning) to improve the analysis tightness while maintaining the efficiency by adaptively determining the most suitable abstraction granularity based on data-flow dependencies.
   This framework allows not only more effective static defect detection than the recent state-of-the-art approach~\cite{zhang2020detecting}~(\textbf{stage \textcircled{1}}), but also efficient abstraction optimization for precondition-fix suggestion~(\textbf{stage \textcircled{3}}).    
    
    \begin{discard}
    We propose the \sysname approach consisting of four novel abstraction and generation techniques.
    \sysname supports the preceding three out of four tasks for testing and debugging numerical failures in DNN architectures: 
    detecting numerical defects~(\textbf{stage \textcircled{1}}),
    generating failure-exhibiting tests~(\textbf{stage \textcircled{2}}),
    and suggesting precondition fixes~(\textbf{stage \textcircled{4}}).
   
    To the best of our knowledge, \sysname is the first approach for automatically generating failure-exhibiting tests and suggesting precondition fixes.
    \end{discard}
    
    \textbf{Two-Step Test Generation.}
    We propose the technique of two-step generation for the task of failure-exhibiting system test generation~(\textbf{stage \textcircled{2}}). 
    First, we generate the pair of model weights and inference example for exhibiting the numerical failure. 
    Second, we generate a training example that leads to these model weights after training.
    We leverage an existing gradient back-propagation approach~\cite{yan2021exposing} to conduct the first step.
    For the second step, we transform the task into finding an input that leads to specific model gradients, where we extend the deep-leakage-from-gradients~(DLG) attack~\cite{zhu2020deep} by incorporating the straight-through gradient estimator~\cite{bengio2013estimating} for second-order optimization to achieve this goal.
    
      
    \textbf{Abstraction Optimization.}
    We propose the novel technique of abstraction optimization for precondition-fix suggestion~(\textbf{stage \textcircled{3}}). 
    We observe that a valid precondition fix is essentially a set of interval domains for some operators such that the input to each  vulnerable operator is always valid.
   
    Therefore, built upon our static analysis framework that provides differentiable interval abstractions, we propose the novel technique of abstraction optimization to ``deviate away'' the input range of a vulnerable operator from the invalid range guided by gradients of abstraction interval endpoints.
    This technique effectively solves the gradient vanishing and explosion problem encountered in abstraction optimization by leveraging the gradient sign instead of gradient value during the optimization.
   
    
    \begin{discard}
    \textbf{Two-Step Test Generation.}
    We formulate the failure-exhibiting system test generation task~(\textbf{stage \textcircled{2}}) as a two-stage problem: 
    First, we generate the pair of model weights and inference example that exhibit the numerical failure;
    Second, we generate training example that leads to these model weights after training.
    We leverage existing gradient back-propagation approach~\cite{yan2021exposing} to complete the first step.
    For the second step, we transform the task into the problem of finding input that leads to specific model gradients, where we extend the deep-leakage-from-gradients~(DLG) attack~\cite{zhu2020deep} to achieve this goal.
    Our extension of DLG attack incorporates the straight-through gradient estimator~\cite{bengio2013estimating} to automatically enable second-order optimization for general DNN architectures.
    
    \textbf{Abstraction Optimization.}
    To produce precondition fixes~(\textbf{stage \textcircled{3}}), we observe that a valid precondition fix is essentially a set of interval domains for some operators such that the input to the vulnerable operator is always valid.
    For example, for DNNs with an \CodeIn{Sqrt} node, a valid precondition fix is interval domains where any input picked from domains always brings a non-negative input to the \CodeIn{Sqrt} node, since negative input to \CodeIn{Sqrt} node triggers numerical failure.
    With this observation, we build a static analysis framework over DNNs to produce \emph{differentiable} interval abstractions for any operator within the DNN architecture and propose the novel technique of abstraction optimization to ``deviate away'' the range of vulnerable operator from the invalid range guided by gradients of abstraction interval endpoints.
    This technique effectively solves the gradient vanishing and explosion problem encountered in abstraction optimization by (1)~iterating over interval span expansion step and interval center adjustment step; (2)~back-propagating the gradient sign and leveraging the gradient sign instead of gradient value to update the interval centers.
    \end{discard}
 
    
    \begin{discard}
    \textbf{Abstraction Techniques for Static Defect Detection.}
    In \sysname, we implement an efficient static analysis framework that can generate tight and sound interval abstractions~\cite{cousot1977static,blanchet2003static} for given DNN architecture.
    The framework includes novel abstraction techniques to reduce precision loss while avoiding efficiency.
   
   
   
   
    
    \textbf{Extended DLG Attack for System Test Generation.}
    We extend the deep-leakage-from-gradients~(DLG) attack~\cite{zhu2020deep} by incorporating the straight-through estimator~\cite{bengio2013estimating} to effectively generate failure-exhibiting system tests.
   
   
   
   
   
    
    \textbf{Bilevel Optimization for Unit Test Generation.}
    Our goal is to generate massive distinct failure-exhibiting unit tests.
    Existing tools, e.g., GRIST~\cite{yan2021exposing}, can only generate a single failure-exhibiting unit test.
    Instead of generating a single test, \sysname generates \emph{failure-exhibiting interval abstraction domains}.
    Therefore, to generate failure-exhibiting unit tests, we only need to sample from these interval domains, which is efficient for massive generation.
   
   
    To generate failure-exhibiting interval domains, we propose a novel bilevel optimization technique which aims to maximize the interval spans of model weights and the inference example to encourage test diversity while preserving the failure-exhibiting property.
   
    
    
    \textbf{Abstraction Optimization for Precondition-Fixes Suggestion.}
    According to the empirical study discussed before, the most common fix for numerical defects is to impose interval preconditions on some specific nodes in the DNN computational graph.
   
    Therefore, generating precondition fixes is essentially searching for interval domains for given locations such that any values picked from the domains cannot trigger numerical failures.
    For example, for DNNs with an \CodeIn{sqrt} node, we need to find interval preconditions for some specific nodes so that the input of \CodeIn{Sqrt} node is within $(0, +\infty)$.
    Based on the DNN analysis framework in \sysname, we propose the novel technique of abstraction optimization to efficiently search for such interval domains~(i.e., generate preconditions) given arbitrary imposing locations.
   
    Our technique solves gradient vanishing and explosion problems as encountered by existing techniques.
   
   
    
   
   
   
   
   
   
   
    
    
   
   
   
   
   
    \end{discard}
    
    We implement a tool for \sysname.\footnote{Code will be made public after the de-anonymization process.}
   
    With the tool, we evaluate \sysname on a benchmark of $85$ real-world DNN architectures with numerical defects collected by \citeauthor{yan2021exposing}~\cite{yan2021exposing}.
    The evaluation results show that \sysname is both effective and efficient in all three stage tasks for DNN reliability assurance. 
    (1)~\sysname has stronger defect detection capability  than the state-of-the-art approach named DEBAR~\cite{zhang2020detecting}, being able to detect numerical defects on $>50\%$ more DNN programs.
    (2)~\sysname can efficiently generate failure-exhibiting system tests for 85.8\% cases while the baseline approach is successful for only 73.8\% cases and spends $2.76\times$ running time.
    \begin{discard}
    (2)~\sysname can efficiently generate massive failure-exhibiting unit tests for 6.6\% more cases compared to the best baseline approach.
    A manual study reveals that almost all failed cases correspond to impractical numerical defects.
    \end{discard}
    (3)~\sysname can efficiently suggest precondition fixes with 100\% success rate while the baseline approaches either require $> 10\times$ running time or fail on $> 25\%$ more cases.
    We also compare precondition fixes generated by \sysname with those from developers' subsequent fixes.
    For 41 cases where fixes are needed, \sysname is equivalent to or better than developers' fixes in 38 cases.
   
    
    This paper makes the following main contributions: 
    \begin{itemize}[leftmargin=*]
        \item 
        We formulate the high reliability assurance problem for DNN architecture against numerical defects and provide explicit definitions for the three-stage tasks within the problem.
    
        \item 
        We propose \sysname---the first automatic approach that solves all the three-stage tasks.
        \sysname contains novel  solutions in two major aspects: an  efficient and effective DNN static analysis framework for the first stage, and two novel techniques~(two-step test generation and abstraction optimization) to provide automatic support for the remaining two stages.
       
        
       
        
       
        
        \item We implement a  tool for \sysname, based on which we conduct a systematic evaluation on $85$ real-world DNN architectures, showing the high effectiveness and efficiency of our approach compared to both the baseline approaches and human artifacts.
       
       
       
       
    \end{itemize}
    
   
   
   
    
   
   
   
    
   
   
   
    
   
    
   
    
   
    
   
   
   
   
    
   
   
   
   
   
   
   
    
   
   
   
   
    
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
    
   
   
   
   
   
   
   
    
   
    
   
    
   
    
   
    
   
    
    
    
   
   
   
   
    
   
        
   
        
   
        
   
   
   
    
   
   

\section{Background and Task Definitions}
    \label{sec:background-task}

    \begin{figure}[t]
        \centering
        \includegraphics[width=0.9\linewidth]{figures/DNN_code_example.pdf}
        \vspace{-1.7em}
        \caption{A DL program snippet from the real-world numerical defect dataset~(Case 2a) collected by \citeauthor{yan2021exposing}~\cite{yan2021exposing}.}
        \label{fig:dnn-code-example}
        \vspace{-1.5em}
    \end{figure}
    
    \begin{figure}[!t]
        \centering
        \includegraphics[width=0.9\linewidth]{figures/DNN_graph_example.pdf}
        \vspace{-1em}
        \definecolor{figbrown}{RGB}{203,97,49}
        \definecolor{figpurple}{RGB}{112,56,153}
        \definecolor{figgreen}{RGB}{120,180,98}
        \definecolor{figyellow}{RGB}{235,163,27}
        \caption{Computational graph encoded by \Cref{fig:dnn-code-example} snippet.
       
       
       
        }
        \label{fig:dnn-graph-example}
        \vspace{-1.9em}
    \end{figure}

    In this section, we introduce the background of DNN numerical defects and failures, and formally define the three-stage tasks for reliability assurance of DNN architecture against numerical defects.

    \subsection{Background}
    
        \label{subsec:background}
    
        In Deep Learning~(DL), the developing paradigm consists of three steps: 
        (1)~Define a specific DNN architecture with code~(model construction);
        (2)~Execute the program to train the weights in the DNN architecture using training data~(model training);
        (3)~Execute the same program with trained weights and inference examples to give predictions~(model inference).
        A DNN model is a DNN architecture along with trained weights.
        
        The DNN architecture can be expressed by a computational graph.
        \Cref{fig:dnn-code-example,fig:dnn-graph-example} shows a real-world example~(a linear regression model).
        DNN architectures in DL programs can be automatically converted to DNN computational graph representation with existing tools.
        Therefore, instead of studying the program, we focus on the computational graph.
        Specifically, our discussion is based on \CodeIn{ONNX}-format computational graph~\cite{onnx} which is supported by popular DL libraries~\cite{tensorflow2015-whitepaper,paszke2019pytorch}.
        
        The computational graph can be viewed as a Directed Acyclic Graph~(DAG): $\gG = \langle \gV, \gE \rangle$, where $\gV$ and $\gE$ are sets of nodes and edges respectively.
        We call nodes with zero in-degree as \emph{initial nodes}, which correspond to input, weight, or constant nodes.
       
        Initial nodes provide concrete data for the DNN model.
        The data from each node is formatted as a \emph{tensor}, i.e., a \emph{multidimensional array}, with a specified data type and array shape annotated alongside the node definition.
        We call nodes with positive in-degree as \emph{internal nodes}, which correspond to concrete operators, such as matrix multiplication~(\CodeIn{MatMul}) and addition~(\CodeIn{Add}).
        During model training, the model weights, i.e., data from weight nodes, are generated by the training algorithm.
        Then, in deployment phase~(i.e., model inference), with these trained weights and an user-specified input called inference example,
        the output of each operator is computed in topological order. 
        The output of some specific node is used as the prediction result.
       
       
        
        We let $\vx$ and $\vw$ to denote the concatenation of data from all input nodes and data from all weight nodes respectively.\footnote{Bolded alphabet stands for vectors or tensors throughout the paper.}
        For example, in \Cref{fig:dnn-graph-example}, $\vx$ concatenates data from nodes 1 and 11; and $\vw$ concatenates data from nodes 2 and 4.
        Therefore, given specific $\vx$ and $\vw$, the input and output for each node are deterministic.\footnote{Some architectures may contain stochastic nodes. We view these nodes as weight nodes whose data are randomly sampled, so as to make the architecture deterministic.}
        We use $\fin_n(\vx; \vw)$ and $\fout_n(\vx; \vw)$ to express the input and output data of node $n$, respectively, given $\vx$ and $\vw$.
        
        \noindent
        \textbf{Numerical Defects in DNN Architecture.}
        We focus on inference-time numerical defects.
        These defects lead to numerical failures when specific operators receive inputs within invalid ranges so that the operators output \NaN or \INF, as formally defined below.
       
       
        
        
        \begin{definition}
            \label{def:numerical-defect}
            For a given computational graph  $\gG=\langle \gV, \gE\rangle$,
            if there is a node $n_0 \in \gV$, such that there exists a valid input and valid weights that can let the input of node $n_0$ be within invalid range, we say there is a \textbf{\emph{numerical defect}} at node $n_0$.
            Formally,
            \vspace{-0.5em}
            \begin{equation}
                \label{eq:numerical-defect}
                \begin{aligned}
                    & \exists \vx_0 \in \gXvalid, \vw_0 \in \gWvalid, 
                    \fin_{n_0}(\vx_0; \vw_0) \in \gI_{n_0,\mathsf{invalid}} \\
                    \Longrightarrow &
                    \text{$\exists$ numerical defect at node $n_0$}.
                \end{aligned}
            \end{equation}
        \end{definition}
        In the definition, $\gXvalid$ and $\gWvalid$ are valid input range and weight range respectively, which are clear given the deployed scenario of the given architecture.
        For example, \CodeIn{Resnet50} model~\cite{torchvision} on ImageNet~\cite{deng2009imagenet} has valid input range $\gXvalid = [0,1]^{3\times 224\times 224}$ since image pixel intensities are within $[0,1]$, and valid weight range $\gWvalid = [-1,1]^{p}$ where $p$ is the number of parameters since weights of well-trained \CodeIn{Resnet50} models are typically within $[-1,1]$.
        The invalid range $\gI_{n_0,\mathsf{invalid}}$ is determined by $n_0$'s operator type, where detailed definitions can be found in \Cref{tab:op_range_loss_form}.
        For example, for node with \CodeIn{Log} operator which computes the natural log of the input tensor, the invalid range $\gI_{n_0,\mathsf{invalid}} = (-\infty, U_{\min})$~($U_{\min}$ is the smallest positive number in input tensor type), since any input within this range would lead the node output to be \NaN or \INF, resulting in a numerical failure.
        
        \paragraph{Example.}
        The DNN architecture shown in \Cref{fig:dnn-graph-example} contains two numerical defects which exist at node 9 and 10 respectively.
        A failure-exhibiting example could be:
        node 1 contains $(-5,5)^\T$, node 2 contains  
        ${\scriptsize\begin{bmatrix} 5 & -5 \\ -5 & 5 \end{bmatrix}}$, 
        and node 4 contains $(0,0)^\T$.
        Then, the input to node 6 is $(-50, 50)^\T$.
        After the \CodeIn{softmax} operation, the output of node 6 becomes $(0,1)^\T$ in 32-bit floating number.
        Therefore, the inputs to nodes 9 and 10 are $(1,0)^\T$ and $(0,1)^\T$ respectively, both leading to failures after applying the \CodeIn{Log} operation because \CodeIn{Log} operation does not allow zero value in its input tensor.
        
        \vspace{-0.5em}
    \subsection{Task Definitions}
        \vspace{-0.3em}
    
        We now illustrate and formally define the three-stage tasks for reliability assurance against numerical defects shown in \Cref{fig:workflow}.
       
    
        \begin{discard}
        We now illustrate and formally define the four-stage tasks for testing and debugging for numerical failures.
        These tasks are shown in \Cref{fig:workflow}.
        Our \sysname provides automatic approaches for three of them: static defect detection, failure-exhibiting test generation, and precondition-fix suggestion.
        \end{discard}
        
        \subsubsection{Static Defect Detection}
            \label{subsec:def-static-detect}
            A defect detection approach should be able to flag out all nodes~(e.g., nodes 9 and 10 in \Cref{fig:dnn-graph-example}) that possibly contain numerical defects.
            Specifically, given a DNN architecture, the approach should output a superset of nodes that contain numerical defects.
            Both existing approaches such as DEBAR~\cite{zhang2020detecting} and our \sysname can complete this task.
            We will describe the \sysname approach for this task in \Cref{subsec:static-analysis} and compare with existing approaches in \Cref{subsec:static-detector-exp}.
           
            
           
           
           
           
            
           
           
           
           
           
        
        \subsubsection{Failure-Exhibiting System Test Generation}
            \label{subsec:def-system-test}
            With the possible numerical defects provided by detection approaches, we need to confirm the feasibility, i.e., whether numerical failures can be triggered in practice.
            \begin{discard}
            Moreover, fault localization techniques need failure-exhibiting tests~\cite{abreu2009practical,abreu2009spectrum,abreu2006evaluation,lou2021boosting}.
            Therefore, we need to generate them.
            We study two types of failure-exhibiting test: system test and unit test.
            \end{discard}
            We confirm such feasibility by generating failure-exhibiting system tests.
            
            \vspace{-0.5em}
            \begin{definition}
                \label{def:training-test}
                A \textbf{\emph{failure-exhibiting system test}} $\tsys$ is a tuple of valid inputs: $\tsys = \langle\vxtrain, \vxinfer\rangle$, where $\vxtrain$, $\vxinfer \in \gXvalid$, 
                such that
                for a given computational graph $\gG = \langle \gV, \gE\rangle$ with randomly initialized weights $\vw_0 \in \gWvalid$, after training with $\vxtrain$, the new weights $\vw_1$ along with inference input $\vxinfer$ trigger numerical defect at some node $n_0$, i.e., $\exists n_0 \in \gV, \fin_{n_0}(\vxinfer; \vw_1) \in \gI_{n_0, \mathsf{invalid}}$.
            \end{definition}
            \vspace{-0.5em}
            
            We call the tuple $\langle\vxtrain, \vxinfer\rangle$ system test because the test involves the whole DL life-cycle including both training and inference.
            The existence of failure-exhibiting system test confirms that the numerical defect is feasible, i.e., the numerical defect could results in numerical failure in practice.
            Our \sysname provides the first approach~(described in \Cref{subsec:system-test-gen}) for effective generation of failure-exhibiting system test to the best of our knowledge.
            
            In contrast to system tests, existing approaches such as GRIST~\cite{yan2021exposing} expose the numerical defect by generating failure-exhibiting unit test which is a tuple of model weights and inference example $\langle\vwinfer, \vxinfer\rangle$.
            We compare system test and unit test as follows:
            \begin{enumerate}[leftmargin=*]
                \item Failure-exhibiting system test $\langle\vxtrain, \vxinfer\rangle$ satisfies \\
                \CodeIn{m = train($\vxtrain$); m.infer($\vxinfer$) \textcolor{darkred}{-> trigger numerical failure}}
                
                \item Failure-exhibiting unit test $\langle\vwinfer, \vxinfer\rangle$ satisfies \\
                \CodeIn{m.weight = $\vwinfer$; \textcolor{darkgreen}{// directly set weights without training} m.infer($\vxinfer$) \textcolor{darkred}{-> trigger numerical failure}}
            \end{enumerate}
            We draw an analogy to traditional software testing:
            system test is a test for the method sequence~(\CodeIn{m = train($\vxtrain$); m.infer($\vxinfer$)}), while unit test directly injects object state to test the inference method~(inject model weights to \CodeIn{m}, then test \CodeIn{m.infer($\vxinfer$)}).
            As a result, the unit test cannot confirm defect feasibility, since  model weights in the unit test may not be accessible via model training.
            
            
            \begin{discard}
            \begin{definition}[failure-exhibiting unit test]
                \label{def:inference-test}
                A failure-exhibiting unit test $\tunit$ is a tuple of valid inputs and weights: $\tunit = \langle\vxinfer, \vwinfer\rangle$, where $\vxinfer \in \gXvalid$, $\vwinfer \in \gWvalid$, such that weights $\vwinfer$ along with inference input $\vxinfer$ trigger numerical defect at some node $n_0$, i.e., $\exists n_0 \in \gV, \fin_{n_0}(\vxinfer; \vwinfer) \in \gI_{n_0, \mathsf{invalid}}$.
            \end{definition}
            
            The unit test only involves the inference phase, which we view as the unit under test.
            The existence of failure-exhibiting unit test is not sufficient to confirm feasibility, since model weights are generated by training process and the model weights in the unit test may not be practical, i.e., may not be achievable via training.
            Although existence of unit test may not confirm feasibility, the failure-exhibiting unit test is helpful for fault localization.
            Existing fault localization approaches, e.g., coverage-based approaches~\cite{abreu2009practical,abreu2009spectrum,abreu2006evaluation,lou2021boosting} , need many distinct failure-exhibiting unit tests.
           
           
           
           
           
            We draw an analogy to traditional software testing:
            system test is similar to the test for the method sequence~(\CodeIn{m = train($\vxtrain$); m.infer($\vxinfer$);}), while unit test is similar to directly injecting object state to test the inference method~(inject model weights $\vwinfer$ to \CodeIn{m}, then run \CodeIn{m.infer($\vxinfer$)}).
            Therefore, system test confirms feasibility, while unit test is more targeted for the unit under test and therefore helpful for fault localization.
           
           
            
            Our \sysname is the first one in the literature that effectively generates failure-exhibiting system test and \emph{many distinct} failure-exhibiting unit test.
            \end{discard}
        
        \begin{discard}
        \subsubsection{Fault Localization}
        
            With failure-exhibiting tests, the next step is fault localization.
            Concretely, we need to locate the concrete operators in the DNN computational graph to impose the fix.
            In other words, the output of this stage is a set of nodes that specify the location to impose the fix.
            In this paper, we do not focus on this task.
            In practice, we can leverage manual inspection or existing fault localization tools, e.g., \cite{wardat2021deeplocalize}, in this stage.
           
           
        \end{discard}
        
        \subsubsection{Precondition-Fix Suggestion}
            \label{subsec:def-precond-fix}
            \begin{discard}
            After fault localization, we obtain a set of nodes $\Vfix \subseteq \gV$ to impose the fix.
            According to the empirical study of developers' fix in \Cref{sec:intro}, 12 out of 13 defects are fixed by imposing interval~(also called box~\cite{gehr2018ai2}) preconditions clipping the inputs of some nodes.
            Therefore, we consider interval precondition fix in this paper which is formally defined as below.
            \end{discard}
            Based on user specifications, we have a set of nodes $\Vfix \subseteq \gV$ to impose the fix.
            For example, if the fix method is clipping the inference input, $\Vfix$ are input nodes~(e.g., nodes 1, 11 in \Cref{fig:dnn-graph-example});
            if the fix method is clipping the model weights during training, $\Vfix$ are weight nodes~(e.g., nodes 2, 4 in \Cref{fig:dnn-graph-example});
            if the fix method is clipping before the defect node \CodeIn{n}, such as calling \CodeIn{torch.clamp($n$.input, l, u)}, $\Vfix$ are defect nodes~(e.g., nodes 9, 10 in \Cref{fig:dnn-graph-example}).
            
            According to the empirical study of developers' fix in \Cref{sec:intro}, 12 out of 13 defects are fixed by imposing interval preconditions clipping the inputs of $\Vfix$.
            Therefore, we consider interval precondition fix in this paper which is formally defined as below.  
        
            \vspace{-0.5em}
            \begin{definition}
                \label{def:fix-precond}
                For given computational graph $\gG = \langle \gV, \gE\rangle$, fix locations $\Vfix \subseteq \gV$, and a given node $n_0 \in \gV$ with numerical defect, a set of tuples $\{ \langle n_i, \vl_i, \vu_i \rangle: n_i \in \Vfix, \vl_i \in \sR^{\text{in\_shape}(n_i)}, \vu_i \in \sR^{\text{out\_shape}(n_i)} \}_{i=1}^m$ is called a \textbf{\emph{precondition fix}}, if 
               
                \vspace{-0.5em}
                \begin{equation}
                    \label{eq:fix-precond}
                    \begin{aligned}
                        & \forall \vx \in \gXvalid, \vw \in \gWvalid, \\
                        & \bigwedge_{1 \le i \le m}
                        \left(\vl_i \le \fin_{n_i}(\vx; \vw) \le \vu_i\right)
                        \Longrightarrow
                        \fin_{n_0}(\vx; \vw) \notin \gI_{n_0, \mathsf{invalid}}.
                    \end{aligned}
                \end{equation}
                In \Cref{eq:fix-precond}, for same-shape tensors $\va$ and $\vb$, $\va \le \vb$ stands for elementwise ``smaller or equal to'', i.e., $\bigwedge_j (\va_j \le \vb_j)$.
               
               
            \end{definition}
            
           
                Given precondition fix $\bigwedge_{1 \le i \le m} (\vl_i \le \fin_{n_i}(\vx; \vw) \le \vu_i)$, the input of node $n_0$ cannot be invalid and according to \Cref{def:numerical-defect} there is no numerical defect at node $n_0$.
               
               
                As described in \Cref{subsec:precond-gen}, given fix locations $\Vfix$ and DNN architecture $\gG$, our \sysname effectively suggest precondition fix.
                
                \begin{discard}
                Different fix locations $\Vfix$ correspond to different fix methods.
                If fix locations are input nodes~(nodes 1, 8 in \Cref{fig:dnn-graph-example}), we preprocess the input with the precondition as the fix;
                if locations are weight nodes~(nodes 2, 3, 7 in \Cref{fig:dnn-graph-example}), we clip the model weights with the precondition after training;
                if locations are some internal nodes~(e.g., node 6 in \Cref{fig:dnn-graph-example}), we modify the architecture to incorporate the clipping operation.
                Our \sysname can effectively suggest precondition fix given specified locations $\Vfix$ and thus applicable for different fix methods.
                \end{discard}
                
               
               
           
    

\section{\sysname}

    \label{sec:method}
    
    In this section, we introduce the \sysname approach in detail.
    We first give an approach overview, then introduce each component.

   
        
       
        
        \begin{figure}[!t]
           
            \centering
            \includegraphics[width=0.7\linewidth]{figures/overview_new.pdf}
            \vspace{-1.2em}
            \caption{Overview of the \sysname approach.
            }
            \label{fig:overview}
            \vspace{-1.85em}
        \end{figure}
        
        In \Cref{fig:overview}, we show the overview structure of the \sysname approach.
        \sysname takes a DNN architecture as the input.
        First, the DNN static analysis framework in \sysname detects all possible numerical defects in the architecture, which is the stage \textcircled{1} in \Cref{fig:workflow}.
        Second, these defects are passed to our two-step system test generalization component including the gradient back-propagation and extended DLG attack techniques to complete the stage \textcircled{2} in \Cref{fig:workflow}.
       
       
       
        Third, the input/output abstractions for all nodes produced by the DNN static analysis framework are supplied to the abstract optimization component.
        Given the fix locations, this component produces valid precondition fixes to fix all numerical defects within the DNN architecture, which completes the stage \textcircled{3} in \Cref{fig:workflow}.
        \begin{discard}
        The static analysis framework in \sysname processes the architecture and generate tight, sound, and differentiable interval abstractions for each node of the architecture graph.
        With these abstractions, we can directly flag out possible numerical defects and output defect detection result, which corresponds to stage \textcircled{1} in \Cref{fig:workflow}.
        Built upon the static analysis framework, we propose a bilevel optimization technique to search for failure-exhibiting interval domains, and then do random sampling within these domains to generate many distinct failure-exhibiting unit tests, which corresponds to stage \textcircled{2}(2) in \Cref{fig:workflow}.
        With a generated failure-exhibiting unit test, we propose extended DLG attack~\cite{zhu2020deep} to generate failure-exhibiting system tests, which corresponds to stage \textcircled{2}(1) in \Cref{fig:workflow}.
        Lastly, for given fix locations~(which can be given by human or fault localization tools), we combine the abstractions from static analysis framework with our abstraction optimization technique to generate precondition fixes, which corresponds to stage \textcircled{4} in \Cref{fig:workflow}. 
        \end{discard}
        
        We will introduce the static analysis framework in \sysname in \Cref{subsec:static-analysis}, 
        the technique of two-step system test generalization in \Cref{subsec:system-test-gen}, 
        \begin{discard}
        the technique of bilevel optimization for unit test generation in \Cref{subsec:unit-test-gen}, 
       
       
        \end{discard}
        and the technique of abstraction optimization for precondition-fix suggestion in \Cref{subsec:precond-gen}.

    \subsection{DNN Static Analysis Framework for Defect Detection}
    
        \label{subsec:static-analysis}
        
        
        \sysname contains a static analysis framework to enable defect detection and support other tasks as shown in \Cref{fig:overview}.
        Given a DNN architecture and the ranges for input and weight nodes, the static analysis framework computes interval abstractions for all possible inputs and outputs of each node. 
        As a result, we can check whether there is overlap between the interval abstraction and invalid input ranges for all nodes in the DNN architecture to statically detect numerical defects.
        The defect nodes are fed into two-step system test generation component to confirm feasibility~(see \Cref{subsec:system-test-gen}); 
        and the differentiable abstractions are fed into abstract optimization component to produce precondition fixes~(see \Cref{subsec:precond-gen}).
        
        The interval abstraction~(also called ``box domain''~\cite{gehr2018ai2}) is a standard technique for analyzing traditional software~\cite{cousot1977static,blanchet2003static} especially for those with numerical components~\cite{gurfinkel2015seahorn,singh2017fast}.
        In our case, formally, for given possible sets of inference input and model weights, namely $\gX$ and $\gW$,
        for each node $n\in \gV$, the framework computes interval abstraction $[\vl_n, \vu_n] := \{\vx: \vl_n \le \vx \le \vu_n\}$ such that $[\vl_n, \vu_n]$ captures all possible inputs~(or outputs) of the node: $[\vl_n, \vu_n] \supseteq \fin_n(\gX; \gW)$~(or $[\vl_n, \vu_n] \supseteq \fout_n(\gX; \gW)$).\footnote{We define $\fin_n(\gX; \gW) := \{\fin_n(\vx, \vw): \vx\in \gX, \vw\in \gW\}$, which is $\fin_n$'s image when the domain is $\gX \times \gW$. Similarly, $\fout_n(\gX; \gW) := \{\fout_n(\vx, \vw): \vx\in \gX, \vw\in \gW\}$.}
        
        
        \begin{figure}[!t]
            \includegraphics[width=0.9\linewidth]{figures/partitioning_example.pdf}
            \vspace{-1.3em}
            \caption{Example of tensor partitioning.
            Tensor partitioning reduces the size of tensors in the abstract domain by sharing one interval bound among all elements inside one subblock.}
            \label{fig:partitioning}
            \vspace{-1.8em}
        \end{figure}
        
        \textbf{Abstract Domain: Interval with Tensor Partitioning.}
        Instead of directly using the interval domain as the abstraction domain, we use the interval with tensor partitioning~\cite{zhang2020detecting} as the abstraction domain, which partitions the tensor into several subblocks and shares the interval abstractions at the block level instead of imposing at the element level, as demonstrated in \Cref{fig:partitioning}.
        Therefore, we can compute the abstraction of a smaller size than the original tensor to improve efficiency.
        Such block-level interval sharing is made possible by the following two observations:
        (1)~The input range of neighboring elements inside one tensor is usually the same in practice.
        For example, for image input, the valid values of all pixels~(i.e., all elements) are the same, such as $[0,1]$.
        (2)~When imposing precondition fixes, human developers usually impose the same intervals for a block of elements.
       
        
        \begin{discard}
        \yuhao{I think the following paragraph can go to appendix. It seems unrelated to the main paper.}
        \linyi{moved}
        Following the notation in abstract interpretation literature~\cite{cousot1977static}, suppose the tensor has $m$ dimensions with shape $n_1 \times n_2 \times \cdots \times n_m$, we define the abstract domain as such:
        \begin{equation}
            \begin{aligned}
                & \sA := \{ \langle \vl, \vu, (S_i)_{i=1}^m \rangle: \vl, \vu \in \sR^{n'_1 \times \cdots \times n'_m} \\
                & \hspace{5em} |S_i| = n'_i, \forall s\in S_i, s\in \sN_+, 0 \le s < n_i \}.
            \end{aligned}
            \label{eq:abs-domain}
        \end{equation}
       
        Each element in $\sA$ is a triplet $a = \langle \vl, \vu, (S_i)_{i=1}^m\rangle$, where the first two elements are subblock-wise interval lower bound and upper bound respectively, and the last element $(S_1, S_2, \dots, S_m)$ contains $m$ sets, where each set $S_i$ corresponds to the \emph{0-indexed} split indices for the $i$-th dimension to form subblocks.
        We use $S_i[j] \in \sN_+$ to represent the $j$-th element of split index set $S_i$ sorted in ascending order and define $S_i[|S_i|] = n_i$.
       
        \end{discard}
        
        \textbf{Abstraction Techniques for Sound, Tight, and Differentiable Abstractions.}
        In this framework, we aim to generate abstractions that are sound, tight, and differentiable.
        The soundness requirement guarantees the defect detection has no false negatives~(i.e., flag out all possible defects) and the validness of generated preconditions-fixes.
        The tightness requirement can reduce the false positives for detecting numerical defects and increase the quality of the generated precondition fix~(by generating wider failure-free intervals). 
        The differentiability requirement is a requisite for our technique of abstraction optimization for precondition-fix generation.
        \begin{discard}
        The soundness requirement guarantees the defect detection has no false negatives, i.e., flag out all possible defects, and the failure-exhibiting and fix properties of generated unit tests and preconditions, respectively.
        The tightness requirement can reduce the false positives for detecting numerical defects, increase the diversity of generated failure-exhibiting unit tests~(by generating wider failure-exhibiting intervals), and increase the quality of the generated precondition fix~(by generating wider failure-free intervals). 
       
        The differentiability requirement is a requisite for our generation techniques to be introduced later.
        \end{discard}
        
        To achieve these goals, we propose two \emph{novel} abstraction techniques:
        (1)~\emph{Backward fine-grained node labelling}:
        Some nodes require fine-grain abstraction, i.e., the finest granularity tensor partitions.
        For example, the \CodeIn{Gather} node takes indices as input to pick elements from another input, and the indices should use abstraction at the element level rather than block level to achieve high abstraction precision.
        We propose to propagate ``fine-grain'' labels from these fine-grain requiring nodes to initial nodes by topologically sorting in the graph with inverted edges, and then apply fine-grained abstractions on labelled initial nodes to satisfy these fine-grain requirements.
        (2)~\emph{Dynamic partitioning}:
        When constructing the abstractions for internal nodes, we propose to construct output abstraction's tensor partitions dynamically based on input partitions to preserve the precision.
        Specifically, we decide the output abstraction's partitions to satisfy both tightness~(bears no precision loss compared to elementwise interval abstraction) and minimality~(using the minimum number of partitions).
        
        With these two novel abstraction techniques, the whole static analysis process consists of three steps: 
        (1)~We determine the tensor partition granularity of all initial nodes in the architecture graph by aforementioned technique of backward fine-grained node labelling.
        (2)~We sort all nodes in the graph in topological order.
        (3)~We apply corresponding abstraction computation algorithms iteratively for all internal nodes in the graph, depending on the node operator types, until the abstractions for all nodes in the graph have been computed.
        Now, we can compare the abstraction intervals with the invalid ranges for each node in the graph, to find out all possible numerical defects, and pass the defect list and abstraction list for system test generation component and precondition-fix suggestion component respectively.
        Due to the page limit, we omit the details of our abstraction techniques, formal definitions, proofs, and discussion with existing approaches to \Cref{adxsec:static-analysis}.
        
    \subsection{Two-Step System Test Generation}
    
        \label{subsec:system-test-gen}
    
        In this section, we generate failure-exhibiting system test for the given DNN to confirm the feasibility of numerical defects.
        Here, we take the DNN architecture as the input.
        From the static analysis framework, we obtain a list of nodes that have potential numerical defects.
        For each node $n_0$ within the list, we apply our technique of two-step system test generation to produce a failure-exhibiting system test $\tunit = \langle \vxtrain, \vxinfer \rangle$ as the output.
        According to \Cref{subsec:def-system-test}, the system test should satisfy that after the architecture is trained with $\vxtrain$, entering $\vxinfer$ would result in numerical failure.
        
        We propose the novel technique of two-step system test generalization:
        First, generate failure-exhibiting unit test $\langle \vwinfer, \vxinfer \rangle$;
        Then, generate training example $\vxtrain$ that leads model weights to be close to $\vwinfer$ after training.
       
        
        \subsubsection{\textbf{Step 1}: Unit Test Generation with Gradient Back-Propagation}
            For this step, we adapt the gradient back-propagation  approach, GRIST~\cite{yan2021exposing}, and further automate it to generate failure-exhibiting unit test.
           
           
            The original GRIST implementation requires the user to manually locate the defect operator and edit the program code to inject gradient-back-propagation-based generalization component.
            To fully automate the process, we reimplement this approach by directly parsing the defect node list from our static analysis framework and imposing the corresponding generation loss functions based on defect operator type.
            The concrete loss expressions are listed in \Cref{tab:op_range_loss_form} under the ``$\gL_{n_0}^{\testgen}$'' column, which satisfies that $\gL_{n_0}^{\testgen}(\vx_{n_0}) < 0$ suffices to trigger numerical failure.
            We minimize the loss with gradient back-propagation following \cite{yan2021exposing}, and extract the unit test $\langle \vwinfer, \vxinfer \rangle$ when the loss value becomes negative.
        
        \subsubsection{\textbf{Step 2}: Training Example Generation with Extended DLG Attack}
            For the second step, we have these as inputs: (1)~the DNN architecture; (2)~the failure-exhibiting unit test $\tunit = \langle \vwinfer, \vxinfer \rangle$ and (3) the randomly initialized weights $\vw_0 \in \gWvalid$.
            Our goal is to generate training example $\vxtrain$, such that the model trained with $\vxtrain$ will contain weights close to $\vwinfer$.
           
            In other words, \emph{given the DNN architecture, initial weights $\vw_0$, and target weights $\vwinfer$, we need to generate $\vxtrain\in \gXvalid$, such that after training with $\vxtrain$, the model weights are updated from $\vw_0$ to (approximately) $\vwinfer$.}
            
            
            DNNs are typically trained with gradient-descent-based algorithms such as stochastic gradient descent~(SGD)~\cite{robbins1951stochastic}.
            In SGD, in each step $t$, we sample a mini-batch of samples from the training dataset to compute their gradients on model weights and use them to update the weights.
            We focus on one-step SGD training with a single training example, since generating a single one-step training example to exhibit failure is more desirable for ease of debugging as discussed in \Cref{sec:intro}.
            In this case, after training, the model weights
           
            \begin{equation}
                \vwinfer = \vw_0 - \gamma \nabla_\vw \gL(\vxtrain; \vw_0),
               
            \end{equation}
            where $\gamma \in \sR_+$ is a predefined learning rate, and $\gL = \fout_{n_l}$ is the loss function, which is the output of some loss node $n_l$ in the DNN architecture.
            Hence, our goal becomes finding $\vxtrain$ that satisfies 
            \vspace{-0.5em}
            \begin{equation}
                \nabla_\vw \gL(\vxtrain; \vw_0) = \dfrac{\vw_0 - \vwinfer}{\gamma}. \label{eq:dlg-attack}
            \end{equation}
            \vspace{-1em}
            
            
            \begin{discard}
            , i.e., a subset of training samples $\gS_t \subseteq \gS$ whose cardinality is $B$, and compute the weight gradients with respect to the loss function which is a specific node $n_l$.
            Then, the model weights are uniformly updated towards the negative direction of the gradients to reduce the loss value.
            Formally, let $\gL = \fout_{n_l}$ be the output function of loss node $n_l$, and let $\gamma \in \sR_+$ be the learning rate, at each step we have
           
            \begin{equation}
               
                \vw_{t+1} = \vw_t - \gamma \cdot \dfrac{1}{B} \sum_{\vx\in \gS_t} \nabla_{\vw} \gL(\vx; \vw_t).
               
            \end{equation}
            \yuhao{Instead of spend space describing widely used SGD. Maybe we should briefly mention SGD and focus on our single-step-and-single-training-example SGD.}
            The SGD can terminate at an arbitrary step as long as the user thinks the model is trained well.
            For simplicity, we consider one-step training, i.e., weights after training $\vwinfer = \vw_1$.
            Extension to multi-step training may further improve the ability to generate failure-exhibiting system tests and we leave as future work. 
            
            We focus on the widely-used SGD training, and generate one valid training example $\vxtrain \in \gXvalid$.
            In this case, any mini-batch contains multiple copies of this training example. 
            \yuhao{I think in case $B=1$ and $|\gXvalid|=1$. Maybe we can add the discussion of single training example in the previous paragraph.}
            Hence, we desire
            \begin{small}
            \begin{equation}
                \vwinfer = \vw_0 - \gamma \nabla_\vw \gL(\vxtrain; \vw_0) 
                \iff
                \nabla_\vw \gL(\vxtrain; \vw_0) = \dfrac{\vw_0 - \vwinfer}{\gamma}. \label{eq:dlg-attack}
            \end{equation}
            \end{small}%
            Our goal is to find $\vxtrain$ that satisfies \Cref{eq:dlg-attack}.
            \end{discard}
            
            The DLG attack~(deep leakage from gradient, \cite{zhu2020deep}) is a method for generating input data that induce specific weight gradients.
            The attack is originally designed and shown effective for recovering training samples from monitored gradient updates.
           
            Since the right hand side~(RHS) of \Cref{eq:dlg-attack} is known, our goal here is also generating input example $\vxtrain$ that induces specific weight gradients~(the RHS of \Cref{eq:dlg-attack}).
            Therefore, we leverage DLG attack to generate training example $\vxtrain$.
           
            
            \paragraph{Extending DLG Attack with Straight-Through Estimator.}
            Directly using DLG attack suffers from an optimization challenge in our scenario, so we propose to leverage the straight-through estimator~\cite{bengio2013estimating} to extend DLG attack to mitigate this issue.
            
            Specifically, in DLG attack, suppose the target weight gradients are $\vdwtarg$, we use gradient descent over the squared error $\| \nabla_\vw \gL(\vx; \vw_0) - \vdwtarg \|_2^2$ to generate $\vx$.
            In this process, we need meaningful gradient information of this squared error loss to perform the optimization.
            However, the gradient of this loss involves second-order derivatives of $\gL(\vx;\vw_0)$, which is zero for many DNNs (e.g., DNNs with $\relu$ as activation function are piecewise linear and have zero second-order derivatives almost everywhere~\cite{serra2018bounding}).
            This optimization challenge is partly addressed in DLG attack by replacing \CodeIn{ReLU} with \CodeIn{Sigmoid}, but it changes the DNN architecture~(i.e., the system under test) and is not suitable for our setting.
            
            Therefore, we leverage the straight-through estimator to mitigate the optimization challenge.
            At a high level, for some operators, such as \CodeIn{ReLU}, we do not change its forward computation, but change its backward gradient to provide second-order derivatives.
            For example, for \CodeIn{ReLU}, in backward computation we use the gradient of $\mathrm{Softplus}$ function, namely $1 - \frac{1}{1+\exp(\vx)}$, as the gradient.
            The $\mathrm{Softplus}$ function is a smoothed approximation of $\relu$ function~\cite{glorot2011deep}.
            We leverage its gradient instead of its value to avoid changing the actual behavior of system under test.
            
            Now we are ready to sketch the whole process of generating failure-exhibiting system test.
            We obtain the list of defect nodes from static analysis.
            For each node, we first generate failure-exhibiting unit test with gradient back-propagation, then extract $\vwinfer$ from the unit test and apply our extended DLG attack to generate corresponding training example.
            In each optimization step of DLG, we clip the current $\vxtrain$ into the valid region $\gXvalid$, and then check whether this $\vxtrain$ can exhibit numerical failures.
            If so, we terminate and output $\tsys = \langle \vxtrain, \vxinfer \rangle$ as the generated failure-exhibiting system test.
            If not, we continue the optimization process.
           
        
        \begin{table}[!t]
            \centering
            \caption{List of representative operators that may contain numerical defects, including their invalid input ranges and the losses in \sysname for failure-exhibiting system test generation and for precondition-fix generation.
            $U_{\min}$ and $U_{\max}$ are minimum and maximum positive floating numbers of node's data type.
            $n_0$ is the node, $x_{n_0}$ is the node's input, and $[l_{n_0}, u_{n_0}]$ is the interval abstraction of node's input.
           
           
            }
            \label{tab:op_range_loss_form}
            \vspace{-1em}
            \resizebox{\linewidth}{!}{
            \begin{tabular}{c|c|c|c}
                \toprule
                \multirow{2}{*}{Op.} & Invalid Input & System Test Generation  & Precondition-Fix Generation  \\
                & Range $\gI_{n_0,\mathsf{invalid}}$ & Step 1 Loss $\gL_{n_0}^{\testgen}$ & Loss $\gL_{n_0}^{\precondgen}$ \\
                \hline\hline
                \CodeIn{Log} & $(-\infty, U_{\min})$ & \multirow{2}{*}{
               
               
               
               
               
               
                $x_{n_0} - U_{\min}$ $^a$
                } & \multirow{2}{*}{$-l_{n_0} + U_{\min}$} \\
                \cline{1-2}
                \CodeIn{Sqrt} & $(-\infty, U_{\min})$ & \\
                \hline
                \CodeIn{Exp} & $[\ln U_{\max} ,\infty)$ & 
                $-x_{n_0} + \ln(U_{\max})$ & $u_{n_0} - \ln(U_{\max})$ \\
                \hline
                \CodeIn{Div} & $(-U_{\min}, U_{\min})$ & $x_{n_0}^2 - U_{\min}^2$ & $(-l_{n_0} + U_{\min})$ \text{ or } $(u_{n_0} + U_{\min})$ $^b$ \\
                \bottomrule
            \end{tabular}
            }
            \begin{flushleft}
                \footnotesize
                $^a$ In implementation, we deploy a more complicated loss to avoid gradient vanishing.\\
                $^b$ Either form can be used, which is a tunable parameter.
                \begin{discard}
                $^b$ Being $< 0$ does not guarantee ``within invalid input range'' for this loss. But minimizing this loss tends to make the guarantee holds empirically, so we adapt this loss as a surrogate.\\
                $^c$ Either form can be used, which is a tunable parameter.
                \end{discard}
            \end{flushleft}
            \vspace{-2em}
        \end{table}
        
    \begin{discard}
    \subsection{Failure-Exhibiting System Test Generation with Extended DLG Attack}
    
        \label{subsec:system-test-gen}
        
        In this section, we aim to generate failure-exhibiting system test for given DNN architecture to confirm the feasibility of numerical defects.
        Here, we take these as inputs: (1)~the DNN architecture; (2)~the failure-exhibiting unit test $\tunit = \langle \vxinfer, \vwinfer \rangle$ (generation process will be described in \Cref{subsec:unit-test-gen}) and (3) the randomly initialized weights $\vw_0 \in \gWvalid$.
        We now describe our extended DLG attack for generating failure-exhibiting system test $\tsys = \langle \vxtrain, \vxinfer \rangle$.
        The definitions of these notations are in \Cref{def:inference-test,def:training-test}.
        
        Since $\vxinfer$ is already provided in inputs as part of $\tunit$, we only need to generate $\vxtrain$.
        Therefore, we reformulate our task as such: \emph{given the DNN architecture, initial weights $\vw_0 \in \gWvalid$, and target weights $\vwinfer \in \gWvalid$, we need to generate $\vxtrain\in \gXvalid$, such that after training with $\vxtrain$, the model weights are updated from $\vw_0$ to $\vwinfer$.}
        
        The DNNs are typically trained with gradient descent based algorithms such as stochastic gradient descent~(SGD)~\cite{robbins1951stochastic} despite their various architectures.
        In SGD, let $\gS$ be the training dataset, in each step $t$, we sample a mini-batch, i.e., a subset of training samples $\gS_t \subseteq \gS$ whose cardinality is $B$, and compute the weight gradients with respect to the loss function which is a specific node in the DNN architecture graph.
        Then, the model weights are uniformly updated towards the negative direction of the gradients to reduce the loss value.
        Formally, let $\gL = \fout_{n_l}$ be the output function of loss node $n_l$, and let $\gamma \in \sR_+$ be the learning rate, at each step we have
        \begin{equation}
            \small
            \vw_{t+1} = \vw_t - \gamma \cdot \dfrac{1}{B} \sum_{\vx\in \gS_t} \nabla_{\vw} \gL(\vx; \vw_t).
        \end{equation}
        The SGD can terminate at an arbitrary step as long as the user thinks the model is trained well.
        For simplicity, we consider one-step training, i.e., weights after training $\vwinfer = \vw_1$.
        Extension to multi-step training may further improve the ability to generate failure-exhibiting system tests and we leave as future work. 
        
        We focus on the widely-used SGD training, and generate a single valid training example $\vxtrain \in \gXvalid$.
        In this case, any mini-batch contains multiple copies of this training example.
        Thus,
        \begin{small}
            \begin{align}
                & \vwinfer = \vw_0 - \gamma \nabla_\vw \gL(\vxtrain; \vw_0) \\
                \iff &
                \nabla_\vw \gL(\vxtrain; \vw_0) = \dfrac{\vw_0 - \vwinfer}{\gamma}. \label{eq:dlg-attack}
            \end{align}
        \end{small}
        Our goal is to find $\vxtrain$ that satisfies \Cref{eq:dlg-attack}.
        
        The DLG attack~(deep leakage from gradient, \cite{zhu2020deep}) is a method for generating input data that induce specific weight gradients.
        The attack is effective for stealing the training samples, since the weight gradients, which are weight updates, are easily accessible in distributed training, and the input that induces given weight gradients is highly similar to the actual training example.
        Since the right hand side~(RHS) of \Cref{eq:dlg-attack} is known, our goal here is also generating input example $\vxtrain$ that induces specific weight gradients~(the RHS of \Cref{eq:dlg-attack}).
        Therefore, we leverage DLG attack to generate training example $\vxtrain$.
        
        \paragraph{Extending DLG Attack with Straight-Through Estimator.}
        Directly using DLG attack suffers from an optimization challenge in our scenario, so we propose to leverage the straight-through estimator~\cite{bengio2013estimating} to extend DLG attack to mitigate this issue.
        
        Concretely, in DLG attack, suppose the target weight gradients are $\vdwtarg$, we use gradient descent over the squared error $\| \nabla_\vw \gL(\vx; \vw_0) - \vdwtarg \|_2^2$ to generate $\vx$.
        In this process, we need meaningful gradient information of this squared error loss to execute the optimization process.
        However, the gradient of this loss involves second-order derivatives of $\gL(\vx;\vw_0)$, which is zero for many DNNs~(e.g., DNNs with $\relu$ as activation function are piecewise linear and have zero second-order derivatives almost everywhere~\cite{serra2018bounding}).
        This optimization challenge is partly addressed in DLG attack by replacing \CodeIn{ReLU} with \CodeIn{Sigmoid}, but it changes the DNN architecture and is not suitable for our setting.
        
        Therefore, we leverage the straight-through estimator to mitigate the optimization challenge.
        At a high level, for some operators, such as \CodeIn{ReLU}, we do not change its forward computation, but change its backward gradient computation process to provide second-order derivatives.
        For example, for \CodeIn{ReLU}, in forward computation we still compute $\relu(\vx)$, but in backward computation we use the gradient of $\mathrm{Softplus}$ function, namely $1 - \frac{1}{1+\exp(\vx)}$, as the gradient.
        The $\mathrm{Softplus}$ function is a smoothed approximation of $\relu$ function~\cite{glorot2011deep}.
        We leverage its gradient instead of its value to avoid changing the actual behavior of system under test.
        
        Now we are ready to sketch the whole process of generating failure-exhibiting system test.
        We start with one failure-exhibiting unit test $\langle \vxinfer, \vwinfer \rangle$, the DNN architecture, and the initial weights $\vw_0$.
        Then, we compute the target weight gradients which is the RHS of \Cref{eq:dlg-attack}.
        Next, we leverage the DLG attack to minimize the squared error with gradient descent, where in gradient computation we use the straight-through estimator for some operators.
        After each optimization step, we clip the current $\vx$ into the valid region $\gXvalid$, and then check whether this $\vx$ can exhibit numerical failures.
        If so, we terminate and output $\tsys = \langle \vxtrain, \vxinfer \rangle$ as the generated failure-exhibiting system test.
        If not, we continue the optimization process until reaching the time limit. 
       
       

    \subsection{Failure-Exhibiting Unit Test Generation with Bilevel Optimization}
    
        \label{subsec:unit-test-gen}
        
        In this section, we aim to generate many distinct failure-exhibiting unit tests for the downstream fault localization task.
        Here, we take the DNN architecture as the input.
        As a preprocessing step, we statically detect the defects in the architecture and specify the node $n_0$ with potential defects for which we would like to trigger numerical failure.
        After preprocessing, we generate failure-exhibiting interval domains for weights and input nodes, and then sample from these intervals to generate many distinct failure-exhibiting unit tests.
        
        We now focus on the task of generating failure-exhibiting interval domains.
       
       
        We formally define this task as an optimization problem.
        We denote the intervals of weight and input nodes by $[\vl_w, \vu_w]$ and $[\vl_x, \vu_x]$ respectively where $\vl_w, \vu_w, \vl_x, \vu_x$ are optimizable variables, and denote the valid ranges of weight and input nodes by $[\vl_w^{\valid}, \vu_w^{\valid}]$ and $[\vl_x^{\valid}, \vu_x^{\valid}]$.
        Then, we define the optimization as below.
       
            \begin{align}
                \allowdisplaybreaks
                & \maximize_{\vl_w, \vu_w, \vl_x, \vu_x} \quad  s 
                \label{eq:test-gen-obj} \\
                \mathrm{s.t.}\quad  & \vu_w = \vl_w + s(\vu_w^{\valid} - \vl_w^{\valid}),
                \vu_x = \vl_x + s(\vu_x^{\valid} - \vl_x^{\valid}),
                \label{eq:test-gen-span-cons} \\
                & \vl_w^{\valid} \le \vl_w \le \vu_w \le \vu_w^{\valid},
                \vl_x^{\valid} \le \vl_x \le \vu_x \le \vu_x^{\valid}, 
                \label{eq:test-gen-valid-cons} \\
                & \gL_{n_0}^{\testgen}(\vl_w, \vu_w, \vl_x, \vu_x) < 0.
                \label{eq:test-gen-effective-cons} 
            \end{align}
       
       
        In this optimization problem, the objective is to maximize $s$, which is the relative span of both weight and input nodes with respect to their valid ranges.
        The larger the span is, the larger range from which we can sample the failure-exhibiting unit tests, and therefore the more diverse the generated unit tests are.
        \Cref{eq:test-gen-span-cons} enforces the span constraints, \Cref{eq:test-gen-valid-cons} ensures that all generated unit tests are valid.
        Then, the last constraint~(\Cref{eq:test-gen-effective-cons}) is more involved.
        In this constraint, we require the node-specific test generation loss $\gL_{n_0}^{\testgen}$ to be negative.
       
       
       
        Loss $\gL_{n_0}^{\testgen}$ is defined over the interval abstraction of the defect node $n_0$ so that if the loss value is negative, the input interval of the node is within the invalid input range.
        This property guarantees that all values from the weight and input intervals can make the node input invalid thus triggering the numerical failure as long as the loss value is negative.
        For example, for both \CodeIn{Log} and \CodeIn{Sqrt} operators, when the loss, i.e., $\vu_{n_0} - U_{\min}$, is negative, then $\vu_{n_0} < U_{\min}$ and hence any input to this operator would be smaller than $U_{\min}$ and induce numerical failure.
        The concrete loss form is decided by the node's operator type. Examples are listed in the third column of \Cref{tab:op_range_loss_form}.
        
        As shown by \Cref{tab:op_range_loss_form}, the loss is defined on the interval abstraction of the defect node $n_0$ instead of weight or input nodes.
        Therefore, we need to compute the interval abstraction of $n_0$ from the intervals of weights and inputs.
        We resort to the static analysis framework introduced in \Cref{subsec:static-analysis} to compute these intervals.
        We use the finest granularity in the analysis framework to precisely reflect the intervals of weight and input nodes.
        
        For a concrete set of $\vl_w$, $\vu_w$, $\vl_x$, and $\vu_x$, we can compute the loss value, and now we need to optimize these variables.
        The optimization problem is challenging to solve due to these complex constraints and complex computation of $n_0$'s interval bounds.
        Therefore, instead of precisely solving the problem, we propose the technique of \textbf{bilevel optimization} to approximately but effectively solve it.
        Based on the property that the feasible objective range of the optimization problem is always a continuous interval starting from $0$, our bilevel optimization grows the span $s$ iteratively and exponentially starting from $s=0$.
        For the current $s$, we only need to optimize $\vl_w$ and $\vl_x$ and then $\vu_w$ and $\vu_x$ are determined by eliminating constraint in \Cref{eq:test-gen-span-cons}.
        Now, we optimize $\vl_w$ and $\vl_x$ for minimizing $\gL_{n_0}^{\testgen}$ subject to the box constraint $\vl_w^{\valid} \le \vl_w \le \vu_w^{\valid} - s(\vu_w^{\valid} - \vl_w^{\valid})$ and $\vl_x^{\valid} \le \vl_x \le \vu_x^{\valid} - s(\vu_x^{\valid} - \vl_x^{\valid})$~(transformed from \Cref{eq:test-gen-valid-cons}).
        Thanks to the \emph{differentiable} interval abstractions, we can apply projected gradient descent to optimize $\vl_w$ and $\vl_x$ for minimizing $\gL_{n_0}^{\testgen}$.
        If the current $s$ is feasible, we then enlarge $s$ and start the next iteration.
        
        Once we figure out the largest possible $s$ for the optimization problem, we extract the solution and construct $[\vl_w, \vu_w]$ and $[\vl_x, \vu_x]$.
        Sampling from these failure-exhibiting weight and input intervals yields many distinct failure-exhibiting unit tests.
    \end{discard}
    
        \vspace{-0.8em}
    \subsection{Abstraction Optimization for Precondition-Fix Suggestion}
    
        \label{subsec:precond-gen}
        
        In this stage, we aim to generate the precondition fix for specified fix locations.
        The inputs are the DNN architecture, the node $n_0$ with numerical defects, and a subset of nodes $\Vfix$ to impose the fix.
        We would like to generate interval preconditions~(typically a single interval for all elements in a node) for $\Vfix$ node inputs, so that after the preconditions are imposed, the defect on $n_0$ is fixed.
        
        Formally, our task is to figure out $\langle l_n, u_n \rangle$ for each $n\in \Vfix$~(where $l_n$ and $u_n$ are scalars to indicate same interval bound for all elements), such that for any $\vx, \vw$ satisfying $\fin_{n}(\vx;\vw) \in [l_n, u_n]$, $\forall n\in \Vfix$, for the defect node $n_0$, we have $\fin_{n_0}(\vx;\vw) \notin \gI_{n_0,\textsf{invalid}}$, where invalid input range $\gI_{n_0, \textsf{invalid}}$ is dependent on $n_0$'s operator type and examples are in \Cref{tab:op_range_loss_form}~(second column).
        We define a surrogate optimization problem for this task as below.
        \vspace{-0.5em}
        \allowdisplaybreaks
        \begin{align}
            & \maximize_{l_n, u_n: n\in \Vfix} \quad s 
            \label{eq:precond-gen-obj} \\
            \mathrm{s.t.}\quad & u_n \ge l_n + s(u_n^{\valid} - l_n^{\valid}), \forall n\in \Vfix \label{eq:precond-span-cons} \\
            & l_n^{\valid} \le l_n \le u_n \le u_n^{\valid}, \forall n\in \Vfix \label{eq:precond-valid-cons} \\
            & \gL_{n_0}^{\precondgen}(\{l_n, u_n\}_{n\in\Vfix}) < 0 \label{eq:precond-effective-cons}
            \vspace{-0.5em}
        \end{align}
        In the optimization, 
        $l_n^{\valid}$ and $u_n^{\valid}$ are the valid ranges of the node's input $n$, which is fixed and determined by the valid ranges of input and weights.
        $\gL_{n_0}^{\precondgen}$ is the node-specific precondition generation loss where example definitions are in the fourth column of \Cref{tab:op_range_loss_form}.
        The optimization variables are the precondition interval endpoints $l_n$ and $u_n$ and the objective is the relative span of these intervals.
        The larger the span is, the looser the precondition constraints are, and the less hurt they are for the model's utility.
        \Cref{eq:precond-span-cons} enforces the interval span requirement.
        \Cref{eq:precond-valid-cons} assures that the precondition interval is in valid range.
        The last constraint~(\Cref{eq:precond-effective-cons}) guarantees the validness of the precondition fix.
        In the constraint, we design the precondition generation loss $\gL_{n_0}^{\precondgen}$ depending on the operator type such that if the loss value is negative, any concrete $\vx$ and $\vw$ value satisfying $\bigwedge_{n\in\Vfix} l_n \le \fin_n(\vx; \vw) \le u_n$ would not trigger numerical defect at node $n_0$.
        For example, for both \CodeIn{Log} and \CodeIn{Sqrt}, we design loss to be $-l_{n_0} + U_{\min}$.
        As a result, when $-l_{n_0} + U_{\min} < 0$, we have $l_{n_0} > U_{\min}$ and the possible input range $(l_{n_0}, u_{n_0}] \cap \gI_{n_0,\textsf{invalid}} = \emptyset$, i.e., the input to $n_0$ will always be valid.
        Hence, we can guarantee that the solution $\{l_n, u_n\}_{n\in \Vfix}$ constitutes a valid precondition fix.
        
        For any $\{l_n, u_n\}_{n\in \Vfix}$, thanks to the static analysis framework introduced in \Cref{subsec:static-analysis}, we can compute induced intervals of defect node $n_0$, and thus compute the loss value $\gL_{n_0}^{\precondgen}$.
        We use the coarsest tensor partitioning granularity in static analysis since interval bounds are shared among elements in one node.
        
        We propose the technique of \textbf{abstraction optimization} to effectively and approximately solve this optimization.
        Our technique works iteratively.
        In the first iteration, we set $s=1$, and in the subsequent iterations, we reduce the span $s$ exponentially: $s \gets \gamma_s s$, where hyperparameter $\gamma_s = 0.9$.
        Inside each iteration, for each node to impose precondition $n\in\Vfix$, we use the interval center $c_n = (l_n + u_n)/2$ as the optimizable variable and compute the \emph{sign} of its gradient: $\sgn(\nabla_{c_n} \gL_{n_0}^{\precondgen})$.
        We use this gradient sign to update each $c_n$ towards the direction of reducing loss value: $c_n \gets c_n - \gamma_c \max\{|c_n|, |\minstep|\} \sgn(\nabla_{c_n} \gL_{n_0}^{\precondgen})$.
        We set hyperparameters $\gamma_c = 0.1$ and $\minstep = 0.1$.
        Then, we use $c_n$ and $s$ to recover the actual interval: $l_n = c_n - s(u_n^{\valid}-l_n^{\valid})/2$, $u_n = c_n + s(u_n^{\valid}-l_n^{\valid})/2$, and then clip $l_n$ and $u_n$ by the valid range $[l_n^{\valid}, u_n^{\valid}]$.
        At the end of this iteration, for updated $l_n$ and $u_n$, we compute $\gL_{n_0}^{\precondgen}(\{l_n, u_n\}_{n\in\Vfix})$ and compare it with zero to check whether the precondition is a valid precondition fix.
        If so, we terminate; otherwise, we proceed to the next iteration.
        
        \begin{remark}
            The key ingredient in the technique is the gradient-sign-based update rule: $c_n \gets c_n - \gamma_c \max\{|c_n|, |\minstep|\} \sgn(\nabla_{c_n} \gL_{n_0}^{\precondgen})$ which is much more effective than the normal gradient descent due to following two reasons:
            (1)~Our update rule can get rid of gradient explosion and vanishing problem:.
            For early optimization iterations, the span $s$ is large and the interval bounds from static analysis are generally coarse, which results in too large or too small gradient magnitude.
            For example, the input range for \CodeIn{Log} could be $[1,10^{10}]$ where derivative of $\log(10^{10})$ is $10^{-10}$ which results in almost negligible gradient updates.
            In contrast, our update rule leverages the sign of gradient which always points to correct gradient direction.
            The update step size in our rule is the maximum of current magnitude $|c_n|$ and $\minstep$ to avoid stagnation.
            (2)~Our update rule mitigates the gradient magnitude discrepancy of different $c_n$.
            At different locations, the nodes in DNNs can have diverse value magnitudes and these magnitudes are not aligned with their gradient magnitudes, which makes the gradient optimization challenging.
           
           
           
           
            Therefore, we use this update rule to solve the challenge, where the update magnitude depends on the value magnitude but not on the gradient magnitude.
            We empirically compare our technique with standard gradient descent in \Cref{subsec:precond-gen-exp}.
        \end{remark}
        
        
\section{Implementation}

    \label{sec:impl}
    
    We implement a tool for our \sysname approach in roughly 10k lines of \CodeIn{Python} code.
    Our tool takes \CodeIn{ONNX}-format DNN architecture as the input, and automatically generates detection result, system tests, and precondition fixes~(fix location needs to be specified).
    The \CodeIn{ONNX} format is supported by popular DL libraries such as \CodeIn{PyTorch} \cite{paszke2019pytorch}, \CodeIn{Tensorflow} \cite{tensorflow2015-whitepaper}, and \CodeIn{Keras} \cite{chollet2015keras}.
    Our approach uses \CodeIn{PyTorch}~\cite{pytorchonnx} to provide tensor computation and auto-differentiation functionality, which has high performance on both CPU and GPU environment.
    More hyperparameters in \sysname are presented in \Cref{adxsec:detail-hyperparam}.
    We will open-source the tool after the de-anonymization process.

    
    \vspace{-0.4em}
\section{Experimental Evaluation}

    \label{sec:exp}

    We conduct systematic experimental evaluation to answer the following research questions.
    
    \vspace{-0.3em}
    \begin{enumerate}[leftmargin=*,label={\textbf{RQ\arabic*}}]
        \item 
        For \emph{generating failure-exhibiting system tests}, how more effectively and efficiently can \sysname trigger failures compared to baseline approaches?
        
        \begin{discard}
        \item 
        For \emph{generating many distinct failure-exhibiting unit tests}, how more effectively and efficiently can \sysname trigger failures compared to baseline approaches?
        \end{discard}
        
        \item 
        For \emph{suggesting precondition fixes}, how more efficiently and effectively is \sysname in terms of the compatibility of guarding against numerical failures compared to baseline approaches and developers' fix, respectively?
       
       
        
        \item
        For \emph{static defect detection}, how more effectively can \sysname detect numerical defects compared to state-of-the-art approaches?
       
        
    \end{enumerate}
    \vspace{-0.3em}
    
    We conducted evaluation on GRIST dataset~\cite{yan2021exposing}, which contains $85$ real-world DNN programs with numerical defects collected from previous studies and GitHub.
    This is the largest real-world DNN numerical defect dataset to our knowledge.
    We leverage existing tools~(\CodeIn{tf2onnx} and \CodeIn{PyTorch}) to extract DNN architectures from these programs to use as the input to both our tool and baseline approaches.
    Each program contains a DNN architecture, and each architecture contains one or more numerical defects.
   
    
    For RQ1, we evaluate at the defect level, where we inspect whether the approaches succeed in generating failure-exhibiting system tests for each defect in the architecture.
    For RQ2 and RQ3, we evaluate at the architecture level, where we inspect whether the approaches succeed in fixing and detecting \emph{all} defects in an architecture.
    The main quantitative results are reported in \Cref{tab:systest}~(RQ1), \Cref{tab:precond-summary}~(RQ2), and \Cref{tab:detect}~(RQ3).
    \begin{discard}
    For RQ1 and RQ2, we evaluate at the defect level, where we inspect whether the approaches succeed in generating failure exhibiting tests for each node that contains a numerical defect.
    For RQ3 and RQ4, we evaluate at the architecture level, where we inspect whether the approaches succeed in fixing and detecting \emph{all} defects in an architecture.
    The main quantitative results are reported in \Cref{tab:systest}~(RQ1), \Cref{tab:unittest}~(RQ2), \Cref{tab:precond-summary}~(RQ3), and \Cref{tab:detect}~(RQ4).
    \end{discard}
    
    
   
   
    
    We evaluate all approaches on the same workstation with 24-core Xeon E5-2650 CPU.
    Now we discuss each RQ in detail respectively.
   
    
    \vspace{-0.5em}
    \subsection{RQ1: System Test Generation}
    
        \label{subsec:system-test-gen-exp}
        
        \input{exp/tables/systest_table}
        
       
        
        We evaluate \sysname and baseline approaches on the effectiveness and efficiency for generating failure-exhibiting system tests.
        
        \paragraph{Baseline}
        Since \sysname is the first approach for this task,
       
        we did not compare with existing literature and propose one random based baseline approach~(called ``Random'' hereinafter)  as the baseline.
        In ``Random'', we first generate failure-exhibiting unit test by randomly sampling $1,000$ times.
        If there is any sample that triggers the failure, we keep the inference example part as the desired $\vxinfer$.
        Then, we generate $\vxtrain$ by randomly sampling $300$ times.
        This sampling number is made consistent with maximum iteration number in \sysname to ensure similar running time.
        If any sample, when used for one-step training, could induce model weights $\vw$ that cause numerical failure when using $\vxinfer$ as the inference input, we keep that sample as $\vxtrain$ and terminate.
        If both $\vxinfer$ and $\vxtrain$ are found, we count the case as a successful one for ``Random''; and otherwise a failed case.
        
        \paragraph{Evaluation Protocol}
        For each defect, due to the randomness of model's initial weights, we rerun both \sysname and ``Random'' for $10$ times and record the number of successful runs and the average execution time per run.
        For fair comparison, we fix both approaches to use the same set of random seeds.
       
       
        
        \paragraph{Main Evaluation Result}
        All results are presented in \Cref{tab:systest}.
        We can observe that \sysname succeeds for $892 / (104\times 10) \approx 85.8\%$ runs while the baseline ``Random'' succeeds for $\approx 73.8\%$ runs.
        Overall, \sysname spends only \SI{24.38}{s} time which is $2.76$X speedup compared to ``Random'' running time \SI{67.32}{s}.
       
       
        Specifically, \sysname spends less time and succeeds more times compared to the baseline in \emph{all} the cases.
        Thus, \textbf{\sysname is significantly more effective and efficient for generating failure-exhibiting system tests than the baseline.}
        
       
        \paragraph{Discussion and Failing Case Analysis}
        The high effectiveness of \sysname mainly comes from the advantage of gradient guided search compared with random search.
        As described in \Cref{subsec:system-test-gen}, \sysname leverages both first-order gradients and second-order derivatives to guide the search of system tests.
       
        In contrast, ``Random'' directly use random sampling hoping that failure-exhibiting training examples can reveal after sufficient sampling.
        As a result, when such training examples are relatively rare in the whole valid input space, ``Random'' is significantly less effective.
        Furthermore, we conduct an ablation study in \Cref{adxsec:system-testgen-by-stage} which shows \sysname improves over ``Random'' in both steps inside \sysname.
        
        We manually studied all $14$ defects where \sysname mostly fails~($S \le 4$ among $10$ runs).
       
        After inspecting the logs, we found that
        (1)~$3$ defects~(Defect ID 1, 23, 83) are because the DNN architecture is challenging for gradient-based optimization, e.g., due to \CodeIn{Min} and \CodeIn{Max} operators which have no gradient information.
        Further improvements on \sysname may help to solve these cases, which may need to dynamically detect operators with vanishing gradients and reconstruct the gradient flow which we leave as future work.
       
       
        (2)~$4$ defects~(Defect ID 63, 65, 88, 104) are because there are other numerical defects triggered when seeking training examples for current defect and the failure prevents \sysname from further execution.
        Since \sysname already exposes other defects, we believe \sysname is indeed effective for these cases.
        (3)~$7$ defects~(Defect ID 42-48) correspond to numerical defects that hardly lead to failure in practice since they are in \CodeIn{Div} operators where common real numbers, when used as a divisor, are too far from zero to trigger the failure.
        \sysname cannot find out failure-exhibiting system tests, which is as expected since \sysname is mainly designed to generate  tests for practical numerical defects instead of these impractical ones.
    
    \begin{discard}
    \subsection{RQ2: Unit Test Generation}
    
        \label{subsec:unit-test-gen-exp}
        
        \input{exp/tables/unittest_table}
        
        For this research question, we evaluate \sysname and baseline approaches on the effectiveness and efficiency for generating massive distinct failure-exhibiting unit tests.
        
       
        
        
        \paragraph{Baseline Approaches}
        \sysname is the first approach that generates massive distinct failure-exhibiting unit tests for DNN architectures on numerical failures, so we did not compare with existing literature and propose two baseline approaches to compare with:
        (1)~``Random'': this approach randomly samples model weights and inference-time inputs for $1000$ times, then keeps the instances that lead to numerical failures.
        (2)~``GD'': this baseline approach generates failure-exhibiting intervals, and samples from the intervals to get unit tests similar to \sysname.
        However, when generating these intervals, ``GD'' directly applies gradient descent on $(\vl_w, \vu_w, \vl_x, \vu_x)$ while keeping $s$ fixed~(defined in \Cref{subsec:unit-test-gen}), instead of using \sysname bilevel optimization technique.
        
        \paragraph{Evaluation Protocol}
        For each defect, due to the randomness in the sampling process, we rerun all approaches for 10 times.
        For each run, we aim to generate $1000$ distinct failure-exhibiting unit tests, and we set the time limit to be \SI{180}{s}.
        The number $1000$ is chosen to reflect the scenario where massive unit tests are needed.
        The time limit \SI{180}{s} reflects the efficiency requirement.
        We report the average number of generated failure-exhibiting distinct unit test and the average running time.
        Note that for ``GD'' and ``\sysname'', once the failure-exhibiting intervals are generated, any sample from the intervals is a failure-exhibiting unit test; 
        whereas for ``Random'', some samples may not trigger the numerical failure and these samples are not counted.
        
        \paragraph{Main Evaluation Result}
        All results are presented in \Cref{tab:unittest}.
        We observe that \sysname can generate $95800$ failure-exhibiting unit tests, giving a success rate of $95800 / (104\times 1000) \approx 92.1\%$, which is $6.6\%$ more than the best baseline ``GD''.
        Compared to the simple ``Random'' baseline, \sysname generates $42.9\%$ more unit tests while saving $21.1\%$ time.
        Thus, \textbf{\sysname is significantly more effective and efficient for generating massive distinct failure-exhibiting unit tests}.
        
        \paragraph{Discussion and Failing Case Analysis}
        When compared with ``GD'', \sysname generates more unit tests thanks to the bilevel optimization technique described in \Cref{subsec:unit-test-gen}.
        This technique allows the $(\vl_w, \vu_w, \vl_x, \vu_x)$ to be reused as a suitable initialization point when $s$ increases, thus increasing the successful rate.
        Though at an expense, compared with direct gradient descent, \sysname is slightly slower.
        When compared with ``Random'', \sysname generates more unit tests because the failure-exhibiting unit tests usually occupy a relatively small portion of the valid input and weight space.
        Therefore, random sampling from the space is ineffective.
        
        We further inspected the cases where \sysname fails.
        We found that except one run, all \sysname failure cases are in Case 17, 28c, 28d, and 28e.
        Among them, we found Case 17 corresponds to a deep and large DNN architecture, namely UNet~\cite{ronneberger2015u}, for image segmentation.
        Therefore, \sysname cannot finish in \SI{180}{s} time limit, nor did baseline approaches.
        Case 28c, 28d, and 28e correspond to DNN architectures where the numerical defects hardly lead to failure---the numerical defect is in a \CodeIn{Div} operator where the real-number divisor could be close to zero, leading to divide-by-zero failure, but with a small probability in practice.
        Therefore, it is reasonable for \sysname and all other baseline approaches to fail to generate failure-exhibiting tests for these cases.
    \end{discard}
    
    \vspace{-0.2em}
    \subsection{RQ2: Precondition-Fix Suggestion}
    
        \label{subsec:precond-gen-exp}
        
        We compare the precondition fixes generated by \sysname with fixes generated by baseline approaches and developers' fixes.
        
        \subsubsection{\textbf{Comparison between \sysname and Baselines}}
            \input{exp/tables/precond_summary_table}
            
           
            
            \paragraph{Baselines}
            \sysname is the first approach for this task,
           
            and we propose two baseline approaches to compare with:
            (1)~``\sysname-E'': this approach follows the optimization of \sysname described in \Cref{subsec:precond-gen}.
            However, \sysname-E uses standard interval domain instead of block-level interval domain introduced in \Cref{subsec:precond-gen}.
            Due to the lack of block-level abstraction, this approach is expected to be less efficient than \sysname.
            (2)~``GD'': this approach uses gradient descent (instead of the technique of abstraction optimization in \sysname) over the block-level interval domain to solve the optimization problem which is expected to be less effective than \sysname.
            
           
            \paragraph{Evaluation Protocol}
            We evaluate whether each approach can generate precondition fixes that eliminate \emph{all} numerical defects for a given DNN architecture given imposing locations.
            We evaluate at the case level instead of defect level, since each case corresponds to a DNN architecture that may contain multiple defects and we care about eliminating all defects.
            We consider three types of imposing locations: imposing on both weight and input nodes, on only weight nodes, and on only input nodes.
            In practice, model providers can impose precondition fixes on weight nodes by clipping weights during training; and model users can impose precondition fixes on input nodes by clipping their inputs before feeding into the model.
            We say that the precondition fix eliminates all numerical defects if and only if (1)~the \sysname static analysis framework, which is guaranteed to detect all possible numerical defects, cannot detect any defects from the fixed architecture; and (2)~random sampling $1,000$ times cannot trigger any numerical failures after imposing the fix.
            Since all the evaluated approaches are deterministic, for each case we run only once.
            Time limit is \SI{180}{s} per case.
            
            \paragraph{Main Evaluation Result}
            Due to the space limit, we report the statistics, including number of successful cases among all $85$ cases and the total running time, in \Cref{tab:precond-summary}.
            Detailed results are in \Cref{adxsec:detail-precond-gen-exp}.
            From the table, we observe that on all the three imposing location settings, \sysname succeeds on the most cases and usually spends much less time.
            The two baselines, \sysname-E and GD, either spends much longer time or has much lower success rate.
            For example, when precondition fixes can be imposed on both weight and input nodes, \sysname succeeds on \emph{all} cases with a total running time \SI{53.86}{s}.
            In contrast, \sysname-E requires $>10 \times$ running time, and GD fails on $25.9\%$ cases.
            Thus, \textbf{\sysname is significantly more effective and efficient for suggesting precondition fixes compared to baseline automated approaches.}
            
           
            \paragraph{Discussion}
            The two baseline approaches can be viewed as ablated versions of \sysname.
            Comparing \sysname and GD, we conclude that the technique of abstraction optimization significantly improves the effectiveness and also helps the efficiency of \sysname.
            Comparing \sysname and \sysname-E, we conclude that the block-level abstraction with tensor partitioning significantly improves the efficiency and also helps the effectiveness of \sysname.
            
            Another finding is that it is much easier to find precondition fix when the imposing locations are weight nodes compared to input nodes.
            For example, when imposing on weights, \sysname can succeed on $77$ cases while on input it succeeds on only $37$ cases.
            Since model providers can impose fixes on weights and users impose on inputs, this finding implies that, fixing numerical defects on providers' side may be more effective than on end-users' side.
        
        
           
        \subsubsection{\textbf{Comparison between \sysname and Developer's Fixes}}
                
            We conduct an empirical study to compare the precondition fixes generated by \sysname and by the developers.
            
           
            \paragraph{Evaluation Protocol}
            We manually locate the GitHub repositories from which the GRIST dataset is constructed.
            Among the $85$ architectures with numerical defects, we find the repositories for $55$ cases on GitHub and we study these cases.
            We locate the developers' fixes of the numerical defects by looking at issues and follow-up pull requests~(PRs).
            Then, two authors independently compare the developers' fixes and \sysname fixes and classify the comparison results.
            Since \sysname suggests different precondition fixes for different imposing locations, for each case we first determine the imposing locations from developer's fix, then compare with \sysname's fix for that  location.
            Lastly, two authors coordinate their results and discuss the misaligned results to reach consensus.
            
            \paragraph{Results}
            We categorize the comparison results between \sysname precondition fixes and developer's fixes for these $55$ cases as below.
            \begin{itemize}[leftmargin=*]
                \item \textbf{(19 cases) Better than developer's fixes or no available developer's fix.}
                For these cases, developers or users either propose no fix or use heuristic fixes, such as reducing the learning rate or using the mean value to reduce the variance.
                These fixes may work well in practice but cannot rigorously guarantee the elimination of the numerical failure for any training or inference data.
                In contrast, \sysname can generate valid precondition fixes that rigorously eliminate numerical failure.
               
               
                
                \item \textbf{(19 cases) Equivalent to developer's fixes.}
                Both developers/users and \sysname suggest equivalent or highly similar fixes.
               
                
                \item \textbf{(14 cases) No need to fix.}
                For these cases, there is no need to fix the numerical defect in the given architecture.
                There are mainly three reasons:
                (1)~The DNN architecture is used in the whole project with fixed weights or test inputs.
                As a result, even though the architecture contains defects, it will not lead to system failure.
                (2)~The architecture deliberately includes the defect as a test case for automatic tools, such as a test architecture in \CodeIn{TensorFuzz}~\cite{odena2019tensorfuzz} repository.
                (3)~The defect can be hardly triggered in practice. 
                For example, the defect is in a \CodeIn{Div} operator where the divisor can be close to zero to trigger divide-by-zero failure, but this hardly happens in practice since the divisor is randomly initialized as a non-zero number.
                
                \item \textbf{(3 cases) Inferior than developer's fix or \sysname fix is impractical.}
                In two cases, \sysname is inferior than developer's fix:
                Human developers observe that the defect operator is \CodeIn{Log}, and its input is non-negative.
                So they propose to add $10^{-6}$ to the input of \CodeIn{Log} as the fix.
                In contrast, \sysname can only generate a clipping-based fix, e.g., clipping the input if it is less than $10^{-6}$, due to the nature of precondition.
                When the input is small, \sysname's precondition fix interrupts the gradient flow from output to input while human's fix maintains it.
                As a result, human's fix does less hurt to model's trainability and is better than \sysname's fixes.
                In one other case, the precondition fix suggested by \sysname imposes a small span for some model weights~(less than $0.1$ for each component of that weight node).
                Such a small weight span strongly limits the model's expressiveness/capacity.
               
            \end{itemize}
           
            
            \vspace{-0.8em}
            \paragraph{Discussion}
            From the comparison results, we can conclude that for the 41 cases where numerical defects are needed to be fixed, the automatic \sysname approach suggests equivalent or better fixes than human developers in 38 cases.
            Therefore, \textbf{\sysname is comparably effective as human developers in terms of suggesting numerical defect fixes, and is much more efficient since \sysname is an automatic approach}.
            
            To further improve \sysname,
            two possible directions are 
            (1)~Considering more forms of fixes: 
            \sysname only supports imposing interval preconditions as the fix.
            From human developers, there are more fix forms such as adding or subtracting a constant number or replacing the operator type.
            For example, for a \CodeIn{Log} operator whose input is non-negative, we can suggest adding a constant positive value to input as a better fix than imposing interval precondition.
            (2)~Tightening the static analysis:
            With a tighter static analysis, we can synthesize a looser precondition as the suggested fix.
           
           
           
        
    
    \vspace{-0.7em}
    \subsection{RQ3: Static Defect Detection}
        \label{subsec:static-detector-exp}
        
        \input{exp/tables/detect_table}
        
        We evaluate \sysname ability on static numerical defect detection and compare it with the state-of-the-art approach DEBAR~\cite{zhang2020detecting}.
        
        The results are shown in \cref{tab:detect}.
        For each test case, if the approach successfully detects all numerical defects, we label it with ``$\checkmark$''; otherwise ``$\times$''.
        For DEBAR, we digest its results from \citeauthor{yan2021exposing}~\cite{yan2021exposing} and replicate on our testing platform.
        For DEBAR, there are some cells with labels ``?''.
        These cells correspond to the cases that exist in GRIST dataset but no results for them are reported in \cite{yan2021exposing}.
        Suppose DEBAR succeeds on all ``?'' cases, then DEBAR detects all defects for 54 cases and \sysname detects defects for all 85 cases, being $57.4\%$ more than DEBAR.
        In terms of efficiency, DEBAR and \sysname have similar running time and both finish in \SI{3}{s} per case.
        Therefore, \textbf{\sysname is significantly more effective, in terms of detecting numerical defects, when compared with the state-of-the-art approach}.
        The effectiveness of \sysname is due to its effective static analysis framework with efficient abstractions for 84 DNN operators, being 17 more than DEBAR.
       
        
        
    
   
    
   
    
   
        
        
   
    
   
    
    \vspace{-0.5em}
\section{Threats to Validity}

    A major external threat to validity is the evaluation subject.
    We evaluate \sysname and baseline approaches on the GRIST dataset~\cite{yan2021exposing}, which collects DNNs with numerical defects from real-world open-source projects and previous studies.
    Though the dataset contains over $80$ DNN architectures and is the largest to the best of our knowledge, it may still not be representative enough for real-world DNNs.
    Another external threat is the randomness from the evaluated approaches.
    To mitigate this threat, we repeat all randomized approaches 10 times with randomly generated seeds.
    To mitigate bias in the empirical study, two authors independently conducted the study and discussed all conflicting cases to reach consensus.
    
    A major internal threat to validity comes from the approach implementation.
    We reduce this threat by conducting code reviews and outputting extensive logs.
    An author independently checked the correctness of code implementation under the help of logs and replicated all major experiments.
    We also verified the soundness of our static abstraction framework with over 50 hand-crafted unit tests and extensive random testing.
    
    \vspace{-0.5em}
\section{Related Work}
    \vspace{-0.3em}

    Our work focuses on reliability assurance of DNNs against numerical defects.
    We outline the existing work on understanding and detecting DNN-based system defects and on analyzing, testing, and debugging DNNs, especially for numerical defects.

    \vspace{-0.5em}
    \paragraph{Understanding and Detecting Defects in DNNs}
    Discovering and mitigating the defects, faults, bugs, and failures in DNN based systems is an important research topic where several empirical studies exist~\cite{zhang2018empirical,pham2020problems,zhang2020empirical,humbatova2020taxonomy}.
    Roughly following the taxonomy in \cite{humbatova2020taxonomy}, DNN defects are at four levels from bottom to top:
    (1)~Platform-level defects: Defects can exist in real-world DL compilers and libraries. 
    Approaches exist for understanding, detecting, and testing against these defects~\cite{wang2020deep,guo2020audee, DBLP:conf/sigsoft/ShenM0TCC21, DBLP:conf/icse/LuoCRWFC21,tambon2021silent,liu2021detecting,xie2021leveraging}.
    (2)~Architecture-level defects: \emph{Our work focuses on numerical defects, which is one type of architecture-level defects.} 
    Automatic detection and localization techniques~\cite{wardat2021deeplocalize,liu2021detecting,rezaalipour2021test,hattori2020semi} exist for other architecture-level defects such as suboptimal structure, activation function, initialization and shape mismatch~\cite{humbatova2020taxonomy}.
   
   
    (3)~Model-level defects: Once a model is trained, its defects can be viewed as violations of desired properties which are discussed in detail in \cite{zhang2020machine}.
    Some example defects are correctness~\cite{tizpaz2020detecting,guerriero2021operation}, robustness~\cite{wang2021robot,tian2020testing}, and fairness~\cite{zhang2020white} defects.
    (4)~Interface-level defects: DNN-based systems, when deployed as services, expose interaction interfaces to users where defects may exist as empirical studies on real-world systems reveal~\cite{humbatova2020taxonomy,wan2021machine,wan2022automated}.
    
   
    \paragraph{Testing and Debugging for DNNs}
    Given DNN defects, a rich body of work exists for testing and debugging DNNs, which often targets at one type of defects.
    \citeauthor{zhang2020white}~\cite{zhang2020white} conduct a systematic survey for these works.
    Some representatives are DeepXplore~\cite{pei2017deepxplore} which generates test that joint maximizes behavior diversity and neuron coverage, and DeepGauge~\cite{ma2018deepgauge} which proposes multi-granularity neuron coverage metrics which are further studied~\cite{sun2018concolic,li2019structural,harel2020neuron} and improved~\cite{yuan2021you}.
    Several recent works enable automatic model debugging and repair via consistency checking~\cite{xiao2021self}, log checking~\cite{zhang2021autotrainer}, spectrum analysis~\cite{qi2021archrepair}, or analyzer-guided synthesis~\cite{sotoudeh2021provable}.
    
   
    \paragraph{Static Analysis for DNNs}
    Besides testing and debugging, another direction towards eliminating defects is conducting static analysis on DNNs to rigorously guarantee the existence/nonexistence of defects.
    Abstraction-based DNN analysis is extensively studied for verifying model robustness, where existing work focuses on either proposing tighter abstractions to reduce false positives~\cite{gehr2018ai2,wang2018formal,wang2021beta,muller2022prima,paulsen2022linsyn} or incorporating abstractions into training to improve robustness~\cite{wong2018provable,li2019robustra,mirman2018differentiable,zhang2020towards,zhang2021robustness}.
    Besides robustness, the abastraction-based analysis has also be applied to rigorously bound model difference~\cite{paulsen2020reludiff,paulsen2020neurodiff}.
    \citeauthor{li2020sok}~\cite{li2020sok} and \citeauthor{albarghouthi2021introduction}~\cite{albarghouthi2021introduction} provide comprehensive survey for this field.
    Our approach includes a static analysis framework customized for numerical defect detection and fixing.
    
    \vspace{-0.5em}
    \paragraph{Detecting and Exposing Numerical Defects in DNNs}
    Despite the widespread existence of numerical defect in real-world DNN-based systems~\cite{zhang2018empirical,humbatova2020taxonomy}, relatively few automatic approaches exist for detecting and exposing these defects.
    \citeauthor{zhang2020detecting}~\cite{zhang2020detecting} propose the first approach, DEBAR, for detecting numerical defects based on static analysis.
    Compared to DEBAR, our analysis framework allows automatic differentiation over the constructed abstractions which enables abstraction optimization for precondition-fix generation and our framework includes abstractions for more types of DNN operators which enable more effective defect detection.
    \citeauthor{yan2021exposing}~\cite{yan2021exposing} propose a gradient back-propagation approach, GRIST, for generating failure-exhibiting \emph{unit test} to expose numerical defects.
    Our \sysname builds upon GRIST to generate \emph{system tests} to confirm the feasibility of numerical defects.
   
   
    \begin{discard}
    \linyim{For generating failure-exhibiting massive distinct unit tests, we did not compare with GRIST because GRIST is not designed and optimized for generating massive distinct unit tests. Even though repeating GRIST execution multiple times may generate massive unit tests, comparing the pairwise difference to ensure the ``distinct'' requirement is costly.}
    \end{discard}
    \citeauthor{yan2021exposing}~\cite{yan2021exposing} also collected the large dataset of real-world DNN numerical defects on which we conduct experimental evaluation.
    To the best of our knowledge, \sysname is the first approach that supports all three-stage tasks for reliability assurance against DNN numerical defects. 
    
   
\section{Conclusion}
    In this paper, we have presented a novel automatic approach named  \sysname for reliability assurance of DNNs against numerical defects.
    \sysname supports static defect detection, system test generation, and precondition-fix suggestion.
    \sysname includes a highly  efficient and effective DNN static analysis framework and two novel techniques~(two-step test generation and abstraction optimization).
    Our extensive evaluation on real-world DNN programs has demonstrated high effectiveness and efficiency of \sysname.
   
   



    
\section{Introduction}
   
    \label{sec:intro}
    Deep Neural Networks~(DNNs) are successfully deployed and show remarkable performance in many challenging applications, including facial recognition~\cite{liu2015deep,zhang2016joint}, game playing~\cite{silver2017mastering}, and code completion~\cite{liu2022unified,bui2021infercode}. To develop and deploy DNNs, one needs to attain a DNN architecture, which is usually encoded by program code as the example shown in \Cref{fig:dnn-code-example}.
    First, for training, the user executes the program with the architecture on the given  training/validation data, attains the model weights, and stores them in a weight file.
    The architecture along with the weights is named a model.
    Then, for inference, the user loads the weight file to CPU/GPU memory or AI chips, executes the same program with the given inference sample and weights as arguments, and gets the model prediction result as the program output.
    With the wide deployment of DNN models (resulted from training DNN architectures), reliability issues of DNN-based systems have become a serious concern, where malfunctioning DNN-based systems have led to serious consequences such as fatal traffic accidents~\cite{drivingfatal}.
    
    To assure the reliability of DNN-based systems, it is highly critical to detect and fix numerical defects
   
   
   
    for two main reasons. 
   
   
   
    First, numerical defects widely exist in DNN-based systems. 
    For example, in the DeepStability database~\cite{kloberdanz2022deepstability}, over 250 defects are identified in deep learning~(DL) algorithms where over 60\% of them are numerical defects. 
   
    Moreover, since numerical defects exist at the architecture level, any model using the architecture naturally inherits these defects. 
    Second, numerical defects can result in serious consequences. Once numerical defects~(such as divide-by-zero) are triggered, the faulty DNN model will output \NaN or \INF instead of producing any meaningful prediction, resulting in numerical failures and system crashes~\cite{zhang2018empirical,zhang2019empirical}. Thus, numerical defects  hinder the application of DNNs in scenarios with high reliability and availability requirements such as threat monitoring in cybersecurity~\cite{problemaisecurity} and cloud systems controlling~\cite{sharma2016machine,jay2019deep}.
    
    \begin{figure}[t]
        \includegraphics[width=\linewidth]{figures/workflow_new.pdf}
        \vspace{-2.25em}
        \caption{Workflow for reliability assurance against numerical defects in DNN architectures.
        The left-hand side shows three tasks and the right-hand side shows corresponding examples.
        \sysname supports all three tasks, and is the first automatic approach for system test generation and fix suggestion.}
       
        \label{fig:workflow}
        \vspace{-0.75em}
    \end{figure}
    
   
   
   
    To address numerical defects in DNN architectures in an actionable manner~\cite{SCIS22}, in this paper, we propose a workflow of reliability assurance (as shown in \Cref{fig:workflow}), consisting of three tasks: 
    potential-defect detection, feasibility confirmation, and fix suggestion, along with our proposed approach to support all three tasks.
    
  
   
   
   
   
   
    
    \emph{Potential-Defect Detection.}
    In this task, we detect all potential numerical defects in a DNN architecture, with a focus on inference-phase numerical defects (which are operators that potentially exhibit numerical failures in the inference phase) for two main reasons, 
   
    following the literature~\cite{zhang2020detecting,yan2021exposing}. First, inference-phase numerical defects can be triggered after the model is deployed and thus are more devastating than training-phase defects~\cite{odena2019tensorfuzz,zhang2020detecting}. Second, training-phase numerical defects can usually be detected during the inference phase.
    For example, the type of training-phase \NaN gradient defect is caused by an operator's input that leads to invalid derivatives, where this input also triggers failure in the inference phase~\cite{yan2021exposing}.
    
    \emph{Feasibility Confirmation}.
    In this task, we confirm the feasibility of these potential numerical defects by generating failure-exhibiting system tests.
    As shown in \Cref{fig:workflow}, a system test is a tuple of training example\footnote{In real settings, multiple training examples are used to train an architecture, but generating a single training example to exhibit failures (targeted by our work) is desirable for ease of debugging while being more challenging than generating multiple training examples to exhibit failures.} $\vxtrain$ and inference example $\vx$ such that after the training example is used to train the architecture under consideration, applying the resulting model on the inference example exhibits a numerical failure.
   
   
   
    
    \begin{discard}
    \textbf{Fault Localization.}
    With failure-exhibiting system and unit tests, we locate the faults.
    A fault corresponds to a program location, which is the output of this task.
    Both manual inspection and automatic approaches~\cite{wardat2021deeplocalize} can conduct fault localization.
    Specifically, for coverage-based fault localization tools~\cite{abreu2009practical,abreu2009spectrum,abreu2006evaluation,lou2021boosting}, many distinct failing unit tests are required.
    In this paper, we do not focus on this task.
   
   
    \end{discard}
    
    \emph{Fix Suggestion.}
     In this task, we fix a feasible numerical defect. 
    To determine the fix form,
    we have inspected the developers' fixes of the numerical defects collected by \citeauthor{zhang2018empirical}~\cite{zhang2018empirical} by looking at follow-up Stack Overflow posts or GitHub commits.
    Among the 13 numerical defects whose fixes can be located,
    12 fixes can be viewed as explicitly or implicitly imposing interval preconditions on different locations, such as after inputs or weights are loaded and before suspicious operators (ones that potential defects are associated with) are invoked.
   
    Thus, imposing an interval precondition, e.g., by clipping (i.e., chopping off the input parts that exceed the specified input range) the input for suspicious operator(s), is an effective and common strategy for fixing a numerical defect.
    Given a location~(i.e., one related to an  operator, input, or weight where users prefer to impose a fix), we suggest a fix for the  numerical defect under consideration.
    
   
   
    
   
   
   
   
    
   
   
   
   
   
   
   
   
   
   
   
   
   
    
    
    To support all three tasks of the \underline{\textbf{r}}ealiability \underline{\textbf{a}}ssurance process against DNN \underline{\textbf{num}}erical defects, we propose the \textbf{\sysname} approach in this paper.
    
    \emph{For task \circled{1} and task \circled{2}a, which are already supported by two existing tools (DEBAR~\cite{zhang2020detecting} and GRIST tool~\cite{yan2021exposing}), \sysname introduces novel extensions and optimizations that substantially improve the effectiveness and efficiency}. 
    (1)~The state-of-the-art tool for potential-defect detection is DEBAR.
    However, DEBAR can handle only static computational graphs and does not support widely used dynamic graphs in PyTorch programs~\cite{paszke2019pytorch}.
    \sysname supports dynamic graphs thanks to our novel technique of \underline{\emph{backward fine-grained node labeling}}.
   
   
   
   
    (2)~The recent GRIST tool generates failure-exhibiting unit tests for confirming potential-defect feasibility by gradient back-propagation but uses the original inference input and weights as the starting point.
   
   
    However, recent studies~\cite{madry2018towards,goodfellow2015explaining} on DNN adversarial attacks suggest that using a randomized input as the starting point leads to stronger attacks than using the original input.
    Taking this observation, we combine gradient back-propagation with random initialization in \sysname.
   
   
    
    \emph{For task \circled{2} and task \circled{3}, which are not supported by any existing tool, \sysname is the first automatic approach for them}. 
    
    For \textbf{feasibility confirmation}, \sysname is  the \textbf{first} approach that generates failure-exhibiting \textbf{system} tests that contain training examples. Doing so is a major step further from the existing GRIST tool, which generates failure-exhibiting unit tests ignoring the practicality of generated model weights.
    Given that in practice model weights are determined by training examples, we propose the technique of \underline{\emph{two-step generation}} for this task. 
    First, we generate failure-exhibiting unit tests. 
    Second, we generate a training example that leads to the model weights in the unit test when used for training.
    For the second step, we transform the task into finding an input that leads to specific model gradients, where we extend the deep-leakage-from-gradient~(DLG) attack~\cite{zhu2020deep} by incorporating the straight-through gradient estimator~\cite{bengio2013estimating}.
    
    For \textbf{fix suggestion}, \sysname is the \textbf{first} automatic approach.
    \sysname is based on the novel technique of \underline{\emph{abstraction}} \underline{\emph{optimization}}. 
    We observe that a defect fix in practice is typically imposing interval clipping on some operators such that each later-executed operator (including those suspicious ones) can never exhibit numerical failures.
   
    Therefore, 
   
    we propose the novel technique of abstraction optimization to ``deviate away'' the input range of a suspicious operator from the invalid range, falling in which can cause numerical failures. 
   
   
    
    
    
    
    
      
    
    
 
    
    
    
    
    
    
    
    
    
For \sysname, we implement a tool\footnote{Open source at \texttt{\url{https://github.com/llylly/RANUM}}.} and evaluate it on the   benchmarks~\cite{yan2021exposing} of 63 real-world DNN architectures containing 79 true numerical defects, which are the largest benchmarks of DNN numerical defects to the best of our knowledge.
    The evaluation results show that \sysname is both effective and efficient in all three tasks   for DNN reliability assurance.  
(1)~For potential-defect detection, \sysname detects $>$60\% more true defects than the state-of-the-art DEBAR approach. 
(2)~For feasibility confirmation, \sysname generates  failure-exhibiting unit tests
to confirm  potential numerical defects in the benchmarks with 100\% success rate, whereas with the much higher time cost (17.32X), the state-of-the-art GRIST approach generates unit tests to confirm defects with 96.96\% success rate. More importantly, for the first time, \sysname generates failure-exhibiting system tests that confirm defects (with 92.78\% success rate). 
(3)~For fix suggestion, \sysname proposes fix suggestions for  numerical defects with 100\% success rate. 
In addition, when the \sysname{}-generated fixes are compared with developers' fixes on open-source projects, in 37 out of 40 cases, \sysname{}-generated  fixes are equivalent to or even better than human fixes.   
 
 
 


   
    
    This paper makes the following main contributions: 
    \begin{itemize}[leftmargin=*]
        \item 
        We formulate the reliability assurance problem for DNN architectures against numerical defects and elaborate on three important tasks for this problem.
    
        \item 
        We propose \sysname---the first automatic approach that solves all these three tasks.
        \sysname includes three novel techniques~(backward fine-grained node labeling, two-step test generation, and abstraction optimization) and solves system test generation and fix suggestion for the first time.
       
        
       
        
       
        
        \item We implement  \sysname and apply it on 63 real-world DNN architectures, showing the high effectiveness and efficiency of \sysname  compared to both the state-of-the-art approaches and developers' fixes.
       
       
       
       
    \end{itemize}
    
   
   
   
    
   
   
   
    
   
   
   
    
   
    
   
    
   
    
   
   
   
   
    
   
   
   
   
   
   
   
    
   
   
   
   
    
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
    
   
   
   
   
   
   
   
    
   
    
   
    
   
    
   
    
   
    
    
    
   
   
   
   
    
   
        
   
        
   
        
   
   
   
    
   
   

    \vspace{-0.25em}
\section{Background and Approach Overview}
   
    \label{sec:background-task}


    In this section, we introduce the background of DNN numerical defects and failures, and then give an overview of the \sysname approach with a running example. 
   

       
    \subsection{Background}
       
        \label{subsec:background}
    
    \begin{figure}[t]
        \centering
        \includegraphics[width=0.95\linewidth]{figures/DNN_code_example.pdf}
        \vspace{-1.2em}
        \caption{A DL program snippet that defines a linear regression model from real-world numerical defect benchmarks~(Case 2a in \cite{yan2021exposing}).}
        \label{fig:dnn-code-example}
        \vspace{-1.0em}
    \end{figure}
    
       
       
       
       
       
        
        DL developers define DNN architecture with code using modern DL libraries such as PyTorch~\cite{paszke2019pytorch} and Tensorflow~\cite{tensorflow2015-whitepaper}.
        The DNN architecture can be expressed by a computational graph.
        \Cref{fig:dnn-code-example,fig:dnn-graph-example} depict a real-world example.
        Specifically, the DNN architectures in DL programs can be automatically converted to ONNX-format computational graph~\cite{onnx}.
        
       
       
       
       
        
        The computational graph can be viewed as a Directed Acyclic Graph~(DAG): $\gG = \langle \gV, \gE \rangle$, where $\gV$ and $\gE$ are sets of nodes and edges respectively.
        We call nodes with zero in-degree as \emph{initial nodes}, which correspond to input, weight, or constant nodes.
       
        Initial nodes provide concrete data for the DNN model.
        The data from each node is formatted as a tensor, i.e., a multidimensional array, with a specified data type and array shape annotated alongside the node definition.
        We call nodes with positive in-degree as \emph{internal nodes}, which correspond to concrete operators, such as matrix multiplication~(\CodeIn{MatMul}) and addition~(\CodeIn{Add}).
        During model training, the model weights, i.e., data from weight nodes, are generated by the training algorithm.
        Then, in the deployment phase~(i.e., model inference), with these trained weights and a user-specified input called inference example,
        the output of each operator is computed in topological order. 
        The output of some specific node is used as the prediction result.
       
       
        
        We let $\vx$ and $\vw$ denote the concatenation of data from all input nodes and data from all weight nodes respectively.\footnote{Bolded alphabet stands for vectors or tensors throughout the paper.}
        For example, in \Cref{fig:dnn-graph-example}, $\vx$ concatenates data from nodes 1 and 11; and $\vw$ concatenates data from nodes 2 and 4.
        Given specific $\vx$ and $\vw$, the input and output for each node are deterministic.\footnote{Some architectures contain stochastic nodes. We view these nodes as nodes with randomly sampled data, so the architecture itself is deterministic.}
        We use $\fin_n(\vx; \vw)$ and $\fout_n(\vx; \vw)$ to express input and output data of node $n$, respectively, given $\vx$ and $\vw$.
        
        \noindent
        \textbf{Numerical Defects in DNN Architecture.}
        We focus on inference-phase numerical defects.
        These defects lead to numerical failures when specific operators receive inputs within invalid ranges so that the operators output \NaN or \INF
       
       
        
        
        \begin{definition}
            \label{def:numerical-defect}
            For a given computational graph  $\gG=\langle \gV, \gE\rangle$,
            if there is a node $n_0 \in \gV$, such that there exists a valid input and valid weights that can let the input of node $n_0$ fall within the invalid range, we say there is a \textbf{\emph{numerical defect}} at node~$n_0$.
            \vspace{-0.5em}
            \begin{equation*}
                \label{eq:numerical-defect}
                \begin{aligned}
                    \text{Formally,}\, & \exists \vx_0 \in \gXvalid, \vw_0 \in \gWvalid, 
                    \fin_{n_0}(\vx_0; \vw_0) \in \gI_{n_0,\mathsf{invalid}} \\
                    \Longrightarrow &
                    \text{$\exists$ numerical defect at node $n_0$}.
                \end{aligned}
            \end{equation*}
        \end{definition}
        In the definition, $\gXvalid$ and $\gWvalid$ are valid input range and weight range respectively, which are clear given the deployed scenario.
        For example, ImageNet \CodeIn{Resnet50} model has valid input range $\gXvalid = [0,1]^{3\times 224\times 224}$ since image pixel intensities are within $[0,1]$, and valid weight range $\gWvalid = [-1,1]^{p}$ where $p$ is the number of parameters since weights of well-trained \CodeIn{Resnet50} models are typically within $[-1,1]$.
        The invalid range $\gI_{n_0,\mathsf{invalid}}$ is determined by $n_0$'s operator type with detailed definitions in \Cref{adxsec:list-numerical-defects}.
        For example, for \CodeIn{Log} operator, 
       
        the invalid range $\gI_{n_0,\mathsf{invalid}} = (-\infty, U_{\min})$ where $U_{\min}$ is the smallest positive number of tensor's data type.

        \vspace{-0.5em}
    \subsection{Approach Overview}
        \label{subsec:approach-overview}
        
        \begin{figure}[!t]
            \centering
            \includegraphics[width=0.9\linewidth]{figures/DNN_graph_example.pdf}
            \vspace{-1em}
            \definecolor{figbrown}{RGB}{203,97,49}
            \definecolor{figpurple}{RGB}{112,56,153}
            \definecolor{figgreen}{RGB}{120,180,98}
            \definecolor{figyellow}{RGB}{235,163,27}
            \caption{Computational graph encoded by \Cref{fig:dnn-code-example} snippet.}
           
           
           

            \label{fig:dnn-graph-example}
            \vspace{-1.0em}
        \end{figure}
    
        \begin{figure}[t]
           
            \centering
            \includegraphics[width=0.90\linewidth]{figures/overview_new.pdf}
            \vspace{-0.8em}
            \caption{Overview of the \sysname approach. 
            The output of \sysname indicates confirmation and manifestation of numerical defects (that can be feasibly triggered at the system level) for a given DNN architecture and effective fixes for the architecture's confirmed defects.
           
            }
           
            \label{fig:overview}
            \vspace{-1em}
        \end{figure}
        
        In \Cref{fig:overview}, we show the overview structure of the \sysname approach.
        \sysname takes a DNN architecture as the input.
        Note that though \sysname is mainly designed and illustrated for DNN architectures, the \sysname approach can also be directly applied to general neural network architectures since they can also be expressed by computational graphs.
        First, the DNN static analysis framework in \sysname detects all potential numerical defects in the architecture, which is the task \circled{1} in \Cref{fig:workflow}.
        Second, the two-step test generation component~(unit test generation + training example generation) confirms the feasibility of these potential numerical defects, which is the task \circled{2} in \Cref{fig:workflow}.
       
       
       
        Third, the abstraction optimization component takes the input/output abstractions produced by the DNN static analysis framework along with the user-specified fix locations and produces preconditions to fix all defects,
        which is the task \circled{3} in \Cref{fig:workflow}.
        
        We now go through the whole process in detail taking the DNN architecture shown in \Cref{fig:dnn-graph-example} as a running example.
        
        \paragraph{Task \circled{1}: Potential-Defect Detection via Static Analysis}
        The DNN static analysis framework within \sysname first computes the numerical intervals of possible input and outputs for all nodes within the given DNN architecture, and then flags any nodes whose input intervals overlap with their invalid ranges as nodes with potential numerical defects.
        
        In \Cref{fig:dnn-graph-example}, suppose that the user-specified input \CodeIn{x-input}~(node 1) is within (elementwise, same below) range $[(-10,-10)^\T, (10,10)^\T]$; \CodeIn{weights}~(node 2) is within range ${\scriptsize[\begin{bmatrix} -10 & -10 \\ -10 & -10 \end{bmatrix}, \begin{bmatrix} 10 & 10 \\ 10 & 10 \end{bmatrix}]}$; and \CodeIn{biases}~(node 4) is within range $[(-10,-10)^\T, (10,10)^\T]$.
        Our DNN static analysis framework will compute these interval abstractions for node inputs:
        \begin{enumerate}[leftmargin=*]
            \item Node 5~(after \CodeIn{MatMul}): $[(-200,-200)^\T, (200,200)^\T]$;
            \item Node 6~(after \CodeIn{Add}): $[(-210,-210)^\T, (210,210)^\T]$;
            \item Node 8~(after \CodeIn{Softmax} in \CodeIn{float32}): $[(0,0)^\T, (1,1)^\T]$;
            \item Node 9~(after \CodeIn{Sub} of $[1,1]$ and node 8), 10: $[(0,0)^\T, (1,1)^\T]$.
        \end{enumerate}
        Since node 9 and 10 use \CodeIn{Log} operator whose invalid input range $(-\infty, U_{\min})$ overlaps their input range $[(0,0)^\T, (1,1)^\T]$, we flag node 9 and 10 as potential numerical defects.
        
        This static analysis process follows the state-of-the-art DEBAR tool~\cite{zhang2020detecting}.
        However, we extend DEBAR with a novel technique, backward fine-grained node labeling. 
        This technique detects all nodes that require fine-grain abstractions, e.g., nodes that determine the control flow in a dynamic graph. 
        For these nodes, we apply interval abstractions with the finest granularity to reduce control flow ambiguity.
        For other nodes, we let some neighboring elements share the same interval abstraction to improve efficiency while preserving tightness.
        As a result, the static analysis in \sysname has high efficiency and supports much more DNN operators including dynamic control-flow operators like \CodeIn{Loop}.
       
       
        
        \paragraph{Task \circled{2}: Feasibility Confirmation via Two-Step Test Generation}
        Given nodes that contain potential numerical defects~(node 9 and 10 in our example), we generate failure-exhibiting system tests to confirm their feasibility. 
        A failure-exhibiting system test is a tuple $\langle\vxtrain, \vxinfer\rangle$, such that after training the architecture with the training example $\vxtrain$,\footnote{More explicitly, if our generation technique outputs $\vxtrain$, the numerical failure can be triggered if the training dataset contains \emph{only} $\vxtrain$ or \emph{only} multiple copies of $\vxtrain$ and inference-time input is $\vxinfer$. Our technique can also be applied for generating a batch of training examples by packing the batch as a single example: $\vxtrain = ({\vxtrain}_1, {\vxtrain}_2, \dots, {\vxtrain}_B)$.} with the trained model weights $\vwinfer$, the inference input $\vxinfer$ triggers numerical failure.
        The name ``system test'' is inspired by traditional software testing, where we test the method sequence~(\CodeIn{m = train($\vxtrain$); m.infer($\vxinfer$)}). 
        In contrast, GRIST~\cite{yan2021exposing} generates model weights $\vwinfer$ and inference input $\vxinfer$ that only tests the inference method \CodeIn{m.infer()} and the weights may be infeasible from training. 
        Hence, we view the GRIST-generated tuple $\langle\vwinfer, \vxinfer\rangle$ as a ``unit test''.
        
        
        We propose a two-step test generation technique to generate failure-exhibiting system tests.
        
        \emph{\textbf{Step a}: Generate failure-exhibiting unit test $\langle\vwinfer, \vxinfer\rangle$.}
        The state-of-the-art GRIST tool supports this step.
        However, GRIST solely relies on gradient back-propagation, which is relatively inefficient.
        In \sysname, we augment GRIST by combining its gradient back-propagation with random initialization inspired by recent research on DNN adversarial attacks~\cite{goodfellow2015explaining,madry2018towards}.
        As a result, \sysname achieves 17.32X speed-up with 100\% success rate.
        Back to the running example in \Cref{fig:dnn-graph-example}, our \sysname can generate ${\scriptsize\begin{bmatrix} 5 & -5 \\ -5 & 5 \end{bmatrix}}$ for node 2 and $(0.9, -0.9)^\T$ for node 4 as model weights $\vwinfer$; and $(10, -10)^\T$ for node 1 and $(1,0)^\T$ for node 11 as the inference input $\vxinfer$.
        Such $\vwinfer$ and  $\vxinfer$ induce input $(0,1)^\T$ and $(1,0)^\T$ for node 9 and 10 respectively.
        Since both node 9 and 10 use \CodeIn{log} operator and $\log 0$ is undefined, both node 9 and 10 trigger numerical failures.
        
        \emph{\textbf{Step b}: Generate training example $\vxtrain$ that achieves model weights $\vwinfer$.}
        To the best of our knowledge, there is no automatic approach for this task yet.
        \sysname provides support for this task based on our extension of DLG attack~\cite{zhu2020deep}.
        The DLG attack is originally designed for recovering the training data from training-phase gradient leakage.
        Here, we figure out the required training gradients to trigger the numerical failure at the inference phase and then leverage the DLG attack to generate $\vxtrain$ that leads to such training gradient.
        Specifically, many DNN architectures contain operators such as \CodeIn{ReLU} on which DLG attack is hard to operate~\cite{serra2018bounding}.
        We combine straight-through estimator~\cite{bengio2013estimating} to provide proxy gradients and bypass this barrier.
        Back to the running example in \Cref{fig:dnn-graph-example}, suppose the initial weights are ${\scriptsize\begin{bmatrix} -0.1 & 0.1 \\ 0.1 & -0.1 \end{bmatrix}}$ for node 2 and $(0,0)^\T$ for node 4, our \sysname can generate training example $\vxtrain$ composed of $(5.635, -5.635)^\T$ for node 1 and $(1,0)^\T$ for node 11, such that one-step training with learning rate 1 on this example leads to $\vwinfer$.
        Combining $\vxtrain$ from this step with $\vxinfer$ from \emph{step a}, we obtain a failure-exhibiting system test that confirms the feasibility of potential defects in node 9 and 10.
        
       
        
        \paragraph{Task \circled{3}: Fix Suggestion via Abstract Optimization}
        In this task, we suggest fixes for the confirmed numerical defects.
        \sysname is the first approach for this task to our knowledge.
        
        The user may prefer different fix locations, which correspond to a user-specified set of nodes $\Vfix \subseteq \gV$ to impose the fix.
        For example, if the fix method is clipping the inference input, $\Vfix$ are input nodes~(e.g., nodes 1, 11 in \Cref{fig:dnn-graph-example});
        if the fix method is clipping the model weights during training, $\Vfix$ are weight nodes~(e.g., nodes 2, 4 in \Cref{fig:dnn-graph-example});
        if the fix method is clipping before the suspicious operator,
       
        $\Vfix$ are nodes with numerical defects~(e.g., nodes 9, 10 in \Cref{fig:dnn-graph-example}).
            
        According to the empirical study of developers' fixes in \Cref{sec:intro}, 12 out of 13 defects are fixed by imposing interval preconditions clipping the inputs of $\Vfix$.
        Hence, we suggest interval precondition as defect fix in this paper, which is interval constraint $\vl_n \le \fin_{n}(\vx; \vw) \le \vu_n$ for nodes $n\in \Vfix$.
        A fix should satisfy that, when these constraints $\bigwedge_{n\in\Vfix} (\vl_n \le \fin_{n}(\vx; \vw) \le \vu_n)$ are imposed, the input of any node in the computational graph should always be valid, i.e., $\fin_{n_0}(\vx; \vw) \notin \gI_{n_0, \mathsf{invalid}}, \forall n_0 \in \gV$.
        
        In \sysname, we formulate the fix suggestion task as a constrained optimization problem, taking the endpoints of interval abstractions for nodes in $\Vfix$ as optimizable variables.
        We then propose the novel technique of abstraction optimization to solve this constrained optimization problem.
       
       
       
        Back to \Cref{fig:dnn-graph-example} example, if users plan to impose a fix on inference input, \sysname can suggest the fix $-1 \le $ \CodeIn{x-input} $\le 1$;
        if users plan to impose a fix on nodes with numerical defects, \sysname can suggest the fix $10^{-38} \le $ node 9 \& node 10\CodeIn{.input} $\le +\infty$.
        
    
   
   
   
    
   
        
   
    
   
   
    
   
   
   
   
   
        
   
   
   
   
   
   
   
   
            
   
   
   
   
            
   
   
   
   
   
        
   
   
   
   
   
   
   
   
   
            
   
   
   
   
   
   
   
   
            
   
   
            
   
            
   
   
   
   
   
                
   
   
   
   
   
   
            
            
   
   
   
   
   
            
   
   
   
   
   
   
   
   
   
   
   
   
   
   
            
   
   
        
   
   
        
   
   
   
   
   
   
   
   
        
   
   
        
   
   
   
   
            
   
   
        
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
            
   
   
   
   
   
                
   
   
   
   
   
   
   
                
   
   
   
    

\section{The \sysname Approach}

    \label{sec:method}
    
    In this section, we introduce the three novel techniques in \sysname: backward fine-grained node labeling in \Cref{subsec:static-analysis};
    two-step test generation in \Cref{subsec:system-test-gen}; 
    and abstraction optimization in \Cref{subsec:precond-gen}.

    \subsection{DNN Static Analysis Framework with Backward Fine-Grained Node Labeling for Potential-Defect Detection}
    
        \label{subsec:static-analysis}
        
        \sysname contains a static analysis framework to enable potential-defect detection and support downstream tasks as shown in \Cref{fig:overview}.
        Given a DNN architecture and valid ranges for input and weight nodes, the static analysis framework computes interval abstractions for possible inputs and outputs of each node. 
        As a result, we can check whether overlap exists between the interval abstraction and invalid input ranges for all nodes in the graph to detect potential numerical defects.
        Then, the defect nodes are fed into the two-step test generation component to confirm the feasibility of potential defects; 
        and the differentiable abstractions are fed into the abstract optimization component to produce fixes.
        
        Formally, for given valid ranges of inference input and model weights, namely $\gX$ and $\gW$, for each node $n\in\gV$, our framework will compute \emph{sound} input interval abstraction $[\vl_n, \vu_n] := \{\vx: \vl_n \le \vx \le \vu_n\}$ such that $[\vl_n, \vu_n]$ always captures all possible inputs of the node: $[\vl_n, \vu_n] \supseteq \{\fin_n(\vx, \vw): \vx\in\gX, \vw\in\gW\}$.
        We also compute output interval abstractions similarly.
        
        Compared with traditional analysis tools for numerical software~\cite{gurfinkel2015seahorn,singh2017fast}, \sysname static analysis framework designs abstractions for DNN primitives operating on multi-dimensional tensors that are not supported by traditional tools.
        Compared with the state-of-the-art DEBAR tool~\cite{zhang2020detecting}, \sysname uses the same abstraction domain~(interval domain with tensor partitioning), but incorporates a novel technique~(backward fine-grained node labeling) to improve abstraction precision and support a wider range of DNN architectures.
        
       
       
       
       
        
        
       
       
       
        
        
        
        \textbf{Abstract Domain: Interval with Tensor Partitioning.}
        Following DEBAR's design, we use the interval with tensor partitioning~\cite{zhang2020detecting} as the abstraction domain, which partitions the tensor into several subblocks and shares the interval abstractions at the block level instead of imposing at the element level.
        Therefore, we can compute the abstraction of a smaller size than the original tensor to improve efficiency.
       
       
       
       
       
        
       
        
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
        
       
       
       
       
       
       
       
       
       
       
       
        
       
       
       
       
       
       
       
       
       
        
        \paragraph{Our Technique: Backward Fine-Grained Node Labeling}
        The interval domain with tensor partitioning provides a degree of freedom in terms of the partition granularity, i.e., we can choose the subblock size for each node's abstraction.
        When choosing the finest granularity, i.e., elementwise abstraction, the abstraction interval is the most concrete.
        When choosing the coarsest granularity~(i.e., one scalar to summarize the node tensor), the abstraction saves the most space and computational cost but loses much precision.

        \emph{Example.}
        Suppose the possible input range of a node is $([-1,0], [0,1], [1,2], [-1,0])$, where each interval $[l,u]$ specifies the range of corresponding elements in the four-dimensional vector.
        If we choose the finest granularity, we will use $[\vl_n,\vu_n] = [(-1,0,1,-1), (0,1,2,0)]$ as the input range abstraction.
        If we choose the coarsest granularity, we will use $[\vl_n, \vu_n] = [-1,2]$ as the abstraction where the same interval is shared for all elements.
        As we can see, finer granularity provides tighter abstraction at the expense of larger computational and space costs.
        
        In DEBAR, the coarsest granularity is used by default for most operators.
        However, we find that using the finest instead of the coarsest granularity for some nodes is more beneficial for overall abstraction preciseness.
        For example, the control-flow operators~(e.g., \CodeIn{Loop}) benefit from concrete execution to determine the exact control flow in the dynamic graph, and the indexing operators~(e.g., \CodeIn{Slice}) and shaping operators~(e.g., \CodeIn{Reshape}) benefit from explicit indexers and shapes to precisely infer the output range.
        Hence, we propose to use the finest granularity for some nodes~(namely fine-grained requiring operators) while the coarsest granularity for other nodes during static analysis.
        
        
        
        To benefit from the finest granularity abstraction for required nodes, typically, all of their preceding nodes also need the finest granularity.
        Otherwise, the over-approximated intervals from preceding nodes will be propagated to the required nodes, making the finest abstraction for the required nodes useless.
        To solve this problem, in \sysname, we back-propagate ``fine-grain'' labels from these fine-grain requiring nodes to initial nodes by topologically sorting the graph with \emph{inverted} edges, and then apply the finest granularity abstractions on all labeled nodes.
        In practice, we find this strategy eliminates the control-flow ambiguity and indexing ambiguity with little loss of efficiency\footnote{Theoretically, using the finest granularity for tensor partitioning cannot fully eliminate the ambiguity, since interval abstraction is intrinsically an over-approximation. Nevertheless, in practice, we found this technique eliminates control-flow and indexing ambiguities on all 63 programs in the benchmark.}.
        As a result, \sysname supports all dynamic graphs~(which are not supported by DEBAR) that compose of 39.2\% of the benchmarks proposed by \citeauthor{yan2021exposing}~\cite{yan2021exposing}.
        
        Furthermore, when preceding nodes use finer-grain abstraction granularity, the subsequent nodes should preserve such fine granularity to preserve the analysis preciseness.
        Principally, the choice of abstraction granularity should satisfy both tightness~(bears no precision loss compared to elementwise interval abstraction) and minimality~(uses the minimum number of partitions for high efficiency).
        To realize these principles, we dynamically determine the node's abstraction granularity based on the granularity of preceding nodes.
        The abstraction design for some operators is non-trivial.
        Omitted details~(formulation, illustration, and proofs) about the static analysis framework are in \Cref{adxsec:static-analysis}.
        
       
       
       
       
       
       
       
        
        In summary, the whole static analysis process consists of three steps: 
        (1)~Determine the tensor partition granularity of all initial nodes by our technique of backward fine-grained node labeling.
        (2)~Sort all nodes in the graph in topological order.
        (3)~Apply corresponding abstraction computation algorithms for each node based on the preceding node's abstractions.
       
       
       
        The key insight behind the design of our static analysis framework is the strategic granularity selection for tensor abstraction, which maintains both high efficiency by selecting the coarse granularity for data-intensive nodes~(which are particularly important for DNNs due to the high-dimensional and uniform range nature of their data nodes) and high precision by selecting the fine granularity for some critical nodes~(nodes with control-flow, indexing, and shaping operators).
        
    \subsection{Two-Step Test Generation for Feasibility Confirmation}
    
        \label{subsec:system-test-gen}
    
        \sysname generates failure-exhibiting system tests for the given DNN to confirm the feasibility of potential numerical defects.
        Here, we take the DNN architecture as the input.
        From the static analysis framework, we obtain a list of nodes that have potential numerical defects.
        For each node $n_0$ within the list, we apply our technique of two-step test generation to produce a failure-exhibiting system test $\tsys = \langle \vxtrain, \vxinfer \rangle$ as the output.
        According to \Cref{subsec:approach-overview}, the test should satisfy that after the architecture is trained with $\vxtrain$, entering $\vxinfer$ in the inference phase results in numerical failure.
        
        We propose the novel technique of two-step test generation:
        First, generate failure-exhibiting unit test $\langle \vwinfer, \vxinfer \rangle$;
        Then, generate training example $\vxtrain$ that leads model weights to be close to $\vwinfer$ after training.
       
        
        \paragraph{\textbf{Step a}: Unit Test Generation}
            As sketched in \Cref{subsec:approach-overview}, we strengthen the state-of-the-art unit test generation approach, GRIST~\cite{yan2021exposing}, by combining with random initialization to complete this step.
            Specifically, GRIST leverages the gradients of the defect node's input with respect to the inference input and weights to iteratively update the inference input and weights to generate failure-exhibiting unit tests.
            However, GRIST always conducts updates from existing inference input and weights, which can suffer from local minima problem~\cite{madry2018towards}.
            Instead, motivated by DNN adversarial attack literature~\cite{madry2018towards,tramer2020adaptive}, a sufficient number of random starts helps to find global minima effectively.
            Hence, in \sysname, we first conduct uniform sampling 100 times for both the inference input and weights to trigger the numerical failure.
            If no failure is triggered, we use the sample that induces the smallest loss as the start point for gradient optimization.
            As \Cref{subsec:static-detector-exp} shows, this strategy substantially boosts the efficiency, achieving 17.32X speed-up.
            
           
           
            
           
           
           
           
           
           
        
        \paragraph{\textbf{Step b}: Training Example Generation}
            For this step, \sysname takes the following inputs: (1)~the DNN architecture; (2)~the failure-exhibiting unit test $\tunit = \langle \vwinfer, \vxinfer \rangle$; and (3) the randomly initialized weights $\vw_0$.
            Our goal is to generate a legal training example $\vxtrain$, such that the model trained with $\vxtrain$ will contain weights close to $\vwinfer$.
           
           
            
            
            DNNs are typically trained with gradient-descent-based algorithms such as stochastic gradient descent~(SGD).
            In SGD, in each step $t$, we sample a mini-batch of samples from the training dataset to compute their gradients on model weights and use them to update the weights.
            We focus on one-step SGD training with a single training example, since generating a single one-step training example to exhibit failure is more desirable for  debugging because, in one-step training, the model weights are updated strictly following the direction of the gradient. 
            Therefore, developers can inspect inappropriate weights, easily trace back to nodes with inappropriate gradients, and then fix these nodes. 
            In contrast, in multi-step training, from inappropriate weights, developers cannot trace back to inappropriate gradients because weights are updated iteratively and interactions between gradients and weights are complex~(even theoretically intractable~\cite{li2017convergence}).
            
            In this one-step training case, after training, the model weights $\vwinfer$ satisfies
            \vspace{-0.5em}
            \begin{equation}
               
                \vwinfer = \vw_0 - \gamma \nabla_\vw \gL(\vxtrain; \vw_0),
                \vspace{-0.5em}
                \label{eq:one-step-sgd}
            \end{equation}
            where $\gamma \in \sR_+$ is a predefined learning rate, and $\gL$ is the predefined loss function in the DNN architecture.
            Hence, our goal becomes finding $\vxtrain$ that satisfies 
            \vspace{-0.3em}
            \begin{equation}
               
                \nabla_\vw \gL(\vxtrain; \vw_0) = (\vw_0 - \vwinfer) / \gamma. \label{eq:dlg-attack}
            \end{equation}
            \vspace{-1.5em}
            
            
           
           
           
           
           
           
           
           
           
           
           
           
           
           
            
           
           
           
           
           
           
           
           
           
           
           
           
           
            
            The DLG attack~\cite{zhu2020deep} is a method for generating input data that induce specific weight gradients.
            The attack is originally designed for recovering training samples from monitored gradient updates.
           
            Since the right-hand side~(RHS) of \Cref{eq:dlg-attack} is known, our goal here is also to generate input example $\vxtrain$ that induces specific weight gradients.
            Therefore, we leverage the DLG attack to generate training example $\vxtrain$.
           
            
            \paragraph{Extending DLG Attack with Straight-Through Estimator}
            Directly using DLG attack suffers from an optimization challenge in our scenario.
            Specifically, in DLG attack, suppose the target weight gradients are $\vdwtarg$, we use gradient descent over the squared error $\| \nabla_\vw \gL(\vx; \vw_0) - \vdwtarg \|_2^2$ to generate $\vx$.
            In this process, we need meaningful gradient information of this squared error loss to perform the optimization.
            However, the gradient of this loss involves second-order derivatives of $\gL(\vx;\vw_0)$, which could be zero. 
            For example, DNNs with $\relu$ as activation function are piecewise linear and have zero second-order derivatives almost everywhere~\cite{serra2018bounding}.
            This optimization challenge is partly addressed in DLG attack by replacing \CodeIn{ReLU} with \CodeIn{Sigmoid}, but it changes the DNN architecture~(i.e., the system under test) and hence is unsuitable.
            
            We leverage the straight-through estimator to mitigate the optimization challenge.
            Specifically, for some operators, such as \CodeIn{ReLU}, we do not change its forward computation but change its backward gradient computation to provide second-order derivatives within the DLG attack process.
            For example, for \CodeIn{ReLU}, in backward computation we use the gradient of $\mathrm{Softplus}$ function, namely $1 - \frac{1}{1+\exp(\vx)}$, 
            because $\mathrm{Softplus}$ is an approximation of \CodeIn{ReLU}~\cite{glorot2011deep} with non-zero second-order derivatives.
           
            Note that we modify the computed gradients only within the DLG attack process.
            After such $\vxtrain$ is generated by the attack, we evaluate whether it triggers numerical failure using the original architecture and gradients in \Cref{eq:one-step-sgd}.
           

            \Cref{adxsec:detail-hyperparam} lists hyperparameters used by our implementation.
            
           
            
           
           
           
           
           
           
           
        
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
       
    
       
    \subsection{Abstraction Optimization for Fix Suggestion}
    
        \label{subsec:precond-gen}
        
        In this task, we aim to generate the precondition fix given imposing locations.
        The inputs are the DNN architecture, the node $n_0$ with numerical defects, and a node set $\Vfix$ to impose the fix.
        We would like to generate interval preconditions for $\Vfix$ node inputs so that after these preconditions are imposed, the defect on $n_0$ is fixed.
        
        Formally, our task is to find $\langle l_n, u_n \rangle$ for each $n\in \Vfix$~($l_n$ and $u_n$ are scalars so the same interval bound applied to all elements of $n$'s tensor), such that for any $\vx, \vw$ satisfying $\fin_{n}(\vx;\vw) \in [l_n, u_n]$, $\forall n\in \Vfix$, for the defect node $n_0$, we have $\fin_{n_0}(\vx;\vw) \notin \gI_{n_0,\textsf{invalid}}$, where the full list of invalid input ranges $\gI_{n_0, \textsf{invalid}}$ is in \Cref{adxsec:list-numerical-defects}.
        There is an infinite number of possible $\langle l_n, u_n \rangle$ interval candidates since $l_n$ and $u_n$ are floating numbers.
        Hence, we need an effective approach to find a valid solution from the exceedingly large search space that incurs a relatively small model utility loss.
        To achieve so, we formulate a surrogate optimization problem for this task.
        \vspace{-0.5em}
        \allowdisplaybreaks
        \begin{small}
        \begin{align}
            \maximize_{l_n, u_n: n\in \Vfix} \quad s \quad
           
            \mathrm{s.t.}\quad & u_n \ge l_n + s(u_n^{\valid} - l_n^{\valid}), \forall n\in \Vfix, \label{eq:precond-span-cons} \\
            & l_n^{\valid} \le l_n \le u_n \le u_n^{\valid}, \forall n\in \Vfix, \label{eq:precond-valid-cons} \\
            & \gL_{n_0}^{\precondgen}(\{l_n, u_n\}_{n\in\Vfix}) < 0. \label{eq:precond-effective-cons}
           
        \end{align}
        \end{small}
        Here, 
        $l_n^{\valid}$ and $u_n^{\valid}$ are the valid ranges of the node's input $n$, which is fixed and determined by the valid ranges of input and weights.
        $\gL_{n_0}^{\precondgen}$ is the node-specific precondition generation loss that is the distance between the furthest endpoint of defect node $n_0$'s interval abstraction and $n_0$'s valid input range.
        Hence, when $\gL_{n_0}^{\precondgen}(\{l_n, u_n\}_{n\in\Vfix})$ becomes negative, the solution $\{l_n, u_n\}_{n\in \Vfix}$ is a valid precondition.
        The optimization variables are the precondition interval endpoints $l_n$ and $u_n$ and the objective is the relative span of these intervals.
        The larger the span is, the looser the precondition constraints are, and the less hurt they are for the model's utility.
        \Cref{eq:precond-span-cons} enforces the interval span requirement.
        \Cref{eq:precond-valid-cons} assures that the precondition interval is in the valid range.
        \Cref{eq:precond-effective-cons} guarantees the validness of the precondition as a fix.
       
       
       
       
        
        For any $\{l_n, u_n\}_{n\in \Vfix}$, thanks to \sysname static analysis framework, we can compute induced intervals of defect node $n_0$, and thus compute the loss value $\gL_{n_0}^{\precondgen}$.
       
        
        \setlength{\textfloatsep}{0pt
        \begin{algorithm}[!t]
            \caption{Abstraction Optimization~(\Cref{subsec:precond-gen})}
            \label{alg:abs-opt}
            \begin{algorithmic}[1]
                \small
                \renewcommand{\algorithmicrequire}{\textbf{Input:}}
                \renewcommand{\algorithmicensure}{\textbf{Output:}}
                \REQUIRE DNN architecture $\gG=\langle \gV, \gE \rangle$, defect node $n_0 \in \gV$, nodes to impose fix $\Vfix \subseteq \gV$
               
               
                \STATE $s\gets 1, \gamma_s\gets 0.9, \gamma_c\gets 0.1, \minstep\gets 0.1, \maxiter\gets 1000$ 
                \STATE $c_n \gets (l_n^{\valid} + u_n^{\valid}) / 2, l_n \gets l_n^{\valid}, u_n \gets u_n^{\valid},\,  \forall n \in \Vfix$
                \FOR {$i = 1$ to $\maxiter$}
                    \FOR {$n\in\Vfix$}
                        \STATE \textsf{loss} $\gets \gL_{n_0}^{\precondgen}(\{l_{n'}, u_{n'}\}_{n'\in\Vfix})$
                        \STATE $c_n \gets c_n - \gamma_c \max\{|c_n|,\minstep\}\sgn(\nabla_{c_n} \mathsf{loss})$
                        \STATE $(l_n, u_n) \gets (c_n - \frac{s(u_n^{\valid} - l_n^{\valid})}{2}, c_n + \frac{s(u_n^{\valid} - l_n^{\valid})}{2})$
                        \STATE $(l_n, u_n) \gets (\max \{l_n, l_n^{\valid}\}, \min \{u_n, u_n^{\valid}\})$ 
                       
                       
                    \ENDFOR
                    \IF {$\gL_{n_0}^{\precondgen}(\{l_n, u_n\}_{n\in\Vfix}) < 0$}
                        \RETURN $\{l_n, u_n\}_{n\in\Vfix}$ \COMMENT{Find precondition fix}
                    \ENDIF
                    \STATE $s \gets \gamma_s \cdot s$
                \ENDFOR
                \RETURN ``failed'' \COMMENT{Failed to find precondition fix}
            \end{algorithmic}
        \end{algorithm}

        As shown in \Cref{alg:abs-opt}, we propose the technique of \textbf{abstraction optimization} to effectively and approximately solve this optimization.
        Our technique works iteratively.
        In the first iteration, we set span $s=1$, and in the subsequent iterations, we reduce the span $s$ exponentially as shown in line 13 where hyperparameter $\gamma_s = 0.9$.
        Inside each iteration, for each node to impose precondition $n\in\Vfix$, we use the interval center $c_n = (l_n + u_n)/2$ as the optimizable variable and compute the \emph{sign} of its gradient: $\sgn(\nabla_{c_n} \mathsf{loss})$.
        We use this gradient sign to update each $c_n$ towards reducing the loss value in line 6.
       
        Then, we use $c_n$ and the span $s$ to recover the actual interval in line 7 and clip $l_n$ and $u_n$ by the valid range $[l_n^{\valid}, u_n^{\valid}]$ in line 8.
        At the end of this iteration, for updated $l_n$ and $u_n$, we compute $\gL_{n_0}^{\precondgen}(\{l_n, u_n\}_{n\in\Vfix})$ to check whether the precondition is a fix.
        If so, we terminate; otherwise, we proceed to the next iteration.
        We note that \emph{if the algorithm finds a precondition, the precondition is guaranteed to be a valid fix} by the soundness nature of our static analysis framework and the definition of $\gL_{n_0}^{\precondgen}$.
        When no feasible precondition is found within $\maxiter = 1000$ iterations, we terminate the algorithm and report ``failed to find the fix''.
       
        
        \begin{remark}
            The key ingredient in the technique is the gradient-sign-based update rule shown in line 6 which is much more effective than normal gradient descent for two reasons. 
            (1)~Our update rule can get rid of gradient explosion and vanishing problems.
            For early optimization iterations, the span $s$ is large and interval bounds are generally coarse, which results in too large or too small gradient magnitude.
            For example, the input range for \CodeIn{Log} could be $[1,10^{10}]$ where gradient can be $10^{-10}$, resulting in almost negligible gradient updates.
            In contrast, our update rule leverages the sign of gradient which always points to the correct gradient direction.
            The update step size in our rule is the maximum of current magnitude $|c_n|$ and $\minstep$ to avoid stagnation.
            (2)~Our update rule mitigates the gradient magnitude discrepancy of different $c_n$.
            At different locations, the nodes in DNNs can have diverse value magnitudes that are not aligned with their gradient magnitudes, which makes gradient optimization challenging.
           
           
           
           
            Therefore, we use this update rule to solve the challenge, where the update magnitude depends on the value magnitude ($|c_n|$) instead of gradient magnitude ($\nabla_{c_n} \mathsf{loss}$).
            We empirically compare our technique with standard gradient descent in \Cref{subsec:precond-gen-exp}.
        \end{remark}
        
    \input{exp/tables/unittest_table}
    \vspace{-1em}
\section{Implementation}

    \label{sec:impl}
    
    We have implemented a tool for \sysname in roughly 10k lines of Python code based on PyTorch.
    Our tool leverages existing modules~(\CodeIn{tf2onnx} for PyTorch and \CodeIn{torch.onnx.export} for Tensorflow and Keras) to extract ONNX-format DNN architectures from DL programs to take as the input, and automatically generates detection result, system tests, and defect fixes~(given imposing location).
   
   
   
   
    
   
\section{Experimental Evaluation}

    \label{sec:exp}

    We conduct a systematic experimental evaluation to answer the following research questions.
    
   
   
    
   
   
        
   
   
    
   
    \begin{enumerate}[leftmargin=*,label={\textbf{RQ\arabic*}}]
        \item
        For tasks already supported by existing state-of-the-art~(SOTA)  tools~(task \circled{1} and \circled{2}a), how much more effective and efficient is \sysname compared to these SOTA  tools?
    
        \item 
        For feasibility confirmation via \emph{generating failure-exhibiting system tests}~(task \circled{2}), how much more effectively and efficiently can \sysname confirm potential numerical defects compared to baseline approaches?
        
        \item 
        For \emph{suggesting fixes}~(task \circled{3}), how much more efficient and effective is \sysname in terms of guarding against numerical failures compared to baseline approaches and developers' fixes, respectively?
       
       
        
    \end{enumerate}
   
    
    For RQ1, we compare \sysname with all SOTA tools.
   
    For RQ2 and RQ3, \sysname is the first approach to the best of our knowledge, so we compare \sysname with baseline approaches (constructed by leaving our novel techniques out of  \sysname) and developers' fixes.
    We conduct the evaluation on the GRIST benchmarks~\cite{yan2021exposing}, being the largest dataset of real-world DNN numerical defects to our knowledge. The benchmarks contain 63 real-world DL programs with numerical defects collected from previous studies and GitHub.
    Each program contains a DNN architecture, and each architecture has one or more numerical defects.
    There are 79 true numerical defects in total.
    
   
   
    
   
   
   
   
   
   
   
   
    
    
   
   
    
    We perform all experiments on a Linux workstation with a 24-core Xeon E5-2650 CPU running at 2.20 GHz.
    Throughout the evaluation, we stop the execution after reaching \SI{30}{min} limit by following the evaluation setup by the most recent related work~\cite{yan2021exposing}.
   
   
    
    
   
    \subsection{RQ1: Comparison with SOTA Tools}
        \label{subsec:static-detector-exp}
        \label{subsec:compare-with-sota}
        
       
        
       
        For two tasks, existing tools can provide automatic support: potential-defect detection~(task \circled{1}) where the SOTA tool is DEBAR~\cite{zhang2020detecting}, and failure-exhibiting unit test generation~(task \circled{2}a) where the SOTA tool is GRIST~\cite{yan2021exposing}.
        We compare \sysname with these tools on their supported tasks respectively.
        
        \paragraph{Comparison with DEBAR}
        \sysname successfully detects all 79 true defects and DEBAR only detects 48 true defects according to both our evaluation and the literature~\cite{yan2021exposing}.
        Hence, \sysname detects 64.58\% more true defects than DEBAR.
        In terms of efficiency, DEBAR and \sysname have similar running time, and both finish in \SI{3}{s} per case.
        
       
       
       
       
       
        
        We manually inspect the cases where DEBAR fails but \sysname succeeds.
        They correspond to DL programs written with the PyTorch library, which generates dynamic computational graphs that DEBAR cannot handle.
        In contrast, \sysname provides effective static analysis support for dynamic computational graphs thanks to our backward fine-grained node labeling technique~(see \Cref{subsec:static-analysis}) that is capable of disambiguating the control flow within dynamic graphs.
        
        
       
       
       
        
        \input{exp/tables/systest_table}
        
        \paragraph{Comparison with GRIST}
        Results are shown in \Cref{tab:unittest}.
        Since both \sysname and GRIST have a randomness component where \sysname uses random initialization and GRIST relies on DNN's randomly initialized weights, we repeat both approaches for 10 runs, record the total number of times where a failure-exhibiting unit test is generated, and the average execution time per run.
       
       
        \sysname succeeds in \emph{all} cases and \emph{all} repeated runs, and GRIST fails to generate unit test in 24 out of 790 runs~(i.e., 96.96\% success rate).
        \sysname has \SI{6.66}{s} average execution time and is 17.32X faster than GRIST.
        
        The superior effectiveness and efficiency of \sysname are largely due to the existence of random initialization as introduced in \Cref{subsec:system-test-gen}.
        We observe that since GRIST always takes initial model weights and inference input as the starting point to update from, the unit test is just slightly changed from initial weights and input, which can be hard to trigger the oblivious numerical defect.
        In contrast, \sysname uses random initialization to explore a much larger space and combines gradient-based optimization to locate the failure-exhibiting instances from the large space.
       
        We also evaluated the pure random strategy that only uses random initialization without gradient-based optimization, and such strategy fails in 30 runs, being inferior to both \sysname and GRIST, which means both random  initialization and gradient-based optimization are important.
        Among all 79 cases, \sysname is slower than GRIST on only one case~(28a).
        For this case, we find the default inference input loaded by the DNN program~(used by GRIST) is not far from the failure-exhibiting one, but a randomly-sampled inference input~(used by \sysname) is usually far from that.
        Hence, \sysname takes more iterations to find the failure-exhibiting inference input by the nature of gradient-based optimization.
       
        
        
    
   
    \subsection{RQ2: Feasibility Confirmation via System Test Generation}
   
    
        \label{subsec:system-test-gen-exp}
        
       
        
        
       
        
        In task \circled{2}, \sysname confirms feasibility of potential numerical defects by generating failure-exhibiting system tests.
        
        \paragraph{Baseline}
        Since \sysname is the first approach for this task,
       
        we did not compare with existing literature and propose one random-based approach~(called ``Random'' hereinafter)  as the baseline.
        In ``Random'', we first generate failure-exhibiting unit test with random sampling.
        If there is any sample that triggers the failure, we stop and keep the inference example part as the desired $\vxinfer$.
        Then, we generate $\vxtrain$ again by random sampling.
        If any sample, when used for training, could induce model weights $\vw$ that cause numerical failure when using $\vxinfer$ as the inference input, we keep that sample as $\vxtrain$ and terminate.
        If and only if both $\vxinfer$ and $\vxtrain$ are found, we count this run as a ``success'' one for ``Random''.
        
        For each defect, due to the randomness of model's initial weights, we repeat both \sysname and ``Random'' for $10$ runs.
        Both approaches use the same set of random seeds.
       
       
        
        \paragraph{Evaluation Result}
        Results are in \Cref{tab:systest}.
        We observe that \sysname succeeds in $733 / (79\times 10) = 92.78\%$ runs and the baseline ``Random'' succeeds in $649 / (79\times 10) = 82.15\%$ runs.
        Moreover, \sysname spends only \SI{17.31}{s} time on average per run which is a 19.30X speed-up compared to ``Random''.
        We also observe that \sysname is more reliable across repeated runs. There are only 6 cases with unsuccessful repeated runs in \sysname, but there are 19 such cases in ``Random''. 
       
        Hence, \sysname is substantially more effective, efficient, and reliable for generating system tests than the baseline.
        
       
        \paragraph{Discussion}
        The high effectiveness of \sysname mainly comes from the advantage of gradient-guided search compared with random search.
        As described in \Cref{subsec:system-test-gen}, \sysname leverages both first-order gradients~(in step a) and second-order derivatives~(in step b) to guide the search of system tests.
       
        In contrast, ``Random'' uses random sampling hoping that failure-exhibiting training examples can reveal after sufficient sampling.
        Hence, when such training examples are relatively rare in the whole valid input space, ``Random'' is less effective.
        We conduct an ablation study in \Cref{adxsec:system-testgen-by-stage} which shows \sysname improves over ``Random'' in both steps inside \sysname.
        
        \paragraph{Failing Case Analysis}
        We studied all six defects where \sysname may fail.
        We found that 
       
       
       
       
       
       
        (1)~four defects~(Case ID 1, 15, 37, and 38) are because the architecture is challenging for gradient-based optimization, e.g., due to \CodeIn{Min}/\CodeIn{Max}/\CodeIn{Softmax} operators which provides little or no gradient information.
        We leave it as future work to solve these cases, which may need to dynamically detect operators with vanishing gradients and reconstruct the gradient flow.
        (2)~Two defects~(Case ID 28a and 28b) correspond to those that caused by \CodeIn{Div} operators where only close-to-zero divisor can trigger the numerical failure. 
        Hence, for operators with narrow invalid ranges, RANUM may fail to generate failure-exhibiting system tests. 
        
       
    
   
   
    
   
        
   
        
   
        
   
        
        
   
   
   
   
   
        
   
   
   
   
   
   
   
   
        
   
   
   
   
   
        
   
   
   
   
   
   
        
   
   
   
   
   
   
   
    
   
    \subsection{RQ3: Fix Suggestion}
    
        \label{subsec:precond-gen-exp}
        
        At task \circled{3}, \sysname suggests fixes for numerical defects.
        We compare \sysname with fixes generated by baseline approaches and developers' fixes.
        
        \subsubsection{\textbf{Comparison between \sysname and Baselines}}
            \input{exp/tables/precond_summary_table}
            
           
            \sysname is the first approach for this task,
           
            and we propose two baseline approaches to compare with. 
            (1)~``\sysname-E'': this approach changes the abstraction domain of \sysname from interval with tensor partitioning to standard interval.
            To some degree, \sysname-E represents the methodology of conventional static analysis tools that use standard interval domain for abstraction and search of effective fixes. 
           
            (2)~``GD'': this approach uses standard gradient descent for optimization instead of the abstraction optimization technique in \sysname.
            
           
            \paragraph{Evaluation Protocol}
            We evaluate whether each approach can generate fixes that eliminate \emph{all} numerical defects for a given DNN architecture given imposing locations.
           
            We consider three types of locations: on both weight and input nodes, on only weight nodes, and on only input nodes.
            In practice, model providers can impose fixes on weight nodes by clipping weights after the model is trained; and users can impose fixes on input nodes by clipping their inputs before loading them into the model.
            Since all approaches are deterministic, for each case we run only once.
            We say that the fix eliminates all numerical defects if and only if (1)~the \sysname static analysis framework cannot detect any defects from the fixed architecture; and (2)~1,000 random samples cannot trigger any numerical failures after imposing the fix.
           
            
            \paragraph{Evaluation Result}
            We report the statistics, including number of successful cases among all 79 cases and the total running time, in \Cref{tab:precond-summary}.
           
            From the table, we observe that on all the three imposing location settings, \sysname always succeeds on the most cases and spends much less time.
           
            For example, when fixes can be imposed on both weights and input nodes, \sysname succeeds on \emph{all} cases with a total running time \SI{54.23}{s}.
            In contrast, \sysname-E requires $>10 \times$ time, and GD succeeds in only $72.15\%$ cases.
            Hence, \sysname is substantially more effective and efficient for suggesting fixes compared to baseline approaches.

            Since \sysname is based on iterative refinement~(see \Cref{alg:abs-opt}), we study the number of iterations needed for finding the fix.
            When fixes can be imposed on both weights and input nodes, where \sysname succeeds on all 79 cases, the average number of iterations is $29.80$, the standard deviation is $14.33$, the maximum is $53$, and the minimum is $2$.
            Hence, when \sysname can find the fix, the number of iterations is small, which coincides with the small total running time \SI{54.23}{s}.
            
           
            \paragraph{Discussion}
            The two baseline approaches can be viewed as ablated versions of \sysname.
            Comparing \sysname and GD, we conclude that the technique of abstraction optimization substantially improves the effectiveness and also helps the efficiency.
            Comparing \sysname and \sysname-E, we conclude that the interval abstraction with tensor partitioning as the abstraction domain substantially improves the efficiency and also helps the effectiveness.
            
            From \Cref{tab:precond-summary}, it is much easier to find the fix when imposing locations are weight nodes compared to input nodes.
           
            Since model providers can impose fixes on weights and users impose on inputs, this finding implies that fixing numerical defects on providers' side may be more effective than on end-users' side.
        
        
           
        \subsubsection{\textbf{Comparison between \sysname and Developers' Fixes}}
                
            We conduct an empirical study to compare the fixes generated by \sysname and by the developers.
            
           
            \paragraph{Evaluation Protocol}
            We manually locate GitHub repositories from which the GRIST benchmarks are constructed.
            Among the 79 cases, we find the repositories for 53 cases on GitHub and we study these cases.
            We locate the developers' fixes of the numerical defects by looking at issues and follow-up pull requests.
           
            Since \sysname suggests different fixes for different imposing locations, for each case we first determine the imposing locations from the developer's fix, then compare with \sysname's fix for that location.

            \sysname fixes are on the computational graph and developers’ fixes are in the source code, so we determine to conduct a code-centered comparison: 
            \sysname fixes are considered feasible only when the fixes can be easily implemented by code~(within 10 lines of code) given that developers’ fixes are typically small, usually in 3-5 lines of code. 
            In particular, our comparison is based on two criteria: 
            (1)~which fix is sound on any valid input; 
            (2)~if both are sound, which fix hurts less to model performance and utility (based on the span of imposed precondition, the larger span the less hurt). 
            Two authors independently classified the comparison results for each case and discussed the results to reach a consensus.
            
            \paragraph{Results}
            We categorize the comparison results 
           
            as below.
            \begin{enumerate}[leftmargin=*,label=\Alph*]
                \item \emph{(30 cases) Better than developers' fixes or no available developer's fix.}
                Developers either propose no fixes or use heuristic fixes, such as reducing the learning rate or using the mean value to reduce the variance.
                These fixes may work in practice but are unsound, i.e., cannot rigorously guarantee the elimination of the numerical defect for any training or inference data.
                In contrast, \sysname generates better fixes since these fixes rigorously eliminate the defect.
               
               
                
                \item \emph{(7 cases) Equivalent to developers' fixes.}
                Developers and \sysname suggest equivalent or highly similar fixes.
               
                
                \item \emph{(13 cases) No need to fix.}
                For these cases, there is no need to fix the numerical defect in the given architecture.
                There are mainly three reasons. 
                (1)~The DNN is used in the whole project with fixed weights or test inputs.
                As a result, though the architecture contains defects, it will not lead to system failure.
                (2)~The architecture is injected a defect as a test case for automatic tools, such as a test architecture in \CodeIn{TensorFuzz}~\cite{odena2019tensorfuzz} repository.
                (3)~The defect can be hardly triggered in practice. 
                For example, the defect is in a \CodeIn{Div} operator where the divisor needs to be very close to zero to trigger divide-by-zero failure, but this hardly happens in practice since the divisor is randomly initialized.
                
                \item \emph{(3 cases) Inferior than developers' fixes or \sysname-generated fixes are impractical.}
                In two cases, \sysname-generated fixes are inferior than developers' fixes.
                Human developers observe that the defect operator is \CodeIn{Log}, and its input is non-negative.
                So they propose to add $10^{-6}$ to the input of \CodeIn{Log} as the fix.
                In contrast, \sysname can only generate a clipping-based fix, e.g., clipping the input if it is less than $10^{-6}$.
                When the input is small, \sysname's fix interrupts the gradient flow from output to input while human's fix maintains it.
                As a result, human's fix does less hurt to model's trainability and is better than \sysname's fixes.
                In one other case, the \sysname-generated fix imposes a small span for some model weights~(less than $0.1$ for each component of that weight node).
                Such a small weight span strongly limits the model's expressiveness and utility.
               
                We leave it as the future work to solve these limitations.
            \end{enumerate}
           
            
           
           
            From the comparison results, we can conclude that for the 40 cases where numerical defects are needed to be fixed~(exclude case C),  \sysname suggests equivalent or better fixes than human developers in 37 cases.
            Therefore, \sysname is comparably effective as human developers in terms of suggesting numerical defect fixes, and is much more efficient since \sysname is an automatic approach.
            
           
           
           
           
           
           
           
           
           
           
           

            \paragraph{Guidelines for Users}
            We discuss two practical questions for \sysname users.
            (1)~\emph{Does \sysname hurt model utility, e.g., inference accuracy?}
            If no training or test data ever exposes a numerical defect, \sysname will not confirm a defect and hence no fix will be generated and no hurt to the utility. 
            If \sysname confirms numerical defects, whether the fix affects the utility depends on the precondition-imposing location.
            If imposing locations can be freely selected, \sysname tends to impose the fix right before the vulnerable operator, and hence the fixes do not reduce inference performance. 
            The reason is that the fix will change~(by clipping) the input only when the input falls in the invalid range of the vulnerable operator. 
            In practice, if the imposing locations cannot be freely selected and need to follow developers’ requirements, our empirical study above shows that, in only 3 out of 40 cases, our fixes incur larger hurt to the inference or training performance of the architecture.
            (2)~\emph{Should we always apply \sysname to fix any architecture?}
            We can always apply \sysname to fix any architecture since \sysname fixes do not visibly alter the utility in most cases. 
            Nonetheless, in deployment, we recommend first using \sysname to confirm defect feasibility. 
            If there is no such failure-exhibiting system test, we may not need to fix the architecture; otherwise, we use \sysname to generate fixes.
        
    
   
    
   
    
   
        
        
   
    
   
    
   
\section{Threats to Validity}

    An external threat to validity is the evaluation subject.
    We evaluate \sysname and baseline approaches on the GRIST benchmarks~\cite{yan2021exposing}.
    Though the dataset contains 63 real-world DL programs and is the largest to our knowledge, it may still not be representative enough.
    Another external threat is the approach randomness.
    To mitigate this threat, we repeat all randomized approaches for 10 runs.
    To mitigate bias in the empirical study, two authors independently conducted the study and discussed all cases to reach a consensus.
    
    A major internal threat to validity comes from the approach implementation.
    We reduce this threat by conducting code reviews and checking logs extensively.
    An author independently checked the correctness of code implementation.
    We also verified the soundness of our static analysis framework with over 50 carefully designed unit tests.
    
   
\section{Related Work}
   

   
   

   
    \paragraph{Understanding and Detecting Defects in DNNs}
    Discovering and mitigating the defects, faults, bugs, and failures in DNN based systems is an important research topic~\cite{zhang2018empirical,pham2020problems,kloberdanz2022deepstability}.
    Following the taxonomy in \cite{humbatova2020taxonomy}, DNN defects are at four levels from bottom to top:
    (1)~Platform-level defects: Defects can exist in real-world DL compilers and libraries. 
    Approaches exist for understanding, detecting, and testing against these defects~(e.g., \cite{wang2020deep, DBLP:conf/sigsoft/ShenM0TCC21,tambon2021silent,xie2021leveraging}).
    (2)~Architecture-level defects: \emph{Our work focuses on numerical defects, which is one type of architecture-level defects.} 
    Automatic detection and localization techniques~(e.g., \cite{wardat2021deeplocalize,liu2021detecting}) exist for other architecture-level defects such as suboptimal structure, activation function, initialization and shape mismatch~\cite{humbatova2020taxonomy}.
   
   
    (3)~Model-level defects: Once a model is trained, its defects can be viewed as violations of desired properties which are discussed in detail in \cite{zhang2020machine}.
    Some example defects are correctness~\cite{tizpaz2020detecting,guerriero2021operation}, robustness~\cite{wang2021robot}, and fairness~\cite{zhang2020white} defects.
    (4)~Interface-level defects: DNN-based systems, when deployed as services, expose interaction interfaces to users where defects may exist as empirical studies on real-world systems reveal~\cite{humbatova2020taxonomy,wan2021machine,wan2022automated}.
    
   
    \paragraph{Testing and Debugging for DNNs}
    A rich body of work exists for testing and debugging DNN defects~\cite{zhang2020machine}.
   
   
    Some representatives are DeepXplore~\cite{pei2017deepxplore} and DeepGauge~\cite{ma2018deepgauge}.
   
    Recent works enable automatic model debugging and repair via consistency checking~\cite{xiao2021self}, log checking~\cite{zhang2021autotrainer}, spectrum analysis~\cite{qi2021archrepair}, or analyzer-guided synthesis~\cite{sotoudeh2021provable}.
    
   
    \paragraph{DNN Static Analysis}
    Another solution for eliminating DNN defects is conducting static analysis to rigorously guarantee the non-existence of defects~\cite{li2020sok,albarghouthi2021introduction}.
    Though DNNs essentially conduct numerical computations, traditional numerical analysis tools~\cite{gurfinkel2015seahorn,singh2017fast} are inefficient for DNN analysis due to lack of support for multi-dimensional tensor computations.
    Recently, static analysis tools customized for DNNs are emerging, which mainly focus on proposing tighter abstractions~\cite{gehr2018ai2,muller2022prima,paulsen2022linsyn} or incorporating abstractions into training to improve robustness~\cite{li2019robustra,mirman2018differentiable,zhang2021robustness}.
    Besides robustness, static analysis has also be applied to rigorously bound model difference~\cite{paulsen2020neurodiff}.
   
   
    Our approach includes a static analysis framework customized for numerical defect detection and fixing.
    
   
    \paragraph{Detecting and Exposing Numerical Defects in DNNs}
    Despite the widespread existence of numerical defect in real-world DNN-based systems~\cite{zhang2018empirical,humbatova2020taxonomy,kloberdanz2022deepstability}, few automatic approaches exist for detecting and exposing these defects.
    To the best of our knowledge, DEBAR~\cite{zhang2020detecting} and GRIST~\cite{yan2021exposing} are only two approaches.
    We discuss and compare with both approaches extensively in \Cref{sec:method,sec:exp}. 
   
   
   
   
   
   
   
   
   
   
   
    
    \vspace{-0.5em}
\section{Conclusion}
    \vspace{-0.3em}
    In this paper, we have presented a novel automatic approach  \sysname for reliability assurance of DNNs against numerical defects.
    \sysname supports detection of potential numerical defects, confirmation
of potential-defect feasibility, and suggestion of defect fixes. 
     \sysname includes multiple novel extensions and optimizations upon existing tools, and includes three novel techniques.
    Our extensive evaluation on real-world DNN architectures has demonstrated high effectiveness and efficiency of \sysname
    compared to both the state-of-the-art approaches and developers' fixes.
   

    \paragraph{Data Availability}
    All artifacts including the tool source code and experiment logs are available and actively maintained at \CodeIn{\url{https://github.com/llylly/RANUM}}.
    \ifnum\arxiv=0
    The supplemental materials containing omitted proofs, hyperparameters, and experiments are available in the arXiv version.
    \fi
    

    \vspace{-0.5em}
\section*{Acknowledgements}
    \vspace{-0.3em}

    This work is sponsored by the National Natural Science Foundation of China under Grant No. 62161146003, the National Key Research and Development Program of China under Grant No. 2019YFE0198100, the Innovation and Technology Commission of HKSAR under Grant No. MHP/055/19, and the Tencent Foundation/XPLORER PRIZE. 
    

    



\subsection{List of Supported Operators}

    \label{adxsec:supported-op}
    
    In this section, we provide a list of the 84 supported operators in the \sysname static analysis framework.
    These operators cover common operators used in DNNs. 
    To the best of our knowledge, \sysname provides abstractions for the largest number of operator types among existing DNN abstraction frameworks.
    We believe that the framework can be further extended to support other types of analysis, testing, and debugging for general DNN architectures or models such as robustness analysis and fairness testing.
    
    In particular, when compared to the state-of-the-art tool DEBAR~\cite{zhang2020detecting}, the underlined 17 operators are newly supported by \sysname.
    
   
   
    
    \noindent
    {\scriptsize \hspace{-0.3cm}
    \begin{tabular}{p{2.1cm}p{3.4cm}p{2.2cm}}
        \SmallCodeIn{Sub} & \SmallCodeIn{Add} & \SmallCodeIn{Mul} \\
        \SmallCodeIn{Div} & \SmallCodeIn{Pow} & \SmallCodeIn{MatMul} \\ 
        \underline{\SmallCodeIn{Gemm}} & \underline{\SmallCodeIn{MaxPool}} & \underline{\SmallCodeIn{GlobalMaxPool}} \\
        \SmallCodeIn{GlobalAveragePool} & \SmallCodeIn{Cos} & \SmallCodeIn{Sin} \\
        \SmallCodeIn{AveragePool} & \SmallCodeIn{Conv} & \SmallCodeIn{ConvTranspose} \\
        \underline{\SmallCodeIn{Pad}} & \underline{\SmallCodeIn{Reciprocal}} & \SmallCodeIn{Sqrt} \\
        \SmallCodeIn{Tanh} & \SmallCodeIn{Relu} & \SmallCodeIn{Softplus} \\
        \SmallCodeIn{LeakyRelu} & \underline{\SmallCodeIn{Softsign}} & \SmallCodeIn{Sigmoid} \\
        \SmallCodeIn{Neg} & \SmallCodeIn{Exp} & \SmallCodeIn{Log} \\
        \SmallCodeIn{Softmax} & \SmallCodeIn{LogSoftmax} & \SmallCodeIn{Abs} \\
        \SmallCodeIn{Ceil} & \SmallCodeIn{Floor} & \SmallCodeIn{Sign} \\
        \SmallCodeIn{Reshape} & \SmallCodeIn{Flatten} & \SmallCodeIn{Transpose} \\
        \SmallCodeIn{Shape} & \SmallCodeIn{Cast} & \SmallCodeIn{Slice} \\
        \underline{\SmallCodeIn{Gather}} & \underline{\SmallCodeIn{GatherND}} & \SmallCodeIn{Squeeze} \\
        \SmallCodeIn{Unsqueeze} & \underline{\SmallCodeIn{ScatterElements}} & \SmallCodeIn{Expand} \\
        \SmallCodeIn{Concat} & \SmallCodeIn{Split} & \SmallCodeIn{Identity} \\
        \SmallCodeIn{ConstantOfShape} & \underline{\SmallCodeIn{RandomNormalLike}} & \underline{\SmallCodeIn{RandomUniformLike}} \\
        \underline{\SmallCodeIn{RandomNormal}} & \underline{\SmallCodeIn{RandomUniform}} & \SmallCodeIn{Range} \\
        \SmallCodeIn{Constant} & \SmallCodeIn{Less} & \SmallCodeIn{LessOrEqual} \\
        \SmallCodeIn{Greater} & \SmallCodeIn{GreaterOrEqual} & \SmallCodeIn{Equal} \\
        \SmallCodeIn{Not} & \SmallCodeIn{Or} & \SmallCodeIn{Min} \\
        \SmallCodeIn{Max} & \SmallCodeIn{Clip} & \SmallCodeIn{Sum} \\
        \SmallCodeIn{ReduceMin} & \SmallCodeIn{ReduceMax} & \SmallCodeIn{ReduceSum} \\
        \SmallCodeIn{ReduceMean} & \SmallCodeIn{ArgMax} & \SmallCodeIn{ArgMin} \\
        \SmallCodeIn{Tile} & \underline{\SmallCodeIn{NegativeLogLikelihoodLoss}} & \underline{\SmallCodeIn{Loop}} \\
        \SmallCodeIn{SequenceInsert} & \SmallCodeIn{BatchNormalization} & \SmallCodeIn{OneHot} \\
        \underline{\SmallCodeIn{NonZero}} & \underline{\SmallCodeIn{Resize}} & \SmallCodeIn{ReduceProd} \\
        \SmallCodeIn{ReduceSumSquare} & \SmallCodeIn{IsInf} & \SmallCodeIn{IsNaN} \\
    \end{tabular}
    }

\subsection{List of Operators with Potential Numerical Defects}

    \label{adxsec:list-numerical-defects}
    
    In this section, we provide a full list of DNN operators that may contain numerical defects, along with their invalid ranges $\gI_{n_0,\mathsf{invalid}}$~(see definition of the numerical defect in \Cref{def:numerical-defect}) respectively.
    In the table, $U_{\min}$ and $U_{\max}$ stand for the minimum and maximum positive number of the input tensor's data type, respectively.
    
    {
    \centering \vspace{0.5em}
    \resizebox{\linewidth}{!}{
    \begin{tabular}{cc}
        \toprule
        Op. Type & $\gI_{n_0,\mathsf{invalid}}$ \\
        \hline
        \CodeIn{Pow} & $[-U_{\min}, U_{\min}] \times (-\infty, -U_{\min}]$ \\
        \CodeIn{Div} & $\sR \times [-U_{\min}, U_{\min}]$ \\
        \CodeIn{Reciprocal} & $[-U_{\min}, U_{\min}]$ \\
        \CodeIn{Sqrt} & $(-\infty, U_{\min}]$ \\
        \CodeIn{Exp} & $[\ln U_{\max}, \infty)$ \\
        \CodeIn{Log} & $(-\infty, U_{\min}]$ \\
        \CodeIn{Range} & $\sR \times \sR \times [-U_{\min}, U_{\min}]$ \\
        \CodeIn{NegativeLogLikelihoodLoss} & $[0,0]$ for number of non-zero cells using mean reduction \\
        \bottomrule
    \end{tabular}
    }
    }
    
\subsection{Detail Description of \sysname Static Analysis Framework}
    
    \label{adxsec:static-analysis}
    
    We present the omitted details in \Cref{subsec:static-analysis}.
    
    \subsubsection{Abstraction Domain and Characteristics}
    
        \label{adxsubsec:abstraction-domain}
    
        \begin{figure}[!t]
            \centering
            \includegraphics[width=0.9\linewidth]{figures/partitioning_example.pdf}
           
            \caption{Example of tensor partitioning.
            Tensor partitioning reduces the size of tensors in the abstract domain by sharing one interval bound among all elements inside one subblock. }
            \label{fig:partitioning}
           
        \end{figure}
    
    We first formally define our abstraction domain: interval with tensor partitioning.
    Following the notation in abstract interpretation literature~\cite{cousot1977static}, suppose the tensor has $m$ dimensions with shape $n_1 \times n_2 \times \cdots \times n_m$, we define the abstract domain as such:
    \begin{equation}
        \begin{aligned}
            & \sA := \{ \langle \vl, \vu, (S_i)_{i=1}^m \rangle: \vl, \vu \in \sR^{n'_1 \times \cdots \times n'_m} \\
            & \hspace{5em} |S_i| = n'_i, \forall s\in S_i, s\in \sN, 0 \le s < n_i \}.
        \end{aligned}
        \label{eq:abs-domain}
    \end{equation}
   
    Each element in $\sA$ is a triplet $a = \langle \vl, \vu, (S_i)_{i=1}^m\rangle$, where the first two elements are subblock-wise interval lower bound and upper bound respectively, and the last element $(S_1, S_2, \dots, S_m)$ contains $m$ sets, where each set $S_i$ corresponds to the \emph{0-indexed} split indices for the $i$-th dimension to form subblocks.
    We use $S_i[j] \in \sN$ to represent the $j$-th element of split index set $S_i$ sorted in ascending order and define $S_i[|S_i|] = n_i$.
   
    \Cref{fig:partitioning} illustrates this abstraction domain.
    
    To define the correspondence between the abstract domain and concrete domain $\sC := 2^{\sR^{n_1\times \cdots \times n_m}}$,
    we form Galois connection $\langle \sC, \subseteq\rangle \underset{\gamma}{\overset{\alpha}{\rightleftharpoons}} \langle \sA, \sqsubseteq \rangle$, where an abstraction function $\alpha: \sC \to \sA$ and the concretization function $\gamma: \sA \to \sC$ are defined as follows:
    \begin{subequations}
    \small
    \begin{equation}
        \begin{aligned}
            & \alpha(\gX) = \langle \vl, \vu, (S_i)_{i=1}^m \rangle \quad \text{where} \\
            & \hspace{2em} \vl_{i_1,i_2,\dots,i_m} = \min_{\vx\in \gX}  \min_{\substack{1\le k\le m\\ S_k[i_k] \le j_k < S_k[i_k + 1]}} \vx_{j_1,j_2,\dots,j_m}, \\
            & \hspace{2em} \vu_{i_1,i_2,\dots,i_m} = \max_{\vx\in \gX} \max_{\substack{1\le k\le m\\ S_k[i_k] \le j_k < S_k[i_k + 1]}} \vx_{j_1,j_2,\dots,j_m}, \\
        \end{aligned}
        \label{eq:alpha}
    \end{equation}
    \vspace{-1em}
    \begin{equation}
        \begin{aligned}
            \gamma(\langle \vl, \vu, (S_i)_{i=1}^m\rangle) = \{ \vx: &\vl_{i_1,i_2,\cdots,i_m} \le \vx_{j_1,j_2,\cdots,j_m} \le \vu_{i_1,i_2,\cdots,i_m}, \\
            & S_k[i_k] \le j_k < S_k[i_k+1], 1\le k\le m \}.
        \end{aligned}
        \label{eq:gamma}
    \end{equation}
    \end{subequations}
    In \Cref{eq:alpha}, the split indices $(S_i)_{i=1}^m$ can be arbitrarily chosen but need to be kept consistent with those in \Cref{eq:gamma}.
    Take \Cref{fig:partitioning} as the example, to abstract a set of tensors $\gX$ with shape $12\times 9 \times 4$, we define split indices for each dimension respectively, then impose interval constraints on each subblock of the tensor.
    For example, $[\vl_{0,0,0}, \vu_{0,0,0}]$ constrain any element in $\vx_{0:2,0:2,0:1}$, $[\vl_{3,2,1}, \vu_{3,2,1}]$ constrain any element in $\vx_{9:11,6:8,2:3}$. 
    
    
    
    \textbf{Abstraction Characteristics.}
    There are multiple ways to compute the abstractions, and we design our particular computational algorithms to achieve soundness, tightness, and differentiability.\\
    (1)~\textbf{Soundness}: We guarantee that all abstractions are sound.
    Formally, suppose $a_\vx$ and $a_\vw$ are the abstractions for input and weights respectively, we guarantee $f_n^*(\gamma(a_\vx); \gamma(a_\vw)) \subseteq \gamma(T_n^*(a_\vx; a_\vw))$ where $* \in \{\mathsf{in}, \mathsf{out}\}$, $n$ is any node in the architecture, and $T_n^*(a_\vx; a_\vw)$ is our computed abstract domain for input or output of node $n$.
    The soundness property theoretically guarantees that our detection approach has no false negatives, i.e., flag all potential numerical defects, and the validness of generated preconditions~(see \Cref{subsec:approach-overview}).\\
    (2)~\textbf{Tightness}: For most operators in DNN, given abstractions for its inputs, we compute the tightest possible interval abstraction for the output.
    Formally, for any atomic DNN operator $op$ and any abstract domain $i \in \sA$ of $op$'s input, if $op(\gamma(i)) \subseteq [\vl, \vu]$, i.e., $[\vl, \vu]$ is an interval abstraction of $op$'s output, then $\gamma(T_{op}(i)) \subseteq [\vl, \vu]$, i.e., our generated abstract domain $T_{op}(i)$ is always tighter or equal to any interval abstraction $[\vl, \vu]$.
    Such tightness can reduce false positives for detecting potential numerical defects and increase the span of generated failure-exhibiting intervals and fix intervals, which improves the quality of generated tests and preconditions.\\
    (3)~\textbf{Differentiability}: We compute differentiable abstractions.
    Concretely, if the input of node $n_2$ is deferentially dependent on output of node $n_1$, we can compute out gradients $\nabla_i o$ for any $i \in \{\vl_{n_1}, \vu_{n_1}\}$ and any $o \in \{\vl_{n_2}, \vu_{n_2}\}$, where $\langle \vl_{n_1}, \vu_{n_1}, (S^{n_1}_i)_{i=1}^{\dim(n_1.\mathsf{out})} \rangle$ and $\langle \vl_{n_2}, \vu_{n_2}, (S^{n_2}_i)_{i=1}^{\dim(n_2.\mathsf{in})} \rangle$ are abstractions of $\fout_{n_1}(\cdot; \cdot)$ and $\fin_{n_2}(\cdot; \cdot)$ respectively.
    When tightness and differentiability cannot be achieved at the same time~(e.g., for \CodeIn{floor} operator, tight abstraction $\langle \vl, \vu\rangle \mapsto \langle \lfloor \vl \rfloor, \lfloor \vu \rfloor \rangle$ is not generally differentiable, and differentiable abstraction $\langle \vl, \vu\rangle \mapsto \langle \vl - \mathbf{1}, \vu \rangle$ is not tight), we implement two abstract algorithms to achieve tightness and differentiability respectively.
   
    
    The existing approach DEBAR~\cite{zhang2020detecting} also proposes a static analysis framework for DNN architectures with tensor partitioning.
    In contrast to our framework, the abstract domain in DEBAR contains affine equalities besides interval domains.
    Therefore, DEBAR can be tighter than ours in some cases and can produce fewer positives, but due to the additional complexity of affine equalities, DEBAR tends to use the coarsest abstraction~(i.e., tensor partitioning) granularity and supports fewer operators than ours.
    At the same time, our algorithms produce the tightest interval abstractions while DEBAR has no tightness guarantee.
    As a result, \sysname detects more true numerical defects and has a comparable number of false positives~(see \Cref{subsec:static-detector-exp}), and we can leverage the feasibility confirmation support in \sysname to filter out false positives.
    
    
   
   
   
   
   
   
   
   
   
   
   
   
    
   
   
   
    
    
    \subsubsection{Initial Abstraction Construction with \textbf{Backward Fine-Grained Node Labeling}}
        \label{adxsubsec:backward-labeling}
        
        We construct abstract domains for initial nodes in two steps:
        First, we determine the tensor partitions, i.e., the split index sets $(S_i)_{i=1}^m$~(see \Cref{eq:abs-domain}), which determine the tightness and efficiency of our static analysis framework, because the tensor partitions of all other nodes will be solely dependent on the partitions of initial nodes as we will show in \Cref{adxsubsec:tight-abstract-dynamic-partitioning}.
        Second, we compute the interval bounds $\vl$ and $\vu$.
        
        We use the following principle to decide the tensor partitions:
        For nodes that are connected to operators requiring fine-grained abstractions~(e.g., the \CodeIn{shape} input for operator \CodeIn{Reshape}) with valid paths which we will specify later, we construct tensor partitions with the finest granularity, i.e., $S_i = \{0,1,\dots,n_i\}$.
        We call these nodes fine-grained initial nodes and starting from the next paragraph we introduce our novel technique of \emph{backward fine-grained node labeling} to find them out.
        For other nodes, we rely on downstream task requirements and user specifications to determine the partitions.
        For example, for fix generation~(see \Cref{subsec:precond-gen}), we use the coarsest granularity by default.
        
        \textbf{Backward Fine-Grained Node Labeling.}
        The fine-grained initial nodes are those starting a valid path in DNN computational graph $\gG$, where a path is valid if and only if the path does not traverse through fine-grained stopping operators and terminates at a fine-grained requiring operator with some specific input indices.
        As discussed in \Cref{subsec:static-analysis}, the fine-grained requiring operators fall into three categories: control-flow operators~(e.g., \CodeIn{Loop}), indexing operators~(e.g., \CodeIn{Slice}), and shaping operators~(e.g., \CodeIn{Reshape}). 
       
        Fine-grained stopping operators are those whose output is independent of input abstraction granularity~(so the granularity of preceding nodes does not matter).
        Detail lists of fine-grained stopping operators and fine-grained requiring operators are provided in \Cref{adxsec:fine-grain-stop-requiring-operators}.
        
        To find out fine-grained initial nodes, we propose backward fine-grained node labeling, which is similar to data dependency analysis in traditional programs:
        First, we invert all edges of the given computational graph $\gG = \langle \gV, \gE \rangle$ to get $\gG' = \langle \gV, \gE' \rangle$.
        Second, we attach a boolean label for each node $n\in \gV$ and initialize them with \CodeIn{False}.
        Third, we do topology sorting on $\gG'$.
        When encountering a node $n$ with \CodeIn{True}, if the node does not contain a fine-grained stopping operator, we propagate this label to all its subsequent nodes; otherwise, we propagate \CodeIn{False}.
        When encountering a node $n$ with \CodeIn{False}, if the node is a fine-grained requiring operator, we propagate \CodeIn{True} to some subsequent nodes corresponding to specific input indices and \CodeIn{False} to others; otherwise, propagate \CodeIn{False} to all subsequent nodes.
        Lastly, we collect all labels for initial nodes. 
        The nodes with \CodeIn{True} label are fine-grained initial nodes.
        
        To this point, we have determined the tensor partitions $(S_i)_{i=1}^m$ for all initial nodes, and we now compute $\vl$ and $\vu$ and thus finish the abstraction domain construction for the valid ranges of initial nodes.
        Initial nodes are further divided into input, weight, and constant nodes.
        In practice, most of weight nodes and all constant nodes have their initial values stored in the ONNX file, and we directly use \Cref{eq:alpha} with these initial values to compute $\vl$ and $\vu$.
        Otherwise, we rely on user specifications and built-in heuristics to determine $\vl$ and $\vu$ for initial nodes.
    
    \subsubsection{Internal Abstraction with \textbf{Dynamic Partitioning}}
        \label{adxsubsec:tight-abstract-dynamic-partitioning}
        
        For all supported operator types~(listed in \Cref{adxsec:supported-op}), we propose concrete algorithms to compute abstract domains for output with dynamic tensor partitioning.
        Formally, for each operator $op$, we construct the computable function $T_{op}: \sA \to \sA$ that satisfies soundness and tentatively satisfies tightness and differentiability, where $\sA$ is the abstract domain defined in \Cref{eq:abs-domain}.
        Therefore, following the computational procedure introduced in the end of \Cref{subsec:static-analysis}, we can compute end-to-end abstractions for all nodes in the given DNN architecture.
        As the showcase of our abstraction algorithm, we describe the algorithms for four representative operators: \CodeIn{MatMul}, \CodeIn{Conv}, \CodeIn{Softmax}, and \CodeIn{Loop}.
        
        \paragraph{\CodeIn{MatMul}}
        The \CodeIn{MatMul} operator computes the matrix multiplication of two operands.
        This operator is widely used in DNNs to express fully-connected layers.
        To simplify the narration, we focus on the two-dimensional case where we compute $op(\mA, \mB) = \mC = \mA \mB$ with $\mA \in \sR^{n \times m}$ and $\mB \in \sR^{m \times l}$. 
        Extensions to other dimensions by abstracting the broadcasting mechanism can be found in our open-source implementation.
        
        We denote the input abstractions of $op$ by $a = \langle \mL_a, \mU_a, (S_A^1, S_A^2) \rangle$ and $b = \langle \mL_b, \mU_b, (S_B^1, S_B^2) \rangle$ respectively.
        First, we compute the union $U = S_A^2 \cup S_B^1$.
        Second, we dynamically partition both $a$ and $b$ with split points $(S_A^1, U)$ and $(U, S_B^2)$ respectively and get
        $a' = \langle \mL'_a, \mU'_a, (S_A^1, U) \rangle$ and $b' = \langle \mL'_b, \mU'_b, (U, S_B^2) \rangle$.
        Note that $a$ and $a'$~(or $b$ and $b'$) correspond to the same concrete domain, but $a'$ and $b'$ have finer or equal partition granularity than $a$ and $b$.
        Third, we compute output abstraction $T_{op}(a,b) = c = \langle \mL_c, \mU_c, (S_A^1, S_B^2) \rangle$
        \begin{equation}
            \small
            \label{eq:mat-abs-1}
            \begin{aligned}
                \text{where} \quad (\mL_c)_{ij} & = \sum_{k=1}^{|U|} \vv_k \min_{\mA \in \{\mL'_a, \mU'_a\}, \mB \in \{\mL'_b, \mU'_b\}} \mA_{ik}\mB_{kj}, \\
                (\mU_c)_{ij} & = \sum_{k=1}^{|U|} \vv_k \max_{\mA \in \{\mL'_a, \mU'_a\}, \mB \in \{\mL'_b, \mU'_b\}} \mA_{ik}\mB_{kj}, \\
                \vv_k & = U[k]-U[k-1].
            \end{aligned}
        \end{equation}
        This formulation can guarantee tightness but is not efficient for tensor computation due to the inner minimum and maximum.
        Therefore, we also implement a fast-mode abstraction which trade-off tightness for efficiency: $T_{op}'(a,b) = c' = \langle \mL_c', \mU_c', (S_A^1, S_B^2)\rangle$ where
        \begin{equation}
            \small
            \begin{aligned}
                \mL_c^{'} = &
                \mU_a^{'-} \vv \mU_b^{'-} + 
                \mU_a^{'0} \vv \mL_b^{'-} + 
                \mU_a^{'+} \vv \mL_b^{'-} + \\
                & \mL_a^{'-} \vv \mU_b^{'0} +
                \mU_a^{'0} \vv \mL_b^{'0} +
                \mL_a^{'0} \vv \mU_b^{'0}
                 +
                \mU_a^{'+} \vv \mL_b^{'0} + \\
                & \mL_a^{'-} \vv \mU_b^{'+} + 
                \mL_a^{'0} \vv \mU_b^{'+} +
                \mL_a^{'+} \vv \mL_b^{'+}, \\
               
                \mU_c^{'} = & 
                \mL_a^{'-} \vv \mL_b^{'-} +
                \mL_a^{'0} \vv \mL_b^{'-} +
                \mL_a^{'+} \vv \mU_b^{'-} + \\
                & \mL_a^{'-} \vv \mL_b^{'0} +
                \mL_a^{'0} \vv \mL_b^{'0} +
                \mU_a^{'0} \vv \mU_b^{'0}
                 +
                \mU_a^{'+} \vv \mU_b^{'0} + \\
                & \mU_a^{'-} \vv \mL_b^{'+} +
                \mU_a^{'0} \vv \mU_b^{'+} +
                \mU_a^{'+} \vv \mU_b^{'+}.
            \end{aligned}
            \label{eq:mat-abs-2}
        \end{equation}
        In the above equation, for any $* \in \{a, b\}$, let $\circ$ be the elementwise~(Hadamard) product,
        \begin{equation}
            \small
            \begin{aligned}
                %
                \mL_*^{'-} = & \mL_*^{'} \circ \1[\mU_*^{'} < 0],
                \mU_*^{'-} = \mU_*^{'} \circ \1[\mU_*^{'} < 0], \\
                \mL_*^{'0} = & \mL_*^{'} \circ \1[\mL_*^{'} \le 0, \mU_*^{'} \ge 0],
                \mU_*^{'0} = \mU_*^{'} \circ \1[\mL_*^{'} \le 0, \mU_*^{'} \ge 0], \\
                \mL_*^{'+} = & \mL_*^{'} \circ \1[\mL_*^{'} > 0],
                \mU_*^{'+} = \mU_*^{'} \circ \1[\mL_*^{'} > 0].
               
               
               
               
               
               
               
               
            \end{aligned}
            \label{eq:mat-abs-2-defs}
        \end{equation}
        From \Cref{eq:mat-abs-2}, we can observe that the abstraction can be easily implemented with tensor computations.
        In \Cref{adxsec:proofs}, we prove the soundness and tightness of these abstractions.
        
        \paragraph{\CodeIn{Conv}}
        The \CodeIn{Conv} operator computes the discrete convolution of two operands.
        This operator is widely used in convolutional neural networks~(CNNs, \cite{krizhevsky2012imagenet}).
        To simplify the narration, we focus on the single-channel single-stride two-dimensional case where we compute $op(\mA, \mW) = \mC$ with $\mA$ being the input matrix and $\mW$ being the convolution kernel. 
        Extensions to general cases are provided in our open-source implementation.
        
        We first dynamically split the kernel abstraction to the finest granularity.
        Second, we compute the receptive field of each output position, which is the sub-region of $\mA$ that decides each output position.
        Third, we inspect the alignment between receptive fields and $\mA$'s partitions.
        If neighboring output positions have their receptive fields partitioned in the same sub-block of $\mA$, it means that these positions can be abstracted by a single interval, i.e., these positions can be partitioned together.
        Fourth, using this principle, we derive the paritions of the output tensor, and repeat the input tensors accordingly so that the abstract computation can be written as convolutional operations.
        Last, we modify the abstraction computation equations in \Cref{eq:mat-abs-2} by replacing matrix multiplication with convolution to compute the abstraction of output $\mC$.
        
        \paragraph{\CodeIn{Softmax}}
        The \CodeIn{Softmax} operator computes normalized exponential values for the given input.
        \CodeIn{Softmax} is widely deployed for classification tasks to yield the normalized confidence.
        In one-dimensional case, for input $\vx \in \sR^n$, a \CodeIn{Softmax} operator outputs $op(\vx) = \frac{\exp(\vx)}{\sum_{i=1}^n \exp(\vx)_i}$.
        The output abstraction of \CodeIn{Softmax} operator $op$ can be thus computed: $T_{op}(\langle \vl, \vu, (S)\rangle) = \langle \vl^o, \vu^o, (S)\rangle$ where
        \begin{equation}
            \small
            \resizebox{0.9\linewidth}{!}{$ \displaystyle
            \begin{aligned}
                \vl^o_i & = \dfrac{\exp(\vl_i)}{\exp(\vu)^\T \vv - \exp(\vu_i) + \exp(\vl_i)}, 
                \vu^o_i = \dfrac{\exp(\vu_i)}{\exp(\vl)^\T \vv - \exp(\vl_i) + \exp(\vu_i)}, \\
                \vv_k & = S[k] - S[k-1].
            \end{aligned}
            $}
            \label{eq:softmax-abs}
        \end{equation}
        The output abstraction's partition is dynamically decided by the input abstraction's partition.
        We prove the soundness and tightness of the abstraction in \Cref{adxsec:proofs}.

        \paragraph{\CodeIn{Loop}}
        The node with a \CodeIn{Loop} operator contains a sub-computational graph with dynamic controlling counters and conditions to represent a runtime loop.
        The \CodeIn{Loop} operator can be used to represent recurrent neural networks~(RNNs) and handle the input sequence of variable length.
        We require the controlling counter and loop termination condition to have the finest abstraction granularity.
        Then, we recursively apply our static analysis framework to the sub-graph representing the loop body and update the interval of the loop counter iteratively.
        When the abstraction interval of the loop counter can explicitly decide whether to terminate the loop, we continue or terminate the loop iterations respectively;
        otherwise, we merge the current abstraction interval with the interval obtained after another iteration.
        We repeat this process until termination.
        Theoretically, this execution process cannot guarantee the termination, i.e., the loop body may execute for infinite times.
        However, in practice, on our dataset, our analysis framework already suffices to guarantee the termination in all cases.
        In the future, we can apply a widening operation if termination cannot be tightly abstracted. 
        
        \begin{remark}
            As we can see, the novel dynamic partition technique is incorporated into the computation process of each operator's abstraction.
            The soundness and tightness of our designed abstractions are immediately achieved by design, and the differentiability of our designed abstractions can be implemented via the auto-differentiation functionality of  popular DL libraries like \CodeIn{PyTorch}~\cite{paszke2019pytorch} and \CodeIn{Tensorflow}~\cite{tensorflow2015-whitepaper} where we use \CodeIn{PyTorch} for implementation. 
            The abstractions for DNNs are also implemented for other applications, e.g., for robustness verification~\cite{singh2019abstract,mirman2018differentiable}.
            However, our abstractions are tailored for tensor partitioned interval domain which is particularly suitable for testing and debugging the numerical defects.
            To the best of our knowledge, these abstractions are the first that achieve soundness, tightness, and differentiability.
        \end{remark}
    
    \subsubsection{List of Fine-Grained Requiring and Stopping Operators}
    
        \label{adxsec:fine-grain-stop-requiring-operators}
        
        We use fine-grained requiring and fine-grained stopping operators in \Cref{adxsubsec:backward-labeling}.
        Among all supported operators, the fine-grained requiring operators are 
        \vspace{0.5em}
        
        \begin{tabular}{p{8.1cm}}
            \CodeIn{Reshape} (input \#2), 
            \CodeIn{Slice} (input \#2, \#3, \#4, \#5), 
            \CodeIn{Squeeze} (input \#2), 
            \CodeIn{Unsqueeze} (input \#2), 
            \CodeIn{Tile} (input \#2, \#3), 
            \CodeIn{Loop} (input \#1, \#2),
            \CodeIn{SequenceInsert} (input \#3),
            \CodeIn{ConstantOfShape} (input \#1),
            \CodeIn{Gather} (input \#2),
            \CodeIn{GatherND} (input \#2),
            \CodeIn{ReduceSum} (input \#2),
            \CodeIn{ScatterElements} (input \#2),
            \CodeIn{Expand} (input \#2),
            \CodeIn{Split} (input \#2),
            \CodeIn{Pad} (input \#2, \#3),
            \CodeIn{NegativeLogLikelihoodLoss} (input \#2),
            \CodeIn{Clip} (input \#2, \#3),
            \CodeIn{OneHot} (input \#2),
            \CodeIn{Resize} (input \#2, \#3, \#4).
        \end{tabular}
        \vspace{0.5em}

        \noindent
        The fine-grained stopping operators are 
        \vspace{0.5em}
        
        \begin{tabular}{p{8.1cm}}
            \CodeIn{Shape}, \CodeIn{RandomNormalLike}, \CodeIn{RandomUniformLike}.
        \end{tabular}
        
        \vspace{0.5em}
            
    \subsubsection{Proofs}
        \label{adxsec:proofs}
        \allowdisplaybreaks
        
        Here we present the omitted proofs in \Cref{adxsubsec:tight-abstract-dynamic-partitioning}.
        
        \paragraph{\CodeIn{MatMul}}
        
        \begin{theorem}[Tightness of Abstraction by \Cref{eq:mat-abs-1} for \CodeIn{MatMul}]
            Suppose $op$ is the \CodeIn{MatMul} operator and $T_{op}$ is as defined in \Cref{eq:mat-abs-1}, then if $op(\gamma(a), \gamma(b)) \subseteq [\mL, \mU]$ for $a, b \in \sA$, $\gamma(T_{op}(a, b)) \subseteq [\mL, \mU]$.
            \label{thm:1}
        \end{theorem}
        
        \begin{proof}
            We use $\mL_a$ and $\mU_a$ to denote the element-wise interval lower and upper bounds of $\gamma(a)$; and use $\mL_b$ and $\mU_b$ to denote these bounds of $\gamma(b)$ respectively, where a formal definition is in \Cref{eq:gamma}.
            Then, there exist $\mA$ and $\mB$ with $\mL_a \le \mA \le \mU_a$ and $\mL_b \le \mB \le \mU_b$, such that
            \begin{small}
                \begin{align}
                    (\mA\mB)_{ij} & = \sum_{k=1}^l \min \{ \mL_{a,ik} \mL_{b,kj}, \mL_{a,ik} \mU_{b,kj}, \nonumber \\
                    & \hspace{5em} \mU_{a,ik} \mL_{b,kj}, \mU_{a,ik} \mU_{b,kj} \}, \label{eq:pf-1} \\
                    \text{or } (\mA\mB)_{ij} & = \sum_{k=1}^l \max \{ \mL_{a,ik} \mL_{b,kj}, \mL_{a,ik} \mU_{b,kj}, \nonumber \\
                    & \hspace{5em} \mU_{a,ik} \mL_{b,kj}, \mU_{a,ik} \mU_{b,kj} \}. \label{eq:pf-2} 
                \end{align}
            \end{small}
            By definition, we have
            \begin{equation}
                \small
                \mL_{ij} \le \text{\Cref{eq:pf-1}}, \mU_{ij} \ge \text{\Cref{eq:pf-2}}.
            \end{equation}
            We let $\mL'_c$ and $\mU'_c$ to denote the element-wise interval lower and upper bounds for $\gamma(T_{op}(a,b))$, 
            let $S_A^1[i'] \le i-1 < S_A^1[i'+1]$ and $S_B^2[j'] \le j-1 < S_B^2[j'+1]$,
            then from \Cref{eq:mat-abs-1},
            \begin{equation}
                \small
                \begin{aligned}
                & (\mL'_c)_{ij} = (\mL_c)_{i'j'} \\
                = & \sum_{k=1}^{|U|} \sum_{k'=U[k-1]+1}^{U[k]} \min \{ \mL_{a,ik'} \mL_{b,k'j}, \mL_{a,ik'} \mU_{b,k'j}, \\
                & \hspace{10em} \mU_{a,ik'} \mL_{b,k'j}, \mU_{a,ik'} \mU_{b,k'j} \} \\
                = & \sum_{k=1}^l \min \{ \mL_{a,ik} \mL_{b,kj}, \mL_{a,ik} \mU_{b,kj}, \mU_{a,ik} \mL_{b,kj}, \mU_{a,ik} \mU_{b,kj} \} \\
                = & \text{\Cref{eq:pf-1}} \ge \mL_{ij}.
                \end{aligned}
                \label{eq:pf-5}
            \end{equation}
            Similarly, $(\mU'_c)_{ij} \le \mU_{ij}$.
            Thus, $\gamma(T_{op}(a,b)) \subseteq [\mL, \mU]$.
        \end{proof}
        
        \begin{theorem}[Soundness of Abstraction by \Cref{eq:mat-abs-1} for \CodeIn{MatMul}]
            Suppose $op$ is the \CodeIn{MatMul} opeartor and $T_{op}$ is as defined in \Cref{eq:mat-abs-1}, then $op(\gamma(a), \gamma(b)) \subseteq \gamma(T_{op}(a,b))$.
            \label{thm:2}
        \end{theorem}
        
        \begin{proof}
            We use $\mL_a$ and $\mU_a$ to denote the element-wise interval lower and upper bounds of $\gamma(a)$; and use $\mL_b$ and $\mU_b$ to denote these bounds of $\gamma(b)$ respectively, where a formal definition is in \Cref{eq:gamma}.
            For any $\mA$ and $\mB$ such that $\mL_a \le \mA \le \mU_a$ and $\mL_b \le \mB \le \mU_b$, 
            \begin{equation}
                \small
                \begin{aligned}
                    (\mA\mB)_{ij} & = \sum_{k=1}^l \mA_{ik} \mB_{kj} \\
                    & \ge \min\{ \mL_{a,ik} \mL_{b,kj}, \mL_{a,ik} \mU_{b,kj}, \mU_{a,ik} \mL_{b,kj}, \mU_{a,ik} \mU_{b,kj} \} \\
                    & =: (\mL_{ab}')_{ij}
                \end{aligned}
                \label{eq:pf-3} 
            \end{equation}
            and
            \begin{equation}
                \small
                \begin{aligned}
                    (\mA\mB)_{ij} & \le \max\{ \mL_{a,ik} \mL_{b,kj}, \mL_{a,ik} \mU_{b,kj}, \mU_{a,ik} \mL_{b,kj}, \mU_{a,ik} \mU_{b,kj} \} \\
                    & =: (\mU_{ab}')_{ij}.
                \end{aligned}
                \label{eq:pf-4} 
            \end{equation}
            Thus, $op(\gamma(a), \gamma(b)) \subseteq [\mL'_{ab}, \mU'_{ab}]$.
            On the other hand, $\gamma(T_{op}(a,b)) = [\mL'_{ab}, \mU'_{ab}]$ as seen from \Cref{eq:pf-5}.
            Therefore, $op(\gamma(a), \gamma(b)) \subseteq \gamma(T_{op}(a,b))$.
        \end{proof}
    
        \begin{theorem}[Soundness of Abstraction by \Cref{eq:mat-abs-2} for \CodeIn{MatMul}]
            Suppose $op$ is the \CodeIn{MatMul} operator and $T'_{op}$ is as defined in \Cref{eq:mat-abs-2,eq:mat-abs-2-defs}, then $op(\gamma(a), \gamma(b)) \subseteq \gamma(T'_{op}(a,b))$.
            \label{thm:3}
        \end{theorem}
        
        \begin{proof}
            We use $\mL_a$ and $\mU_a$ to denote the element-wise interval lower and upper bounds of $\gamma(a)$; and use $\mL_b$ and $\mU_b$ to denote these bounds of $\gamma(b)$ respectively, where a formal definition is in \Cref{eq:gamma}.
            We define $\mL'_{ab}$ and $\mU'_{ab}$ by \Cref{eq:pf-3,eq:pf-4}.
            From the proof of \Cref{thm:2}, we have $op(\gamma(a), \gamma(b)) \subseteq [\mL'_{ab}, \mU'_{ab}]$.
            We now only need to show that $[\mL'_{ab}, \mU'_{ab}] \subseteq \gamma(T'_{op}(a,b))$.
            
            $\gamma(T'_{op}(a,b))$ imposes independent interval abstractions element-wise, therefore, we study each element independently.
            For the element $(i,j)$, from \Cref{eq:mat-abs-2,eq:mat-abs-2-defs}, the interval lower bound of $\gamma(T'_{op}(a,b))$, namely $(\mL_c)_{ij}$, satisfies
            \begin{equation}
                \small
                \begin{aligned}
                    & (\mL_c)_{ij} =  \sum_{k=1}^l  
                    (\mU_a^-)_{ik} (\mU_b^-)_{kj} +
                    (\mU_a^{0})_{ik} (\mL_b^{-})_{kj} + 
                    (\mU_a^{+})_{ik} (\mL_b^{-})_{kj} + \\
                    & (\mL_a^{-})_{ik} (\mU_b^{0})_{kj} +
                    (\mU_a^{0})_{ik} (\mL_b^{0})_{kj} +
                    (\mL_a^{0})_{ik} (\mU_b^{0})_{kj} +
                    (\mU_a^{+})_{ik} (\mL_b^{0})_{kj} + \\
                    & (\mL_a^{-})_{ik} (\mU_b^{+})_{kj} + 
                    (\mL_a^{0})_{ik} (\mU_b^{+})_{kj} +
                    (\mL_a^{+})_{ik} (\mL_b^{+})_{kj}.
                \end{aligned}
                \label{eq:pf-6}
            \end{equation}
            By \Cref{eq:mat-abs-2-defs},
            \begin{itemize}
                \item when $\mL_{a,ik} \le \mU_{a,ik} < 0$, \\
                $(\mL_a^-)_{ik} = \mL_{a,ik}, (\mU_a^-)_{ik} = \mU_{a,ik}, 
                (\mL_a^0)_{ik} = 0, \\ (\mU_a^0)_{ik} = 0, 
                (\mL_a^+)_{ik} = 0, (\mU_a^+)_{ik} = 0$;
                
                \item when $\mL_{a,ik} \le 0 \le \mU_{a,ik}$, \\
                $(\mL_a^-)_{ik} = 0, (\mU_a^-)_{ik} = 0, 
                (\mL_a^0)_{ik} = \mL_{a,ik}, \\ (\mU_a^0)_{ik} = \mU_{a,ik}, 
                (\mL_a^+)_{ik} = 0, (\mU_a^+)_{ik} = 0$;
                
                \item when $0 < \mL_{a,ik} \le \mU_{a,ik}$, \\
                $(\mL_a^-)_{ik} = 0, (\mU_a^-)_{ik} = 0, 
                (\mL_a^0)_{ik} = 0, \\ (\mU_a^0)_{ik} = 0, 
                (\mL_a^+)_{ik} = \mL_{a,ik}, (\mU_a^+)_{ik} = \mU_{a,ik}$.
            \end{itemize}
            Similarly for $(\mL_b^-)_{kj}$, $(\mU_b^-)_{kj}$, $(\mL_b^0)_{kj}$, $(\mU_b^0)_{kj}$, $(\mL_b^+)_{kj}$ and $(\mU_b^+)_{kj}$.
            Thus, by enumerating all cases, we have
            \begin{equation}
                \small
                \begin{aligned}
                \text{\Cref{eq:pf-6}} \le \min\{ & \mL_{a,ik} \mL_{b,kj},\mL_{a,ik} \mU_{b,kj},  \\
                & \mU_{a,ik} \mL_{b,kj}, \mU_{a,ik} \mU_{b,kj} \} = (\mL'_{ab})_{ij}.
                \end{aligned}
            \end{equation}
            Similarly, the interval lower bound of $\gamma(T'_{op}(a,b))$, namely $(\mU_c)_{ij}$,
            \begin{equation}
                \small
                \begin{aligned}
                (\mU_c)_{ij} \ge \max\{ & \mL_{a,ik} \mL_{b,kj}, \mL_{a,ik} \mU_{b,kj}, \\
                & \mU_{a,ik} \mL_{b,kj}, \mU_{a,ik} \mU_{b,kj} \} = (\mU'_{ab})_{ij}.
                \end{aligned}
            \end{equation}
            Thus, $[\mL'_{ab}, \mU'_{ab}] \subseteq \gamma(T'_{op}(a,b))$.
        \end{proof}
        
        \paragraph{\CodeIn{Softmax}}
        
        \begin{theorem}[Tightness of Abstraction for \CodeIn{Softmax}]
            Suppose $op: \sR^n \to \sR^n$ is the \CodeIn{Softmax} operator and $T_{op}$ is as defined in \Cref{eq:softmax-abs}, if $op(\gamma(\langle \vl, \vu, (S)\rangle)) \subseteq [\vl^r, \vu^r]$, then $\gamma(T_{op}(\langle \vl, \vu, (S)\rangle)) \subseteq [\vl^r, \vu^r]$.
            \label{thm:4}
        \end{theorem}
        
        \begin{proof}
            We use $\vl'$ and $\vu'$ to denote the element-wise interval lower and upper bounds of $\gamma(\langle \vl, \vu, (S)\rangle)$.
            Formally, $\vl'_i = \vl_{i'}$ and $\vu'_i = \vu_{i'}$ where $S[i'] \le i - 1 < S[i'+1]$.
            Then, 
           
            for each $i$, by setting $\vx_i = \vl_i$ and $\vx_j = \vu_j$ for all $j \neq i$, we get
            \begin{equation}
                \small
                \begin{aligned}
                    \vl_i^r & \le \dfrac{\exp(\vl'_i)}{\exp(\vl'_i) + \sum_{k=1, k\neq i}^n \exp(\vu'_k)} \\
                    & = \dfrac{\exp(\vl_{i'})}{\exp(\vl_{i'}) + \sum_{k=1}^n \exp(\vu'_k) - \exp(\vu'_i)} \\
                    & = \dfrac{\exp(\vl_{i'})}{\exp(\vl_{i'}) + \vv^\T \exp(\vu_k) - \exp(\vu_{i'})} = \vl_{i'}^o.
                \end{aligned}
            \end{equation}
            In addition, by setting $\vx_i = \vu_i$ and $\vx_j = \vl_j$ for all $j \neq i$, we get $\vu_i^r \ge \vu_{i'}^o$.
            Here, $\vl^o$ and $\vu^o$ are as defined in \Cref{eq:softmax-abs}.
            Combining these two arguments, we get 
            $$\gamma(T_{op}(\langle \vl, \vu, (S)\rangle)) \subseteq [\vl^r, \vu^r].$$
        \end{proof}
    
        \begin{theorem}[Soundness of Abstraction for \CodeIn{Softmax}]
            Suppose $op: \sR^n \to \sR^n$ is the \CodeIn{Softmax} operator and $T_{op}$ is as defined in \Cref{eq:softmax-abs}, then $op(\gamma(\langle \vl, \vu, (S)\rangle)) \subseteq \gamma(T_{op}(\langle \vl, \vu, (S)\rangle))$.
            \label{thm:5}
        \end{theorem}
        
        \begin{proof}
            Leveraging the fact that softmax function is monotinically increasing, i.e., $\frac{\d op(\vx)_i}{\d \vx_j} > 0$, we have
            $\gamma(\langle \vl, \vu, (S)\rangle)_i \in [\vl^o_{i'}, \vu^o_{i'}]$,
            where $S[i'] \le i-1 < S[i'+1]$.
            Since $[\vl^o_{i'}, \vu^o_{i'}] = \gamma(T_{op}(\langle \vl, \vu, (S)\rangle))_i$ by our definition of $T_{op}$, $op(\gamma(\langle \vl, \vu, (S)\rangle)) \subseteq \gamma(T_{op}(\langle \vl, \vu, (S)\rangle))$ follows.
        \end{proof}

\balance
\subsection{Hyperparameters}

    \label{adxsec:detail-hyperparam}

    In this section, we listed hyperparameters in our two-step system test generation technique:
   
   
   
    For the unit test generation, we use the Adam optimizer where the learning rate is $1$ and the maximum iteration number is $100$.
    For the training example generation, we target for training example under learning rate $\gamma=1$ and the approach has similar performance under other learning rates.
    We follow the convention in the DLG attack, where we use the L-BFGS method as the optimizer for gradient-based minimization.
    We terminate the method and return ``\CodeIn{failed}'' if either the running time exceeds \SI{1800}{s}~(universal execution time limit for all experimental evaluations), or a failure-exhibiting example training is not found after $300$ iterations of L-BFGS optimization.
    
   
   
   
   
    
   
    
\subsection{Ablation Study of System Test Generation}

    \label{adxsec:system-testgen-by-stage}
    
   
    \sysname uses the two-step test generation technique to
    produce failure-exhibiting system tests:
    it first generates failure-exhibiting unit tests with gradient back-propagation, then generates failure-exhibiting training example via the extended DLG attack.
    To isolate the impact of \sysname at each step, we replace either step with random sampling:
    ``Random + \sysname'' which first generates failure-exhibiting unit tests via random sampling, then generates training example via \sysname;
    ``\sysname + Random'' which first generates failure-exhibiting unit tests via \sysname, then generates training example via random sampling.
    We follow the same evaluation protocol as in \Cref{subsec:system-test-gen-exp} in this ablation study.
    We find that ``\sysname + Random'' takes 9.38X running time than \sysname and fails for 68 runs~(\sysname only fails for 57 runs); and ``Random + \sysname'' fails for 113 runs~(roughly 2X failed runs compared to \sysname).
    This study implies that \sysname's technique helps to improve effectiveness and efficiency at both steps of failure-exhibiting system test generation compared to the pure random baseline.
    The improvement for the first step is mainly from the effectiveness perspective, and the improvement for the second step is mainly from the efficiency perspective.
   
   
    






\section{List of Supported Operators}

    \label{adxsec:supported-op}

    
    \hspace{10em}
    
    In this appendix, we provide a list of the 84 supported operators in the \sysname static analysis framework.
    These operators cover common operators used in DNNs. 
    To the best of our knowledge, \sysname provides abstractions for the largest number of operator types among existing DNN abstraction frameworks.
    We believe the framework can be further extended to support other types of analysis, testing, and debugging for general DNN architectures or models such as robustness analysis and fairness testing.
    
   
   
    
    \noindent
    {\scriptsize
    \begin{tabular}{p{1.3cm}p{1.3cm}p{1.3cm}p{1.3cm}p{1.3cm}}
        \CodeIn{Sub} & \CodeIn{Add} & \CodeIn{Mul} & \CodeIn{Div} & \CodeIn{Pow} \\
        \CodeIn{MatMul} & \CodeIn{Gemm} & \CodeIn{MaxPool} & \CodeIn{GlobalMaxPool} & \CodeIn{GlobalAveragePool} \\
        \CodeIn{Cos} & \CodeIn{Sin} & \CodeIn{AveragePool} & \CodeIn{Conv} & \CodeIn{ConvTranspose} \\
        \CodeIn{Pad} & \CodeIn{Reciprocal} & \CodeIn{Sqrt} & \CodeIn{Tanh} & \CodeIn{Relu} \\
        \CodeIn{Softplus} & \CodeIn{LeakyRelu} & \CodeIn{Softsign} & \CodeIn{Sigmoid} & \CodeIn{Neg} \\
        \CodeIn{Exp} & \CodeIn{Log} & \CodeIn{Softmax} & \CodeIn{LogSoftmax} & \CodeIn{Abs} \\
        \CodeIn{Ceil} & \CodeIn{Floor} & \CodeIn{Sign} & \CodeIn{Reshape} & \CodeIn{Flatten} \\
        \CodeIn{Transpose} & \CodeIn{Shape} & \CodeIn{Cast} & \CodeIn{Slice} & \CodeIn{Gather} \\
        \CodeIn{GatherND} & \CodeIn{Squeeze} & \CodeIn{Unsqueeze} & \CodeIn{ScatterElements} & \CodeIn{Expand} \\
        \CodeIn{Concat} & \CodeIn{Split} & \CodeIn{Identity} & \CodeIn{ConstantOfShape} & \CodeIn{RandomNormalLike} \\
        \CodeIn{RandomUniform} \CodeIn{Like} & \CodeIn{RandomNormal} & \CodeIn{RandomUniform} & \CodeIn{Range} & \CodeIn{Constant} \\
        \CodeIn{Less} & \CodeIn{LessOrEqual} & \CodeIn{Greater} & \CodeIn{GreaterOrEqual} & \CodeIn{Equal} \\
        \CodeIn{Not} & \CodeIn{Or} & \CodeIn{Min} & \CodeIn{Max} & \CodeIn{Clip} \\
        \CodeIn{Sum} & \CodeIn{ReduceMin} & \CodeIn{ReduceMax} & \CodeIn{ReduceSum} & \CodeIn{ReduceMean} \\
        \CodeIn{ArgMax} & \CodeIn{ArgMin} & \CodeIn{Tile} & \CodeIn{NegativeLog} \CodeIn{LikelihoodLoss} & \CodeIn{Loop} \\
        \CodeIn{SequenceInsert} & \CodeIn{Batch} \CodeIn{Normalization} & \CodeIn{OneHot} & \CodeIn{NonZero} & \CodeIn{Resize} \\
        \CodeIn{ReduceProd} & \CodeIn{ReduceSumSquare} & \CodeIn{IsInf} & \CodeIn{IsNaN} \\
    \end{tabular}
    }
    
    \hspace{2em}
    
    
\section{Detail Description of Static Analysis Framework}
    
    \label{adxsec:static-analysis}
    
    In this appendix, we present the omitted details in \Cref{subsec:static-analysis}.
    
    We first formally define our abstraction domain: interval with tensor partitioning.
    Following the notation in abstract interpretation literature~\cite{cousot1977static}, suppose the tensor has $m$ dimensions with shape $n_1 \times n_2 \times \cdots \times n_m$, we define the abstract domain as such:
    \begin{equation}
        \begin{aligned}
            & \sA := \{ \langle \vl, \vu, (S_i)_{i=1}^m \rangle: \vl, \vu \in \sR^{n'_1 \times \cdots \times n'_m} \\
            & \hspace{5em} |S_i| = n'_i, \forall s\in S_i, s\in \sN_+, 0 \le s < n_i \}.
        \end{aligned}
        \label{eq:abs-domain}
    \end{equation}
   
    Each element in $\sA$ is a triplet $a = \langle \vl, \vu, (S_i)_{i=1}^m\rangle$, where the first two elements are subblock-wise interval lower bound and upper bound respectively, and the last element $(S_1, S_2, \dots, S_m)$ contains $m$ sets, where each set $S_i$ corresponds to the \emph{0-indexed} split indices for the $i$-th dimension to form subblocks.
    We use $S_i[j] \in \sN_+$ to represent the $j$-th element of split index set $S_i$ sorted in ascending order and define $S_i[|S_i|] = n_i$.
   
    
    To define the correspondence between the abstract domain and concrete domain $\sC := 2^{\sR^{n_1\times \cdots \times n_m}}$,
    we form Galois connection $\langle \sC, \subseteq~\rangle \underset{\gamma}{\overset{\alpha}{\rightleftharpoons}} \langle \sA, \sqsubseteq \rangle$, where an abstraction function $\alpha: \sC \to \sA$ and the concretization function $\gamma: \sA \to \sC$ are defined as follows:
    \begin{subequations}
    \small
    \begin{equation}
        \begin{aligned}
            & \alpha(\gX) = \langle \vl, \vu, (S_i)_{i=1}^m \rangle \quad \text{where} \\
            & \hspace{2em} \vl_{i_1,i_2,\dots,i_m} = \min_{\vx\in \gX}  \min_{\substack{1\le k\le m\\ S_k[i_k] \le j_k < S_k[i_k + 1]}} \vx_{j_1,j_2,\dots,j_m}, \\
            & \hspace{2em} \vu_{i_1,i_2,\dots,i_m} = \max_{\vx\in \gX} \max_{\substack{1\le k\le m\\ S_k[i_k] \le j_k < S_k[i_k + 1]}} \vx_{j_1,j_2,\dots,j_m}, \\
        \end{aligned}
        \label{eq:alpha}
    \end{equation}
    \vspace{-1em}
    \begin{equation}
        \begin{aligned}
            \gamma(\langle \vl, \vu, (S_i)_{i=1}^m\rangle) = \{ \vx: &\vl_{i_1,i_2,\cdots,i_m} \le \vx_{j_1,j_2,\cdots,j_m} \le \vu_{i_1,i_2,\cdots,i_m}, \\
            & S_k[i_k] \le j_k < S_k[i_k+1], 1\le k\le m \}.
        \end{aligned}
        \label{eq:gamma}
    \end{equation}
    \end{subequations}
    In \Cref{eq:alpha}, the split indices $(S_i)_{i=1}^m$ can be arbitrarily chosen but need to be kept consistent with those in \Cref{eq:gamma}.
    Take \Cref{fig:partitioning} as the example, to abstract a set of tensors $\gX$ with shape $12\times 9 \times 4$, we define split indices for each dimension respectively, then impose interval constraints on each subblock of the tensor.
    For example, $[\vl_{0,0,0}, \vu_{0,0,0}]$ constrain any element in $\vx_{0:2,0:2,0:1}$, $[\vl_{3,2,1}, \vu_{3,2,1}]$ constrain any element in $\vx_{9:11,6:8,2:3}$. 
    
    
    
    \textbf{Abstraction Characteristics.}
    There are multiple ways to compute the abstractions, and we design our particular computational algorithms to achieve soundness, tightness, and differentiability.\\
    (1)~\textbf{Soundness}: We guarantee that all abstractions are sound.
    Formally, suppose $a_\vx$ and $a_\vw$ are the abstraction for input and weights respectively, we guarantee $f_n^*(\gamma(a_\vx); \gamma(a_\vw)) \subseteq \gamma(T_n^*(a_\vx; a_\vw))$ where $* \in \{\mathsf{in}, \mathsf{out}\}$, $n$ is any node in the architecture, and $T_n^*(a_\vx; a_\vw)$ is our computed abstract domain for input or output of node $n$.
    The soundness property theoretically guarantees that our defect detection approach has no false negatives, i.e., flag out all possible numerical defects, and the validness~(\Cref{def:fix-precond}) of generated precondition fixes.\\
    (2)~\textbf{Tightness}: For most operators in DNN, given abstractions for its inputs, we compute the tightest possible interval abstraction for the output.
    Formally, for any atomic DNN operator $op$ and any abstract domain $i \in \sA$ of $op$'s input, if $op(\gamma(i)) \subseteq [\vl, \vu]$, i.e., $[\vl, \vu]$ is an interval abstraction of $op$'s output, then $\gamma(T_{op}(i)) \subseteq [\vl, \vu]$, i.e., our generated abstract domain $T_{op}(i)$ is always tighter or equal to any interval abstraction $[\vl, \vu]$.
    Such tightness can reduce false positives for detecting numerical defects and increase the span of generated failure-exhibiting intervals and precondition-fix intervals, which improves the quality of generated tests and preconditions.\\
    (3)~\textbf{Differentiability}: We compute differentiable abstractions.
    Concretely, if the input of node $n_2$ is deferentially dependent on output of node $n_1$, we can compute out gradients $\nabla_{\vl_{n_1}/\vu_{n_1}} \vl_{n_2}/\vu_{n_2}$, where $\langle \vl_{n_1}, \vu_{n_1}, (S^{n_1}_i)_{i=1}^{\dim(n_1.\mathsf{out})} \rangle$ and $\langle \vl_{n_2}, \vu_{n_2}, (S^{n_2}_i)_{i=1}^{\dim(n_2.\mathsf{in})} \rangle$ are abstractions of $\fout_{n_1}(\cdot; \cdot)$ and $\fin_{n_2}(\cdot; \cdot)$ respectively.
    When tightness and differentiability cannot be achieved at the same time~(e.g., for \CodeIn{floor} operator, tight abstraction $\langle \vl, \vu\rangle \mapsto \langle \lfloor \vl \rfloor, \lfloor \vu \rfloor \rangle$ is not generally differentiable, and differentiable abstraction $\langle \vl, \vu\rangle \mapsto \langle \vl - \mathbf{1}, \vu \rangle$ is not tight), we implement two abstract algorithms to achieve tightness and differentiability respectively.
   
    
    The existing approach DEBAR~\cite{zhang2020detecting} also proposes a static analysis framework for DNN architectures with tensor partitioning.
    In contrast to our framework, the abstract domain in DEBAR contains affine equalities besides interval domains.
    Therefore, DEBAR can be tighter than ours in some cases and can produce fewer positives, but due to the additional complexity of affine equalities, DEBAR supports fewer operators than ours.
    At the same time, our algorithms produce tightest interval abstractions while DEBAR has no tightness guarantee.
    As a result, \sysname usually detects more numerical defects and has a comparable number of false positives~(see \Cref{subsec:static-detector-exp}), and we can leverage whether failure-exhibiting system tests exist or not in \sysname to filter out false positives.
    
    
    \textbf{Computational Procedure.}
    We now discuss the computing procedure of the static analyis framework, and introduce our two novel techniques in detail in \Cref{adxsubsec:backward-labelling,adxsubsec:tight-abstract-dynamic-partitioning}.
   
    The input to the framework is a DNN architecture given as a computational graph $\gG = \langle \gV, \gE \rangle$ where the graph is directed and acyclic.
    First, we compute the abstractions for all initial nodes which corresponds to input, weight, and constant nodes.
    In this step, we construct abstractions based on their valid ranges~(for input and weight nodes) and concrete values~(for constant nodes)~(details in \Cref{adxsubsec:backward-labelling}).
    Then, we sort all nodes in topological order~(which is possible for a directed acyclic graph).
    Third, we execute our corresponding algorithms to compute the output abstractions for internal nodes based on operator types.
    Detail algorithms will be discussed in \Cref{adxsubsec:tight-abstract-dynamic-partitioning}.
    We enforce the soundness, tightness, and differentiability properties in these algorithms.
   
    By computing in topological order, we can assure that for the current node, all its input nodes have been processed and the corresponding input abstractions are available, so that computing output abstraction is feasible.
    
    Some DNN architectures may contain loops, such as recursive neural networks~(RNNs, \cite{rumelhart1986learning,hochreiter1997long}).
    In \CodeIn{ONNX} computational graph, RNNs are expressed by \CodeIn{Loop} nodes whose annotations contain the loop body, loop condition, and returned tensors.
    In \sysname, we rely on backward labelling~(see \Cref{adxsubsec:backward-labelling}) to conduct concrete execution of loop controlling variables and then unfold the loop to compute out the abstractions for returned tensors.
    
    \subsection{Initial Abstraction Construction with \textbf{Backward Fine-Grained Node Labelling}}
        \label{adxsubsec:backward-labelling}
        
        We construct abstract domains for initial nodes in two steps:
        First, we determine the tensor partitions, i.e., the split index sets $(S_i)_{i=1}^m$~(see \Cref{eq:abs-domain}), which decide the tightness and efficiency of our static analysis framework, because the tensor partitions of all other nodes will be solely dependent on the partitions of initial nodes as shown in \Cref{adxsubsec:tight-abstract-dynamic-partitioning}.
        Second, we compute the interval bounds $\vl$ and $\vu$.
        
        We use the following principle to decide their tensor partitions:
        For nodes that are connected to operators requiring fine-grained abstractions~(e.g., the \CodeIn{shape} input for operator \CodeIn{Reshape}) with valid paths which we will specify later, we construct tensor partitions with the finest granularity, i.e., $S_i = \{0,1,\dots,n_i-1\}$.
        We call these nodes as fine-grained initial nodes, and starting from next paragraph we introduce our novel technique of \emph{backward fine-grained node labelling} to find them out.
        For other nodes, we rely on downstream task requirements and user specifications to determine the partitions.
        For example, for precondition-fix generation~(see \Cref{subsec:precond-gen}) use coarsest granularity by default.
        
        \textbf{Backward Fine-Grained Node Labelling.}
        The fine-grained initial nodes are those starting a valid path in DNN computational graph $\gG$, where a path is valid if and only if the path does not traverse through fine-grained stopping operators and terminates at a fine-grained requiring operator with some specific input index.
        Fine-grained requiring operators are those taking indices as inputs~(so each index value is important and finest partition granularity is needed) and those controlling the loop execution~(so we can need finest granularity to know looping times to unfold).
        Fine-grained stopping operators are those which output is independent of input abstraction granularity~(so granularity of preceding nodes do not matter).
        Detail lists of fine-grained stopping operators and fine-grained requiring operators are provided in \Cref{adxsec:fine-grain-stop-requiring-operators}.
        
        To find out fine-grained initial nodes, we propose backward fine-grained node labelling, which is similar to data dependency analysis in traditional programs~\cite{hennessy2011computer}:
        First, we invert all edges of the given computational graph $\gG = \langle \gV, \gE \rangle$ to get $\gG' = \langle \gV, \gE' \rangle$.
        Second, we attach a boolean label for each node $n\in \gV$ and initialize them with \CodeIn{False}.
        Third, we do topology sort on $\gG'$.
        When encountering a node $n$ with \CodeIn{True}, if the node does not contain a fine-grained stopping operator, we propagate this label to all its subsequent nodes; otherwise, we propagate \CodeIn{False}.
        When encountering a node $n$ with \CodeIn{False}, if the node is a fine-grained requiring operator, we propagate \CodeIn{True} to corresponding subsequent nodes and \CodeIn{False} to others; otherwise propagate \CodeIn{False} to all subsequent nodes.
        Fourth, we collect all labels for initial nodes. 
        The nodes with \CodeIn{True} label are fine-grained initial nodes.
        
        To this point, we have determined the tensor partitions $(S_i)_{i=1}^m$ for all initial nodes, and we now compute $\vl$ and $\vu$ and thus finish the abstraction domain construction for the valid ranges of initial nodes.
        Initial nodes are further divided into input, weight, and constant nodes.
        In practice, most of weight nodes and all constant nodes have their initial values stored in the \CodeIn{ONNX} file, and we directly use \Cref{eq:alpha} with these initial values to compute $\vl$ and $\vu$.
        Otherwise, we rely on user specifications and built-in heuristics to determine $\vl$ and $\vu$.
    
    \subsection{Internal Abstraction with \textbf{Dynamic Partitioning}}
        \label{adxsubsec:tight-abstract-dynamic-partitioning}
        
        For all supported operator types~(listed in \Cref{adxsec:supported-op}), we propose concrete algorithms to compute abstract domains for output with our dynamic partitioning technique.
        Formally, for each operator $op$, we construct the computable function $T_{op}: \sA \to \sA$ that satisfies soundness and tentatively satisfies tightness and differentiability, where $\sA$ is the abstract domain defined in \Cref{eq:abs-domain}.
        Therefore, following the computational procedure introduced before, we can compute end-to-end abstractions for all nodes in the given DNN architecture.
        Given the limited space, we only describe the algorithms for these representative operators: \CodeIn{MatMul}, \CodeIn{Conv}, and \CodeIn{Softmax}.
        
        \paragraph{\CodeIn{MatMul}}
        The \CodeIn{MatMul} operator computes the matrix multiplication of two operands.
        This operator is widely used in DNNs to express fully-connected layers.
        To simplify the narration, we focus on two-dimensional case where we compute $op(\mA, \mB) = \mC = \mA \mB$ with $\mA \in \sR^{n \times m}$ and $\mB \in \sR^{m \times l}$, and extensions to other dimensions by abstracting the broadcasting mechanism can be found in our open-source implementation.
        
        We denote the input abstractions of $op$ by $a = \langle \mL_a, \mU_a, (S_A^1, S_A^2) \rangle$ and $b = \langle \mL_b, \mU_b, (S_B^1, S_B^2) \rangle$ respectively.
        First, we compute the union $U = S_A^2 \cup S_B^1$.
        Second, we dynamically partition both $a$ and $b$ with split points $(S_A^1, U)$ and $(U, S_B^2)$ respectively and get
        $a' = \langle \mL'_a, \mU'_a, (S_A, U) \rangle$ and $b' = \langle \mL'_b, \mU'_b, (U, S_B^2) \rangle$.
        Note that $a$ and $a'$~(or $b$ and $b'$) correspond to the same concrete domain, but $a'$ and $b'$ have finer or equal partition granularity than $a$ and $b$.
        Third, we compute output abstraction $T_{op}(a,b) = c = \langle \mL_c, \mU_c, (S_A^1, S_B^2) \rangle$
        \begin{equation}
            \small
            \label{eq:mat-abs-1}
            \begin{aligned}
                \text{where} \quad (\mL_c)_{ij} & = \sum_{k=1}^{|U|} \vv_k \min_{\mA \in \{\mL'_a, \mU'_a\}, \mB \in \{\mL'_b, \mU'_b\}} \mA_{ik}\mB_{kj}, \\
                (\mU_c)_{ij} & = \sum_{k=1}^{|U|} \vv_k \max_{\mA \in \{\mL'_a, \mU'_a\}, \mB \in \{\mL'_b, \mU'_b\}} \mA_{ik}\mB_{kj}, \\
                \vv_k & = U[k]-U[k-1].
            \end{aligned}
        \end{equation}
        This formulation can guarantee the tightness but is not efficient for tensor computation due to inner minimum and maximum.
        Therefore, we also implement a fast-mode abstraction which trade-off tightness for efficiency: $T_{op}'(a,b) = c' = \langle \mL_c', \mU_c', (S_A^1, S_B^2)\rangle$ where
        \begin{equation}
            \small
            \begin{aligned}
                \mL_c^{'} = &
                \mU_a^{'-} \vv \mU_b^{'-} + 
                \mU_a^{'0} \vv \mL_b^{'-} + 
                \mU_a^{'+} \vv \mL_b^{'-} + \\
                & \mL_a^{'-} \vv \mU_b^{'0} +
                \mU_a^{'0} \vv \mL_b^{'0} +
                \mL_a^{'0} \vv \mU_b^{'0}
                 +
                \mU_a^{'+} \vv \mL_b^{'0} + \\
                & \mL_a^{'-} \vv \mU_b^{'+} + 
                \mL_a^{'0} \vv \mU_b^{'+} +
                \mL_a^{'+} \vv \mL_b^{'+}, \\
               
                \mU_c^{'} = & 
                \mL_a^{'-} \vv \mL_b^{'-} +
                \mL_a^{'0} \vv \mL_b^{'-} +
                \mL_a^{'+} \vv \mU_b^{'-} + \\
                & \mL_a^{'-} \vv \mL_b^{'0} +
                \mL_a^{'0} \vv \mL_b^{'0} +
                \mU_a^{'0} \vv \mU_b^{'0}
                 +
                \mU_a^{'+} \vv \mU_b^{'0} + \\
                & \mU_a^{'-} \vv \mL_b^{'+} +
                \mU_a^{'0} \vv \mU_b^{'+} +
                \mU_a^{'+} \vv \mU_b^{'+}.
            \end{aligned}
            \label{eq:mat-abs-2}
        \end{equation}
        In the above equation, for any $* \in \{a, b\}$, let $\circ$ be the elementwise~(Hadamard) product,
        \begin{equation}
            \small
            \begin{aligned}
                %
                \mL_*^{'-} = & \mL_*^{'} \circ \1[\mU_*^{'} < 0],
                \mU_*^{'-} = \mU_*^{'} \circ \1[\mU_*^{'} < 0], \\
                \mL_*^{'0} = & \mL_*^{'} \circ \1[\mL_*^{'} \le 0, \mU_*^{'} \ge 0],
                \mU_*^{'0} = \mU_*^{'} \circ \1[\mL_*^{'} \le 0, \mU_*^{'} \ge 0], \\
                \mL_*^{'+} = & \mL_*^{'} \circ \1[\mL_*^{'} > 0],
                \mU_*^{'+} = \mU_*^{'} \circ \1[\mL_*^{'} > 0].
               
               
               
               
               
               
               
               
            \end{aligned}
            \label{eq:mat-abs-2-defs}
        \end{equation}
        From \Cref{eq:mat-abs-2}, we can observe that the abstraction can be easily implemented with tensor computations.
        In \Cref{adxsec:proofs}, we prove the soundness and tightness of these abstractions.
        
        \paragraph{\CodeIn{Conv}}
        The \CodeIn{Conv} operator computes the discrete convolution of two operands.
        This operator is widely used in convolutional neural networks~(CNNs, \cite{krizhevsky2012imagenet}).
        To simplify the narration, we focus on single-channel single-stride two-dimensional case where we compute $op(\mA, \mW) = \mC$ with $\mA$ being the input matrix and $\mW$ being the convolution kernel. 
        Extensions to general cases are provided in our open-source code.
        
        We first dynamically split the kernel abstraction to the finest granularity.
        Second, we compute the receptive field of each output position, which is the sub-region of $\mA$ that decides each output position.
        Third, we inspect the alignment between receptive fields and $\mA$'s partitions.
        If neighboring output positions have their receptive fields partitioned in the same sub-block of $\mA$, it means that these positions can be abstracted by a single interval, i.e., these positions can be partitioned together.
        Fourth, using this principle, we derive the paritions of the output tensor, and repeat the input tensors accordingly so that the abstract computation can be written as convolutional operations.
        Last, we modify the abstraction computation equations in \Cref{eq:mat-abs-2} by replacing matrix multiplication with convolution to compute the abstraction of output $\mC$.
        
        \paragraph{\CodeIn{Softmax}}
        The \CodeIn{Softmax} operator computes normalized exponential values for the given input.
        \CodeIn{Softmax} is widely deployed for classification tasks to yield the normalized confidence.
        In one-dimensional case, for input $\vx \in \sR^n$, a \CodeIn{Softmax} operator outputs $op(\vx) = \frac{\exp(\vx)}{\sum_{i=1}^n \exp(\vx)_i}$.
        The output abstraction of \CodeIn{Softmax} operator $op$ can be thus computed: $T_{op}(\langle \vl, \vu, (S)\rangle) = \langle \vl^o, \vu^o, (S)\rangle$ where
        \begin{equation}
            \small
            \resizebox{0.9\linewidth}{!}{$ \displaystyle
            \begin{aligned}
                \vl^o_i & = \dfrac{\exp(\vl_i)}{\exp(\vu)^\T \vv - \exp(\vu_i) + \exp(\vl_i)}, 
                \vu^o_i = \dfrac{\exp(\vu_i)}{\exp(\vl)^\T \vv - \exp(\vl_i) + \exp(\vu_i)}, \\
                \vv_k & = S[k] - S[k-1].
            \end{aligned}
            $}
            \label{eq:softmax-abs}
        \end{equation}
        The output abstraction's partition is dynamically decided by the input abstraction's partition.
        We prove the soundness and tightness of the abstraction in \Cref{adxsec:proofs}.
        
        \begin{remark}
            As we can see, the novel dynamic partition technique is incorporated in the computation process of each operator's abstraction.
            The soundness and tightness of our designed abstractions are immediately achieved by design, and the differentiability of our designed abstractions relies on the auto-differentiation functionality of  popular DL libraries like \CodeIn{PyTorch}~\cite{paszke2019pytorch} and \CodeIn{Tensorflow}~\cite{tensorflow2015-whitepaper} where we use \CodeIn{PyTorch} for implementation. 
            The abstractions for DNNs are also implemented for other applications, e.g., for robustness verification~\cite{singh2019abstract,mirman2018differentiable,xu2020automatic}.
            However, our abstractions are tailored for tensor partitioned interval domain which is particularly suitable for testing and debugging for numerical failures as discussed before.
            To the best of our knowledge, these abstractions are the first that achieve soundness, tightness, and differentiability.
        \end{remark}
    
    \subsection{List of Fine-Grained Requiring and Stopping Operators}
    
        \label{adxsec:fine-grain-stop-requiring-operators}
        
        We use fine-grained requiring and fine-grained stopping operators in \Cref{adxsubsec:backward-labelling}.
        Among all supported operators, the fine-grained requiring operators are \\
        \begin{tabular}{p{8.1cm}}
            \CodeIn{Reshape} (input \#2), 
            \CodeIn{Slice} (input \#2, \#3, \#4, \#5), 
            \CodeIn{Squeeze} (input \#2), 
            \CodeIn{Unsqueeze} (input \#2), 
            \CodeIn{Tile} (input \#2, \#3), 
            \CodeIn{Loop} (input \#1, \#2),
            \CodeIn{SequenceInsert} (input \#3),
            \CodeIn{ConstantOfShape} (input \#1),
            \CodeIn{Gather} (input \#2),
            \CodeIn{GatherND} (input \#2),
            \CodeIn{ReduceSum} (input \#2),
            \CodeIn{ScatterElements} (input \#2),
            \CodeIn{Expand} (input \#2),
            \CodeIn{Split} (input \#2),
            \CodeIn{Pad} (input \#2, \#3),
            \CodeIn{NegativeLogLikelihoodLoss} (input \#2),
            \CodeIn{Clip} (input \#2, \#3),
            \CodeIn{OneHot} (input \#2),
            \CodeIn{Resize} (input \#2, \#3, \#4).
        \end{tabular}\\
        The fine-grained stopping operators are \CodeIn{Shape}, \CodeIn{RandomNormalLike}, and \CodeIn{RandomUniformLike}.
            
    \subsection{Proofs}
        \label{adxsec:proofs}
        \allowdisplaybreaks
        
        Here we present the omitted proofs in \Cref{adxsubsec:tight-abstract-dynamic-partitioning}.
        
        \subsubsection{\CodeIn{MatMul}}
        
        \begin{theorem}[tightness of first abstraction]
            Suppose $op$ is the \CodeIn{MatMul} operator and $T_{op}$ is as defined in \Cref{eq:mat-abs-1}, then if $op(\gamma(a), \gamma(b)) \subseteq [\mL, \mU]$ for $a, b \in \sA$, $\gamma(T_{op}(a, b)) \subseteq [\mL, \mU]$.
            \label{thm:1}
        \end{theorem}
        
        \begin{proof}
            We use $\mL_a$ and $\mU_a$ to denote the element-wise interval lower and upper bounds of $\gamma(a)$; and use $\mL_b$ and $\mU_b$ to denote these bounds of $\gamma(b)$ respectively, where a formal definition is in \Cref{eq:gamma}.
            Then, there exist $\mA$ and $\mB$ with $\mL_a \le \mA \le \mU_a$ and $\mL_b \le \mB \le \mU_b$, such that
            \begin{small}
                \begin{align}
                    (\mA\mB)_{ij} & = \sum_{k=1}^l \min \{ \mL_{a,ik} \mL_{b,kj}, \mL_{a,ik} \mU_{b,kj}, \mU_{a,ik} \mL_{b,kj}, \mU_{a,ik} \mU_{b,kj} \}, \label{eq:pf-1} \\
                    \text{or } (\mA\mB)_{ij} & = \sum_{k=1}^l \max \{ \mL_{a,ik} \mL_{b,kj}, \mL_{a,ik} \mU_{b,kj}, \mU_{a,ik} \mL_{b,kj}, \mU_{a,ik} \mU_{b,kj} \}. \label{eq:pf-2} 
                \end{align}
            \end{small}
            By definition, we have
            \begin{equation}
                \small
                \mL_{ij} \le \text{\Cref{eq:pf-1}}, \mU_{ij} \ge \text{\Cref{eq:pf-2}}.
            \end{equation}
            We let $\mL'_c$ and $\mU'_c$ to denote the element-wise interval lower and upper bounds for $\gamma(T_{op}(a,b))$, 
            let $S_A^1[i'] \le i-1 < S_A^1[i'+1]$ and $S_B^2[j'] \le j-1 < S_B^2[j'+1]$,
            then from \Cref{eq:mat-abs-1},
            \begin{equation}
                \small
                \begin{aligned}
                & (\mL'_c)_{ij} = (\mL_c)_{i'j'} \\
                = & \sum_{k=1}^{|U|} \sum_{k'=U[k-1]+1}^{U[k]} \min \{ \mL_{a,ik'} \mL_{b,k'j}, \mL_{a,ik'} \mU_{b,k'j}, \mU_{a,ik'} \mL_{b,k'j}, \mU_{a,ik'} \mU_{b,k'j} \} \\
                = & \sum_{k=1}^l \min \{ \mL_{a,ik} \mL_{b,kj}, \mL_{a,ik} \mU_{b,kj}, \mU_{a,ik} \mL_{b,kj}, \mU_{a,ik} \mU_{b,kj} \} \\
                = & \text{\Cref{eq:pf-1}} \ge \mL_{ij}.
                \end{aligned}
                \label{eq:pf-5}
            \end{equation}
            Similarly, $(\mU'_c)_{ij} \le \mU_{ij}$.
            Thus, $\gamma(T_{op}(a,b)) \subseteq [\mL, \mU]$.
        \end{proof}
        
        \begin{theorem}[soundness of first abstraction]
            Suppose $op$ is the \CodeIn{MatMul} opeartor and $T_{op}$ is as defined in \Cref{eq:mat-abs-1}, then $op(\gamma(a), \gamma(b)) \subseteq \gamma(T_{op}(a,b))$.
            \label{thm:2}
        \end{theorem}
        
        \begin{proof}
            We use $\mL_a$ and $\mU_a$ to denote the element-wise interval lower and upper bounds of $\gamma(a)$; and use $\mL_b$ and $\mU_b$ to denote these bounds of $\gamma(b)$ respectively, where a formal definition is in \Cref{eq:gamma}.
            For any $\mA$ and $\mB$ such that $\mL_a \le \mA \le \mU_a$ and $\mL_b \le \mB \le \mU_b$, 
            \begin{equation}
                \small
                \begin{aligned}
                    (\mA\mB)_{ij} & = \sum_{k=1}^l \mA_{ik} \mB_{kj} \\
                    & \ge \min\{ \mL_{a,ik} \mL_{b,kj}, \mL_{a,ik} \mU_{b,kj}, \mU_{a,ik} \mL_{b,kj}, \mU_{a,ik} \mU_{b,kj} \} =: (\mL_{ab}')_{ij}
                \end{aligned}
                \label{eq:pf-3} 
            \end{equation}
            and
            \begin{equation}
                \small
                \begin{aligned}
                    (\mA\mB)_{ij} & \le \max\{ \mL_{a,ik} \mL_{b,kj}, \mL_{a,ik} \mU_{b,kj}, \mU_{a,ik} \mL_{b,kj}, \mU_{a,ik} \mU_{b,kj} \} =: (\mU_{ab}')_{ij}.
                \end{aligned}
                \label{eq:pf-4} 
            \end{equation}
            Thus, $op(\gamma(a), \gamma(b)) \subseteq [\mL'_{ab}, \mU'_{ab}]$.
            On the other hand, $\gamma(T_{op}(a,b)) = [\mL'_{ab}, \mU'_{ab}]$ as seen from \Cref{eq:pf-5}.
            Therefore, $op(\gamma(a), \gamma(b)) \subseteq \gamma(T_{op}(a,b))$.
        \end{proof}
    
        \begin{theorem}[soundness of second abstraction]
            Suppose $op$ is the \CodeIn{MatMul} operator and $T'_{op}$ is as defined in \Cref{eq:mat-abs-2,eq:mat-abs-2-defs}, then $op(\gamma(a), \gamma(b)) \subseteq \gamma(T'_{op}(a,b))$.
            \label{thm:3}
        \end{theorem}
        
        \begin{proof}
            We use $\mL_a$ and $\mU_a$ to denote the element-wise interval lower and upper bounds of $\gamma(a)$; and use $\mL_b$ and $\mU_b$ to denote these bounds of $\gamma(b)$ respectively, where a formal definition is in \Cref{eq:gamma}.
            We define $\mL'_{ab}$ and $\mU'_{ab}$ by \Cref{eq:pf-3,eq:pf-4}.
            From the proof of \Cref{thm:2}, we have $op(\gamma(a), \gamma(b)) \subseteq [\mL'_{ab}, \mU'_{ab}]$.
            We now only need to show that $[\mL'_{ab}, \mU'_{ab}] \subseteq \gamma(T'_{op}(a,b))$.
            
            $\gamma(T'_{op}(a,b))$ imposes independent interval abstractions element-wise, therefore, we study each element independently.
            For the element $(i,j)$, from \Cref{eq:mat-abs-2,eq:mat-abs-2-defs}, the interval lower bound of $\gamma(T'_{op}(a,b))$, namely $(\mL_c)_{ij}$, satisfies
            \begin{equation}
                \small
                \begin{aligned}
                    (\mL_c)_{ij} = \sum_{k=1}^l & 
                    (\mU_a^-)_{ik} (\mU_b^-)_{kj} +
                    (\mU_a^{0})_{ik} (\mL_b^{-})_{kj} + 
                    (\mU_a^{+})_{ik} (\mL_b^{-})_{kj} + \\
                    & (\mL_a^{-})_{ik} (\mU_b^{0})_{kj} +
                    (\mU_a^{0})_{ik} (\mL_b^{0})_{kj} +
                    (\mL_a^{0})_{ik} (\mU_b^{0})_{kj} +
                    (\mU_a^{+})_{ik} (\mL_b^{0})_{kj} + \\
                    & (\mL_a^{-})_{ik} (\mU_b^{+})_{kj} + 
                    (\mL_a^{0})_{ik} (\mU_b^{+})_{kj} +
                    (\mL_a^{+})_{ik} (\mL_b^{+})_{kj}.
                \end{aligned}
                \label{eq:pf-6}
            \end{equation}
            By \Cref{eq:mat-abs-2-defs},
            \begin{itemize}
                \item when $\mL_{a,ik} \le \mU_{a,ik} < 0$, \\
                $(\mL_a^-)_{ik} = \mL_{a,ik}, (\mU_a^-)_{ik} = \mU_{a,ik}, 
                (\mL_a^0)_{ik} = 0, \\ (\mU_a^0)_{ik} = 0, 
                (\mL_a^+)_{ik} = 0, (\mU_a^+)_{ik} = 0$;
                
                \item when $\mL_{a,ik} \le 0 \le \mU_{a,ik}$, \\
                $(\mL_a^-)_{ik} = 0, (\mU_a^-)_{ik} = 0, 
                (\mL_a^0)_{ik} = \mL_{a,ik}, \\ (\mU_a^0)_{ik} = \mU_{a,ik}, 
                (\mL_a^+)_{ik} = 0, (\mU_a^+)_{ik} = 0$;
                
                \item when $0 < \mL_{a,ik} \le \mU_{a,ik}$, \\
                $(\mL_a^-)_{ik} = 0, (\mU_a^-)_{ik} = 0, 
                (\mL_a^0)_{ik} = 0, \\ (\mU_a^0)_{ik} = 0, 
                (\mL_a^+)_{ik} = \mL_{a,ik}, (\mU_a^+)_{ik} = \mU_{a,ik}$.
            \end{itemize}
            Similarly for $(\mL_b^-)_{kj}$, $(\mU_b^-)_{kj}$, $(\mL_b^0)_{kj}$, $(\mU_b^0)_{kj}$, $(\mL_b^+)_{kj}$ and $(\mU_b^+)_{kj}$.
            Thus, by enumerating all cases, we have
            \begin{equation}
                \small
                \text{\Cref{eq:pf-6}} \le \min\{ \mL_{a,ik} \mL_{b,kj}, \mL_{a,ik} \mU_{b,kj}, \mU_{a,ik} \mL_{b,kj}, \mU_{a,ik} \mU_{b,kj} \} = (\mL'_{ab})_{ij}.
            \end{equation}
            Similarly, the interval lower bound of $\gamma(T'_{op}(a,b))$, namely $(\mU_c)_{ij}$,
            \begin{equation}
                \small
                (\mU_c)_{ij} \ge \max\{ \mL_{a,ik} \mL_{b,kj}, \mL_{a,ik} \mU_{b,kj}, \mU_{a,ik} \mL_{b,kj}, \mU_{a,ik} \mU_{b,kj} \} = (\mU'_{ab})_{ij}.
            \end{equation}
            Thus, $[\mL'_{ab}, \mU'_{ab}] \subseteq \gamma(T'_{op}(a,b))$.
        \end{proof}
        
        \subsubsection{\CodeIn{Softmax}}
        
        \begin{theorem}[tightness]
            Suppose $op: \sR^n \to \sR^n$ is the \CodeIn{Softmax} operator and $T_{op}$ is as defined in \Cref{eq:softmax-abs}, then if $op(\gamma(\langle \vl, \vu, (S)\rangle)) \subseteq [\vl^r, \vu^r]$, then $\gamma(T_{op}(\langle \vl, \vu, (S)\rangle)) \subseteq [\vl^r, \vu^r]$.
            \label{thm:4}
        \end{theorem}
        
        \begin{proof}
            We use $\vl'$ and $\vu'$ to denote the element-wise interval lower and upper bounds of $\gamma(\langle \vl, \vu, (S)\rangle)$.
            Formally, $\vl'_i = \vl_{i'}$ and $\vu'_i = \vu_{i'}$ where $S[i'] \le i - 1 < S[i'+1]$.
            Then, 
           
            for each $i$, by setting $\vx_i = \vl_i$ and $\vx_j = \vu_j$ for all $j \neq i$, we get
            \begin{equation}
                \small
                \begin{aligned}
                    \vl_i^r & \le \dfrac{\exp(\vl'_i)}{\exp(\vl'_i) + \sum_{k=1, k\neq i}^n \exp(\vu'_k)}
                    = \dfrac{\exp(\vl_{i'})}{\exp(\vl_{i'}) + \sum_{k=1}^n \exp(\vu'_k) - \exp(\vu'_i)} \\
                    & = \dfrac{\exp(\vl_{i'})}{\exp(\vl_{i'}) + \vv^\T \exp(\vu_k) - \exp(\vu_{i'})} = \vl_{i'}^o.
                \end{aligned}
            \end{equation}
            In addition, by setting $\vx_i = \vu_i$ and $\vx_j = \vl_j$ for all $j \neq i$, we get $\vu_i^r \ge \vu_{i'}^o$.
            Here, $\vl^o$ and $\vu^o$ are as defined in \Cref{eq:softmax-abs}.
            Combining these two arguments, we get 
            $$\gamma(T_{op}(\langle \vl, \vu, (S)\rangle)) \subseteq [\vl^r, \vu^r].$$
        \end{proof}
    
        \begin{theorem}[soundness]
            Suppose $op: \sR^n \to \sR^n$ is the \CodeIn{Softmax} operator and $T_{op}$ is as defined in \Cref{eq:softmax-abs}, then $op(\gamma(\langle \vl, \vu, (S)\rangle)) \subseteq \gamma(T_{op}(\langle \vl, \vu, (S)\rangle))$.
            \label{thm:5}
        \end{theorem}
        
        \begin{proof}
            Leveraging the fact that softmax function is monotinically increasing, i.e., $\frac{\d op(\vx)_i}{\d \vx_j} > 0$, we have
            $\gamma(\langle \vl, \vu, (S)\rangle)_i \in [\vl^o_{i'}, \vu^o_{i'}]$,
            where $S[i'] \le i-1 < S[i'+1]$.
            Since $[\vl^o_{i'}, \vu^o_{i'}] = \gamma(T_{op}(\langle \vl, \vu, (S)\rangle))_i$ by our definition of $T_{op}$, $op(\gamma(\langle \vl, \vu, (S)\rangle)) \subseteq \gamma(T_{op}(\langle \vl, \vu, (S)\rangle))$ follows.
        \end{proof}






\section{Hyperparameters}

    \label{adxsec:detail-hyperparam}

    In this appendix, we listed the omitted hyperparameters in our implementation.
    
    \paragraph{System Test Generation}
    The system test generation approach in \sysname is introduced in \Cref{subsec:system-test-gen}.
    The approach is composed of two steps:  generation of unit test $\langle\vwinfer, \vxinfer\rangle$, and generation of training example $\vxtrain$.
    For the unit test generation, we use the Adam optimizer~\cite{kingma2015adam} with learning rate $1$ and maximum iteration $100$.
    For the training test generation, we target for training example under learning rate $\gamma=1$ and the approach has similar performance under other learning rates.
    We follow the convention in DLG attack, where we use the L-BFGS~\cite{liu1989limited} method as the optimizer for gradient-based minimization.
    We terminate the method and return ``\CodeIn{failed}'' if either the running time exceeds \SI{180}{s}, or a failure-exhibiting training is not found after $300$ iterations of L-BFGS optimization.
    
    \paragraph{Precondition-Fix Suggestion}
    Key hyperparameters for \sysname precondition-fix suggestion are presented in \Cref{subsec:precond-gen}.
    We execute the optimization for at most $1000$ iterations and stop it when the running exceeds \SI{180}{s}.
    If no valid precondition fix is found, we return ``\CodeIn{failed}''.
    
    More hyperparameters and implementation details can be found in our attached code.
    
\section{ablation Study of System Test Generation}

    \label{adxsec:system-testgen-by-stage}
    
    
   
    \sysname contains the novel technique of two-step generation for generating failure-exhibiting system tests:
    it first generates failure-exhibiting unit tests with gradient back-propagation, then generates failure-exhibiting training example via the extended DLG attack.
    To isolate the impact of \sysname at each step, we replace either step with random sampling:
    ``Random + \sysname'' which first generates failure-exhibiting unit tests via random sampling, then generates training example via \sysname at the second stage;
    ``\sysname + Random'' which first generates failure-exhibiting unit tests via \sysname, then generates training example via random sampling.
    We follow the same evaluation protocol in this ablation study.
    We find that ``Random + \sysname'' takes $6\times$ more running time than \sysname but succeeds for only $818$ runs~($8.3\%$ fewer than \sysname); and similarly ``\sysname + Random'' succeeds for only $848$ runs~($5.2\%$ fewer than \sysname).
    This ablation study implies that, \sysname's technique helps to improve effectiveness and efficiency at both steps, and the improvement for the first step is slightly more significant than that for the second step.
   
   
    
\section{Detail Evaluation Results on Precondition-Fix Suggestion}

    \label{adxsec:detail-precond-gen-exp}

    \input{exp/tables/precond_table_weight_input}
    \input{exp/tables/precond_table_weight}
    \input{exp/tables/precond_table_input}

    In \Cref{subsec:precond-gen-exp}, due to the page limit, we only presented the summary statistics of the comparison between \sysname and baseline approaches~(in \Cref{tab:precond-summary}).
    In this appendix, we present detail comparison results for each case.
    \Cref{tab:precond_weight_input} shows results when precondition fixes can be imposed on both weight and inputs nodes.
    \Cref{tab:precond_weight} shows results when precondition fixes can be imposed on only weight nodes.
    \Cref{tab:precond_input} shows results when precondition fixes can be imposed on only input nodes.
    Discussion of these results can be found in \Cref{subsec:precond-gen-exp}.














