SparseProp: Efficient Sparse Backpropagation for Faster Training of
  Neural Networks
We provide a new efficient version of the backpropagation algorithm,
specialized to the case where the weights of the neural network being trained
are sparse. Our algorithm is general, as it applies to arbitrary (unstructured)
sparsity and common layer types (e.g., convolutional or linear). We provide a
fast vectorized implementation on commodity CPUs, and show that it can yield
speedups in end-to-end runtime experiments, both in transfer learning using
already-sparsified networks, and in training sparse networks from scratch.
Thus, our results provide the first support for sparse training on commodity
hardware.
Mahdi Nikdan, Tommaso Pegolotti, Eugenia Iofinova, Eldar Kurtic, Dan Alistarh
# SparseProp: Efficient Sparse Backpropagation

###### Abstract

We provide a new efficient version of the backpropagation algorithm, specialized to the case where the weights of the neural network being trained are _sparse_. Our algorithm is general, as it applies to arbitrary (unstructured) sparsity and common layer types (e.g., convolutional or linear). We provide a fast vectorized implementation on commodity CPUs, and show that it can yield speedups in end-to-end runtime experiments, both in transfer learning using already-sparsified networks, and in training sparse networks from scratch. Thus, our results provide the first support for sparse training on commodity hardware.

Machine Learning, SparseProp, Efficient Sparse Backpropagation

## 1 Introduction

The significant computational costs of deep learning have led to massive interest in approaches for leveraging _sparsity_ in neural networks, which have been investigated in great breadth and depth (Hoefler et al., 2021). On the inference side, there is already emerging algorithmic and system support for sparsity on both GPUs (Mishra et al., 2021; Gale et al., 2020) and CPUs (Elsen et al., 2020; NeuralMagic, 2022), as well as a wide range of methods for obtaining models which are both highly-sparse and highly-accurate.

A new frontier in the area is _accurate and efficient sparse training_. On the algorithmic side, there are several interesting proposals for _sparse training_ algorithms (Dettmers and Zettlemoyer, 2019; Kusupati et al., 2020; Evci et al., 2020; Jayakumar et al., 2020; Schwarz et al., 2021), i.e. variants of stochastic gradient descent (SGD) which aim to keep as many weights as possible sparse during training. Another interesting related approach is _sparse transfer_(Zafrir et al., 2021; Chen et al., 2021; Iofinova et al., 2022; Kurtic et al., 2022), by which models sparsified on a large pretraining corpus are then used for _transfer learning_ on different tasks, while preserving the sparsity mask.

Despite this progress on the optimization side, the vast majority of these approaches lack _system support for fast training_, in that they do not provide any practical speedups. This is because the weight sparsity introduced is _unstructured_, which is notoriously hard to leverage for computational gains. Specifically, there is no general implementation of backpropagation that can leverage unstructured weight sparsity for practical speedup on common hardware. At the same time, approaches leveraging _block sparsity_(Mishra et al., 2021; Gray et al., 2017) can only reach lower sparsity without significant accuracy drops, and require specialized training algorithms (Lagunas et al., 2021; Jiang et al., 2022). As such, unstructured weight sparsity is often dismissed as a practical way of accelerating model training.

**Contribution.** We contradict this conventional wisdom by presenting a new vectorized implementation of backpropagation (Rumelhart et al., 1986), designed to be efficient in the case where the weights of the neural network are _sparse_, i.e. contain a significant fraction of zero values, and show its potential for practical speedups in common edge training scenarios, for both vision and language tasks.

More precisely, our algorithm, called SparseProp, is general in the sense that 1) it applies to arbitrary sparsity patterns, 2) general layer types, and 3) can be efficiently vectorized using standard CPU-supported approaches. The asymptotic complexity of the algorithm is _linear_ in the _layer density_, i.e. the number of non-zero weights in the layer, providing proportional runtime improvements to the weight sparsity, for both linear and convolutional layers.

To illustrate practical efficiency, we provide a fast vectorized implementation of SparseProp aimed at general-purpose Intel and AMD CPU architectures. Specifically, our implementation provides drop-in replacement implementations for standard layer types, and only relies on widely-supported AVX2 instructions. We show that SparseProp can lead to practical runtime improvements both on single sparse layers, validating our linear sparsity scaling claims, as well as on end-to-end training of sparse models.

Specifically, we provide results for a preliminary integration with Pytorch (Paszke et al., 2019), which can run sparse backpropagation for linear and convolutional layers,
covering most popular model families. As such, SparseProp can provide direct support for methods like Gradual Pruning (Zhu & Gupta, 2017), RigL (Evci et al., 2020) or AC/DC (Peste et al., 2021), which assume a fixed sparsity mask for any fixed forward and backward pass, and can be modified to support more complex methods (Jayakumar et al., 2020), which specify different sparsities for weights and gradients. We believe it is the first implementation to do so on commodity hardware.

Our end-to-end experiments aim to make the case that sparsity can be a viable option for DNN training _at the edge_. That is, we explore settings where a device with moderate computational power (e.g., a CPU with a limited number of cores) performs either _sparse transfer_ or _from-scratch sparse training_ over a specialized task. We investigate sparsity-versus-accuracy trade-offs in two model/task combinations: 1) ResNets (He et al., 2016) applied to twelve popular vision tasks (Kornblith et al., 2019), and 2) a standard BERT-base model (Devlin et al., 2019) applied to GLUE language modelling tasks (Wang et al., 2018).

