Year: 2019

Abstract:
Deep Neural Networks (DNNs) are increasingly deployed in highly energyconstrained environments such as autonomous drones and wearable devices while at the same time must operate in real-time. Therefore, reducing the energy consumption has become a major design consideration in DNN training. This paper proposes the first end-to-end DNN training framework that provides quantitative energy consumption guarantees via weighted sparse projection and input masking. The key idea is to formulate the DNN training as an optimization problem in which the energy budget imposes a previously unconsidered optimization constraint. We integrate the quantitative DNN energy estimation into the DNN training process to assist the constrained optimization. We prove that an approximate algorithm can be used to efficiently solve the optimization problem. Compared to the best prior energy-saving methods, our framework trains DNNs that provide higher accuracies under same or lower energy budgets.

heading: 1 INTRODUCTION
text: Deep Neural Networks (DNNs) have become the fundamental building blocks of many emerging application domains such as computer vision (Krizhevsky et al., 2012; Simonyan & Zisserman, 2014), speech recognition (Hinton et al., 2012), and natural language processing (Goldberg, 2016). Many of these applications have to operate in highly energy-constrained environments. For instance, autonomous drones have to continuously perform computer vision tasks (e.g., object detection) without a constant power supply. Designing DNNs that can meet severe energy budgets has increasingly become a major design objective.
The state-of-the-art model compression algorithms adopt indirect techniques to restrict the energy consumption, such as pruning (or sparsification) (He et al., 2018; Han et al., 2015a; Liu et al., 2015; Zhou et al., 2016; Li et al., 2016; Wen et al., 2016) and quantization (Gong et al., 2014; Wu et al., 2016; Han et al., 2015a; Courbariaux et al., 2015; Rastegari et al., 2016). These techniques are agonistic to energy consumption; rather they are designed to reduce the amount of computations and the amount of model parameters in a DNN, which do not truly reflect the energy consumption of a DNN. As a result, these indirect approaches only indirectly reduce the total energy consumption. Recently, Energy-Aware Pruning (EAP) (Yang et al., 2017) proposes a more direct manner to reduce the energy consumption of DNN inferences by guiding weight pruning using DNN energy estimation, which achieves higher energy savings compared to the indirect techniques.
However, a fundamental limitation of all existing methods is that they do not provide quantitative energy guarantees, i.e., ensuring that the energy consumption is below a user-specified energy budget. In this paper, we aspire to answer the following key question: how to design DNN models that satisfy a given energy budget while maximizing the accuracy? This work provides a solution to this question through an end-to-end training framework. By end-to-end, we refer to an approach that directly meets the energy budget without relying heuristics such as selectively restoring pruned weights and layer by layer fine-tuning (Han et al., 2015b; Yang et al., 2017). These heuristics are effective in practice but also have many hyper-parameters that must be carefully tuned.
Our learning algorithm directly trains a DNN model that meets a given energy budget while maximiz-
ing model accuracy without incremental hyper-parameter tuning. The key idea is to formulate the DNN training process as an optimization problem in which the energy budget imposes a previously unconsidered optimization constraint. We integrate the quantitative DNN energy estimation into the DNN training process to assist the constrained optimization. In this way, a DNN model, once is trained, by design meets the energy budget while maximizing the accuracy.
Without losing generality, we model the DNN energy consumption after the popular systolic array hardware architecture (Kung, 1982) that is increasingly adopted in today’s DNN hardware chips such as Google’s Tensor Processing Unit (TPU) (Jouppi et al., 2017), NVidia’s Tensor Cores, and ARM’s ML Processor. The systolic array architecture embodies key design principles of DNN hardware that is already available in today’s consumer devices. We specifically focus on pruning, i.e., controlling the DNN sparsity, as the main energy reduction technique. Overall, the energy model models the DNN inference energy as a function of the sparsity of the layer parameters and the layer input.
Given the DNN energy estimation, we formulate DNN training as an optimization problem that minimizes the accuracy loss under the constraint of a certain energy budget. The key difference between our optimization formulation and the formulation in a conventional DNN training is two-fold. First, our optimization problem considers the energy constraint, which is not present in conventional training. Second, layer inputs are non-trainable parameters in conventional DNN training since they depend on the initial network input. We introduce a new concept, called input mask, that enables the input sparsity to be controlled by a trainable parameter, and thus increases the energy reduction opportunities. This lets us further reduce energy in scenarios with known input data pattern.
We propose an iterative algorithm to solve the above optimization problem. A key step in optimization is the projection operation onto the energy constraint, i.e., finding a model which is closest to the given (dense) model and satisfies the energy constraint. We prove that this projection can be casted into a 0/1 knapsack problem and show that it can be solved very efficiently. Evaluation results show that our proposed training framework can achieve higher accuracy under the same or lower energy compared to the state-of-the-art energy-saving methods.
In summary, we make the following contributions in this paper:
• To the best of our knowledge, this is the first end-to-end DNN training framework that provides quantitative energy guarantees;
• We propose a quantitative model to estimate the energy consumption of DNN inference on TPU-like hardware. The model can be extended to model other forms of DNN hardware;
• We formulate a new optimization problem for energy-constrained DNN training and present a general optimization algorithm that solves the problem.

heading: 2 RELATED WORK
text: Energy-Agnostic Optimizations Most existing DNN optimizations indirectly optimize DNN energy through reducing the model complexity. They are agonistic to the energy consumption, and therefore cannot provide any quantitative energy guarantees.
Pruning, otherwise known as sparsification, is perhaps the most widely used technique to reduce DNN model complexity by reducing computation as well as hardware memory access. It is based on the intuition that DNN model parameters that have low-magnitude have little impact on the final prediction, and thus can be zeroed-out. The classic magnitude-based pruning (Han et al., 2015b) removes weights whose magnitudes are lower than a threshold. Subsequent work guides pruning using special structures (Liu et al., 2015; Zhou et al., 2016; Li et al., 2016; Wen et al., 2016; He et al., 2017), such as removing an entire channel, to better retain accuracy after pruning.
Quantization reduces the number of bits used to encode model parameters, and thus reducing computation energy and data access energy (Gong et al., 2014; Wu et al., 2016; Han et al., 2015a). The extreme case of quantization is using 1-bit to represent model parameters (Courbariaux et al., 2015; Rastegari et al., 2016). Such binary quantization methods are usually trained from scratch instead of quantizing a pre-trained DNN.
Energy-Aware Optimizations Recently, energy-aware pruning (EAP) (Yang et al., 2017) proposes to use a quantitative energy model to guide model pruning. Different from pure magnitude-based pruning methods, EAP selectively prunes the DNN layer that contributes the most to the total energy consumption. It then applies a sequence of fine-tuning techniques to retain model accuracy. The pruning step and fine-tuning step are alternated until the accuracy loss exceeds a given threshold.
Although EAP a promising first-step toward energy-aware optimizations, its key limitation is that it does not provide quantitative energy guarantees because it does not explicitly consider energy budget as a constraint. Our work integrates the energy budget as an optimization constraint in model training.
Latency-Guaranteed Compression Lately, model compression research has started providing guarantees in execution (inference) latency, which theoretically could be extended to providing energy guarantees as well. However, these methods are primarily search-based through either reinforcement learning (He et al., 2018) or greedy-search (Yang et al., 2018). They search the sparsity setting for every single layer to meet the given budget. Thus, they may require a large number of trials to achieve a good performance, and may not ensure that the resulting model accuracy is maximized.

