Sound Mixed Fixed-Point Quantization of Neural Networks

1 INTRODUCTION

Neural networks (NN) are increasingly being considered as components in safety-critical systems. For example, neural network controllers have been shown to be effective for a variety of closed-loop systems, including simple car and airplane models, adaptive cruise control, and aircraft collision avoidance [15, 16, 34]. Ensuring their correctness for safety-critical applications is nontrivial, but several recent approaches and tools can already automatically verify the safety of limited-size but interesting systems [27, 31, 49, 51].

To make such neural networks practical, they need to be safe but also efficient. While NNs are typically trained in high-precision floating-point arithmetic on server-like machines with graphics processing units (GPUs), this high precision can be prohibitively expensive for resource-constrained embedded systems. To obtain efficiency, e.g., in terms of area, latency, or memory usage, the trained NNs are quantized to use, for example, low-precision fixed-point arithmetic [21, 23].

The increased efficiency due to low precision comes at the cost of reduced accuracy of the NN computations since each operation potentially incurs a (larger) roundoff error. To ensure the overall correctness of systems with NNs, we thus need to choose the precision for quantization such that the safety of the overall system can still be guaranteed, i.e., the roundoff error is within some application-specific bound. While the safety proofs of NN controllers, for instance, do not take into account finite-precision roundoff errors directly and assume exact real-valued arithmetic for performance reasons, they typically allow to account for bounded errors stemming from a noisy environment or from inaccurate implementations.

Numerous approaches have been proposed for the quantization of neural networks, demonstrating promising results on standard machine learning benchmarks [21, 23, 37, 44, 47, 48]. However, these techniques are not applicable to safety-critical closed-loop control systems for two primary reasons. First, they are specifically tailored for neural network classifiers, focusing on dynamically comparing classification accuracy on specific test datasets. As a result, they cannot handle neural network controllers that do not perform classifications tasks and rather implement regression tasks and compute (continuous) control values. Secondly, they are not sound which means they cannot guarantee (classification) accuracy for all possible inputs, which is essential for safety-critical systems.

On the other hand, the available tools that do focus on sound quantization or precision tuning are primarily designed for general-purpose arithmetic code [10, 11]. When applied to neural networks, they often provide overly conservative results or are inefficient, as they do not leverage the unique characteristics, structures and optimizations specific to neural networks.

In this paper, we present the first sound and fully automated mixed-precision quantizer for fixed-point arithmetic that specifically targets fully connected feed-forward deep neural networks. We support an arbitrary number of layers with ReLU or linear activation functions and efficiently minimize the number of bits needed to implement a network while guaranteeing a provided error bound wrt. an idealized real-valued implementation. While our implementation targets these networks as they are sufficient for many circumstances in neural network controllers (all our benchmarks are real-world controllers that were used in the context of sound verification), the proposed technique is not fundamentally limited to feed-forward networks; it handles networks with (sparse) matrix multiplications and linearized activations.

We focus on mixed-precision fixed-point arithmetic because it allows for efficient implementations on resource-constrained systems. Fixed-point arithmetic can be implemented with standard integer operations (in particular it avoids specialized floating-point hardware), and when implemented on configurable hardware such as FPGAs allows operations to use an arbitrary number of bits. Thus, unlike floating-point arithmetic, which is limited to typically 16, 32, or 64 bits of precision, fixed-point arithmetic allows for a larger scope for optimization.

Optimizing mixed-precision fixed-point arithmetic efficiently and accurately is challenging for two main reasons. First, fixed-point arithmetic is fundamentally discrete, and the continuous abstractions that allow to use efficient continuous optimization techniques for floating-point arithmetic [10] are not applicable. Secondly, due to the large choice (of combinations) of different precisions for individual operations, the search space for quantization is enormous and heuristic search techniques can explore this space only limitedly.

To overcome these challenges, we phrase sound fixed-point quantization of feed-forward neural networks as a mixed integer linear programming (MILP) problem. A naive formulation, however, results in non-linear constraints that are computationally intractable. We show how to use over-approximations to relax the problem to linear constraints that are quickly solvable. Furthermore, we leverage the special structure of NNs to efficiently, yet accurately to optimize the dot product operations by building on an existing technique for their correctly rounded implementation [13].

We implement our approach in a prototype tool called Aster that we will release as open source. Aster takes as input a trained neural network written in a small real-valued domain-specific language, a specification of the possible inputs, and an error bound on the result, and outputs a mixed fixed-point precision assignment that is guaranteed to satisfy the error bound (wrt. the real-valued input) and to minimize a (customizable) cost function.

We evaluate Aster on NN embedded controller benchmarks used for verification from the literature [22, 34, 40] and compare it with an existing sound quantizer [11]. Our results show that Aster can generate feasible implementations for significantly more benchmarks. We also show that for most of our benchmarks, Aster’s implementations take fewer machine cycles when compiled by Xilinx Vivado HLS [52] for an FPGA. In addition, Aster is substantially faster, especially for larger networks with thousands of parameters, with improvements in optimization time on average of \(\sim\) 67%.

Contributions. To summarize, in this paper we present:

(i) a novel MILP-based mixed fixed-point quantization approach that guarantees a given roundoff error bound,

(ii) an experimental comparison against the state-of-the-art to demonstrate the effectiveness of our approach on a benchmark set of neural network controllers, and

(iii) a prototype tool called Aster available on Zenodo: https://doi.org/10.5281/zenodo.8123416.

2 OVERVIEW AND BACKGROUND