In the _sparse transfer_ scenario, we are provided an already-sparse model pretrained on a large corpus, e.g. ImageNet (Russakovsky et al., 2015) respectively WikiText (Merity et al., 2016), and wish to finetune the corresponding sparse weights on a (usually smaller) target dataset. This application has gained significant popularity (Zafrir et al., 2021; Chen et al., 2021; Iofinova et al., 2022; Kurtic et al., 2022), and pretrained sparse models are available for several standard tasks and models (Wolf et al., 2019; SparseZoo, 2022). In this context, we show that, for both vision and language tasks, SparseProp can lead to end-to-end sparse transfer speedups of up to 1.85x, at similar accuracies, relative to CPU-based finetuning of _dense models_ in the same environment. Measured only over backward-pass operations--and thus omitting framework-level overheads-our algorithm provides speedups of 3.6x at 95% model sparsity.

In the second scenario, we examine the ability of SparseProp to provide speedups for _sparse training from scratch_, on the same series of tasks, adapting variants of sparse training (Zhu & Gupta, 2017) to our setting. Experiments show that, in this scenario, SparseProp leads to end-to-end speedups of up to 1.4x, with moderate accuracy loss.

In sum, our results show that SparseProp can efficiently provide system support for CPU-based unstructured sparse training, ensuring speedups for both from-scratch training and sparse transfer. We believe our approach could lead to additional practical impact for research on sparsity, especially given that our end-to-end runtime numbers can still be improved via additional optimizations, and by mitigating external, framework-specific overheads.

## 2 Related Work

**Sparse Inference.** One of the key motivations behind sparsity in DNNs is reducing inference costs. For this, an impressive number of weight pruning techniques have been introduced, e.g. (LeCun et al., 1990; Hagiwara, 1994; Han et al., 2016; Singh & Alistarh, 2020; Sanh et al., 2020). Complementing this work, there have been a number of algorithmic proposals for efficient sparse inference algorithms over DNNs, e.g. (Park et al., 2016; Han et al., 2016; Gale et al., 2020; Elsen et al., 2020), although it is known that layer-wise gains can be difficult to translate into end-to-end speedups (Wang, 2020). Nevertheless, sparse inference support is now available on both CPUs, e.g. (NeuralMagic, 2022) and GPUs (Mishra et al., 2021).

Hubara et al. (2021) proposed a theoretically-justified approach for identifying sparse transposable masks matching the NVIDIA 2:4 sparsity pattern, which could be leveraged for faster training on GPUs. However, they do not provide an implementation, and, currently, GPU-based 2:4 sparsity speedups tend to be minimal (NVIDIA, 2021).

**Sparse SGD-Based Training.** As noted, there has been a significant amount of work on SGD-like algorithms for _sparse training_ of DNNs, balancing accuracy while trying to maximize sparsity in the models' internal representations (Mocanu et al., 2016; Bellec et al., 2018; Zhu & Gupta, 2017; Mostafa & Wang, 2019; Lis et al., 2019; Dettmers & Zettlemoyer, 2019; Zhang et al., 2020; Wiedemann et al., 2020; Kusupati et al., 2020; Evci et al., 2020; Jayakumar et al., 2020; Peste et al., 2021; Schwarz et al., 2021). Unfortunately, a precise comparison is quite difficult, since each makes different assumptions regarding the degree of sparsity in the network's internal representations, potentially even varying the amount of sparsity between weights and gradients, e.g. (Jayakumar et al., 2020).

**Sparse Training for Speedup.** Leveraging sparsity for practical speedups has been a major goal in model compression (Hoefler et al., 2021). Yang et al. (2020) proposed a specialized hardware accelerator which is specialized to the DropBack pruning algorithm (Lis et al., 2019). SWAT (Raihan & Aamodt, 2020) proposed a sparsity-aware algorithm for the case where both weights and activations have high sparsity, and showed speedups in a simulated environment. Their approach works only for specific networks, and can lose significant accuracy. More recently, Jiang et al. (2022) proposed an algorithm-hardware co-design approach, for the case of GPU-based training. Specifically, their approach imposes _block sparsity_ in GPU-friendly patterns, and leverages it for speedup.

By contrast to this work, our approach considers efficient support for backpropagation for _unstructured sparse weights_, implements this efficiently for commodity CPUs, and shows that this can be leveraged for _end-to-end speedup
during training. Specifically, this provides support to the vast amount of existing work on unstructured sparse training algorithms, on commodity hardware.

**System Support.** Pytorch (Paszke et al., 2019) introduced partial sparse tensor support recently, while the STen (Ivanov et al., 2022) provides a general interface for such representations. Our work is complementary to this direction, as our implementation can be interfaced with Pytorch or specifically STen to provide training speedups.

## 3 The Sparse Backpropagation Algorithm

### Background

**SIMD Instructions.** Fast and efficient numerical code heavily relies on Single Instruction Multiple Data (SIMD) instructions to improve performance. These instructions operate on specialized machine registers (xmm, ymm, and zmm) that contain multiple values.

Our implementations currently support x86 machines that provide the standard AVX2 instruction set, which uses \(256\) bit registers, or \(8\) single precision floating point values. Table 1 provides an overview of the instructions employed by our library. Our SIMD implementation structure follows the Load-Compute-Store paradigm, where data is explicitly transferred to registers via the **loadv** and **broadcastv** instructions. Computation is performed on the data in the registers using fused multiply-add instructions **vfmadd** (\(r=a\cdot b+c\)), and the results are subsequently moved back to memory with the **vstore** instruction. Following this structure, we significantly increase performance since the data loaded into the registers can be used for multiple operations before being stored back in memory.