heading: 3 MODELING DNN INFERENCE ENERGY CONSUMPTION
text: This section introduces the model of estimating energy consumption of a single DNN inference. We consider the widely-used feed-forward DNNs. Note that our proposed methodology can be easily extended to other network architectures as well. In this section, we first provide an overview of our energy modeling methodology (Section 3.1). We then present the detailed per-layer energy modeling (Section 3.2 and Section 3.3), which allow us to then derive the overall DNN energy consumption (Section 3.4). Our energy modeling results are validated against the industry-strength DNN hardware simulator ScaleSim (Samajdar et al., 2018).
DNN model sparsity (via pruning) is well recognized to significantly affect the execution efficiency and thus affect the energy consumption of a DNN model (He et al., 2018; Yang et al., 2017; Han et al., 2015a; Liu et al., 2015; Zhou et al., 2016). We thus use pruning as the mechanism to reduce energy consumption1. Note, however, that model sparsity is not the end goal of our paper; rather we focus on reducing the energy consumption directly. Many dedicated DNN hardware chips (a.k.a., Neural Processing Units, NPUs) (Jouppi et al., 2017; Chen et al., 2016; Han et al., 2016; Parashar et al., 2017) have been developed to directly benefit from model sparsity, and are already widely available in today’s consumer devices such as Apple iPhoneX, Huawei Mate 10, and Microsoft HoloLens. Our paper focuses on this general class of popular, widely-used DNN chips.

heading: 3.1 ENERGY MODELING OVERVIEW
text: A DNN typically consists of a sequence of convolution (CONV) layers and fully connected (FC) layers interleaved with a few other layer types such as Rectified Linear Unit (ReLU) and batch normalization. We focus mainly on modeling the energy consumption of the CONV and FC layers. This is because CONV and FC layers comprise more than 90% of the total execution time during a DNN inference (Chen et al., 2016) and are the major energy consumers (Han et al., 2015a; Yang et al., 2017). Energy consumed by other layer types is insignificant and can be taken away from the energy budget as a constant factor.
A DNN inference’s energy consumption is tied to the underlying hardware that performs the inference. In particular, we assume a systolic-array-based DNN hardware architecture. Systolic array (Kung, 1982) has long been know as an effective approach for matrix multiplication. Many DNN hardware architectures adopt the systolic array, most notably the Google Tensor Processing Unit (TPU) (Jouppi
1Quantization is another useful mechanism to reduce energy consumption. It is orthogonal to the pruning mechanism and they could be combined. This paper specifically focuses on the pruning mechanism.
et al., 2017), Nvidia’s Tensor Cores in their most recent Volta GPUs, and ARM’s ML Processor. Targeting systolic-array-based DNN hardware ensures that our approach has a wide applicability. However, our modeling and training strategies can generally be applied to other DNN architectures.
Figure 1 shows the overall hardware architecture. The systolic array comprises of several compute units that perform the Multiply-and-Accumulate (MAC) operation, which conducts the following computation: a← a+ (b× c), where b and c are the two scalar inputs and a is the scalar intermediate result called “partial sum.” MAC operation is the building block for matrix multiplication. The MAC units are organized in a 2-D fashion. The data is fed from the edges, both horizontally and vertically, which then propagate to the MAC units within the same row and columns.
We decompose the energy cost into two parts: computation energy Ecomp and data access energy Edata. Ecomp denotes the energy consumed by computation units, and Edata denotes the energy consumed when accessing data from the hardware memory. Since we mainly use pruning as the energy reduction technique, we now model how Ecomp and Edata are affected by DNN sparsity.

heading: 3.2 ENERGY CONSUMPTION FOR COMPUTATION
text: CONV layers perform convolution and FC layer perform matrix-vector multiplication. Both operations can be generalized to matrix-matrix multiplication, which involves only the MAC operation (Chetlur et al., 2014; Jouppi et al., 2017). Figure 1 illustrates how a matrix-matrix multiplication is carried out on the systolic array hardware. Given X and W , the systolic array computes XW by passing each row of X to each row in the systolic array and passing each column of W to each column in the systolic array. If the width of the systolic array, denoted by sw, is less than the width of W , the hardware will fold W column-wise in strides of sw. Similarly, if the height of X is greater than the height of the systolic array (sh), X is folded row-size in strides of sh. Figure 1 illustrates a 2× 2 systolic array multiplying two 4× 4 matrices. Both matrices are folded twice in strides of 2. Critically, if either inputs of a MAC operation is zero, we can skip the MAC operation entirely and thus save the computation energy. At a high-level, the total computation energy, Ecomp, can be modeled as eMACNMAC, where eMAC denotes the energy consumption of one MAC operation whereas NMAC denotes the total number of MAC operations that are actually performed. The challenge is to identify NMAC for CONV and FC layers, which we discuss below.
Fully connected layer LetX(v) ∈ R1×c be the input vector andW (v) ∈ Rc×d be the weight matrix of the FC layer v. The FC layer performs matrix-vector multiplication X(v)W (v). The number of MAC operations NMAC is sum(supp(X)supp(W )), where supp(T ) returns a binary tensor indicating the nonzero positions of tensor T . So the computation energy for a fully connected layer v:
E(v)comp = eMACsum(supp(X (v))supp(W (v))) ≤ eMAC‖W (v)‖0, (1)
where the equality is reached when the input is dense.
Convolution layer The CONV layer performs the convolution operation between a 4-D weight (also referred to as kernel or filter) tensor and a 3-D input tensor. Let W (u) ∈ Rd×c×r×r be the weight tensor, where d, c, and r are tensor dimension parameters. Let X(u) ∈ Rc×h×w be the input tensor, where h and w are the input height and width. The convolution operation in the CONV layer u generates a 3-dimensional tensor:
(X(u) ∗W (u))j,y,x = c∑ i=1 r−1∑ r′,r′′=0 X (u) i,y+r′,x+r′′W (u) j,i,r′,r′′ , (2)
where x, y indicate the position of the output tensor, which has height h′ = b(h+ 2p− r)/sc+ 1 and width w′ = b(w + 2p− r)/sc+ 1 (p is the convolution padding and s is the convolution stride). Tensor convolution (2) can be seen as a special matrix-matrix multiplication (Chellapilla et al., 2006; Chetlur et al., 2014). Specifically, we would unfold the tensor X(u) to a matrix X̄(u) ∈ Rh′w′×cr2 , and unfold the tensor W (u) to a matrix W̄ (u) ∈ Rcr2×d. X̄(u) and W̄ (u) are then multiplied together in the systolic array to compute the equivalent convolution result between X(u) and W (u).
Nonzero elements inX(u) andW (u) incur actual MAC operations. Thus,NMAC = sum(supp(X(u))∗ supp(W (u))) ≤ h′w′‖W (u)‖0 (the equality means the input is dense), resulting in the following computation energy of a CONV layer u:
E(u)comp = eMACsum(supp(X (u)) ∗ supp(W (u))) ≤ eMACh′w′‖W (u)‖0. (3)