Consider a unicycle model of a car [34] that models the dynamics of a car with 4 parameter variables—the Cartesian coordinates (x, y), the speed, and the steering angle. We use a neural network that was trained as a controller for this model as our initial example. It is a fully connected feed-forward neural network with 1 hidden layer. The network is fed as input a 4-parameter vector, denoted by \(\overline{x}_0 = [\!\!\!\begin{array}{llll} x_0^1 & x_0^2 & x_0^3 & x_0^4 \end{array}\!\!\!]\) where each parameter is constrained by real-valued intervals: \(x_0^1 \in [-0.6, 9.55]\), \(x_0^2 \in [-4.5, 0.2]\), \(x_0^3 \in [-0.6, 2.11]\), and \(x_0^4 \in [-0.3, 1.51]\), each specifying valid input values to the network. The inputs are propagated through each layer by successive application of the dot product operations, bias additions and activation functions: (1) \(\begin{align} \text{layer 1:} \quad \overline{x}_1 = ReLU(W_1 \cdot \overline{x}_0 + \overline{b}_1) \quad \quad \text{output:} \quad \overline{x}_2 = Linear(W_2 \cdot \overline{x}_1 + \overline{b}_2) \end{align}\) \(\overline{x}_1\) is the output of the first layer, which is then fed as input to the next, output layer. The output of the overall network is \(\overline{x}_2\). Each layer is parameterized by a weight matrix \(W_i\) and a bias vector \(\overline{b_i}\): \(\begin{align*} W_1 &= \begin{bmatrix} -0.037 & \cdots & 0.129 \\ -0.003 & \cdots & 0.099 \\ -0.128 & \cdots & 0.047\\ 0.045 & \cdots & -0.166 \end{bmatrix},\; \overline{b}_1 = [\!\!\!\begin{array} {ccc} 0.028 & \cdots & 0.342 \end{array}\!\!\!],\quad W_2 = \begin{bmatrix} 0.052 & 0.137\\ \cdots & \cdots \\ -0.200 & 0.154 \\ \end{bmatrix}, \overline{b}_2 = [\!\!\!\begin{array} {cc} 0.273 & 0.304 \end{array}\!\!\!] \end{align*}\) While the network has only a single hidden layer, it has 500 neurons (elided above) and thus results in a non-trivial size of the overall network.

To be deployed in a safety-critical system, such a neural network controller needs to be proven safe before implementation, for instance, one may need to prove that the system reaches a set of safe states within a given time window [34]. Since exact reasoning about finite-precision arithmetic does not scale, existing verification techniques assume real-valued parameters and arithmetic operations for the network [27, 31, 49, 51], but can typically deal with bounded uncertainties, from the implementation or measurements. We will thus assume that a bound on the output error \(\epsilon\) is given.

As a controller is primarily executed on resource-constrained hardware, using floating-point arithmetic may be overly expensive, as it requires either additional floating-point processor support or slow software emulations. The alternative is to quantize the NN controller in fixed-point arithmetic.

Fixed-point Quantization. In a fixed-point implementation, all program variables and constants are implemented using integers, and have an (implicit) representation < \(Q,\pi\)>, consisting of the total word length \(Q \in \mathbb {N}\) (overall number of bits including a sign bit), and the position of the binary point \(\pi \in \mathbb {N}\) counting from the least significant bit. Arithmetic operations on fixed-point variables can be implemented efficiently using only integer arithmetic and bit-shifting [3], or can use equivalent efficient hardware implementations, e.g., on FPGAs. While fixed-point arithmetic is not standardized, we employ the most commonly used fixed-point representation [55] and review the parts relevant for our paper here.

This representation effectively divides the number of overall bits into an integer part \(I = Q-\pi - 1\), and a fractional part \(\pi\). The integer part needs enough bits to ensure that it can hold large enough values to not overflow, i.e., the range of representable numbers is [\(-2^{I}, 2^{I}\)]. The fractional part controls the precision of a variable or an operation—the larger the number of fractional bits, the more precisely can the value be represented. Assuming that each operation uses truncation as the rounding mode, the maximum roundoff error with \(\pi\) fractional bits is \(2^{-\pi }\). In this paper, we use truncation as the rounding mode as it is the most commonly used one (the default with circuit design and synthesis compilers like Xilinx), and more efficient at runtime than e.g., rounding to nearest.

In our example (Equation (1)), the input \(x^{1}_0\) is in the range [\(-0.6, 9.55\)]. The number of integer bits required to hold this range, i.e., to represent the maximum absolute value 9.55 without overflow, is 5 (including the sign bit). Assuming 32 bits are available for wordlength (\(Q = 32\)), then we have \((32-5) = 27\) bits remaining for the fractional part. Hence, the maximum roundoff error for the input is \(2^{-27}\).

We need to assign each operation enough integer and fractional bits to guarantee overflow-freedom and sufficient overall accuracy of the final result. To ensure that the overall roundoff error does not exceed a given bound \(\epsilon\), we need to propagate and accumulate the errors of individual operations, which in general happens in nontrivial ways and is challenging to do manually. At the same time, we want to assign only as many bits as are actually needed to not waste resources.

Mixed-Precision Tuning. Using uniform fixed-point precision, i.e., the same word length for all operations, can be suboptimal; if one precision is not enough even at a single point in the program, we need to upgrade all operations to the next higher precision. However, not all layers necessarily have the same effect on the overall accuracy. Hence, it can be more resource-efficient to assign different precisions (word lengths) to different operations to achieve a target error bound, and thus implement the model in mixed precision.

Existing sound techniques applicable to fixed-point arithmetic rely on a heuristic search that repeatedly evaluates roundoff errors for different precision assignments [11]. While this technique works well for small programs, the repeated global roundoff error analysis quickly becomes expensive as we show in our evaluation in Section 5. Moreover, these techniques are designed for general-purpose straight-line programs, and do not take into account the structure of neural networks. To apply existing tools, one needs to assign all weight matrix and bias vector elements to individual scalar variables, and unroll all loops over these, resulting in an enormous straight-line expression.