**Backpropagation.** Let \(f(\mathbf{X};\mathbf{W})\) represent a layer (fully-connected or convolution) in a neural network \(\mathcal{N}\); \(\mathbf{W}\) represents the parameters of this layer and \(\mathbf{X}\) represents a batch of inputs. Let \(B\) be the batch size. Additionally, denote the output of this layer by \(\mathbf{O}=f(\mathbf{X};\mathbf{W})\). Let \(L\) be the loss of the whole network \(\mathcal{N}\) for this batch of inputs. Backpropagating through this layer involves calculating the gradients \(\partial L/\partial\mathbf{W}\) and \(\partial L/\partial\mathbf{O}\), given \(\partial L/\partial\mathbf{O}\).

Consider the situation where we have a highly sparsified matrix \(\mathbf{W}\) that is stored as a sparse matrix. During the backpropagation process, it is necessary to calculate the gradient of this matrix. However, in practice, the full gradient of the dense matrix is often calculated, even though the pruned elements are not updated and their gradients are discarded. This can be inefficient, as it consumes a significant amount of computation and time without providing any benefits.

### The Case of Fully-connected Layers

We now focus on the case where \(f(.)\) is a fully-connected layer. Assume \(\mathbf{X}\) and \(\mathbf{W}\) are \(B\times M\) and \(M\times N\) matrices, respectively. Consequently, \(\mathbf{O}=f(\mathbf{X};\mathbf{W})=\mathbf{X}\mathbf{W}\) will be a \(B\times N\) matrix. The gradients of \(L\) with respect to \(\mathbf{X}\) and \(\mathbf{W}\) are calculated as follows:

\[\frac{\partial L}{\partial\mathbf{X}}=\frac{\partial L}{\partial \mathbf{O}}\mathbf{W}^{T} \tag{1}\] \[\frac{\partial L}{\partial\mathbf{W}}=\mathbf{X}^{T}\frac{ \partial L}{\partial\mathbf{O}} \tag{2}\]

If we examine equations (1) and (2), we can see that the former is a General Sparse Matrix-Matrix Multiplication (SpGEMM) operation, while the latter is a Sampled Dense Dense Matrix Multiplication (SDDMM) operation.

**Sparse Representation.** The matrix \(\mathbf{W}\) is stored in a compressed sparse row (CSR) format, a standard representation for sparse matrices. The non-zero values of the matrix are stored in the arrays \(W_{\mathrm{vals}}\) and \(W_{\mathrm{cols}}\), which correspond to the values and column indices of the non-zero elements, respectively. The array \(W_{\mathrm{rows}}\) encodes each row's start and end indices. For example, the non-zero values and column indices of a row \(i\) of \(\mathbf{W}\) are contained between positions \(W_{\mathrm{rows}}[i]\) and \(W_{\mathrm{rows}}[i+1]\).

**Algorithm.** In Algorithm 1, we present high-level pseudocode for backpropagation in our linear layer. The calculations for (1) and (2) are performed in a single pass by utilizing the sparsity pattern of \(\mathbf{W}\), which is identical to \(\partial L/\partial\mathbf{W}\). Specifically, the result of \((\partial L/\partial\mathbf{O})\mathbf{W}^{T}\) is computed as a sparse matrix-matrix multiplication. Whereas \(\mathbf{X}^{T}(\partial L/\partial\mathbf{O})\) is computed as an SDDMM, with \(\mathrm{nnz}\) dot-products, where \(\mathrm{nnz}\) is the number of non-zero elements of \(\mathbf{W}\). In more detail, the computation is divided into \(3\) loops. The innermost loop contains the core of the computation. It computes at each iteration \(16\) floating point operations using \(2\)**fmadd** instructions: the first **fmadd** computes \(8\) entries of \(\partial L/\partial\mathbf{X}\) and the second accumulate a dot-product in a register acc.

**Implementation Details.** We operate on the transposed matrices in our linear layer implementation to improve cache utilization. Specifically, in both the forward and backward passes, we operate on the transposed version of the input matrix, \(\mathbf{X}^{T}\), which is a column-major representation of \(\mathbf{X}\). By doing so, we achieve a streaming access pattern over the rows of \(\mathbf{X}^{T}\), as we see from the innermost loop Algorithm 1.

\begin{table}
\begin{tabular}{l l} \hline \hline
**vload**(address) & load from memory address \\
**vstore**(address, a) & store a at memory address \\
**vbroadcast**(a) & fill a register with a \\
**vfmadd**(a, b, c) & return \(a\cdot b+c\) \\
**vaddreduce**(a) & return sum elements of a \\ \hline \hline \end{tabular}
\end{table}
Table 1: List of vector instructions used in the implementation and their semantics.
Additionally, we leverage that the transpose of a CSR matrix is none other than the same matrix in Compressed-Sparse-Column (CSC) format to avoid expensive sparse transpose operations. An example computation is given in Figure 1.

### Sparse Backpropagation for Convolutional Layers

