Title: DECOUPLING WEIGHT REGULARIZATION FROM BATCH SIZE FOR MODEL COMPRESSION

Abstract: Batch size selection affects the generalization ability of gradient descent, and small batch size is usually preferred. Conventionally, weight regularization is performed for every mini-batch without considering when would be the right time to regularize weights during optimization steps. For model compression, which regularizes weights to follow compressed forms, compression-aware training also performs weight compression for every mini-batch to compute the impact of compression on the loss function. In this paper, we propose a new hyper-parameter called "Non-Regularization period" or NR period during which weights are not regularized (or compressed). We first investigate the influence of NR period on regularization using weight decay and weight random noise insertion. Throughout various experiments, we show that stronger weight regularization demands longer NR period (regardless of batch size) to best utilize regularization effects. From our empirical evidence, we argue that weight regularization for every mini-batch allows small weight decay/noise only and limited regularization effects such that there is a need to search for right NR period and weight regularization strength to enhance model accuracy. Consequently, NR period becomes especially crucial for model compression where strong weight regularization is necessary to increase compression ratio. Using various models, we show that simple weight regularization to comply with compression formats along with long NR period is enough to achieve high compression ratio and model accuracy.

1. INTRODUCTION:

For Deep Neural Networks (DNNs), a common training method updates weights by gradient descent and explicit weight manipulation through regularization. To overcome some practical issues of gradient descent on non-convex optimization problems, there have been several enhancements such as learning rate scheduling and adaptive update schemes using momentum and update history (Ruder, 2016) . Optimizing batch size is another way to yield efficient gradient descent. Note that large batch size has the advantage of enhancing parallelism of the training system in order to speed up training, critical for DNN research (Dean et al., 2012) . Despite such advantages, small batch size is preferred because it improves generalization associated with flat minima search (Keskar et al., 2016) and other hyper-parameter explorations are more convenient (Masters & Luschi, 2018) . Small batch size also affects weight regularization if weight updates for gradient descent and weight regularization are supposed to happen for every mini-batch. For example, for weight decay conducted for every minibatch, if batch size is modified, then the weight decay factor should also be adjusted accordingly (Loshchilov & Hutter, 2017) .

Weight regularization is a process that adds information to the model as a way to avoid overfitting (Goodfellow et al., 2016; van Laarhoven, 2017) . In this paper, we explore weight compression as a form of weight regularization as it severely restricts the search space of weights (i.e., regularized by compression forms). Moreover, model compression shrinks the effective model size, which is an important regularization principle (Goodfellow et al., 2016) (note that improved model accuracy by model compression is reported (Frankle et al., 2019) ). Weights are regularized in numerous ways by model compression. For example, each weight can be pruned (e.g., Han et al. (2015) ) or quantized (e.g., Guo et al. (2017) ) to yield a sparse model representation or to reduce the number of bits to represent each weight. While weight regularization for model compression can simply be performed once after training, compression-aware training can improve model accuracy by reflecting the impact of model compression on the loss function for every mini-batch update (Courbariaux et al., 2015; Zhu & Gupta, 2017; Zhu et al., 2017; Guo et al., 2017) . Now, our question in this paper is the following: Should weight regularization follow batch size that is considered for gradient descent? If the answer is "No," then we would need to consider a new hyper-parameter search space with the right weight regularization timing for existing weight regularization methods and model compression schemes.

In this paper, we report the following observations: 1) Model compression tends to induce weight decay and random weight noise; 2) For weight decay and weight noise insertion, less frequent weight regularization allows stronger regularization to best utilize regularization effects regardless of batch size; and 3) Similarly, if weight regularization for compression is performed less frequently, training for model compression permits stronger weight regularization with improved model accuracy, and thus, achieves higher compression ratios. We confirm our above observations through various experiments and propose an occasional regularization method to decouple weight regularization from batch size selection that is optimized for gradient descent. We verify that our simple model compression techniques (without modifying the underlying training procedures) based on occasional weight regularization (with longer nonregularization period) can achieve higher compression ratio and higher model accuracy compared to previous techniques that demand substantial modifications to the training process.

2. NOISE MODEL ON WEIGHT COMPRESSION:

We first study the relationship between model compression ratio and weight regularization strength using quantization and singular-value decomposition (SVD) as model compression techniques. We assume a popular quantization method based on binary codes for which a weight vector w is approximated to be q i=1 α i b i for q-bit quantization, where α is a scaling factor and b(= {-1, +1} n ) is a binary vector, and n is the vector size. The quantization error ||wi α i b i || 2 is minimized by a method proposed by Xu et al. (2018) to compute α and b. For SVD, a weight matrix W ∈ R m×n is approximated to be W ∈ R m×n by minimizing ||W -W || subject to rank(W ) ≤ R, where R is the target rank.

For our experiments, we use a synthetic (2048 × 2048) weight matrix where each element is randomly generated from the Gaussian distribution N (µ = 0, σ 2 = 1). Then, we are interested in the amount of change of each weight after quantization and SVD. Assuming that weight noise through compression is expressed as in the form of w = w(1 + ), Figure 1 shows the distribution of with various quantization bits or target ranks. From distributions skewed to be negative, it is clear that weights tend to decay more with higher compression ratio, along with a wider range of random noise. Reasonable explanations of Figure 1 would include: 1) weights generated from the Gaussian distribution are uncorrelated such that an approximation step (by compression) using multiple weights would result in noise for each weight, 2) in the case of SVD, elements associated with small eigenvalues are eliminated, 3) averaging effects in quantization reduce the magnitude of large weights. For weight pruning, becomes -1 or 0 (i.e., weight decay for selected weights). Correspondingly, we study on weight decay and weight noise insertion in the next two sections as an effort to gain a part of basic knowledge on improved training for model compression, even though actual model compression would demand much more complicated weight noise models.

Step 1

Step 2

Step 3

∆W ∆L:

Step 4

Step 5 

Loss:



Non-Regularization Period=pNR:

Step 1

Step 2

Step 3

Step 5

Step 4 Weight Updates w/o Regularization Weight Regularization

Figure 2 : Gradient descent and weight regularization when NR period is given as a multiple of batches. Depending on the loss surface and/or strength of regularization, regularization would lead to step 2 (escaping from a local minimum) or step 5 (returning to a local minimum).

3. NON-REGULARIZATION PERIOD:

Since weight regularization cannot precede updates for gradient descent, in order to control the frequency of weight regularization, an available option is to skip a few batches without regularization.