For the unicycle example, only computing the roundoff error for a uniform 32-bit implementation using the state-of-the-art tool Daisy [11] takes 5.91 minutes, and mixed-precision tuning of this example takes more than 2.7 hours.

Our Approach. We encode fixed-point precision tuning as a mixed integer linear programming (MILP) problem, and perform several over-approximations to generate a linearized problem from the fundamentally nonlinear constraints (see Section 3). Our MILP constraints optimize the number of fractional bits such that overflow-freedom is ensured even in the presence of errors and a cost function is minimized. Our approach is parametric in the cost function to be optimized, but for the purpose of evaluation we follow Daisy’s cost function and count the total number of bits needed.

Assuming \(\epsilon = 0.001\) as the error bound for the unicycle example, Aster with our MILP-based mixed-precision tuning assigns different precisions (using 18, 19, 20, 21, 24, 30, and 34 bits) to different variables, constants, and operations in just under 50 seconds (i.e. significantly less than the 2.7 hours taken by Daisy). When compiled for an FPGA architecture with Xilinx’ Vivado HLS tool, Aster’s generated code takes only 27 machine cycles to execute, whereas Daisy’s generated code, both uniform and mixed-precision, take 178 cycles.

3 MILP-BASED MIXED-PRECISION TUNING

In this work we specifically focus on feed-forward neural networks with ‘relu’ and ‘linear’ activation functions that solve regression problems and compute continuous outputs. Such networks may, for instance, be utilized in closed-loop control systems, where they can be proven safe with different techniques [27, 31, 49, 51] that typically assume real-valued arithmetic operations, inputs and constants. Our objective is to minimize the resource usage of the implemented networks, while preserving safety by ensuring that computed (control) outputs remain within a specified error bound (that e.g., arises from a safety proof).

Problem Definition. Given a real-valued neural network architecture, ranges of inputs, and a roundoff error bound at the output, the goal is to generate a fixed-point mixed-precision neural network implementation that minimizes the precision of the variables and constants while ensuring that the roundoff errors at the output remain within the specified bound.

Our approach is inspired by successful mixed-precision optimization techniques for floating-point programs [10] that phrase precision tuning as an optimization problem (which it ultimately is). They rely on the dynamic ranges of floating-point arithmetic that allows to bound floating-point roundoff errors by nonlinear continuous abstractions, which, in turn, enables continuous, purely real-valued optimization techniques for precision tuning.

However, such continuous techniques cannot be applied to fixed-point arithmetic programs, because the ranges of individual operations need to be fixed at compile time, i.e., the number of integer bits for each operation has to be determined for all possible inputs up-front. Hence, a continuous abstraction is not possible. Additionally, while floating-point arithmetic supports only a small number of precisions (typically 16, 32, 64 bits), fixed-point arithmetic allows any number of bits to be used for any operation and a more precise encoding is necessary to capture this.

We can capture different numbers of bits precisely by using integer constraints. However at the same time, we also need to guarantee that a target error bound is satisfied, and this error is a non-integer. Hence, we choose to encode fixed-point precision tuning as an Mixed-Integer Linear Programming (MILP) problem [4]. In MILP, some decision variables are constrained to be integers, and other variables can be real-valued. Integers allow us to directly encode the discrete decisions about how many bits to use for operations, and real-valued constraints effectively express error constraints.

The primary constraint of our optimization problem guarantees that the total roundoff error remains inside the target error bound. At the same time, we also need to ensure no overflow. Encoding these two conditions together would result in non-linearity because of the dot operations in a neural network that perform multiplications. Unfortunately, only linear constraints can be efficiently handled by state-of-the-art MILP solvers; nonlinear mixed integer arithmetic is in general NP-hard, with non-convex problems being undecidable. Hence, we soundly over-approximate and linearize constraints as and when necessary to make our solution efficient and feasible.

There are two primary sources of non-linearity in the optimization problem: (1) computation of ranges of all arithmetic operations in order to ensure that there is no overflow, and (2) optimization of fractional bits for all variables and constants to guarantee the error bound.

To avoid the first case, we pre-compute a sound over-approximation of real-valued ranges for all operations efficiently using interval arithmetic (Section 3.1). From these we compute the integer bits needed to represent the real-valued ranges. However, we still need to ensure that the finite-precision ranges (real-valued ranges + roundoff errors) do not overflow. We encode this as a linear range constraint in our optimization problem (Section 3.2.4), thus ensuring no overflow even in the presence of errors.

For the latter, instead of optimizing for all intermediate variables and constants individually, we treat the dot product as a single operation and encode only the assignment of the fractional bits of the dot product results. However, to generate the implementation, we still need to assign precisions to all intermediate variables and constants such that the roundoff errors in intermediate computations do not affect the result. For this, we utilize existing techniques [6, 13] (discussed in more detail in Section 3.3) to determine the number of fractional bits of the intermediate variables in a dot product while guaranteeing an overall error bound.

Provided that the cost function is also linear, we can then encode the precision tuning problem for neural networks with ReLU and linear activation functions with purely linear constraints that can be solved efficiently by state-of-the-art MILP solvers. Other activation functions like sigmoid need to be linearized first before encoding them into our MILP based quantization method. As an additional performance optimization, we consider uniform bit lengths for the operations within each individual layer. That is, within one layer, all dot products will be assigned the same bit length (and similarly for bias and activations), but the bit lengths will vary from layer to layer. This choice is to limit the number of the constraints, but it is not a fundamental limitation of our approach.

Overall, our technique consists of three steps:

Our tool Aster performs all these steps fully automatically. We will explain each of these steps in detail next. Though Aster is primarily designed for continuous feed-forward neural network controllers, our proposed approach can be extended to other types of neural networks that have sparse matrices and activation functions that can be piece-wise linearized [54].