We now examine convolutional layers. Denote the number of input channels by \(IC\), and the number of output channels by \(OC\). The input width and height are represented by \(M\) and \(N\), respectively, and let the kernel size be \(K\times K\). The input tensor, \(\mathbf{X}\), and the weights tensor, \(\mathbf{W}\), have dimensions \(B\times IC\times M\times N\) and \(OC\times IC\times K\times K\), respectively. For simplicity, here we only consider the case where padding is \(0\) and stride is \(1\). For larger padding and stride, the generalization is not difficult. The output tensor, \(\mathbf{O}\), will be of size \(B\times OC\times OM\times ON\), where \(OM=M-K+1\) and \(ON=N-K+1\). For \(0\leq b<B,0\leq oc<OC,0\leq p<OM,0\leq q<ON\) we have:

\[\mathbf{O}[b,oc,p,q]=\sum_{ic=0}^{IC-1}\sum_{i=0}^{K-1}\sum_{j=0}^ {K-1}\mathbf{W}[oc,ic,i,j] \tag{3}\] \[.\mathbf{X}[b,ic,p+i,q+j].\]

Using the chain rule, it is easy to check that for \(0\leq b<B,0\leq ic<IC,0\leq m<M,0\leq n<N\):

\[\frac{\partial L}{\partial\mathbf{X}}[b,ic,m,n]=\sum_{oc=0}^{OC-1} \sum_{p=p_{s}}^{m}\sum_{q=q_{s}}^{n}\frac{\partial L}{\partial\mathbf{O}}[b, oc,p,q] \tag{4}\] \[.\mathbf{W}[oc,ic,m-p,n-q],\]

with \(p_{s}=m-K+1\), and \(q_{s}=n-K+1\). And for \(0\leq oc<OC,0\leq ic<IC,0\leq i<K,0\leq j<K\), we have:

\[\frac{\partial L}{\partial\mathbf{W}}[oc,ic,i,j]=\sum_{b=0}^{B-1} \sum_{p=0}^{M-K}\sum_{q=0}^{N-K}\frac{\partial L}{\partial\mathbf{O}}[b,oc,p,q] \tag{5}\] \[.\mathbf{X}[b,ic,p+i,q+j].\]

It is assumed that the weight matrix \(\mathbf{W}\) is sparse. In accordance with equation (4), when a weight \(\mathbf{W}[oc,ic,m-p,n-q]\) is pruned, the multiplication and corresponding addition operations can be skipped. Furthermore, when a weight \(\mathbf{W}[oc,ic,i,j]\) is pruned, the calculation of the gradient for this parameter is not necessary, as it will not be updated, and therefore the computation outlined in equation (5) can be skipped.

```
for\(i=0\) to \(N-1\)do for\(j=\mathbf{W}_{\mathrm{coh}}[i]\)do \(\prime\)\(\prime\) repeat one weight entry 8 times*/ \(v\leftarrow\textbf{broadcast}(\mathbf{W}_{\mathrm{vals}}[j]\) /* initialize \(acc\) to zero*/ \(acc\leftarrow\textbf{vbroadcast}(0)\) \(r\leftarrow\textbf{W}_{\mathrm{rows}}[j]\) for\(k=0\) to \(B-1\)with \(k\) += 8do /* load * values*/ \(dx\leftarrow\textbf{vload}((\partial L/\partial\mathbf{X})_{r,k})\) \(do\leftarrow\textbf{vload}((\partial L/\partial\mathbf{O})_{r,k})\) \(x\leftarrow\textbf{vload}(\mathbf{X}_{i,k})\) /* compute \(8\)\(dx\leftarrow\textbf{vfindact}(do,v,dx)\) /* compute \(8\)\(acc\gets do\cdot x+acc*/\) \(acc\leftarrow\textbf{vfindact}(do,x,acc)\) /* store updated \(dx\) back*/ \(\textbf{vstore}((\partial L/\partial\mathbf{X})_{r,k},dx)\) endfor /* sum the 8 values in \(acc*/\) \((\partial L/\partial\mathbf{W})_{\mathrm{vals}}[j]\leftarrow\textbf{vaddreduce }(acc)\) endfor endfor
```

**Algorithm 1** AVX2 Linear Backward Pass

**Sparse Representation.** To efficiently represent a sparse tensor, we employ a representation akin to the compressed sparse row (CSR) format used for sparse matrices. Four arrays, \(W_{\mathrm{och}}\), \(W_{\mathrm{ich}}\), \(W_{\mathrm{x}}\), and \(W_{\mathrm{y}}\), are used to store the indices, and an array \(W_{\mathrm{vals}}\) is used to store the non-zero values. Specifically, \(W_{\mathrm{x}}\) and \(W_{\mathrm{y}}\) are arrays of size \(\mathrm{nnz}\), which contain the coordinates of the non-zero values of each filter. \(W_{\mathrm{och}}\) is an array of size \(OC+1\), which encodes the start of each output channel's entries in \(W_{\mathrm{ich}}\). Finally, \(W_{\mathrm{ich}}\) is an array of size \(OC\times(IC+1)\), which encodes the indices in \(W_{\mathrm{x}}\), \(W_{\mathrm{y}}\), and \(W_{\mathrm{vals}}\) of each input channel. For example, for an output channel \(oc\), the non-zero elements of the input channel \(ic\) are stored between indices \(W_{\mathrm{och}}[oc]+W_{\mathrm{ich}}[oc\cdot(IC+1)+ic]\) and \(W_{\mathrm{och}}[oc]+W_{\mathrm{ich}}[oc\cdot(IC+1)+ic+1]\). As example, consider the following sparse tensor of dimensions \((3,2,2,3)\)