In this paper, we propose a new hyper-parameter, called "Non-Regularization period" or NR period, to enable occasional regularization and to define the interval of two consecutive regularization events as shown in Figure 2 . NR period is expressed as a multiple of batches.

Strong weight regularization facilitates the chance of escaping a local minimum (depicted as step 2 in Figure 2 ) or require longer NR period to return to a local minimum (described as step 5 in Figure 2 ). Let us estimate NR period (pNR) for a returning case. Given a parameter set w (that is assumed to be close enough to a local minimum) and a learning rate γ, the loss function of a model L(w) can be approximated as

EQUATION

using a local quadratic approximation where H is the Hessian of L and w 0 is a set of parameters at a local minimum. After regularization is performed at step t, w can be updated by gradient descent as follows:

EQUATION

Thus, after pNR, we obtain w t+pNR = w 0 + (I -γH(w 0 )) pNR (w t -w 0 ), where I is an identity matrix. Suppose that H is positive semi-definite and all elements of I -γH(w 0 ) are less than 1.0, w t+pNR can converge to w 0 with long pNR which should be longer with larger (w t -w 0 ) (i.e., stronger weight regularization) or smaller γH(w 0 ). In the next sections, we focus on the relationship between pNR and the strength of weight regularization.

4. NR PERIOD STUDY ON WEIGHT DECAY AND WEIGHT NOISE INSERTION:

Weight decay is one of the most well-known regularization techniques (Zhang et al., 2018) and different from L 2 regularization in a sense that weight decay is separated from the loss function calculation (Loshchilov & Hutter, 2017) . Weight decay is performed as

EQUATION

where θ is a constant weight decay factor. Weight noise insertion is another regularization technique aiminig at reaching flat minima (Goodfellow et al., 2016; Hochreiter & Schmidhuber, 1995) . Suppose that random Gaussian noise is added to weights such that w = w + when ∼ N (0, ηI). Then, L(w ) = E[f w+ (x) -y] 2 where x, y, f are input, target, and prediction function, respectively. Using Taylor-series expansion to second-order terms, we obtain f w+ (x) ≈ f w (x) + ∇f (x) + ∇ 2 f (x) /2. Correspondingly, the loss function can also be approximated as

EQUATION

)

where the second term disappears near a local minimum and the third term induces flat minima. Random noise insertion with other distribution models can be explained in a similar fashion (Goodfellow et al., 2016) . We study the impact of NR period on weight decay and weight noise insertion using ResNet-32 on CIFAR-10 model (He et al., 2016) and a long short-term memory (LSTM) model on PTB dataset (Zaremba et al., 2014) . For LSTM model, we use 2 layers with 200 hidden units and the hyperparameter set introduced by Zaremba et al. (2014) . For weight noise model, we plug θ ∼ U(-a, +a) (uniform distribution) into Eq. (3) to simplify the experiments. Note that compared with a conventional weight decay factor selection (i.e., θ in Eq. (3) when pNR = 1), weight decay factor can be approximately 1,000 times larger with pNR ≈ 1, 000 in Figure 3 and Figure 4 . If e w is the mean absolute difference given as (= E[|w -w |]), where w is a weight vector after regularization, then Figure 3 and Figure 4 imply that optimal e w /pNR seems to be constant. Consequently, optimal e w depends on pNR, and hence, we argue that a wide exploration of various NR period (for occasional regularization) and different weight regularization strength is necessary in order to best utilize regularization effects on a DNN model while previous attempts employ pNR = 1 only. In addition, the observation that stronger weight regularization is enabled by longer pNR is our basic training principle for model compression. 1 0 2 0 3 0 4 0 5 0 1 0 0 1 5 4 2 9 0 4 6 4 5 8 0 7 7 4  1 1 6 1  2 3 2 3  2 7 8 7  3 2 5 2  3 7 1 6  4 1 8 1  4 6 4 6  5 1 1 0  5 5 7 5  6 0 3 9  6 5 0 4  6 9 6 9 Non-Regularization Period (a) Perplexity when the weights are compressed by SVD.

5 1 0 2 0 3 0 4 0 5 0 1 0 0 1 5 4 2 9 0 4 6 4 5 8 0 7 7 4 1 1 6 1 2 3 2 3 2 7 8 7 3 2 5 2 3 7 1 6 4 1 8 1 4 6 4 6 5 1 1 0 5 5 7 5 6 0 3 9 6 5 0 4 6 9 6 9

Non-Regularization Period 

5. NR PERIOD FOR MODEL COMPRESSION:

As discussed, weight compression incurs a much more complicated weight regularization model than weight decay or uniform weight noise insertion because 1) as shown in Figure 1 , diversified noise models need to be combined to describe weight regularization after model compression and 2) compression-aware training methods would reduce the strength of weight regularization as training is performed with more epochs and weights converge to a compressed form. Nonetheless, we can conjecture that the best training scheme for model compression may require the condition of pNR = 1 that can be empirically justified.

We apply weight quantization, low-rank approximation (SVD), and pruning to an LSTM model on PTB that we selected for the previous section. We do not modify underlying training principles and use the following simple strategy:

1) Train the model during NR period (as if model compression is not being considered.) 2) Then, perform weight compression in the form of w = h(w).

3) With new full-precision weight w , repeat the above two steps.

h(w) can be a magnitude-based pruning (i.e., h(w)=w if |w| is larger than a certain threshold, or h(w)=0, otherwise), αb for quantization, SVD function, or even as-yet undiscovered functions.

Figure 5 shows model accuracy associated with a number of different sets of pNR and model compression strength (i.e., target rank for low-rank approximation and the number of quantization bits).

Optimal pNR for the best model accuracy is definitely much larger than 1. To explain how Figure 5 is aligned with Section 4, we investigate the relationship between model accuracy and the average of e w /pNR

(e w = E[|w -w |],

where w is a weight vector after weight decay, SVD, or pruning) during entire training. We first determine optimal e w /pNR (= e opt ) for the best model accuracy.

As discussed in Section 4, e opt is assumed to be constant regardless of compression ratio or decay

1 5 2 0 4 0 1 0 0 2 9 0 5 8 0 1 1 6 1

Non-Regularization Period Error(|ew/pNR -eopt|)

PPL(pruning rate=0.8) Error(pruning rate=0.8) PPL(pruning rate=0.9)

Error(pruning rate=0.9) (c) Pruning (eopt=0.035)