Running Example. We will illustrate our approach using the following small (artificial) neural network in Figure 1 as our running example: The input to the network is a vector \(\overline{x_0}\) consisting of 2 elements, whose ranges are provided by the user: \(x_0^1 = [-10, 10]\) and \(x_0^2 = [-5, 5]\). The neural network produces a single output \(x_2\). To generate an implementation of the network, the user needs to provide the network architecture (i.e., the weight matrices and bias vectors, as well as activation functions for each layer) as input to Aster. Additionally, the user specifies the precision of the input and if the input is represented exactly, which incurs no roundoff error in the input. The precision of inputs is typically known from the specification of sensors or similar in the context of embedded systems, but our approach also supports inputs with initial roundoff errors.

Fig. 1.

View Figure

Assume, for the running example, the inputs are exactly represented with 10 bits. The input specification implies that 5 bits are required to represent the maximum absolute value of the ranges, which is 10, leaving 5 bits for the fractional part of the inputs, with no initial roundoff error. The goal is to generate a quantized mixed-precision implementation of the NN such that its cost is minimized and the output error is bounded by a user-specified target error \(\epsilon _{target} = 0.1\).

3.2 Step 2: Optimizing Fractional Bits

Next, we reduce the problem of computing the fractional bits of dot product results, the addition of bias, and the activation functions to a mixed integer linear programming problem. We first provide an overview of the relevant variables that we will use to formulate our MILP problem in Table 1. The variables specific to a layer i are referenced along with the subscript i, and op can denote the dot product, the bias addition, or the activation function.

Table 1.

type	variables	definitions
user inputs	\(k_i\)	number of neurons at layer i
	\(Q_{W_i}\)	word length of maximum weight at layer i
	\(Q_{x_0}\)	maximum word length of inputs
	lo	lower bound on fractional bits \(\pi\)
	hi	upper bound on fractional bits \(\pi\)
pre-computed inputs	\(R^{op}_i\)	real-valued ranges of operation op at layer i
	\(I^{op}_i\)	number of integer bits of the finite-precision range of operation op at layer i
	\(A^{op}_i\)	maximum representable ranges of operation op without overflow at layer i
decision variables	\(\epsilon ^{prop}_i\)	propagation error at layer i
	\(\epsilon ^{new}_i\)	new roundoff error at layer i
	\(\gamma ^{op}_i\)	cost of an operation op at layer i
	\(\pi ^{op}_i\)	number of fractional bits of op at layer i
other variables	\(b_i^{op}\)	indicator variable for operation op at layer i
	\(opt^{op}_i\)	optimal variable for \(\pi ^{op}\) at layer i

• op: Dot product/bias addition/activation function.

View Table

Table 1. Variable notations

• op: Dot product/bias addition/activation function.

User inputs: The variable \(k_i\), that denotes the number of neurons at layer i, is deduced by Aster from the neural network given by the user. The remaining values are directly provided by the user. The variables \(Q_{W_i}\) and \(Q_{x_0}\) denote the word length of the maximum weight and the largest input from the input vector \(\overline{x}\). The input word length is used as is, however, the word lengths of maximum weights are required in the beginning only to make the optimization problem linear. We later deduce these word lengths automatically considering the result’s precision of the dot product, which we will explain in Section 3.3. The user additionally specifies a range for the fractional bits \(\left[lo, hi\right]\) (one for all operations and layers) that will be considered during optimization.

Pre-computed inputs: The values of these variables are pre-computed by the first step. \(R^{op}_i\) and \(I^{op}_i\) denote the real-valued range and the number of integer bits required for the finite-precision range of each operation op at layer i. \(A^{op}_i\) denotes the maximum representable range of each operation op without overflow at layer i.

Decision variables: The decision variables are those for which the MILP problem will be solved. The variables \(\epsilon ^{prop}_i\) and \(\epsilon ^{new}_i\) are the error variables representing the propagation error and newly introduced error at layer i. These variables constitute the error constraint. \(\gamma ^{op}_i\) represents the cost of each operation op at layer i where op ranges over dot, bias and activation function, and is used to formulate the objective function. The variables \(\pi ^{op}_i\) denote the fractional bits of all operations at layer i.

Other variables: The other variables are used internally by Aster to linearize non-linear constraints (without over-approximations). The indicator variables are used to select the optimized fractional bit length and the optimal variables are used to store the corresponding errors. We explain these variables in more detail in Section 3.2.5.

3.2.1 Problem Formulation.

Our objective is to minimize a cost (function \(\gamma\)) subject to error and range constraints denoted by \(C_1\) and \(C_2\), where the constraints ensure that the user-provided error bound is respected and that no overflow occurs, respectively. Figure 2 presents the overall formulation of the MILP constraints for optimizing fractional bits. We provide a detailed explanation of each constraint; we emphasize the references to individual constraints in Figure 2 using underlines.

Fig. 2.

View Figure

3.2.2 Cost Function.

In each layer, there are three operations: first the dot product is computed, then the bias vector is added, and finally, an activation function is applied. We compute the costs of these operations individually and add them up to compute the total cost of each layer. The total cost over all layers is then computed by adding up the costs of n layers (see \(\gamma\) in Figure 2).

This cost function closely follows previous work [33] and computes the total number of bits needed to implement the neural network. Our approach supports other cost functions, e.g., capturing performance or energy, as long as they can be expressed as linear expressions. Since the latter are highly hardware dependent, we implement the cost function used in previous work.

The cost of the dot operation at layer i, denoted by \(\gamma _i^{dot}\), depends on the number of inputs \(k_{i-1}\), the number of neurons \(k_i\), and the weight matrix \(W_{i}\) of the layer. The dot product at layer i is defined by expanding it into multiplications and additions of the weights and inputs and adding them for all \(k_i\) neurons of the layer: \({dot_i} = \sum _{j=1}^{k_{i}} (\sum _{l=1}^{k_{i-1}} W_{i,jl}\times \overline{x}_{i,l})\).