heading: 3.3 ENERGY CONSUMPTION FOR DATA ACCESS
text: Accessing data happens in every layer. The challenge in modeling the data access energy is that modern hardware is equipped with a multi-level memory hierarchy in order to improve speed and save energy (Hennessy & Patterson, 2011). Specifically, the data is originally stored in a large memory, which is slow and energy-hungry. When the data is needed to perform certain computation, the hardware will load it from the large memory into a smaller memory that is faster and consume less energy. If the data is reused often, it will mostly live in the small memory. Thus, such a multi-level memory hierarchy saves overall energy and improves overall speed by exploiting data reuse.
Without losing generality, we model a common, three-level memory hierarchy composed of a Dynamic Random Access Memory (DRAM), a Cache, and a Register File (RF). The cache is split into two halves: one for holding X (i.e., the feature map in a CONV layer and the feature vector in a FC layer) and the other for holding W (i.e., the convolution kernel in a CONV layer and the weight matrix in an FC layer). This is by far the most common memory hierarchy in DNN hardware such as Google’s TPU (Jouppi et al., 2017; Chen et al., 2016; Zhu et al., 2018; Han et al., 2016). Data is always loaded from DRAM into cache, and then from cache to RFs.
In many today’s DNN hardwares, the activations and weights are compressed in the dense form, and thus only non-zero values will be accessed. This is done in prior works (Chen et al., 2016; Parashar et al., 2017). Therefore, if the value of the data that is being loaded is zero, the hardware can skip the data access and thereby save energy. There is a negligible amount of overhead to “unpack” and “pack” compressed data, which we simply take away from the energy budget as a constant factor. This is also the same modeling assumption used by Energy-Aware Pruning (Yang et al., 2017).
To compute Edata, we must calculate the number of data accesses at each memory level, i.e., NDRAM, Ncache, NRF. Let the unit energy costs of different memory hierarchies be eDRAM, ecache, and eRF, respectively, the total data access energy consumptionEdata will be eDRAMNDRAM +ecacheNcache + eRFNRF. We count the number of data accesses for both the weights and input, then combine them together. The detailed derivation of data access energy is included in the Appendix.

heading: 3.4 THE OVERALL ENERGY ESTIMATION FORMULATION
text: Let U and V be the sets of convolutional layers and fully connected layers in a DNN respectively. The superscript (u) and (v) indicate the energy consumption of layer u ∈ U and v ∈ V , respectively. Then the overall energy consumption of a DNN inference can be modeled by
E(X,W ) := ∑ u∈U (E(u)comp + E (u) data) + ∑ v∈V (E(v)comp + E (v) data), (4)
where X stacks input vectors/tensors at all layers and W stacks weight matrices/tensors at all layers.

heading: 4 ENERGY-CONSTRAINED DNN MODEL
text: Given the energy model presented in Section 3, we propose a new energy-constrained DNN model that bounds the energy consumption of a DNN’s inference. Different from prior work on model pruning in which energy reduction is a byproduct of model sparsity, our goal is to directly bound the energy consumption of a DNN while sparsity is just used as a means to reduce energy.
This section formulates training an energy-constrained DNN as an optimization problem. We first formulate the optimization constraint by introducing a trainable mask variable into the energy modeling to enforce layer input sparsity. We then define a new loss function by introducing the knowledge distillation regularizer that helps improve training convergence and reduce overfitting.
Controlling Input Sparsity Using Input Mask The objective of training an energy-constrained DNN is to minimize the accuracy loss while ensuring that the DNN inference energy is below a given budget, Ebudget. Since the total energy consumption is a function of ‖X(u)‖0 and ‖W (u)‖0, it is natural to think that the trainable parameters are X and W . In reality, however, X depends on the input to the DNN (e.g., an input image to an object recognition DNN), and thus is unknown during training time. Therefore, in conventional DNN training frameworks X is never trainable.
To include the sparsity of X in our training framework, we introduce a trainable binary mask M that is of the same shape of X , and is multiplied with X before X is fed into CONV or FC layers, or equivalently, at the end of the previous layer. For example, if the input to a standard CONV layer is X(u), the input would now be X(u) M (u), where denotes the element-wise multiplication. In practice, we do not really do this multiplication but only read X(u) on the nonzero positions of M (u).
Algorithm 1: Energy-Constrained DNN Training. Input: Energy budget Ebudget, learning rates η1, η2, mask sparsity decay step ∆q. Result: DNN weights W ∗, input mask M∗.
1 Initialize W = Wdense,M = 1, q = ‖M‖0 −∆q; 2 while True do
// Update DNN weights 3 while W has not converged do 4 W = W − η1∇̂W L̄(M,W ) ; // SGD step 5 W = PΩ(Ebudget)(W ) ; // Energy constraint projection for weights W 6 end 7 If previous_accuracy > current_accuracy, exit loop with previous W and M ;
// Update input mask 8 while M has not converged do 9 M = M − η2∇̂M L̄(M,W ) ; // SGD step
10 Clamp values of M into [0, 1]: assign 1 (or 0) to the values if they exceeds 1 (or negative); 11 M = P‖M‖0≤q(M) ; // L0 constraint projection for input mask M 12 end 13 Round values of M into {0, 1}; 14 Decay the sparsity constraint q = q −∆q; 15 end 16 W ∗ = W,M∗ = M .
With the trainable mask M , we can ensure that ‖X(u) M (u)‖0 ≤ ‖M (u)‖0, and thereby bound the sparsity of the input at training time. In this way, the optimization constraint during training becomes E(M,W ) ≤ Ebudget, where E(M,W ) denotes the total DNN inference energy consumption, which is a function of X and W (as shown in Equation (4)), and thus a function of M and W .
Knowledge Distillation as a Regularizer Directly optimizing over the constraint would likely lead to a local optimum because the energy model is highly non-convex. Recent works (Mishra & Marr, 2017; Tschannen et al., 2017; Zhuang et al., 2018) notice that knowledge distillation is helpful in training compact DNN models. To improve the training performance, we apply the knowledge distillation loss (Ba & Caruana, 2014) as a regularization to the conventional loss function. Intuitively, the regularization uses a pre-trained dense model to guide the training of a sparse model. Specifically, our regularized loss function is:
L̄λ,Wdense(M,W ) := (1− λ)L(M,W ) + λEX [‖φ(X;W )− φ(X;Wdense)‖2/|φ(·;W )|], (5) where Wdense is the original dense model, and L(M,W ) is the original loss, e.g., cross-entropy loss for classification task. φ(X;W ) is the network’s output (we use the output before the last activation layer as in Ba & Caruana (2014)), |φ(·;W )| is the network output dimensionality and 0 ≤ λ ≤ 1 is a hyper parameter similar to other standard regularizations.
Thus, training an energy-constrained DNN model is formulated as an optimization problem:
min M,W L̄λ,Wdense(M,W ) s.t. E(M,W ) ≤ Ebudget. (6)