Figure 1: Visual representation of the core computation of Algorithm 1 using vector registers of size 4. We represent elementwise multiplication with \(\times\) and **vaddreduce** with \(\bigoplus\).
\[\left\{\begin{bmatrix}\cdot&a&\cdot\end{bmatrix}\begin{bmatrix}\cdot&\cdot&\cdot \\ \cdot&\cdot&\cdot\end{bmatrix}\begin{bmatrix}\cdot&\cdot&\cdot\\ b&\cdot&c\end{bmatrix}\right\},\] \[\left\{\begin{bmatrix}\cdot&\cdot&\cdot\\ \cdot&\cdot\end{bmatrix}\begin{bmatrix}\cdot&d&\cdot\\ e&e\end{bmatrix}\begin{bmatrix}\cdot&\cdot&f\\ \cdot&\cdot&\cdot\end{bmatrix}\right\},\]

where \(\cdot\) represents a zero value. Its sparse representation is given by

\[W_{\mathrm{och}} =\begin{pmatrix}0&1&3&6\end{pmatrix},\] \[W_{\mathrm{ich}} =\begin{pmatrix}0&1&1&0&0&2&0&2&3\end{pmatrix},\] \[W_{\mathrm{x}} =\begin{pmatrix}0&0&1&1&1&0\end{pmatrix},\] \[W_{\mathrm{y}} =\begin{pmatrix}1&1&1&0&2&2\end{pmatrix},\] \[W_{\mathrm{vals}} =\begin{pmatrix}a&b&c&d&e&f\end{pmatrix}.\]

Although the usage of a \(W_{\mathrm{och}}\) may seem superfluous, its inclusion allows us to reduce total memory usage. Indeed, assuming \(K<256\) and \(IC<8192\), we can store entries in \(W_{\mathrm{x}}\) and \(W_{\mathrm{y}}\) using uint8_t and \(W_{\mathrm{ich}}\) using int16_t. Therefore, using \(W_{\mathrm{och}}\) lowers memory usage from \(4B\times(IC+1)\times OC\) bytes to \(2B\times(IC+1)\times OC+4\times OC\) bytes.

**Algorithm.** In Algorithm 2, we present an overview of our backpropagation algorithm for a convolutional layer. If we ignore at first the pointer arithmetic needed to traverse the structures, the main structure remains similar to that of Algorithm 1 as both innermost loops use the same instructions.

**Implementation Details.** We developed two kernels for fast 2D convolutions based on the input dimensions. For larger \(M\) and \(N\), we found that no preprocessing was needed. Keeping the input dimensions as \(B\times IC\times M\times N\) offers both computing batches in parallel on multiple threads and high single-core parallelization using AVX2 instructions. On the other hand, for small \(M\) and \(N\), which often occur for the last layers of a network, we found it more efficient to permute the tensors to have the batch as the last index. In particular, \(\mathbf{X}\) and \(\partial L/\partial\mathbf{X}\) became of size \(IC\times M\times N\times B\), and \(\mathbf{O}\) and \(\partial L/\partial\mathbf{O}\) became of size \(OC\times OM\times ON\times B\). Indeed, setting \(B\) as the last dimension allows the usage of SIMD instructions even for small \(M\) and \(N\).

```
for\(ic=0\)to\(IC\)do for\(oc=0\)to\(OC\)do \(si_{i}=\mathbf{W}_{\mathrm{och}}[oc]+\mathbf{W}_{\mathrm{ich}}[ic]\) \(si_{e}=\mathbf{W}_{\mathrm{och}}[oc]+\mathbf{W}_{\mathrm{ich}}[ic+1]\) for\(si=si_{e}\)to /*repeat 8 times */ \(v\leftarrow\)vbroadcast(\(\mathbf{W}_{\mathrm{vals}}[j]\)) * initialize \(acc\) */ \(acc\leftarrow\)vbroadcast(\(0\)) \(p_{*}\leftarrow\max(0,pad-\mathbf{W}_{\mathrm{x}}[si])\) \(p_{e}\leftarrow\min(pad-\mathbf{W}_{\mathrm{x}}[si]+M,OM)\) \(q_{*}\leftarrow\max(0,pad-\mathbf{W}_{\mathrm{y}}[si])\) \(q_{e}\leftarrow\min(pad-\mathbf{W}_{\mathrm{y}}[si]+N,ON)\) for\(p=p_{*}\)do for\(q=q_{*}\)to\(q_{*}\)do for\(k=0\)to\(B-1\)with\(k==8\)do /* load 8 values*/ \(do\leftarrow\)vload(\((\partial L/\partial\mathbf{O})_{oc,p,q,k}\)) \(dx\leftarrow\)vload(\((\partial L/\partial\mathbf{X})_{ic,p,q,k}\)) \(x\leftarrow\)vload(\(\mathbf{X}_{ic,p,q,k}\)) /* compute 8 \(dx\) - \(do\cdot v+dx\)*/ \(dx\leftarrow\)vfmadd(\(do,v,dx\)) /* compute 8 \(acc=do\cdot x+acc\)*/ \(acc\leftarrow\)vfmadd(\(do,acc\)) /* store the updated \(dx\)*/ \(\mathbf{vstore}((\partial L/\partial\mathbf{X})_{ic,p,q,k},dx)\) endfor endfor endfor endfor endfor endfor endfor /* sum the 8 values in \(acc\)*/ \((\partial L/\partial\mathbf{W})_{\mathrm{vals}}(j)\leftarrow\)vaddreduce(\(acc\)) endfor endfor
```