Recall that our approach takes as input the maximum word length of the weights \(Q_{W_i}\) for each layer. With the \(Q_{W_i}\) and the fractional bits of the dot product \(\pi ^{dot}_i\) (decision variable), we can over-approximate the total cost of the dot operation at layer i by multiplying it with \(Q_{W_i}\), the number of neurons of the previous layer (\(k_{i-1}\)) and the current layer (\(k_i\)) (\(\gamma _i^{dot}\) in Figure 2). For the first layer, instead of \(k_{i-1}\) we consider the number of input variables. Note that we only use an estimate of the maximum word length of W, as the actual word lengths of W are not known beforehand, and defining them as variables makes the cost non-linear. We later assign precisions for W from the solutions of the MILP problem, which we explain in Section 3.3.

As the bias vector is only added with the result of the dot product, the cost of this addition denoted as \(\gamma _i^{bias}\) is the maximum of the fractional bits of the result of the addition and fractional bits of the dot operation (\(\gamma _i^{bias}\) in Figure 2). Technically, the ‘max’ function is also non-linear. We show how to linearize this function in Section 3.2.5.

Finally, the cost of the activation is a function of \(\pi _i^{bias}\). We assume ReLU and linear activation functions for the networks. As these two activation functions are linear functions of \(\pi _i^{bias}\), the cost is the same as \(\pi _i^{bias}\) (\(\gamma _i^{\alpha }\) in Figure 2).

Let us assume for our running example in Figure 1 that \(Q_{W} = 8\) for both layers. The network has 2 input variables, 2 neurons in the first layer and 1 neuron in the output layer. With our formulation, the cost of the whole network is defined as: \(\begin{equation*} \underbrace{(2 * 2 * 8 * \pi ^{dot}_1)+ \text{max}(\pi ^{dot}_1, \pi ^{bias}_1) + \pi ^{bias}_1}_{\text{layer}1} + \underbrace{(2 * 1 * 8 * \pi ^{dot}_2) + \text{max}(\pi ^{dot}_2, \pi ^{bias}_2) + \pi ^{bias}_2}_{\text{output}} \end{equation*}\)

3.2.3 Error Constraint.

The error constraint \(C_1\) in Figure 2 states that the overall roundoff error of the full network \(\epsilon _n\) needs to be bounded by the user specified error \(\epsilon _{target}\). For this constraint, we need to express the roundoff error as a function of the precisions of individual operations.

To compute the roundoff error at each layer i, we need to track the propagated error \(\epsilon _i^{prop}\) from the previous layer and compute the new roundoff error \(\epsilon _i^{new}\) committed by the operations at layer i. The total roundoff error \(\epsilon _n\) is then defined as the sum of all the errors at all layers.

The initial error \(\epsilon _0\) is considered as the propagated error at layer 1. This error \(\epsilon _0\) is bounded by the function of fractional bits of inputs (\(\epsilon _{1}^{prop}\) in Figure 2) that is determined by the word length of the inputs \(Q_{x_0}\) and the integer bits needed to represent the range of the inputs \(R_{x_0}\) computed using the intBits function: \(\pi _0 = Q_{x_0} - {\texttt {intBits}}(R_{x_0})\).

The propagation error at layer > 1 depends on the errors from previous layers as well as the absolute magnitudes of the weights (a weight that is bigger than one will magnify an existing error). Thus, the propagation error \(\epsilon _i^{prop}\) for layer > 1 is defined as the error at the previous layer \(\epsilon _{i-1}\) multiplied by the sums of the absolute weights of each neurons. In practice, we sum the absolute weights of each neuron (absolute sum over the number of inputs \(k_{i-1}\)) and take the maximum result as a factor to amplify the propagation error (maximum over the number of neurons \(k_i\)) (see \(\epsilon _i^{prop}\) in Figure 2).

This is a sound over-approximation of the total propagated error as we assume the maximum magnification of errors for all neurons in the previous layer. For this constraint, Aster computes the maximum values of the weights over the neurons for each layer beforehand and uses them as constants in the optimization problem, thus preserving linearity.

The new roundoff error at layer i is defined as the sum of the errors for the activation function, for the dot computation and the bias addition. The ReLU activation function (\(\alpha (x) = \max (0, x)\)) and the Linear activation function (\(\alpha (x) = x\)) do not affect the error. Thus the new error is bounded by \(2^{-\pi _{i}^{dot}}\) and \(2^{-\pi _i^{bias}}\), where \(\pi _{i}^{dot} = Q_i^{dot} - I_i^{dot}\) and \(\pi _{i}^{bias} = Q_i^{bias} - I_i^{bias}\). For the dot and bias operation, the error is computed considering the fractional bits of these operations.

This error constraint is non-linear (\(2^{-\pi _{i}^{op}}\)). We linearize the error constraint for the optimization exactly by considering the user-provided range \([lo, hi]\) for \(\pi _i^{op}\). We explain the linearization process in Section 3.2.5. We also assume a roundoff error on the result of the dot product only, and do not account for the roundoff errors of the individual operations of the dot product in our MILP constraints. We do this to avoid nonlinearity in the constraints, and rely on the fact that the dot product can be computed with faithful rounding up to the chosen format using the technique from de Dinechin et al. [13]. We explain the details of computing the intermediate formats in Section 3.3.

Let us derive the error constraint for our running example in Figure 1. We considered no input error here: \(\epsilon _0 = 0.0\). We first compute the total error in layer 1. This error is then considered as the propagated error for the output layer by magnifying it with the maximum absolute sum of weights: \(\max \left((0.1+0.2),(0.2 + 0.15)\right) = 0.35\). Thus the overall error constraint is (assuming \(\epsilon _{target} = 0.1\)): \((2^{-\pi _{1}^{dot}} + 2^{-\pi _1^{bias}}) * 0.35 + (2^{-\pi _{2}^{dot}} + 2^{-\pi _2^{bias}}) \le 0.1\).