Figure 6 : Relationship between model accuracy and error (defined as the difference between e w /pNR and e opt ) using PTB LSTM model when weights are regularized by weight decay, SVD, or pruning.

factor, and obtained by finding hyper-parameter sets associated with maximum model accuracy in Figure 4 and Figure 5 and by taking the average of corresponding e w values. When error is defined to be |e w /pNR -e opt |, Figure 6 shows test perplexity and error of PTB LSTM model with different pNR. Such defined error is affected by pNR as shown in Figure 6 , and indeed, when error approaches to an optimal value we gain improved model accuracy. Unlike weight decay where e w is directly computed by decay factors, for model compression techniques, e w is not directly related to compression-related hyper-parameters (such as ranks and pruning rates). As a result, while Figure 4 shows clear correlation between decay factors and pNR for best model accuracy, Figure 5 suggests that compression ratio and pNR are weakly correlated. Hence, pNR is a hyper-parameter to be determined empirically for model compression. Nonetheless, optimal pNR is definitely much larger than 1, as shown in Figure 6 , and decoupled from batch size selection. That means weight regularization for model compression needs to be conducted much less frequently compared with gradient descent since batch size selection considers generalization ability of gradient descent, not regularization effects.

Note that periodic compression has been introduced in the literature to gradually improve compression ratio or automate hyper-parameter search process. DropPruning repeats dropping weights randomly and retraining the model while some previously dropped weights are unpruned until pruning rate reaches a target number (Jia et al., 2018) . Weights are incrementally quantized to improve model accuracy (Zhou et al., 2017) or the number of quantization bits can be controlled differently for each layer by a loop based on reinforcement learning (Elthakeb et al., 2018) . Structured pruning and fine-tuning process can be iterated to increase pruning rate (Molchanov et al., 2016; Liu et al., 2017) . All of these previous works assume pNR = 1 (i.e., performing compression for every mini batch) while the goal is increasing compression ratio slowly or finding a set of hyper-parameters through iterative fine-tuning stages. Our proposed compression technique can be combined with such periodic compression methods (incremental compression or automatic hyper-parameter selection are also applicable to our proposed method). In the work by He et al. (2018) , soft filter pruning is conducted with pNR = 1 epoch without analysis of why such occasional pruning improves model accuracy.

6. COMPARISON WITH PREVIOUS MODEL COMPRESSION TECHNIQUES:

In this section, we compare some of previous model compression techniques with our compression scheme that introduces pNR and obviates special training algorithm modifications.

6.1. FINE-GRAINED WEIGHT PRUNING:

The initial attempt of pruning weights was to locate redundant weights by computing the Hessian to calculate the sensitivity of weights to the loss function (LeCun et al., 1990) . However, such (Han et al., 2015) has become popular because one can quickly find redundant weights by simply measuring the magnitude of weights. Since then, numerous researchers have realized higher compression ratio largely by introducing Bayesian inference modeling of weights accompanying supplementary hyper-parameters.

For example, dynamic network surgery (DNS) (Guo et al., 2016) permits weight splicing when a separately stored full-precision weight becomes larger than a certain threshold. Optimizing splicing threshold values, however, necessitates extensive search space exploration, and thus, longer training time. Variational dropout method (Molchanov et al., 2017) introduces an explicit Bayesian inference model for a prior distribution of weights, which also induces various hyper-parameters and increased computational complexity.

We perform magnitude-based pruning at every pNR step. As a result, even though weights are pruned and replaced with zero at pNR steps, pruned weights are still updated in full precision during NR period. If the amount of updates of a pruned weight grows large enough between two consecutive regularization steps, then the weight pruned at last pNR step may not be pruned at the next pNR step. Such a feature (i.e., pruning decisions are not fixed) is also utilized for weight splicing in DNS (Guo et al., 2016) . Weight splicing in DNS relies on a hysteresis function (demanding sophisticated finetuning process with associated hyper-parameters) to switch pruning decisions. Pruning decisions through our scheme, on the other hand, are newly determined at every pNR step.

We present experimental results with LeNet-5 and LeNet-300-100 models on MNIST dataset which are also reported by Guo et al. (2016) ; Molchanov et al. (2017) . LeNet-5 consists of 2 convolutional layers and 2 fully connected layers while 3 fully connected layers construct LeNet-300-100. We train both models for 20000 steps using Adam optimizer where batch size is 50. All the layers are pruned at the same time and the pruning rate increases gradually following the equation introduced in Zhu & Gupta (2017):

EQUATION

)

where E is a constant, p f is the target pruning rate, p i is the initial pruning rate, t is the current step, and the pruning starts at training step t i and reaches p f at training step t f . After t f steps, pruning rate is maintained to be p f . For LeNet-5 and LeNet-300-100, t i , p i , E are 8000 (step), 25(%), and 7, respectively. t f is 12000 (step) for LeNet-5 and 13000 (step) for LeNet-300-100. Note that these choices are not highly sensitive to test accuracy as discussed in Zhu & Gupta (2017) . We exclude dropout to improve the accuracy of LeNet-300-100 and LeNet-5 since pruning already works as a regularizer (Han et al., 2015; Wan et al., 2013) . We keep the original learning schedule and the total number of training steps (no additional training time for model compression).

Table 1 presents the comparison on pruning rates (see Appendix for test accuracy). Despite the simplicity, our pruning scheme produces higher pruning rate compared with DNS and similar com-pared with variational dropout technique which involves much higher computational complexity.

For Table 1 , we use pNR=10 for LeNet-5 and pNR=5 for LeNet-300-100.

6.2. LOW-RANK APPROXIMATION:

We apply our proposed occasional regularization algorithm integrated with Tucker decomposition (Tucker, 1966) 

where C i,j,rs,rt is the reduced kernel tensor, R s is the rank for input feature map dimension, R t is the rank for output feature map dimension, and P S and P T are 2D filter matrices to map C i,j,rs,rt to K i,j,s,t . Each component is obtained to minimize the Frobenius norm of ( Ki,j,s,t -K i,j,s,t ). As a result, one convolution layer is divided into three convolution layers, specifically, (1×1) convolution for P S , (d × d) convolution for C i,j,rs,rt , and (1 × 1) convolution for P T (Kim et al., 2016) .

In prior tensor decomposition schemes, model training is performed as a fine-tuning procedure after the model is restructured and fixed (Lebedev et al., 2015; Kim et al., 2016) . On the other hand, our training algorithm is conducted for Tucker decomposition as follows:

Step 1: Perform normal training for pNR (batches) without considering Tucker decomposition

Step 2: Calculate C, P S , and P T using Tucker decomposition to obtain K

Step 3: Replace K with K

Step 4: Go to Step 1 with updated K

After repeating a number of the above steps towards convergence, the entire training process should stop at Step 2, and then the final decomposed structure is extracted for inference. Because the model is not restructured except in the last step, Steps 2 and 3 can be regarded as special steps to encourage wide search space exploration so as to find a compression-friendly local minimum where weight noise by decomposition does not noticeably degrade the loss function.

Using the pre-trained ResNet-32 model with CIFAR-10 dataset (He et al., 2016; Kossaifi et al., 2019) , we compare two training methods for Tucker decomposition: 1) typical training with a decomposed model and 2) pNR-based training, which maintains the original model structure and occasionally injects weight noise through decomposition. Using an SGD optimizer, both training methods follow the same learning schedule: learning rate is 0.1 for the first 100 epochs, 0.01 for the next 50 epochs, and 0.001 for the last 50 epochs. Except for the first layer, which is much smaller than the other layers, all convolution layers are compressed by Tucker decomposition with rank R s and R t selected to be S and T multiplied by a constant number R c (0.3 ≤ R c ≤ 0.7 in this experiment). Then, the compression ratio of a convolution layer is

d 2 ST /(SR s + d 2 R s R t + T R t ) = d 2 ST /(S 2 R c +d 2 R 2 c ST +T 2 R c ), which can be approximated to be 1/R 2 c if S = T and d R c

. pNR is chosen to be 200.

Figure 7 shows test accuracy after Tucker decomposition1 by two different training methods. Note that test accuracy results are evaluated only at Step 3 where the training process can stop to generate a decomposed structure. In Figure 7 , across a wide range of compression ratios (determined by R c ), the proposed scheme yields higher model accuracy compared to typical training. Note that even higher model accuracy than that of the pre-trained model can be achieved by our method if the compression ratio is small enough. In fact, Figure 8 shows that our technique improves training loss and test accuracy throughout the entire training process. Initially, the gap of training loss and test accuracy between pre-regularization and post-regularization is large. Such a gap, however, is quickly reduced through training epochs. Overall, ResNet-32 converges successfully through the 

7. CONCLUSION:

In this paper, we show that weight regularization should be decoupled from batch size to enhance regularization effects. To accomplish such decoupling, we introduce a new hyper-parameter called non-regularization period or NR period during which weights are updated only for gradient computations. A larger amount of weight decay and weight noise insertion need to be supported by longer NR period to maximize the regularization effect. We find that model compression with higher compression ratio is also better trained by longer NR period. We demonstrate that various DNN models can be compressed by pruning, quantization, and low-rank approximation successfully with longer NR period while the underlying training algorithm do not need to be modified.

A APPENDIX:

A.1 SUPPLEMENTARY EXPERIMENTS FOR WEIGHT DECAY AND WEIGHT NOISE 1 2 5 1 0 2 0 5 0 1 0 0 2 0 0 5 0 0 1 0 0 0 2 0 0 0 Non-Regularization Period 8 114.1 114.6 114.3 114.3 114.3 114.2 114.6 114.5 114.4 114.7 114.7 114.6 114.8 114.2 114.4 114.2 114.4 114.3 114.6 114.5 114.1 114.4 114.1 115.1 114.6 113.9 114.1 114.5 114.1 113.7 114.8 114.6 114.2 114.7 114.5 113.9 114.1 114.3 114.3 113.9 114.1 113.9 114.1 114.0 114.0 114.2 113.7 114.7 115.0 114.5 114.9 114.2 114.3 114.3 114.4 114.5 114.2 114.3 114.0 114.7 114.8 114.1 114.1 113.6 114.3 114.1 114.7 114.2 114.3 114.4 114.4 114.2 114.5 114.2 114.0 114.4 114.4 112.8 114.1 114.4 114.2 114.5 114.3 114.7 114.6 114.8 114.9 114.7 114.5 114.3 114.1 114.5 114.5 139.4 114.6 113.2 114.4 114.2 114.2 114.2 115.0 114.2 114.8 114.5 114.5 114.3 114.1 114.5 114.0 499.9 153.9 114.2 113.2 113.7 113.8 114.3 113.8 114.2 113.9 114.2 113.8 114.4 114.7 114.6 114.8 215.2 418.1 170.8 118.7 112.9 113.1 113.7 113.9 113.9 114.2 114.0 114.7 114.3 113.9 114.6 114.0 200.7 209.9 521.9 164.9 118.6 113.5 113.0 112.4 114.0 113.9 114.3 113.8 113.9 114.2 114.9 114.6 189.1 198.9 232.3 544.8 168.7 127.5 119.1 114.8 112.9 113.5 113.8 114.7 114.4 113.8 114.2 114.3 186.8 196.3 205.9 272.3 450.8 173.7 138.3 125.0 113.4 112.6 113.6 114.0 114.0 113.5 114.3 114.4 182.3 190.4 6 114.6 114.9 114.5 114.7 114.7 114.6 114.2 114.9 114.3 114.2 114.2 114.3 114.8 115.0 114.3 114.6 114.4 114.9 114.5 114.5 114.6 114.6 114.6 115.0 114.8 114.4 114.6 114.8 114.6 114.9 114.6 114.5 114.9 114.5 115.0 114.8 114.3 115.1 114.9 114.8 114.6 114.4 114.9 114.5 114.5 114.5 115.1 114.7 114.6 115.0 114.2 114.9 114.5 114.6 114.7 114.3 114.6 114.5 114.7 114.3 114.7 114.6 114.8 114.9 114.8 114.8 114.4 114.7 114.5 114.7 114.7 114.6 115.1 115.0 114.8 114.7 114.6 115.0 114.5 114.0 115.2 114.6 114.8 114.5 114.5 114.5 114.5 114.8 114.8 114.6 114.7 114.3 114.9 114.7 114.2 114.2 114.3 114.5 114.6 114.6 114.5 114.8 114.9 114.7 115.0 114.6 114.9 114.6 114.8 114.9 119.2 113.8 114.1 114.2 114.5 115.0 114.4 114.7 114.6 114.8 114.3 115.1 114.4 114.6 114.6 114.6 540.9 173.2 115.6 113.7 113.8 114.3 114.1 114.6 114.8 114.8 114.7 114.8 114.9 114.8 114.9 114.3 217.8 216.7 455.4 274.0 126.4 116.6 114.5 113.8 113.8 113.7 114.5 114.9 114.8 114.4 114.5 114.5 202.9 212.5 219.8 225.3 538.8 530.4 268.7 174.9 118.9 114.6 113.3 113.9 114.0 114.0 114.0 114.5 192.0 200.2 214.5 220.2 222.7 264.3 376.2 591.7 210.5 134.9 115.2 113.8 113.8 113.6 113.9 114.2 190.8 192.6 207.2 216.7 224.0 221.9 232.3 253.9 693.3 428.5 132.5 117.1 114.5 113.8 113.7 113.9 189.7 190.7 199.3 212.4 We apply pNR-based pruning to an RNN model to verify the effectiveness of pNR. We choose an LSTM model (Zaremba et al., 2014) on the PTB dataset (Marcus et al., 1993) . Following the model structure given in Zaremba et al. (2014) , our model consists of an embedding layer, 2 LSTM layers, and a softmax layer. The number of LSTM units in a layer can be 200, 650, or 1500, depending on the model configurations (referred as small, medium, and large model, respectively). The accuracy is measured by Perplexity Per Word (PPW), denoted simply by perplexity in this paper. We apply gradual pruning with E = 3, t i = 0, p i = 0, t f = 3 rd epoch (for medium) or 5 th epoch (for large) to the pre-trained PTB models. pNR-based pruning for the PTB models is performed using pNR = 100 and the initial learning rate is 2.0 for the medium model (1.0 for pre-training) and 1.0 for the large model (1.0 for pre-training) while the learning policy remains to be the same as in Zaremba et al. (2014) . For all of the pruning rates selected, Table 4 shows that our compression scheme improves perplexity better than the technique in Zhu & Gupta (2017) which is based on Han et al. (2015) . The superiority of pNR-based pruning is partly supported by the observation that non-zero weights successfully avoid to be small through retraining while the conventional pruning still keeps near-zero (unmasked) weights as depicted in Figure 15 . 6 114.3 114.1 113.2 114.3 114.3 113.1 114.2 113.3 112.5 113.7 111.7 112.2 111.9 111.7 111.6 112.5 112.4 113.7 114.1 115.4 116.3 116.5 116.8 117.2 118.9 114.1 113.7 113.8 114.5 114.1 113.1 113.4 113.4 113.6 113.0 112.7 111.7 111.7 111.4 111.0 111.2 111.2 111.5 112.4 112.5 113.8 114.1 115.0 114.2 116.1 116.2 114.0 113.9 113.6 114.4 113.7 114.6 113.9 113.7 113.5 113.3 112.9 112.1 112.0 112.4 111.3 111.1 111.8 111.4 111.8 111.5 113.2 113.9 113.5 114.8 114.5 114.8 114.0 114.5 113.5 113.5 113.5 114.4 113.8 113.9 113.4 113.7 113.5 112.3 112.6 112.3 111.8 111.6 112.0 111.2 111.7 111.6 112.9 113.7 113.2 113.6 114.0 115.0 115.4 114.3 114.6 114.0 114.2 113.9 114.5 113.8 113.8 112.7 113.4 113.2 113.3 112.7 112.6 111.4 111.7 112.0 111.7 112.1 112.9 112.3 113.4 114.1 113.9 113.9 114.4 115.2 113.8 114.3 114.8 114.0 114.1 114.3 114.2 114.6 113.5 113.5 113.7 112.8 112.5 111.6 111.3 111.9 112.0 112.4 112.9 112.5 112.3 112.9 113.4 113.6 115.1 115.7 115.7 114.7 114.1 113.9 114.8 115.2 114.7 113.7 113.8 113.1 113.5 113.6 113.2 112.1 112.7 113.0 112.3 112.5 112.9 113.0 113.4 113.6 113.3 113.4 114.9 115.6 115.3 115.2 114.5 115.0 114.3 114.2 114.3 113.8 114.3 114.0 114.4 114.3 113.7 112.5 113.0 112.5 112.7 113.1 113.5 113.5 113.9 114.4 112.7 114.0 111.2 112.0 112.8 113.6 114.4 115.2 [Quantization] 1 128.6 123.1 120.9 118.5 117.7 114.1 113.2 112.5 111.7 111.6 112.6 112.8 115.0 113.8 114.3 114.5 115.6 116.0 115.8 116.8 116.4 117.7 119.1 118.0 114.3 112.5 112.1 111.9 111.8 112.0 113.1 113.6 113.8 114.3 114.5 115.6 116.0 115.3 115.6 116.2 116.8 117.8 117.5 119.0 117.6 117.6 120.7 111.1 111.1 112.1 112.3 111.2 112.1 113.5 113.4 114.0 114.4 114.7 115.3 115.7 116.6 116.2 116.4 116.5 116.2 117.0 116.9 117.7 117.4 116.5 117.4 117.9 117.6 117.1 116.6 116.2 115.2 116.1 115.5 115.6 114.7 115.0 114.0 114.9 114.7 114.5 114.2 113.7 114.4 114.3 114.3 114.4 114.9 114.8 114.3 114.8 115.2 118.0 116.6 116.2 116.5 116.2 116.5 116.1 116.1 115.2 115.5 114.0 114.0 113.9 114.1 113.7 113.5 113.7 113.9 115.4 115.6 114.1 114.3 114.1 114.2 114.6 114.0 118.3 117.4 116.6 117.3 115.5 116.5 116.2 115.5 115.8 114.3 114.7 114.0 114.2 113.9 112.9 112.6 113.5 113.4 113.3 113.3 114.5 114.2 114.0 112.9 114.2 115.1 117.6 117.3 117.0 116.3 116.3 115.8 115.6 115.4 114.9 114.7 114.2 113.1 113.1 112.7 113.0 113.2 113.1 113.0 113.3 114.2 114.7 114.3 114.2 114.3 114.6 115.5 117.2 118.4 117.2 116.8 115.6 115.8 116.0 115.2 114.7 113.8 114.0 113.4 113.2 113.4 113.3 114.0 114.3 114.3 115.7 116.2 117.1 116.8 116.4 117.4 117.6 Convolution can be performed by matrix multiplication if an input matrix is transformed into a Toeplitz matrix with redundancy and a weight kernel is reshaped into a T × (S × d × d) matrix (i.e., a lowered matrix) Chetlur et al. (2014) . Then, commodity computing systems (such as CPUs and GPUs) can use libraries such as Basic Linear Algebra Subroutines (BLAS) without dedicated hardware resources for convolution Cho & Brand (2017) . Some recently developed DNN accelerators, such as Google's Tensor Processing Unit (TPU) Jouppi et al. (2017) , are also focused on matrix multiplication acceleration (usually with reduced precision).