**Algorithm 2** AVX2 Convolutional Backward Pass

## 4 Experiments

**Setup and Goals.** We now experimentally validate our approach. First, we perform an in-depth exploration of our algorithm's runtime relative to weight sparsity in a synthetic scenario, i.e. for standard layer and input shapes. Then, we examine performance for two _end-to-end training_ scenarios, as part of a Pytorch integration. Specifically, we examine performance for _sparse transfer_, i.e. fine-tuning of already-sparsified accurate models on a different "transfer" dataset, and _from-scratch sparse training_ on some specialized tasks.

**Pytorch Modules.** We provide two Pytorch modules, one for linear and one for convolution layers, facilitating the integration of our algorithm to different applications. To this end, we employ the pybind11 library (Jakob et al., 2016) for fast communication between Pytorch and the C** backend.

### Synthetic Performance Evaluation

This section analyzes the runtime scaling for our linear and convolutional layers implementations of Algorithms 1 and 2. To validate our linear-scaling claims, we specifically examine the runtime dependence of our implementations relative to sparsity. These implementations are compared to their dense counterparts available in Pytorch. Additionally, for the linear layer implementation, we compare it against the sparse implementation offered by Pytorch. We present the improvement in our implementations' performance as a function of a layer's sparsity ranging from \(80\%\) to \(99\%\). In particular, we give the runtime for our multithreaded implementations run on \(8\)-threads in Figures 4 and 3, and the runtime for the single-core implementations in Figure 2. To highlight the speedup between implementations, we give all our measurements in _log-lin_ plots.
**Linear layers.** We evaluate the performance of our linear layer by measuring the runtime of forward and backward pass through a large layer of dimensions \((M,N)=(768,3072)\) and an input of size \((B,M)=(902,768)\). We report the single-core results for the backward pass in Figure 2(a), and we compare them to a Pytorch-based sparse implementation, which is currently only single-core on CPUs. Our sequential performance increases linearly with the sparsity of \(\mathbf{W}\), we match a dense implementation around \(90\%\) sparsity, and we obtain a \(5\times\) speedup at \(99\%\). In Figure 4, we plot the results for the multithreaded forward and backward passes. For our multithreaded forward pass, we observe similar behavior to the single-core variant beating the dense implementation after \(95\%\) sparsity. It should be noted that the performance of our multithreaded implementation of the backward pass is currently limited by synchronization overheads, which can be removed with additional optimizations. As of now, it matches the dense implementation only at very high levels of sparsity.

**Convolutional layers.** In the case of the convolutional layers, we present the runtime performance for two distinct scenarios: small input tensor and large input tensor. The dimensions of the layer are set to \((OC,IC,K,K)=(256,128,3,3)\), and the input tensor's dimensions are set to \((B,IC,M,N)=(8,128,7,7)\) and \((B,IC,M,N)=(32,256,244,244)\). These scenarios highlight the difference between the two sparse convolution kernels we developed: SparseConv and SparseConvOverON. The former is presented in Algorithm 2, and the latter is its permuted variant that vectorizes over the \(ON\) dimension of the tensor.

We show that by permuting the input, we substantially improve the performance of our algorithm. The single-core runtime performance over a small input is presented in Figure 2(b), where we observe a significant speedup compared to a dense implementation, with \(19\)x speedup at \(99\%\) sparsity. The parallel runtime for forward and backward passes are in Figure 4. Timings for small input sizes are given in Figures 4(a) and 4(b), and for large inputs in Figures 4(c) and 4(d). We see how permuting impacts our performance for small \(ON\) values achieving a speedup over the dense implementation of up to \(5\)x for the forward and backward pass at \(99\%\) sparsity. On the other hand, for larger \(M\) and \(N\), our results indicate that vectorizing over \(ON\) yields the best performance, resulting in a speedup of \(3.35\)x for the forward pass and \(9\)x for the backward.

### End-to-End Training Experiments

We now evaluate SparseProp on _sparse transfer learning_, and _sparse training from scratch_. In each case we examine

Figure 4: Runtime measurement of our parallel sparse convolutional layers (base and vectorized over the \(ON\) dimension) on \(8\) threads compared against a dense implementation. The layer has dimension \((OC,IC,K,K)=(256,128,3,3)\), and the input’s dimensions are set to \((B,IC,M,N)=(8,128,7,7)\) in (a) and (b) and to \((B,IC,M,N)=(32,256,244,244)\) in (c) and (d).

Figure 3: Runtime measurement of our parallel sparse linear layer compared against a dense implementation. The layer has dimension \((M,N)=(768,3072)\) the input is of size \((B,M)=(902,768)\).

Figure 2: Runtime of our single-core backward propagation of sparse linear and convolutional layers. The former has dimension \((M,N)=(768,3072)\) and the input is of size \((B,M)=(902,768)\), while the latter is of size \((OC,IC,K,K)=(256,128,3,3)\), and the input is \((B,IC,M,N)=(8,128,7,7)\).
sparsity settings which, while maintaining reasonable accuracy, can achieve non-trivial speedups. Since we executed over more than 15 different tasks, on multiple models, the full experiments used to determine accuracy are executed on GPU. At the same time, we computed CPU speedups using proportionally-shortened versions of the training schedules on CPU, and validated correctness in a few end-to-end runs on different tasks.