3.2.4 Range Constraint.

As we have only considered the real-valued ranges to compute the integer bits of the variables and constants in the first phase, we need to ensure that the finite ranges after each operation in each layer do not overflow. Our hypothesis is that the roundoff errors are small enough to keep the integer bits unaffected. Accordingly, in our range constraint \(C_2\) in Figure 2, we want to ensure that the number of integer bits required to represent the finite-precision range after each operation \(I_i^{op}\) is enough to store the integer bits of the finite-precision result (real-valued ranges together with the roundoff errors).

However, directly implementing this constraint would result in applying the function intBits on error decision variables, which is non-linear. Hence, we reduce the problem of computing the integer bits to ensuring the ranges of the finite-precision result remains inside the maximum representable range (\(A_i^ {op}\)) with the integer bits of the real-ranges after every operation at layer i.

We pre-compute the integer bits required for the real-valued result, and use them to generate the maximum representable range after each operation. For the dot operation, we have the propagation errors \(\epsilon _i^{prop}\) from the previous layer and the operation itself introduces the new roundoff error \(\epsilon _i^{dot}\). Hence, after the dot operation, the finite range includes these two errors. For bias addition, however, we have the error from the dot operation as propagated error along with the new roundoff error \(\epsilon _i^{bias}\) introduced by the addition. The total error after the bias addition is essentially the total error \(\epsilon _i\) of the layer i as the activation error is zero for ReLU and Linear functions.

Let us consider our running example in Figure 1. As there is no initial error, for the first layer we simply add the roundoff error committed by the dot operation. For the bias, we add roundoff error for both the dot and the bias operation. The range constraints for the first layer is thus: \([-4.0, 4.0] \ge [-3.0, 3.0] + \epsilon _1^{dot}\) and \([-8.0, 8.0] \ge [-3.0, 6.0] + \epsilon _1^{dot} + \epsilon _1^{bias}\). Similarly, we generate range constraints for all the layers. Given the user inputs, our prototype tool Aster automatically encodes the objective cost function, the error and range constraints for the optimization problem.

3.2.5 Linearization of Constraints.

However, we still have nonlinearity in the objective cost function which includes the nonlinear function max. The error constraints are also nonlinear in terms of the number of unknown fractional bits (\(\pi\)).

Hence, we linearize these constraints exactly without introducing any approximation. Figure 3 presents our linearization constraints for the optimization. Linearization of the cost objective is straight-forward. The max function of two variables \(\pi _i^{dot}\) and \(\pi _i^{bias}\) (as shown in Figure 2) in the form \(\gamma _{i}^{bias} = \text{max}(\pi _i^{dot}, \pi _i^{bias})\) can always be divided into a set of two linear constraints such that \(\gamma _{i}^{bias}\) is always bigger than or equal to all values of both \(\pi _i^{dot}\) and \(\pi _i^{bias}\). We utilize this argument to linearize our max function as \(C_{0(1-2)}\) in Figure 3.

Fig. 3.

View Figure

Linearizing the error constraint is however more tricky. To linearize the error constraint, we use indicator constraints [7]. The nonlinear error constraint is defined in Figure 2 as: \(\epsilon _i^{new} \le 2^{-\pi _i^{dot}} + 2^{-\pi _i^{bias}}\) where \(\pi _i^{dot}\) and \(\pi _i^{bias}\) are unknown integers in an user-provided integer range \([lo, hi]\). Recall that, we want to compute the values of \(\pi _i^{dot}\) and \(\pi _i^{bias}\) that minimize our objective function.

As we already know the integer ranges of the variables, the possible values for \(\pi _i^{dot}\) and \(\pi _i^{bias}\) become finite ( m) as they are exactly the same as the number of integers in the given range. With these m values, we define a set of m discrete reals that represent the set of possible values of \(2^{-\pi _i^{dot}}\) and \(2^{-\pi _i^{bias}}\): \(T = [2^{-hi}, 2^{-(hi-1)}, \ldots , 2^{-lo}]\).

We introduce m binary indicator variables \(b^{dot}\) and \(b^{bias}\) for each valuation of \(2^{-\pi _i^{dot}}\) and \(2^{-\pi _i^{bias}}\) within the specified range. Intuitively, the indicator variables select only one specific value from the list T. We formulate two indicator constraints \(C_{1(1-2)}\) for each layer i as presented in Figure 3.

The constraints \(C_{1(1-2)}\) select the values of \(\pi _i^{dot}\) and \(\pi _i^{bias}\) that are optimal (\(opt^{dot}_i\) and \(opt^{bias}_i\)) respectively. They also state that if \(b_j^{dot}\) (or \(b_j^{bias}\)) is true, we select the value \(2^{-(hi-j-1)}\) from the list T for \(opt^{dot}_i\) (or \(opt^{bias}_i\)). Obviously, we want only one of the \(b_j^{dot}\)(or \(b_j^{bias}\)) to be true. Hence, we add another two constraints \(C_{1(3-4)}\) in Figure 3 to enforce that only one of these binary indicator variables is true. With these new indicator variables and constraints, finally, we linearize the original nonlinear error constraint \(\epsilon _i^{new}\) in Figure 2 as \(C_{1(5)}\) in Figure 3. Our tool Aster encodes these linearization constraints for all layers fully automatically.