For BLAS-based CNN inference, reshaping a 4D tensor K and performing SVD is preferred for low-rank approximation rather than relatively inefficient Tucker decomposition followed by a lowering technique. However, a critical problem with SVD (with a lowered matrix) for convolution layers is that two decomposed matrices by SVD do not present corresponding (decomposed) convolution layers, because of intermediate lowering steps. As a result, fine-tuning methods requiring a structurally modified model for training are not available for convolution layers to be compressed by SVD. On the other hand, occasional regularization does not alter the model structure for training. For occasional regularization, SVD can be performed as a way to feed noise into a weight kernel K for every regularization step. Once training stops at a regularization step, the final weight values can be decomposed by SVD and used for inference with reduced memory footprint and computations. In other words, occasional regularization enables SVD-aware training for CNNs. 

A.4.2 TILING-BASED SVD FOR SKEWED WEIGHT MATRICES:

A reshaped kernel matrix K ∈ R T ×(S×d×d) is usually a skewed matrix where row-wise dimension (n) is smaller than column-wise dimension (m) as shown in Figure 17 (i.e., n m). A range of available rank r for SVD, then, is constrained by small n and the compression ratio is approximated to be n/r. If such a skewed matrix is divided into four tiles as shown in Figure 17 and the four tiles do not share much common chateracteristics, then tiling-based SVD can be a better approximator and rank r can be further reduced without increasing approximation error. Moreover, fast matrix multiplication is usually implemented by a tiling technique in hardware to improve the weight reuse rate Fatahalian et al. (2004) . Hence, tiling could be a natural choice not only for high-quality SVD but also for high-performance hardware operations.

To investigate the impact of tiling on weight distributions after SVD, we tested a (1024 × 1024) random weight matrix in which elements follow a Gaussian distribution. A weight matrix is divided by (1 × 1), (16 × 16), or (128 × 128) tiles (then, each tile is a submatrix of (1024 × 1024), (64 × 64), or (8×8) size). Each tile is compressed by SVD to achieve the same overall compression ratio of 4× for all of the three cases. As described in Figure 17 (on the right side), increasing the number of tiles tends to increase the count of near-zero and large weights (i.e., variance of weight values increases). Figure 17 can be explained by sampling theory where decreasing the number of random samples (of small tile size) increases the variance of sample mean. In short, tiling affects the variance of weights after SVD (while the impact of such variance on model accuracy should be empirically studied). We applied the tiling technique and SVD to the 9 largest convolution layers of ResNet-32 using the CIFAR-10 dataset. Weights of selected layers are reshaped into 64 × (64 × 3 × 3) matrices with the tiling configurations described in Table 5 . We perform training with the same learning schedule and pNR(=200) used in Section 3. Compared to the test accuracy of the pre-trained model (=92.63%), all of the compressed models in 18 use the same learning rate scheduling and the number of training epochs as described in Section 3. Results show that tiled SVD yields the best test accuracy and test accuracy is not highly sensitive to tile configuration. SVD presents competitive model accuracy for small compression ratios. As compression ratio increases, however, model accuracy using SVD significantly degrades. From Figure 18 , tiled SVD associated with occasional regularization is clearly the best low-rank approximation scheme. We compare Tucker decomposition trained by a typical fine-tuning process and tiled SVD trained by occasional regularization using the VGG19 model2 with CIFAR-10. Since this work mainly discusses compression on convolution layers, fully-connected layers of VGG19 are compressed and fixed before compression of convolution layers (refer to Appendix for details on the structure of VGG19). Except for small layers with S < 128 (that presents small compression ratio as well), all convolution layers are compressed with the same compression ratio. During 300 epochs to train convolution layers, learning rate is initially 0.01 and is then halved every 50 epochs. In the case of tiled SVD, pNR is 300 and tile size is fixed to be 64×64 (recall that the choice of pNR and tile size do not affect model accuracy significantly). As described in Table 6 , while Tucker decomposition with conventional fine-tuning shows degraded model accuracy through various R c , occasional-regularization-assisted tiled SVD presents noticeably higher model accuracy. We also test our proposed low-rank approximation training technique with the ResNet-34 model 3 He et al. (2016) using the ImageNet dataset Russakovsky et al. (2015) . A pre-trained ResNet-34 is fine-tuned for Tucker decomposition (with conventional training) or tiled SVD (with occasional regularization) using the learning rate of 0.01 for the first 20 epochs, 0.001 for the next 30 epochs, and 0.0001 for the remaining 30 epochs. Similar to our previous experiments, the same compression ratio is applied to all layers except the layers with S < 128 (such exceptional layers consist of 1.4% of the entire model). In the case of Tucker decomposition, selected convolution layers are compressed with R c = 0.46 to achieve an overall compression of 3.1×. For tiled SVD, lowered matrices are tiled and each tile of (64×64) size is decomposed with r=10 to match an overall compression of 3.1×. As shown in Figure 19 , occasional-regularization-based tiled SVD yields better training loss and test accuracy compared to Tucker decomposition with typical training. At the end of the training epoch in Figure 19 , tiled SVD and Tucker decomposition achieves 73.00% and 72.31% for top-1 test accuracy, and 91.12% and 90.73% for top-5 test accuracy, while the pre-trained model shows 73.26% (top-1) and 91.24% (top-5).

A.6 LOWERING TECHNIQUE FOR CNNS:

Figure 20 describes a kernel matrix reshaped from a 4D kernel tensor and an input feature map matrix in the form of a Toeplitz matrix. At the cost of redundant memory usage to create a Toeplitz matrix, lowering enables matrix multiplication which can be efficiently implemented by BLAS libraries. A kernal matrix can be decomposed by 2D SVD. 

A.7 MODEL DESCRIPTIONS FOR LOW-RANK APPROXIMATION EXPERIMENTS:

In this section, we describe model structures and layers selected for low-rank approximation experiments. Small layers close to the input are not compressed because both weight size and compression rate are too small.

:

https://github.com/larry0123du/Decompose-CNN

https://github.com/chengyangfu/pytorch-vgg-cifar10

https://pytorch.org/docs/stable/torchvision/models.html

Figures and Tables:

Figure 1: Distribution of weight noise after low-rank approximation via SVD (Left) or quantization (Right) when R is the rank and Q is the number of quantization bits.

Figure 3: Model accuracy of ResNet-32 on CIFAR-10 using various NR period and amount of weight decay or noise for regularization (original model accuracy without regularization is 92.6%). (Left): Weight decay. (Right): Uniform weight noise insertion.

Figure 4: Perplexity of LSTM model on PTB dataset using various NR period and amounts of weight decay or noise. (Left): Weight decay. (Right): Uniform weight noise insertion.

Figure 3 shows model accuracy of ResNet-32 given different NR period and weight decay factors. For both weight decay and weight noise insertion, the choice of θ (representing the amount of each weight regularization) has a clear correlation with NR period (refer to Appendix for training and test accuracy graphs). If we wish to apply stronger weight regularization, then such weight regularization should be conducted less frequently (i.e., larger weight decay factor requires longer NR period) to maximize the regularization effect and achieve high model accuracy. Similar observations are discovered by LSTM model on PTB as shown in Figure 4. Lower perplexity (indicating better generalization) is obtained when NR period increases as weight regularization becomes stronger for each regularization event.

when the weights are compressed by quantization based on binary codes. when the weights are compressed by magnitude-based pruning.

Figure 5: Model accuracy of an LSTM model on PTB compressed by quantization, low-rank approximation or pruning. Original perplexity without model compression is 114.6. For more details, refer to Figure 14.

Figure 7: Test accuracy comparison on ResNet-32 using CIFAR-10 trained by typical training method and the proposed training method with various compression ratios. For the proposed scheme, test accuracy is measured only at Step 3 that allows to extract a decomposed structure, and pNR is 200.

Figure 9: Model accuracy of ResNet-32 on CIFAR-10 using various NR period and amount of weight regularization (original model accuracy without regularization is 92.6%). (Left): Weight decay. (Right): Uniform weight noise insertion.

Figure 13: Training and test accuracy of ResNet-32 on CIFAR-10 for variuos NR period and the amount of weight regularization. (Left): Weight decay. (Right): Uniform weight noise.

115.4 115.3 114.4 115.1 115.0 115.0 114.4 113.6 115.0 114.0 113.2 113.2 113.2 112.7 114.1 116.1 117.3 118.6 119.0 120.6 122.6 123.6 124.5 127.0 128.0 114.

138.8 135.2 131.1 128.0 126.8 125.0 120.3 118.3 115.9 113.2 113.0 112.7 113.1 115.3 116.6 119.4 120.3 123.6 126.2 128.0 133.2 134.1 135.1 139.4 133.

Figure 14: Relationship between test perplexity and pNR using PTB LSTM model.

Figure 17: Skewed matrix and a tiling technique are illustrated on the left side, while the right side presents distributions of weights after SVD with different tiling schemes (only positive weights are included).

Tiled SVD(64x64,r=10) Typical Scheme, Tucker Decomp.

Figure 19: Comparison of two compression schemes on training loss and (top-1) test accuracy of ResNet-34 model using ImageNet. pNR=500.

Figure 20: An example of lowering technique using im2col.

Pruning rate comparison using LeNet-300-100 and LeNet-5 models on MNIST dataset. DC (Deep Compression) and Sparse VD represent a magnitude-based technique(Han et al., 2016) and variational dropout method(Molchanov et al., 2017), respectively.

Table content: <table><tr><td>Model</td><td>Layer</td><td>Weight Size</td><td colspan="4">Pruning Rate (%) DC DNS Sparse VD Ours</td></tr><tr><td/><td>FC1</td><td colspan="2">235.2K 92</td><td>98.2</td><td>98.9</td><td>98.9</td></tr><tr><td>LeNet-300-100</td><td>FC2 FC3</td><td>30K 1K</td><td>91 74</td><td>98.2 94.5</td><td>97.2 62.0</td><td>96.0 62.0</td></tr><tr><td/><td colspan="3">Total 266.2K 92</td><td>98.2</td><td>98.6</td><td>98.4</td></tr><tr><td/><td>Conv1</td><td>0.5K</td><td>34</td><td>85.8</td><td>67</td><td>60.0</td></tr><tr><td/><td>Conv2</td><td>25K</td><td>88</td><td>96.9</td><td>98</td><td>97.0</td></tr><tr><td>LeNet-5</td><td>FC1</td><td>400K</td><td>92</td><td>99.3</td><td>99.8</td><td>99.8</td></tr><tr><td/><td>FC2</td><td>5K</td><td>81</td><td>95.7</td><td>95</td><td>95.0</td></tr><tr><td/><td colspan="3">Total 430.5K 92</td><td>99.1</td><td>99.6</td><td>99.5</td></tr></table>

to convolutional neural network (CNN) models and demonstrate superiority of pNRbased scheme over conventional training methods. In CNNs, the convolution operation requires a 4D kernel tensor K = R d×d×S×T where each kernel has d × d dimension, S is the input feature map size, and T is the output feature map size. Then, following the Tucker decomposition algorithm, K is decomposed into three components as

Table content: <table><tr><td>Rs</td><td>Rt</td></tr><tr><td>Ki,j,s,t =</td><td>C i,j,rs,rt P S s,rs P T t,rt ,</td></tr><tr><td>rs=1</td><td>rt=1</td></tr></table>

Model accuracy comparison using LeNet-300-100 and LeNet-5 models on MNIST dataset. DC (Deep Compression) and Sparse VD represent a magnitude-based technique(Han et al., 2016) and variational dropout method(Molchanov et al., 2017), respectively.

Table content: <table><tr><td>Model</td><td colspan="3">Accuracy (%) DC DNS Sparse VD Ours</td></tr><tr><td colspan="2">LeNet-300-100 98.4 98.0</td><td>98.1</td><td>98.1</td></tr><tr><td>LeNet-5</td><td>99.2 99.1</td><td>99.2</td><td>99.1</td></tr></table>

Comparison on perplexity using various pruning rates. p f is the target pruning rates for the embedded layer, LSTM layer, and softmax layer.