In all the experiments, dense modules (linear or convolution) are executed dense as long as they are less than \(80\%\) sparse. Once a module reaches at least \(80\%\) sparsity (which may happen during training for gradual pruning scenarios), we measure the time to run one batch through both dense and sparse versions, and choose the fastest version based on forward+backward time (in the case of convolution, both sparse implementations are considered). Notice that in all experiments, the sparsity patterns change only a few times during the whole run, meaning the overhead of these few extra batches is negligible. (An alternative approach would be to generate a static database of the best implementation choices for each layer type and size, and greedily adopt the implementation for each layer in turn.)

#### 4.2.1 Application 1: Sparse Transfer

**Image Classification.** We first consider a standard transfer learning setup for image classification using CNNs, in which a model pretrained on ImageNet-1K has its last (FC) layer resized and re-initialized, and then is further finetuned on a smaller target dataset. As targets, we focus on twelve datasets that are conventionally used as benchmarks for transfer learning, e.g. in (Kornblith et al., 2019; Salman et al., 2020; Iofinova et al., 2022). See Table 5 for a summary of the tasks. Importantly, input images are scaled to standard ImageNet size, i.e. \(224\times 224\times 3\), resulting in proportional computational costs.

We consider dense and sparse ResNet50 models pre-trained (and sparsified) on the ImageNet-1K dataset. Models are pruned using the AC/DC method (Peste et al., 2021), which shows high accuracy on ImageNet-1K, and produces transferable sparse features (Iofinova et al., 2022). (We adopt their publicly-available models.) To explore the accuracy-vs-speedup trade-off, we consider both a _Uniform pruning_ scenario, in which all convolutional layers except for the input are pruned to a uniform sparsity (90% and 97%), and a _Global pruning_ scenario, in which all convolutional layers are pruned jointly, to a target average sparsity (95%), using the global magnitude pruning criterion. The former two models have been trained for 200 epochs each, and have 76.01% and 74.12% Top-1 accuracy, whereas the latter is trained for 100 epochs and has 73.1% Top-1 accuracy. The dense model has 76.8% Top-1 accuracy. For transfer learning, we maintain the sparsity pattern of the models, reinitialize the final fully-connected layer, and train for 150 epochs on each of the 12 "downstream" tasks.

In Figure 5, we aggregated results across all tasks, in terms of mean and variance of the _accuracy drop relative to transferring the dense model_ (the full per-task accuracy results are presented in Table 6, and speedups are presented in Table 2). As expected, the aggregated sparse test accuracy drops relative to the dense baseline, proportionally to the ImageNet Top-1 accuracy. The Uniform-90 model shows the smallest drops (1% on average), but also the lowest end-to-end speedup (25%), while the Uniform-97 and Global-95 models have slightly worse average drops (around 2%). Remarkably, due to higher initial accuracy, the Uniform-97 model has similar accuracy to Global-95, but much higher end-to-end speedup of 1.75x.

**Case Study: 95% Uniformly-Pruned ResNet18.** We now analyze in detail both the accuracy drops and the per-layer and global speedups for a 95% uniformly-pruned ResNet18 model. On ImageNet, the AC/DC pruned model ResNet18 model has 68.9% Top-1 accuracy, a relative 1% drop from the Torchvision baseline. Figure 7 depicts the transfer performance of the respective models on twelve target datasets, using exactly the same transfer recipe for both sparse and dense models. The accuracy loss is moderate (2.85% average, with a maximum gap of 4.5% Top-1, across all tasks).

Figure 5: Accuracy vs. speedup for _transfer learning_ experiments on the ResNet50 architecture. The boxplots represent aggregated performance across all twelve target tasks.

Figure 6: Accuracy vs. estimated speedup for _from-scratch learning_ experiments on the ResNet18 architecture. The boxplots represent aggregated performance across all twelve target tasks.
Figure 8 depicts layer-wise backward speedups with respect to Pytorch's dense implementation. The results show that the overall end-to-end speedup is 1.58x, with a 1.26x speedup end-to-end for forward computations and a 2.11x speedup end-to-end for backward computations. A closer examination reveals that if we only measure the time spent in convolution and linear modules' forward and backward functions we get a 2.53x speedup, suggesting the presence of significant overheads outside of the convolution and linear computations (such as batch normalization layers and ReLU) for the Pytorch CPU implementation. More precisely, our implementations provide speedups of 1.57x and 3.60x, for the forward and backward multiplications, respectively. This highlights the efficiency of our backward algorithms, but also the overheads of Pytorch's current CPU training implementation. (Specifically, in our sparse implementation, the batch normalization and ReLU computations, executed via Pytorch, take approximately 25% of the total training time.)

**Language Modelling.** Next, we replicate the setup of Kurtic et al. (2022), where a BERT-base (Devlin et al., 2019) model is pruned in the pre-training stage on BookCorpus and English Wikipedia (Lhoest et al., 2021) with the state-of-the-art unstructured pruner oBERT (Kurtic et al., 2022). After that, the remaining weights are fine-tuned on several downstream tasks with fixed sparsity masks. We consider a model with 97% global sparsity, and maintain its masks throughout finetuning on the MNLI and QQP tasks from the GLUE benchmark (Wang et al., 2018). Both accuracy and speedup results are shown in Table 3, and show 37% speedup on a single core for inference on this model, at the price of \(\sim\) 1-3.5% accuracy drop.