Let us assume for our running example in Figure 1 that the range is provided as \([lo, hi] = [4, 8]\) which makes the possible values for \(\pi ^{dot}_i\) and \(\pi ^{bias}_i\): \(T = [2^{-8}, 2^{-7}, 2^{-6}, 2^{-5}, 2^{-4}]\). Next, we define 5 binary indicator variables for \(b^{dot}\) and \(b^{bias}\) each. The indicator constraints for dot product in the 1st layer are as follows: \(opt^{dot}_1 = \sum _{j = 4}^{8} 2^{-j} \times b^{dot}_j\) and \(\sum _{j = 4}^{8} b_j^{dot} = 1\). Similarly we have indicator constraints for bias. If the solver picks \(b^{dot}_7 =1\) and \(b^{bias}_8 =1\), the corresponding new error is then bounded by \(2^{-7} + 2^{-8}\), and the optimized fractional bit lengths of dot and bias are 7 and 8, respectively.

4 IMPLEMENTATION

Our tool Aster takes as input the network architecture and specification written in a small domain-specific language. The input corresponding to our running example from Figure 1 is shown in Figure 4. The ‘require’ clause specifies the input ranges, and the ‘ensuring’ clause specifies the overall error bound. In addition, Aster takes the maximum fractional bit length of the input vector, the word length of maximum weights, and the range of the fractional bits for the optimization as command-line user inputs. As explained in Section 3.2, we need these values to make the optimization problem linear.

Fig. 4.

View Figure

Note that we have implemented Aster to handle feed-forward NNs with relu and linear activations fully automatically; however, Aster can straightforwardly be extended to handle neural networks with sparse matrix multiplication, and piece-wise linear activation functions with bounded domains.

Aster generates fully quantized code including (more) accurate precisions for the weights written in C using the ap_fixed library. This code can directly be compiled for an FPGA with the state-of-the-art HLS compiler by Xilinx [52]. Aster’s generated code either expresses the matrix operations as for-loops, or alternatively it can fully unroll these loops.

Integer Bit Computation. We observed that a straight-forward interval analysis computation of the real-valued ranges quickly becomes expensive with increasing complexity of the network. We thus leverage the structure of the neural network for a more efficient, though over-approximate range analysis. Specifically, we abstract the input variables to each layer by a single range that soundly covers all individual variables, and do so similarly for the weights and biases at each layer. Doing so, we can compute the output range of a layer by computing a single dot product with one addition. This over-approximation makes the interval analysis scalable even for large networks with thousands of parameters, while we have not observed it to significantly effect the optimization results.

Choice of MILP Solver. For our fixed-point precision optimization, the values we encounter can be very small, e.g., for fixed 32 bit precision the roundoff errors are on the order of 1e-9. State-of-the-art MILP solvers use finite precision internally as well, and it is thus crucial to choose a solver that is precise enough to be able to distinguish its own internal roundoff errors from the values in our constraints. We integrate Aster with the SCIP optimization Suite [20] with the underlying SoPlex solver to solve our mixed-integer linear problem. SCIP internally uses extended-precision which goes beyond the limits of floating-point arithmetic, and thus allows us to deal with values as small as 1e-15. We found that other widely-used industrial solvers (such as CPLEX [29] and Gurobi [24]) have tolerances that are bounded by 1e-6. These are unfortunately too coarse for our fixed-point precision optimization. Note that SCIP’s precision is not unlimited either. If the roundoff error goes beyond the tolerance limits (that is 1e-15), Aster cannot optimize and will report an error (i.e., it will not return an incorrect result).

5 EVALUATION

In this section we evaluate Aster based on the following research questions:

Benchmarks. We have collected a set of 18 neural network controllers, consisting of 15 models from the competition at the Applied Verification for Continuous and Hybrid Systems (ARCH) workshop from the years 2019 [40] and 2020 [34], where the networks were provided by academia as well as industry. Additionally, we included 3 controller models from the VNN-LIB standard benchmark set [22], which is widely used for the verification of neural networks. Verification of these controllers is becoming increasingly important due to their usage in safety and operational critical systems. We took all benchmarks for which we could extract all values of weights and biases from the repository. Table 2 provides details of the network architectures. We present here a brief description of the benchmarks; a more detailed discussion can be found in [14, 22, 34].

The InvPendulum, SinglePend, DoublePend, Acrobot (originally from [32]) are the neural network controllers for pendulums in different settings. The Unicycle and MountainCar, originally taken from [15] and [32] respectively, model cars to control two separate control objectives. We also have four different versions of ACC, an Adaptive Cruise Control system of a car, with 3, 5, 7, and 10 hidden layers. The MPC benchmark models a quadrotor (introduced in [32]) and the ControllerTora benchmark controls translational oscillations [16]. The Airplane and VCAS (originally from [35]) model a dynamical system and collision avoidance system ACAS X of an aircraft. The AC benchmarks are derived from a baseline case study involving a Drone controller [22] that regulates a single RPM value across all motors to facilitate takeoff and maintain a fixed hovering position.

We have extracted the safe input ranges, the weight matrices and bias vectors of these controllers from the competition’s repository [1, 14], as well as the VNN-LIB repository [22]. The provided MATLAB files, hierarchical data format (HDF) files, and open format for ML models (ONNX) files were converted into the input format of Aster. The networks do not come with specified target error bounds, so we choose two target absolute error bounds uniformly for all networks, 1e-3 and 1e-5, that we believe to be in a reasonable range for embedded controllers, which are typically implemented in lower precision (16 bits or 32 bits).

Experimental Setup. All experiments are done on an Intel Core i5 Debian system with 3.3 GHz and 16 GB RAM. Aster uses SCIP Optimization Suite 7.0.3 as the external MILP solver. We set 1e-15 as the zero and feasibility tolerance limits for the SCIP as these tolerances are precise and efficient enough for our purpose. We further assume 32 as the initial guess for the word length of the maximum weights to start the optimization. We choose to compare Aster with the state-of-the-art precision tuner Daisy [11], as this is the only available sound mixed-precision tuner for fixed-point programs. For the comparison we use Daisy’s version downloaded on 2^nd March, 2021 (there were no major commits since). We set a 5-hour time budget for each optimization run in our experiments. We believe that a total of 5-hour is a reasonable time for an analysis to generate a sound implementation once that can directly be synthesized on FPGAs. For synthesis of FPGA designs, we use Xilinx’s Vivado HLS tool [52] (version v2020.1) downloaded on 27^th May, 2020, and set a timeout of 5 hours for our experiments.