Table content: <table><tr><td>Model Size</td><td>Pruning Method</td><td>p f =</td><td>0%</td><td>80%</td><td>Perplexity 85% 90%</td><td>95% 97.5%</td></tr><tr><td>Medium</td><td>(Zhu &amp; Gupta, 2017)</td><td/><td colspan="4">83.37 83.87 85.17 87.86 96.30 113.6</td></tr><tr><td>(19.8M)</td><td>Proposed Scheme</td><td/><td colspan="4">83.78 81.54 82.62 84.64 93.39 110.4</td></tr><tr><td>Large</td><td>(Zhu &amp; Gupta, 2017)</td><td/><td colspan="4">78.45 77.52 78.31 80.24 87.83 103.20</td></tr><tr><td>(66M)</td><td>Proposed Scheme</td><td/><td colspan="4">78.07 77.39 77.73 78.28 84.69 99.69</td></tr></table>

Test accuracy(%) of ResNet-32 model using CIFAR-10 dataset while the 9 largest convolution layers (T =S=64, d=3) are compressed by SVD using different tiling configurations. For each tile size, rank r is selected to achieve compression ratio of 2× or 4×. pNR=200 is used for occasional regularization.

Table content: <table><tr><td>Pre-Trained</td><td>Compression Ratio</td><td>64×64</td><td>Size of Each Tile 32×32 16×16</td><td>8×8</td></tr><tr><td>92.63</td><td>2× 4×</td><td colspan="3">93.34 (r=16) 93.11 (r=8) 93.01 (r=4) 93.23 (r=2) 92.94 (r=8) 92.97 (r=4) 93.00 (r=2) 92.81 (r=1)</td></tr></table>

Table5achieves higher model accuracy due to the regularization effect of our compression scheme. Note that for each target compression ratio, the relationship between tile size and model accuracy is not clear. Hence, various configurations of tile size need to be explored to enhance model accuracy, even though variation of model accuracy for different tile size is small.Figure18: Test accuracy of ResNet-32 model using CIFAR-10 with various target compression ratio and decomposition methods. Except the first small convolution layer, all layers are compressed by the same compression ratio. Convolution layers can be grouped according to 3 different S values (=16, 32, or 64). For tiled SVD, three groups (of different S) are tiled in (k 1 × k 1 ), (k 2 × k 2 ), or (k 3 × k 3 ) tile size. (k 1 , k 2 , k 3 ) configuration is described in legends.

Table content: <table><tr><td/><td>93</td><td/><td/><td/><td/><td/><td/></tr><tr><td/><td>93</td><td>93.03(r=5) 92.98 92.92</td><td>93.11(r=2) 92.96(r=4) 92.78(r=8)</td><td/><td>92.78(r=3)</td><td/><td/></tr><tr><td/><td/><td/><td>92.78 92.62</td><td/><td>92.78(r=6)</td><td/><td colspan="2">Pre-Trained(92.63)</td></tr><tr><td>Accuracy(%)</td><td>91 92 92</td><td colspan="3">92.20 Typical Scheme (pN R=1) 92.09 Tucker Decomp.(pN R=200)) Tiled SVD with (16,32,32)(pN R=200)</td><td>92.07 92.43 91.75</td><td/><td>91.95(r=4) 91.60(r=1) 91.49 91.41(r=2)</td></tr><tr><td/><td/><td colspan="3">Tiled SVD with (16,16,16)(pN R=200)</td><td/><td/><td/></tr><tr><td/><td>91</td><td colspan="3">Tiled SVD with (8,8,8)(pN R=200) SVD(pN R=200)</td><td/><td/><td>90.96</td></tr><tr><td/><td/><td/><td/><td/><td/><td/><td>90.84</td></tr><tr><td/><td>1</td><td>1</td><td>2</td><td>2</td><td>3</td><td>3</td><td>4</td><td>4</td></tr><tr><td/><td/><td/><td colspan="3">Compression Ratio</td><td/><td/></tr></table>

Comparison on various low-rank approximation schemes of VGG19 (using CIFAR-10 dataset). To focus on convolution layers only, fully-connected layers are compressed by 8× and trained by occasional regularization. Then, fully-connected layers are frozen and convolution layers are compressed (except small layers of S < 128) by Tucker decomposition or tiled SVD.

Table content: <table><tr><td>Comp. Scheme</td><td>Parameter</td><td>Weight Size</td><td>FLOPs</td><td>Accuracy(%)</td></tr><tr><td>Pre-Trained</td><td>-</td><td>18.98M</td><td>647.87M</td><td>92.37</td></tr><tr><td>Tucker</td><td>R c =0.6</td><td colspan="2">9.14M (2.08×) 319.99M (2.02×)</td><td>91.97</td></tr><tr><td>Decomposition</td><td>R c =0.5</td><td colspan="2">6.71M (2.83×) 235.74M (2.75×)</td><td>91.79</td></tr><tr><td>(Typical</td><td>R c =0.45</td><td colspan="2">5.49M (3.45×) 191.77M (3.38×)</td><td>91.36</td></tr><tr><td>Scheme)</td><td>R c =0.4</td><td colspan="2">4.61M (4.11×) 161.60M (4.01×)</td><td>91.11</td></tr><tr><td>Tiled</td><td colspan="3">64×64 (r=16) 9.49M (2.00×) 316.28M (2.04×)</td><td>92.42</td></tr><tr><td>SVD</td><td colspan="3">64×64 (r=11) 6.52M (2.91×) 214.25M (3.02×)</td><td>92.33</td></tr><tr><td>(Occasional</td><td colspan="3">64×64 (r=10) 5.93M (3.20×) 193.85M (3.34×)</td><td>92.23</td></tr><tr><td>Regularization,</td><td>64×64 (r=9)</td><td colspan="2">5.55M (3.41×) 173.44M (3.73×)</td><td>92.22</td></tr><tr><td>pNR=300)</td><td>64×64 (r=8)</td><td colspan="2">4.74M (4.00×) 153.04M (4.33×)</td><td>92.07</td></tr></table>

Peer review (iclr2020):
I think that the paper is very well written, I like it. The authors localized a phenomenon and demonstrated how to exploit it. I trust the results because I performed exactly the same experiments for CIFAR-10 with longer non-regularization periods and found that there is no effect (this is also that the authors show in the paper) but I didn*t test on other datasets and obviously didn*t think about potential benefits for compression. Since model compression is not my field, I would just make a general suggestion to include and compare to state-of-the-art techniques. --- Update: The authors updated the paper and replied to all reviewers in detail. I still think that the paper should be accepted.