#### 4.2.2 Application 2: Sparse Training

**Image Classification.** Finally, we evaluate SparseProp in the from-scratch sparse training scenario. We first consider the same 12 specialized datasets as in Section 4.2.1, and train a ResNet18 architecture _from scratch_. We apply 90% and 95% sparsity using Gradual Magnitude Pruning (Zhu and Gupta, 2017), in two scenarios--Uniform sparsity, in which all convolutional layers, except the first, are pruned to the same target, and the Global scenario, in which the magnitudes of the weights are aggregated across all layers, and the smallest-magnitude weights are pruned to reach the target sparsity. (For the latter scenario, we consider only 95% sparse models.) Note that in the Uniform scenario, we do not prune the initial convolution, nor the final FC layer, while in the Global scenario, we do. In both cases, we train the model dense for ten epochs before pruning the lowest 5% of the weights (either uniformly or globally) and then prune gradually every 10 epochs until epoch 80, at which point we fine-tune for a further 20 epochs.

The results are given in Figure 6. In terms of accuracy, Global-GMP outperforms Uniform-GMP at the same target sparsity, with Global-95% showing similar accuracy to Uniform-90%, though in all cases performance is inferior to when ImageNet weights are used for pre-training. The highest speedup is of 1.31x, for 95% Uniform sparsity.

Additionally, we evaluate SparseProp on the CelebA dataset (Liu et al., 2015), which consists of a training set

\begin{table}
\begin{tabular}{l|l l l|l l l} \hline \hline Backtrack & Forward & Backward & End-to-End & Normal & Backward & End-to-End \\ \hline Dense & \(18.0\pm 0.14\) & \(18.6\pm 0.16\) & \(28.18\pm 0.20\) & \(1.14\pm 0.00\) & \(3.04\pm 0.00\) & \(3.00\pm 0.00\) \\ Uniform 90\% & 1.11\% & 1.31\% & 1.10\% & 1.00\% & 1.36\% & 1.16\% \\ Global 95\% & 1.26\% & 1.28\% & 1.27\% & 1.12\% & 1.26\% & 1.15\% \\ Uniform 95\% & 1.30\% & 1.86\% & 1.76\% & 1.40\% & 1.36\% & 1.80\% \\ Global 95\% & 1.00\% & 1.71\% & 1.86\% & 1.92\% & 1.63\% & 1.46\% \\ \hline \hline \end{tabular}
\end{table}
Table 2: Single-core (left) and parallel (right) relative speedup for transfer learning on sparse ResNet50 models pretrained on ImageNet1k over a dense implementation.

Figure 8: Layer-wise back-propagation time comparison between Dense and 95% uniformly-pruned Resnet18. Note that first and last layers are always dense and are hence removed from the comparison.

\begin{table}
\begin{tabular}{l|l l|l l l} \hline \hline BERT-base & Forward & Backward & End-to-End & MNLI & QQP \\ \hline Dense & \(0.72\pm 0.01s\) & \(1.38\pm 0.00s\) & \(\mathbf{2.67s\pm 0.01s}\) & \(84.54\%\) & \(91.06\%\) \\ Global 97\% & \(1.27\times\) & \(1.40\times\) & \(\mathbf{1.27\times}\) & \(80.91\%\) & \(90.33\%\) \\ \hline \hline \end{tabular}
\end{table}
Table 3: Accuracies and single-core relative speedup for transfer learning sparse BERT-base models.

Figure 7: Top-1 validation transfer accuracy of dense and 95% sparse ResNet18 models pre-trained on ImageNet-1K.
of 162'770 images, and a validation set of 19'962 images from 10 000 well-known individuals, each annotated with forty binary attributes, such as "Smiling", "Male", "Wearing Necklace", etc. We consider the task of jointly predicting all forty attributes, in a similar setup as above, and show the resulting AUC in Table 4. We observe that AUC stays fairly constant even at high sparsities, even as the speed of training increases.

**Language Modelling.** For sparse fine-tuning from scratch on language models, we start from the pre-trained BERTbase (Devlin et al., 2019) model, which we fine-tune for 3 epochs on the target downstream task, and then prune in _one-shot_ with the state-of-the-art unstructured pruner oBERT (Kurtic et al., 2022), uniformly to 90%, 95% or 97% per-layer sparsity. After one-shot pruning, we fine-tune the remaining weights for 5 epochs and examine accuracy and speedup versus the dense variant. The results are presented in Table 4 (right), and show end-to-end speedups of up to 30%, at an accuracy loss between 1 and 3.5%.

## 5 Discussion

We have provided an efficient vectorized algorithm for sparse backpropagation, with linear runtime dependency in the density of the layer weights. We have also provided an efficient CPU-based implementation of this algorithm, and integrated it with the popular Pytorch framework. Experimental evidence validates the runtime scaling of our algorithm on various layer shapes and types. We complemented this algorithmic contribution with an extensive study of the feasibility of sparse transfer learning and from-scratch training in edge scenarios. We observed consistent speedups across scenarios, at the cost of moderate accuracy loss. Our results should serve as motivation for further research into accurate sparse training in this setting, in particular for leveraging sparsity on highly-specialized tasks, which is an under-studied area.