6 RELATED WORK

Since there is a substantial body of work on efficient neural network inference, we focus here on a representative selection of the most immediately relevant techniques, and point to existing surveys [43, 53] for an overview of other complementary approaches (e.g., model and weight pruning). Verification techniques of other aspects of neural networks (e.g., adversarial examples) are, for instance, covered in the survey by Huang et al. [28]. A more recent work [41] introduces a MILP formulation for verifying robustness of quantized neural networks. However, this work is orthogonal to ours as it focuses on post-quantization verification and does not consider the impact of roundoff errors. Another work [54] focuses on computing output ranges of feed-forward neural networks with non-linear activations (e.g., sigmoid and exponential linear units) in safety-critical systems, where it encodes the non-linear activations with linear constraints. It could be an interesting future direction for Aster to extend it for non-linear activations by combining their technique.

Quantization is a frequently applied technique to reduce the memory footprint of neural networks, but most techniques have been applied on neural network classifiers outside of safety-critical applications and unlike Aster do not provide any accuracy guarantees. We present a brief overview of these approaches in this paper, but note that they are fundamentally not comparable with our sound mixed-precision tuning for regression problems. These methods typically focus on selecting a particular uniform (custom) floating-point [36, 45] or fixed-point precision [21, 23, 39] and showing empirically that it performs well on a particular data/benchmark set. Shiftry [37] automatically chooses a mixed fixed-point precision by iteratively reducing the precision of variables in recurrent neural networks from 16 to 8 bits to run on memory restricted hardware. Another approach to mixed-precision is to dynamically adapt to outliers, i.e., particularly large values, in inputs or weights by providing specialized hardware architectures [44, 47, 48]. Alternative number representations have also been considered, such as stochastic computing [17], posits [9] and floating-points with different, tunable exponent and mantissa widths [19, 50].

Recent work [26] on quantization-aware training for embedded controllers offers new tools for automatic parameter compression or quantization by retraining the network and using a multiplier-less approach for dense layers. This work does not analyze the rounding error propagation, but it is complementary to our approach, and combining these two is a potential future work. There exist different techniques to build custom operators for sums of products by constants with an accuracy requirement. The first work of Dinechin et al. [13] presents an approach that permits to use only a necessary amount of bits for intermediate products, based on a worst-case error-analysis, and we use this approach in this paper. Another technique is a multiplierless approach: multiplications by constants are replaced by bit-shifts and additions, and intermediate terms are shared when several constants are multiplied [38]. A new ILP-based solution for sound truncated multipliers was recently proposed [46], and it takes the result’s format as its input. Aster could be extended to use this approach, instead of [13] for the Step 3, but the scalability of [46] is a major concern for word-lengths beyond 16 bits. For the best fine-tuned FPGA-based dot products the standard HLS tools have to be replaced by custom hardware code generators, such as FloPoCo tool [12]. But then, a whole new problem of a custom NN accelerator arises.

The tool POP finds mixed-precision assignments for floating-point arithmetic code by phrasing the problem as an ILP problem [2]. POP uses a dynamic analysis to infer variable ranges and thus does not guarantee complete soundness. A similar technique has been applied to tuning the floating-point precision of neural networks [30]. FPTuner [10] provides sound mixed-precision tuning for floating-point arithmetic by solving a nonlinear interval-valued optimization problem. Unfortunately, approaches for floating-point arithmetic are generally not applicable to fixed-point arithmetic due to its non-dynamic range.

The tool Daisy [11, 33] includes a sound mixed-precision tuning procedure for both floating-point and fixed-point arithmetic, based on an iterative search using delta-debugging. Other iterative approaches have been explored for mixed-precision tuning for fixed-point arithmetic, for example Min+1 or Max-1 [8], and a combination of Bayesian optimization and Min+1 [25]. Iterative approaches have been proposed to overcome difficulties in phrasing and more importantly solving sound global optimization problems. Alternatives are relaxing the integer to real-valued optimization [18], but work for Digital signal processors (DSPs) only. Generally, these techniques (with the exception of Daisy) treat errors as uncorrelated and additive or evaluate them dynamically with simulation, and thus do not guarantee soundness. We have shown that for the restricted class of neural networks of limited size and with appropriate approximations, a sound optimization is indeed feasible.

Fixed-point arithmetic was also recently formalized as a theory in the context of Satisfiability Modulo Theory (SMT) solvers [5]. However, to compute or bound roundoff errors, such a theory must be combined with the real-valued theory (representing the ideal program), which is in general inefficient.

7 CONCLUSION

This paper presents a novel sound quantization approach that assigns fixed-point mixed precision to neural networks that solve regression tasks, guaranteeing a user-provided roundoff error bound. We reduce the problem of precision optimization to an MILP problem and show how to effectively introduce sound over-approximations to be able to solve it efficiently using a state-of-the-art solver. Our experiments show that our method is fast, even for large networks with thousands of parameters. This efficiency indicates that constraint optimization is not a bottleneck here. Our proposed technique is viable for soundly quantizing neural networks that appear as verified embedded controllers. It is able to handle more and larger benchmarks than existing fixed-point tuners, and mostly generates more efficient code for custom hardware such as FPGAs. While this paper primarily focuses on neural network controllers which are typically regression models, there is potential for extending Aster’s applicability to quantizing classifiers with more work on bridging the gap between numerical and classification accuracy.