heading: 5 OPTIMIZATION
text: This section introduces an algorithm to solve the optimization problem formulated in (6). The overall algorithm is shown in Algorithm 1. Specifically, the algorithm includes three key parts:
• Initialization by training a dense model. That is, Wdense := arg min
W L(M,W ) (7)
• Fix M and optimize W via approximately solving (using Wdense initialization): min W L̄(M,W ) s.t. E(M,W ) ≤ Ebudget (8)
• Fix W and optimize M by approximately solving : min M L̄(M,W ) s.t. ‖M‖0 ≤ q,M ∈ [0,1] (9)
After the initialization step (Line 1 in Algorithm 1), the training algorithm iteratively alternates between the second (Line 3-6 in Algorithm 1) and the third step (Line 8-13 in Algorithm 1) while gradually reducing the sparsity constraint q (Line 14 in Algorithm 1) until the training accuracy converges. Note that Equation (7) is the classic DNN training process, and solving Equation (9) involves only the well-known L0 norm projection P‖M‖0≤q(Q) := arg min‖M‖0≤q ‖M −Q‖
2. We thus focus on how Equation (8) is solved.
Optimizing Weight Matrix W To solve (8), one can use either projected gradient descent or projected stochastic gradient descent. The key difficulty in optimization lies on the projection step
PΩ(Ebudget)(Z) := arg min W∈Ω(Ebudget)
‖W − Z‖2 (10)
where Z could be W − η∇W L̄(W,M) or replacing ∇W L̄(W,M) by a stochastic gradient ∇̂W L̄(W,M). To solve the projection step, let us take a closer look at the constraint Equation (4). We rearrange the energy constraint Ω(Ebudget) into the following form with respect to W :{
W ∣∣ ∑ u∈U∪V α (u) 1 min(k, ‖W (u)‖0) + α(u)2 max(0, ‖W (u)‖0 − k) + α(u)3 ‖W (u)‖0 + α(u)4 ≤ Ebudget
} ,
(11)
whereW stacks all the variable {W (u)}u∈U∪V , and α(u)1 , α (u) 2 , α (u) 3 , α (u) 4 and k are properly defined
nonnegative constants. Note that α(u)1 ≤ α (u) 2 and k is a positive integer. Theorem 1 casts the energyconstrained projection problem to a 0/1 knapsack problem. The proof is included in the Appendix. Theorem 1. The projection problem in (10) is equivalent to the following 0/1 knapsack problem:
max ξ is binary
〈Z Z, ξ〉, s.t. 〈A, ξ〉 ≤ Ebudget − ∑
u∈U∪V α
(u) 4 , (12)
where Z stacks all the variables {Z(u)}u∈U∪V , A and ξ are of the same shape as Z, and the j-th element of A(u) for any u ∈ U ∪ V is defined by
A (u) j =
{ α
(u) 1 + α (u) 3 , if Z (u) j is among the top k elements of Z (u) in term of magnitude; α
(u) 2 + α (u) 3 , otherwise.
(13)
The optimal solution of (10) is Z ξ∗, where ξ∗ is the optimal solution to the knapsack problem (12).
The knapsack problem is NP hard. But it is possible to find approximate solution efficiently. There exists an approximate algorithm (Chan, 2018) that can find an −accurate solution in O(n log(1/ ) + −2.4) computational complexity. However, due to some special structure in our problem, there exists an algorithm that can find an an −accurate solution much faster. In the Appendix, we show that an (1 + )-approximate solution of problem (10) can be obtained in Õ(n + 1 2 ) time complexity (Õ omits logarithm), though the implementation of the algorithm is complicated. Here we propose an efficient approximate algorithm based on the “profit density.” The profit density of item j is defined as Z2j /Aj . We sort all items based on the “profit density” and iteratively select a group of largest items until the constraint boundary is reached. The detailed algorithm description is shown in the Appendix (Algorithm 2). This greedy approximation algorithm also admits nice property as shown in Theorem 2. Theorem 2. For the projection problem (10), the approximate solution W ′′ ∈ Ω(Ebudget) to the greedy approximation algorithm admits
‖W ′′−Z‖2 ≤ ‖PΩ(Ebudget)(Z)−Z‖ 2+Top‖W ′′‖0+1((Z Z) A)·min((max(A)−gcd(A)), R(W ′′)) (14) where max(A) is the maximal element of A, which is a nonnegative matrix defined in (13); Topk(·) returns the k-th largest element of ·; denotes the element-wise division. gcd(·) is the largest positive rational number that divides every argument, e.g., gcd(0, 1/3, 2/3) = 1/3. In (14), gcd(A) denotes the greatest common divisor of all elements in A2, and R(W ′′) denotes the remaining budget
R(W ′′) = ( Ebudget −
∑ u∈U∪V α (u) 4 − 〈A, supp(W ′′)〉
) .
2Here we assume A only contains rational numbers since gcd is used.
The formal proof is in the Appendix. W ′′ is the optimal projection solution to (10) if either of the following conditions holds:
1. (The remaining budget is 0.) It means that the greedy Algorithm 2 runs out of budget; 2. (The matrix A satisfies max(A) = gcd(A).) It implies that all elements in A have the same
value. In other words, the weights for all items are identical.

heading: 6 EVALUATION
text: The evaluations are performed on ImageNet (Deng et al., 2009), MNIST, and MS-Celeb-1M (Guo et al., 2016) datasets. For the MS-Celeb-1M, we follow the baseline setting reported in the original paper (Guo et al., 2016), which selects 500 people who have the most face images. We randomly sample 20% images as the validation set. We use both classic DNNs, including AlexNet (Krizhevsky et al., 2012) and LeNet-5 (LeCun et al., 1998), as well as recently proposed SqueezeNet (Iandola et al., 2016) and MobileNetV2 (Sandler et al., 2018).
We compare our method mainly with five state-of-art pruning methods: magnitude-based pruning (MP) (Han et al., 2015b;a), structured sparsity learning (SSL) (Wen et al., 2016), structured bayesian pruning (SBP) (Neklyudov et al., 2017), bayesian compression (BC) (Louizos et al., 2017) and energy-aware pruning (EAP) (Yang et al., 2017). Filter pruning methods (Li et al., 2016; He et al., 2017) require a sparsity ratio to be set for each layer, and these sparsity hyper-parameters will determine the energy cost of the DNN. Considering manually setting all these hyper-parameters in energy-constrained compression is not trivial, we directly compare against NetAdapt (Yang et al., 2018) which automatically searches such sparsity ratios and use filter pruning to compress DNN models. We implement an energy-constrained version of NetAdapt, which is originally designed to restrict the inference latency. Note that MobileNetv2 and SqueezeNet have special structures (e.g. residual block) that are not fully supported by NetAdapt. Thus, we show the results of NetAdapt only for AlexNet and LeNet-5.
A cc
ur ac
y D
ro p
(% )
Hyper-parameters In the experiment, we observe that knowledge distillation can improve the performance of MP and SSL, so we apply knowledge distillation to all methods including the baseline for a fair comparison. The results of removing knowledge distillation on MP and SSL are included in the Appendix. We choose the distillation weight λ = 0.5. EAP proposes an alternative way to solve the overfitting issue, so we directly use their results. For all the DNNs, we turn off the dropout layers since we find the knowledge distillation regularization will perform better. In all the experiments, we choose ∆q = 0.1|M | where |M | is the number of all mask elements. For optimizing W , we use a pre-trained dense initialization and update W by SGD with the learning rate η1 = 0.001 and weight decay 10−4. For optimizing input mask parameters M , we use the Adam optimizer (Kingma & Ba, 2014) with η2 = 0.0001 and weight decay 10−5(MNIST)/10−6(MS-Celeb-1M). To
stabilize the training process, we exponentially decay the energy budget to the target budget, and also use this trick in MP training (i.e. decaying the sparsity budget) for fair comparisons.

heading: 6.1 IMAGENET
text: We set an energy budget to be less than the minimal energy consumption among the three baseline methods. We use the same performance metric (i.e. top-5 test accuracy) and hardware parameters, i.e., eMAC, eDRAM, ecache, eRF, sh, sw, kW , kX , as described in the EAP paper (Yang et al., 2017). We initialize all the DNNs by a pre-trained dense model, which is also used to set up the knowledge distillation regularization. The top-5 test accuracies on the dense models are 79.1% (AlexNet), 80.5% (SqueezeNet), and 90.5% (MobileNetV2). We use batch size 128 and train all the methods with 30 epochs. For SSL and NetAdapt, we apply 20 additional epochs to achieve comparable results. We implement the projection operation PΩ(Ebudget) on GPU, and it takes < 0.2s to perform it in our experiments. The detailed wall-clock result is included in the Appendix.
Table 1 shows the top-5 test accuracy drop and energy consumption of various methods compared to the dense model. Our training framework consistently achieves a higher accuracy with a lower
energy consumption under the same energy budget. For instance on AlexNet, under a smaller energy budget (26% < 27%), our method achieves lower accuracy drop over EAP (0.5% vs. 0.8%). The advantage is also evident in SqueezeNet and MobileNetV2 that are already light-weight by design. EAP does not report data on MobileNetV2. We observe that weight sparsity is not a good proxy for energy consumption. Our method achieves lower energy consumption despite having higher density.
Figure 2 comprehensively compares our method with prior work. Solid markers represent DNNs trained from our framework under different energy budgets (x-axis). Empty markers represent DNNs produced from previous techniques. DNNs trained by our method have lower energies with higher accuracies (i.e., solid markers are closer to the bottom-left corner than empty markers). For instance on SqueezeNet, our most energy-consuming DNN still reduces energy by 23% while improves accuracy by 0.2% compared to EAP.

heading: 6.2 MNIST AND MS-CELEB-1M
text: 0
10
20
30 0.0
1.0
1.0
MNIST and MS-Celeb-1M (Guo et al., 2016) represent datasets where inputs have regular patterns that are amenable to input masking. For instance, MS-Celeb-1M is a face image dataset and we use its aligned face images where most of the facial features are located in the center of an image. In such scenarios, training input masks lets us control the sparsity of the layer inputs and thus further reduce energy than merely pruning model parameters as in conventional methods. We do not claim that applying input mask is a general technique; rather, we demonstrate its effectiveness when applicable. We compare our method with MP and SSL using
LeNet-5 and MobileNetV2 for these two datasets, respectively. The pre-trained dense LeNet-5 has 99.3% top-1 test accuracy on MNIST, and the dense MobileNetV2 has 65.6% top-5 test accuracy on MS-Celeb-1M. EAP does not report data on these two datasets. Similar to the evaluation on ImageNet, we set the energy budget to be lower than the energy consumptions of MP and SSL. We use batch size 32 on MNIST and 128 on MS-Celeb-1M, and number of epochs is set the same as the ImageNet experiments. Table 2 compares the energy consumption and accuracy drop. Our method consistently achieves higher accuracy with lower energy under the same or even smaller energy budget. We visualize the sparsity of the learned input masks in Figure 3.

heading: 7 CONCLUSION
text: This paper demonstrates that it is possible to train DNNs with quantitative energy guarantees in an end-to-end fashion. The enabler is an energy model that relates the DNN inference energy to the DNN parameters. Leveraging the energy model, we augment the conventional DNN training with an energy-constrained optimization process, which minimizes the accuracy loss under the constraint of a given energy budget. Using an efficient algorithm, our training framework generates DNNs with higher accuracies under the same or lower energy budgets compared to prior art.

heading: APPENDICES

heading: DETAIL OF ENERGY CONSUMPTION FOR DATA ACCESS

heading: FULLY CONNECTED LAYER
text: To multiply X(v) ∈ Rc and W (v) ∈ Rc×d, each nonzero element of W (v) is used once but loaded three times, once each from DRAM, cache and RF, respectively. Thus, the number of DRAM, cache, and RF accesses for weight matrix W (v) is:
NweightsDRAM = N weights cache = N weights RF = ‖W (v)‖0. (15)
Input X(v) is fed into the systolic array dd/swe times, where sw denotes the the systolic array width. Thus, the number of cache accesses for X(v) is:
N inputcache = dd/swe‖X (v)‖0. (16)
Let kX be the cache size for inputX(v). If kX is less than ‖X(v)‖0, there are ‖X(v)‖0−kX elements that must be reloaded from DRAM every time. The rest kX elements need to load from only DRAM once as they will always reside in low-level memories. Thus, there are dd/swe(‖X(v)‖0− kX) + kX DRAM accesses for X(v). In addition, the output vector of the FC layer (result of X(v)W (v)) needs to be written back to DRAM, which further incurs d DRAM accesses. Thus, The total number of DRAM accesses to retrieve X(v) is:
N inputDRAM = dd/swemax(0, ‖X (v)‖0 − kX) + min(kX , ‖X(v)‖0) + d. (17)
Each input element is loaded from RF once for each MAC operation, and there are two RF accesses incurred by accumulation for each MAC operation (one read and one write). Thus, the total number of RF accesses related to X(v) is:
N inputRF = d‖X (v)‖0 + 2‖W (v)‖0. (18)
In summary, the data access energy of a fully connected layer v is expressed as follows, in which each component follows the derivations in Equation (15) through Equation (18):
E (v) data = eDRAM(N input DRAM +N weights DRAM ) + ecache(N input cache +N weights cache ) + eRF(N input RF +N weights RF ). (19)

heading: CONVOLUTION LAYER
text: Similar to a FC layer, the data access energy of a CONV layer u is modeled as:
E (u) data = eDRAM(N input DRAM +N weights DRAM ) + ecache(N input cache +N weights cache ) + eRF(N input RF +N weights RF ). (20)
The notations are the same as in FC layer. We now show how the different components are modeled.
To convolve W (u) ∈ Rd×c×r×r with X(u) ∈ Rc×h×w, each nonzero element in the weight tensor W (u) is fed into the systolic array dh′w′/she times, where sh denotes the height of the systolic array and h′ and w′ are dimension parameters of X(u). Thus,
Nweightscache = dh ′w′/she‖W (u)‖0. (21)
Similar to the FC layer, the number of RF accesses for W (u) during all the MAC operations is:
NweightsRF = h ′w′‖W (u)‖0. (22)
Let kW be the cache size for the weight matrix W (u). If ‖W (u)‖0 > kW , there are kW nonzero elements of W (u) that would be accessed from DRAM only once as they would reside in the cache, and the rest ‖W (u)‖0 − kW elements would be accessed from DRAM by dh′w′/she times. Thus,
NweightsDRAM = dh ′w′/shemax(0, ‖W (u)‖0 − kW ) + min(kW , ‖W (u)‖0). (23)
Let kX be the cache size for input X(u). If every nonzero element in X(u) is loaded from DRAM to cache only once, N inputDRAM would simply be ‖X(u)‖0. In practice, however, the cache size kX is much
smaller than ‖X(u)‖0. Therefore, some portion of X(u) would need to be re-loaded. To calculate the amount of re-loaded DRAM access, we observe that in real hardware X(u) is loaded from DRAM to the cache at a row-granularity.
When the input X(u) is dense, there are at least cw elements loaded at once. In this way, the cache would first load bkX/(cw)c rows from DRAM, and after the convolutions related to these rows have finished, the cache would load the next bkX/(cw)c rows in X(u) for further processing. The rows loaded in the above two rounds have overlaps due to the natural of the convolution operation. The number of overlaps Roverlap is dh/(bkX/cwc− r+ s)e− 1, and each overlap has cw(r− s) elements. Thus, Roverlap × cw(r − s) elements would need to be reloaded from DRAM. Finally, storing the outputs of the convolution incurs an additional dh′w′ DRAM writes. Summing the different parts together, the upper bound (X(u) is dense) number of DRAM accesses for X(u) is:
N inputDRAM = ‖X (u)‖0 + (dh/(bkX/cwc − r + s)e − 1)cw(r − s) + dh′w′. (24)
When the input X(u) is not dense, we can still count the exact number of elements in the overlaps Noverlap of the consecutive loading rounds, so we have:
N inputDRAM = ‖X (u)‖0 +Noverlap + dh′w′. (25)
Every nonzero element in the unfolded input X̄(u) would be fed into the systolic array dd/swe times (for grouped convolution, this number is divided by the number of groups). Each MAC operation introduces 2 RF accesses. Thus,
N inputcache = dd/swe‖X̄‖0, N input RF = d‖X̄ (u)‖0 + 2h′w′‖W (u)‖0. (26)

heading: PROOF TO THEOREM 1
text: Proof. First, it is easy to see that (10) is equivalent to the following problem
max ξ is binary
〈Z Z, ξ〉, s.t. ξ ∈ Ω(Ebudget). (27)
Note that if the optimal solution to problem (27) is ξ̄, the solution to problem (10) can be obtained by Z ξ̄; given the solution to (10), the solution to (27) can be obtained similarly. Therefore, we only need to prove that (27) is equivalent to (12). Meeting the following two conditions guarantees that (27) and (12) are equivalent since they have identical objective functions:
1. Any optimal solution of problem (27) is in the constraint set of problem (12);
2. Any optimal solution of problem (12) is in the constraint set of problem (27).
Let us prove the first condition. Let ξ̂ be the optimal solution to (27). Then for any u ∈ U ∪ V , the elements of Z(u) selected by ξ̂(u) are the largest (in terms of magnitude) ‖ξ̂(u)‖0 elements of Z(u); otherwise there would exist at least one element that can be replaced by another element with a larger magnitude, which would increase the objective value in (27). Since ξ̂ ∈ Ω(Ebudget), according to the definition of A in (13), ξ̂ satisfies the constraint of (12).
Let us now prove the second condition. The definition of A in (13) show that there could at most be two different A(u) values for each element u, and the largest k elements in Z(u) always have the smaller value, i.e., α(u)1 + α (u) 3 . Let ξ̄ be the optimal solution to the knapsack problem (12). For any u ∈ U ∪ V , the elements selected by ξ̄(u) are also the largest elements in Z(u) in terms of magnitude; otherwise there would exist an element Z(u)j that has a larger magnitude but corresponds to a smaller A (u) j ((13) shows that A (u) i ≥ A (u) j when |Z (u) i | ≤ |Z (u) j |). This would contradict the fact that ξ̄ is optimal. In addition, ξ̄ meets the constraint in problem (12). Therefore, ξ̄ ∈ Ω(Ebudget). It completes the proof.
AN (1 + )-APPROXIMATE SOLUTION FOR PROBLEM (10)
Theorem 3. For the projection problem (10), there exists an efficient approximation algorithm that has a computational complexity of O (( n+ (|U |+|V |) 3
2
) log nmax(A)min(A+) ) and generates a solution

heading: W ′ ∈ Ω(Ebudget) that admits
text: ‖W ′ − Z‖2 ≤ ∥∥∥∥PΩ( Ebudget
1+O( )
)(Z)− Z ∥∥∥∥2 , (28)
where min(A+) is the minimum of the positive elements in A.
|U | and |V | denote the number of CONV and FC layers, respectively. They are very small numbers that can be treated as constants here. Thus, the computational complexity for our problem is reduced to Õ(n+ 1 2 ), where Õ omits the logarithm term. In the following, we will prove this theorem by construction.

heading: PROBLEM FORMULATION
text: Definition 1. Inverted knapsack problem. Given n objects I := {(vi, wi)}ni=1 each with weight wi > 0, and value vi ≥ 0, define hI(x) to be the smallest weight budget to have the total value x:
hI(x) := min ξ∈{0,1}n n∑ i=1 wiξi (29)
s. t. n∑ i=1 viξi ≥ x
We are more interested in the case that the weights of n objects are in m clusters, i.e. there are only m distinct weights,
|{wi}ni=1| = m.
In our case, m is proportional to the number of layers in DNN, and n is the number of all the learnable weights in W , so m n. Definition 2. Inverse of step function. The inverse of the step function f is defined as the maximal x having the function value y:
f−1(y) := max f(x)≤y x (30)
Observation The inverse of the step function h−1I (y) is just the maximal value we can get given the weight budget, i.e. the original knapsack problem:
h−1I (y) = max ξ∈{0,1}n n∑ i=1 viξi, s. t. n∑ i=1 wiξi ≤ y. (31)
Observation Given a step function with l breakpoints, its inverse can be generated with O(l) time complexity, and vice versa.
Thus, given the step function of hI in (29) which has l breakpoints, we can get h−1I (i.e. the original knapsack problem) within O(l) time complexity. Definition 3. w-uniform. Step function f is w-uniform if the ranges of f is from −∞, 0, w, 2w, ..., lw.
Observation If all the objects in I have the same weight w, i.e. m = 1, then the function hI(x) is nondecreasing and w-uniform. Moreover, its breakpoints are:
(0, 0), (v1, w), (v1 + v2, 2w), ..., ( n∑ i=1 vi, nw ) ,
if the objects’ indices follows the decreasing order in terms of the values, i.e. v1 ≥ v2 ≥ ... ≥ vn. Thus we can get all possible function values of hI(x):
hI(x) = kw, ∀x ∈ ( k−1∑ i=1 vi, k∑ i=1 vi ] .
Definition 4. (min, +)-convolution. For functions f, g, the (min, +)-convolution is:
(f ⊕ g)(x) = min x′
(f(x′) + g(x− x′)).
Observation If object sets I1 ∩ I2 = ∅, then
fI1∪I2 = fI1 ⊕ fI2 .
Observation The inverse of (min, +)-convolution between w-uniform function f and w-uniform function g is the (max, +)-convolution between f−1 and g−1:
(f ⊕ g)−1(y) = max y′∈{0,1w,...,lw} (f−1(y′) + g−1(y − y′)). (32)
Lemma 4. For any f and g nonnegative step functions, given an arbitrary number b, we always have
min{f ⊕ g, b} = min{min{f, b} ⊕min{g, b}, b} (33)
Proof. Given any x, let z ∈ Arg minx′ f(x′) + g(x − x′) and z̄ ∈ Arg minx′ min(f(x′), b) + min(g(x− x′), b), so we have (f ⊕ g)(x) = f(z) + g(x− z) and (min{f, b} ⊕min{g, b})(x) = min(f(z̄), b) + min(g(x− z̄), b). Consider the following cases:
1. (f ⊕ g)(x) ≥ b. In this case, we claim that (min{f, b} ⊕min{g, b})(x) ≥ b. We prove it by contradiction. Suppose (min{f, b} ⊕min{g, b})(x) < b which implies min(f(z̄), b) + min(g(x − z̄), b) < b. Because both f and g are nonnegative, we have f(z̄) < b and g(x − z̄) < b which imply min(f(z̄), b) + min(g(x − z̄), b) = f(z̄) + g(x − z̄) < b, However, this contradicts (f ⊕ g)(x) ≥ b. Therefore, we have min((f ⊕ g)(x), b) = min((min{f, b} ⊕min{g, b})(x), b) = b.
2. (f ⊕ g)(x) < b. In this case, we have f(z) < b and g(x − z) < b, so min(f(z̄), b) + min(g(x− z̄), b) ≤ min(f(z), b)+min(g(x−z), b) = f(z)+g(x−z) = (f ⊕g)(x) < b. Since both f and g are nonnegative, we have f(z̄) < b and g(x − z̄) < b which imply min(f(z̄), b) + min(g(x − z̄), b) = f(z̄) + g(x − z̄) ≥ (f ⊕ g)(x). Therefore, we have min(f(z̄), b) + min(g(x − z̄), b) = f(z) + g(x − z) ⇔ (min{f, b} ⊕ min{g, b})(x) = (f ⊕ g)(x).

heading: EFFICIENCY OF (MIN, +)-CONVOLUTION
text: Lemma 5. Let f and g be nondecreasing w-uniform functions with O(l) breakpoints, the (min, +)-convolution f ⊕ g (having O(l) breakpoints) can be generated with O(l2) time complexity.
Proof. Firstly, we compute the inverse representation of f and g, i.e. compute f−1 and g−1 from Equation (30). The inverse representation can be computed in O(l) time (proportional to the number of breakpoints). From Equation (32), we can compute the inverse of f ⊕ g. For each y ∈ {0, 1w, ..., 2lw}, function (f ⊕ g)−1(y) can be computed in O(l) time by brute force. Thus a total O(l2) is enough to get (f ⊕ g)−1 which has O(l) breakpoints. We can get f ⊕ g via (f ⊕ g)−1 by the inverse definition (30) in O(l) time.
Lemma 6. Let f and g be nondecreasing step functions with l breakpoints in total, min{f ⊕ g, b} can be approximated by a step function φb with O(l + 1 2 ) complexity and 2 b additive error, i.e. min{f ⊕ g, b} ≤ φb ≤ min{f ⊕ g, b}+ 2 b. The resultant function φb has O(1/ ) breakpoints.
Proof. We can construct ( b)-uniform functions f ′b, g ′ b which have d1/ e breakpoints:
f ′b(x) =
⌈ min(b, f(x))
b
⌉ b, g′b(x) = ⌈ min(b, g(x))
b
⌉ b.
This needsO(l) computational complexity. From Lemma 5, we can compute f ′b⊕g′b withO( 1 2 ) time complexity and φb = min{f ′b ⊕ g′b, b} has O(1/ ) breakpoints. Because f ′b and g′b are constructed by ceiling min{f, b} and min{g, b}, we have:
min{f, b} ⊕min{g, b} ≤ f ′b ⊕ g′b ≤ min{f, b} ⊕min{g, b}+ 2 b,
which implies
min{min{f, b} ⊕min{g, b}, b} ≤ min{f ′b ⊕ g′b, b} ≤ min{min{f, b} ⊕min{g, b}, b}+ 2 b.
From Lemma 4, we know that min{min{f, b}⊕min{g, b}, b} = min{f ⊕ g, b}, so it completes the proof.
Lemma 7. Let f1, f2, ..., fm be nondecreasing step functions with l breakpoints in total, min{f1 ⊕ f2 ⊕ ... ⊕ fm, b} can be approximated by a step function ψb with O(l + m/ 2) computational complexity and m b additive error. The resultant function ψb has O(1/ ) breakpoints.
Proof. From Lemma 6, we have shown the case m = 2. For general m > 2, we can construct a binary tree to approximate pairs of functions, e.g., if m = 4, we can firstly approximate ψ(1) ≈ min{f1 ⊕ f2, b}, and ψ(2) ≈ min{f3 ⊕ f4, b}, then approximate ψ(3)b ≈ min{ψ(1) ⊕ ψ(2), b}. By this way, we construct a binary tree which has O(logm) depth and O(m) nodes. In the beginning, we use ceil function to construct m new b-uniform functions:
f ′i,b(x) =
⌈ min(b, fi(x))
b
⌉ b,∀i ∈ {1, 2, ...,m}.
Then we can use the binary tree to “merge” all the m functions in pairs, via O(logm) iterations. Without loss of generality, we assume m is a power of two. We can recursively merge t functions into t/2 functions:
1. Initialize t = m, g′i,b = f ′ i,b,∀i ∈ {1, ..., t}.
2. Reassign g′i,b = min{g′2i−1,b ⊕ g′2i,b, b},∀i ∈ {1, ..., t/2}. According to Lemma 6, the number of break points of min{g′2i−1,b ⊕ g′2i,b, b} is still O(1/ ).
3. t = t/2. If t > 1, go back to Step 2.
4. Return ψb := min{g′1,b, b}.
For this binary tree, functions of the bottom leaf nodes have b additive error, and every (min, +)- convolution f ′ ⊕ g′ will accumulate the additive error from the two functions f ′ and g′. The root node of the binary tree will accumulate the additive errors from all the m leaf nodes, thus the resultant function ψb ≤ min{f1⊕ ...⊕ fm, b}+m b. For the computational complexity, initializing f ′i,b takes O(l), Step 1 takes O(l), Step 2 and 3 take O(m/ 2) (since there are O(m) nodes in the binary tree), and Step 4 takes O(m/ ). Therefore, there is O(l +m/ 2) in total.
Lemma 8. For the inverted knapsack problem defined in Equation (29), if all the n objects can be separated into m groups I1, ..., Im which have m distinct weights, there exists an approximate algorithm with computational complexity O((n+ m 3
2 ) log nmax(w) min(w) ) which can approximate hI by
h̃I : hI(x) ≤ h̃I(x) ≤ (1 +O( ))hI(x), ∀x.
Proof. Firstly, the step function hIi ,∀i ∈ {1, 2, ...,m} can be easily generated within O(n log n) by sorting the objects of each group according to their values (in descending order). From the definition of (min, +)-convolution, we know that hI = hI1 ⊕ ... ⊕ hIm . Let us construct an algorithm to approximate hI :
1. Construct a set B := {2inmax(w) ∈ [min(w), nmax(w)]; i ∈ Z≤0}, where min(w) and max(w) are the minimum and maximum weight of items respectively, and Z≤0 is the nonpositive integer set. We have |B| = O(log nmax(w)min(w) ).
2. For every b ∈ B, construct ψb to approximate min{hI1 ⊕ ...⊕ hIm , b} based on Lemma 7.
3. Construct function h̃−1I :
h̃−1I (y) = { ψ−1b (y), if b/2 < y ≤ b and y > min(B); ψ−1min(B)(y), if y ≤ min(B).
where min(B) is the minimum element in B. The resultant function h̃−1I (or h̃I ) has at most O( 1 log
nmax(w) min(w) ) breakpoints.
4. Compute the original function h̃I from h̃−1I .
According to the above procedure, for any hI(x) ∈ (b/2, b], h̃I(x) approximate hI(x) with additive error O(m b), so we have hI(x) ≤ h̃I(x) ≤ (1 + O(m ))hI(x). The algorithm takes O((n + m/ 2) log nmax(w)min(w) ), if we require the approximation factor to be 1 +O( ), i.e.,
hI(x) ≤ h̃I(x) ≤ (1 +O( ))hI(x), ∀x, we need
O ( (n+m3/ 2) log nmax(w)
min(w) ) time complexity.
Theorem 9. For the knapsack problem defined in Equation (12), if all the n objects have m distinct weights, there exists an approximate algorithm with computational complexity O((n + m3
2 ) log nmax(w) min(w) ) to generate a function h̃ −1 I satisfying:
h−1I
( y
1 +O( )
) ≤ h̃−1I (y) ≤ h −1 I (y), ∀y.
Proof. From Lemma 8, we have h̃I(x) ≤ (1 +O( ))hI(x) which implies
{x | (1 +O( ))hI(x) ≤ y} ⊆ {x | h̃I(x) ≤ y}. So
max hI(x)≤y/(1+O( )) x ≤ max h̃I(x)≤y
x⇔ h−1I (
y
1 +O( )
) ≤ h̃−1I (y).
Similarly, we can get {x | h̃I(x) ≤ y} ⊆ {x | hI(x) ≤ y} from Lemma 8, so we have
max hI(x)≤y x ≥ max h̃I(x)≤y x⇔ h−1I (y) ≥ h̃ −1 I (y).
Let I be the set of objects whose weights are nonzero elements in A and values are the corresponding elements in Z Z, i.e. I+ = {(Z2i , Ai) | ∀i ∈ {1, 2, ..., |A|} and Ai > 0}, ξ̃+ be the solution corresponding to h̃−1I (Ebudget− ∑ u∈U∪V α (u) 4 ). Let ξ̃Ic+ = 1 and ξ̃I+ = ξ̃
+, where Ic+ = {(Z2i , Ai) | ∀i ∈ {1, 2, ..., |A|} and Ai = 0} is the complement of I+. Here we have m ≤ 2|U |+ |V | distinct values in A. According to Theorem 9, we have 〈Z Z, ξ̃〉 ≥ maxξ〈Z Z, ξ〉, s.t. 〈A, ξ〉 ≤ Ebudget− ∑ u∈U∪V α (u) 4
1+O( ) , which implies
〈Z Z, ξ̃〉 ≥ max ξ∈Ω(Ebudget/(1+O( ))) 〈Z Z, ξ〉.
From Theorem 9, we can directly get Theorem 3.
Algorithm 2: Greedy Algorithm to Solve Problem (12). Input: Z,A,Ebudget, {α(u)}u∈U∪V as in (12). Result: Greedy solution ξ̃ for problem (12).
1 Initialize b = 0, ξ = 0. 2 Generate the profit density δ:
δj =
{ (Zj)
2/Aj , if Aj > 0; ∞, if Aj = 0.
3 Sort δ, let I be the indices list of the sorted δ (in descending order). 4 foreach index j ∈ I do 5 b = b+Aj ; 6 If b > Ebudget − ∑ u∈U∪V α (u) 4 , exit loop; 7 ξj = 1; 8 end 9 ξ̃ = ξ.

heading: PROOF TO THEOREM 2
text: Proof. From Theorem 1, we know the original projection problem (10) is equivalent to the knapsack problem (12). So proving the inequality (14) is equivalent to proving
〈Z Z, ξ̃〉 ≥ 〈Z Z, ξ∗〉 − Top‖ξ̃‖0+1((Z Z) A) ·R(ξ̃) (34) and
〈Z Z, ξ̃〉 ≥ 〈Z Z, ξ∗〉 − Top‖ξ̃‖0+1((Z Z) A) · (max(A)− gcd(A)), (35)
where ξ̃ is the greedy solution of knapsack problem corresponding to W ′′, and ξ∗ is the exact solution of knapsack problem corresponding to PΩ(Ebudget)(Z), i.e.,
W ′′ = Z ξ̃, PΩ(Ebudget)(Z) = Z ξ ∗.
Firstly, let us prove the inequality (35). If we relax the values of ξ to be in the range [0, 1] instead of {0, 1}, the discrete constraint is removed so that the constraint set becomes
∆ = { ξ | 0 ≤ ξ ≤ 1 and 〈A, ξ〉 ≤ Ebudget −
∑ u∈U∪V α (u) 4
} .
So the 0/1 knapsack problem is relaxed as a linear programming. This relaxed problem is called fractional knapsack problem, and there is a greedy algorithm (Dantzig, 1957) which can exactly solve the fractional knapsack problem. Slightly different from our Algorithm 2, the greedy algorithm for the fractional knapsack can select a fraction of the item, so its remaining budget is always zero. The optimal objective value of the fractional knapsack is
max ξ∈∆ 〈Z Z, ξ〉 = 〈Z Z, ξ̃〉+ Top‖ξ̃‖0+1((Z Z) A) ·R(ξ̃).
Since the constraint set of the fractional knapsack problem is a superset of the constraint of the original knapsack problem, we have 〈Z Z, ξ∗〉 ≤ max0≤ξ≤1〈Z Z, ξ〉, that leads to inequality (34). Secondly, we show that the inequality (35) is also true. Since all the coefficients in A are multiples of gcd(A), we can relax the original 0/1 knapsack problem in this way: for each item, split them to several items whose coefficients in the constraint are gcd(A), and the coefficients in the objective function are split equally. For the j-th item, the coefficient in the constraint is Aj and the coefficient in the objective function is (Z Z)j . It will be split into Aj/ gcd(A) items, and the j-th item is associated with coefficient (Z2j /Aj) · gcd(A) in the objective function. This relaxation gives us a new 0/1 knapsack problem, where all the items have the same coefficient in the constraint, so the optimal solution is just selecting the ones with the largest coefficients in the objective function. We can formulate this problem as a relaxed knapsack problem by replacing the constraint of ξ into ξ ∈ Γ, where
Γ = { ξ | for all j, ξj is a multiple of gcd(A)
Aj , 0 ≤ ξj ≤ 1, and 〈A, ξ〉 ≤ Ebudget − ∑ u∈U∪V α (u) 4
} .
All the elements of the solution are either 0 or 1 except the last picked one which corresponds to Top‖ξ̃‖0+1((Z Z) A). Let the (‖ξ̃‖0 + 1)-th largest element in (Z Z) A be indexed by t. We have 0 ≤ ξ̃t ≤ 1− gcd(A)/At. Therefore, comparing with the original 0/1 knapsack problem, we have
max ξ∈Γ 〈Z Z, ξ〉 ≤ 〈Z Z, ξ̃〉+ (Z Z)t · (1− gcd(A)/At)
= 〈Z Z, ξ̃〉+ Top‖ξ̃‖0+1((Z Z) A) ·At · (1− gcd(A)/At)
= 〈Z Z, ξ̃〉+ Top‖ξ̃‖0+1((Z Z) A) · (At − gcd(A))
≤ 〈Z Z, ξ̃〉+ Top‖ξ̃‖0+1((Z Z) A) · (max(A)− gcd(A))
Since {ξ | ξ is binary} ⊆ Γ, we have 〈Z Z, ξ∗〉 ≤ maxξ∈Γ〈Z Z, ξ〉. So we have the inequality (35).

heading: SUPPLEMENTARY EXPERIMENT RESULTS

heading: RESULTS OF BASELINE WITHOUT KNOWLEDGE DISTILLATION
text: Table 3 shows the energy and accuracy drop results of the baseline methods MP and SSL when the knowledge distillation is removed from their loss function. By using knowledge distillation, the results in Table 1 are much better. Therefore, we use knowledge distillation in all the experiments when it is applicable.

heading: ENERGY-CONSTRAINED PROJECTION EFFICIENCY
text: The projection operation PΩ(Ebudget) in Algorithm 1 can be implemented on GPU. We measured its wall-clock time on a GPU server (CPU: Xeon E3 1231-v3, GPU: GTX 1080 Ti), and the result is shown in Table 4 (the time is averaged over 100 iterations).

Review summaries:
The paper proposes a method for neural network training under a hard energy constraint (i.e. the method guarantees the energy consumption to be upper bounded). Based on a systolic array hardware architecture the authors model the energy consumption of transferring the weights and activations into different levels of memory (DRAM, Cache, register file) during inference. The energy consumption is therefore determined by the number of nonzero elements in the weight and activation tensors. To minimize the network loss under an energy constraint, the authors develop a training framework including a novel greedy algorithm to compute the projection of the weight tensors to the energy constraint.

The paper is dedicated to energy-based compression of deep neural networks. While most works on compression are dedicated to decreasing the number of parameters or decreasing the number of operations to speed-up or reducing of memory footprint, these approaches do not provide any guarantees on energy consumption. In this work the authors derived a loss for training NN with energy constraints and provided an optimization algorithm for it. The authors showed that the proposed method achieves higher accuracy with lower energy consumption given the same energy budget. The experimental results are quite interesting and include even highly optimized network MobileNetV2.