Sliding Window Sum Algorithms for Deep Neural Networks
Sliding window sums are widely used for string indexing, hashing and time
series analysis. We have developed a family of the generic vectorized sliding
sum algorithms that provide speedup of O(P/w) for window size $w$ and number of
processors P. For a sum with a commutative operator the speedup is improved to
O(P/log(w)). Even more important, our algorithms exhibit efficient memory
access patterns. In this paper we study the application of the sliding sum
algorithms to the training and inference of the Deep Neural Networks. We
demonstrate how both pooling and convolution primitives could be expressed as
sliding sums and evaluated by the compute kernels with the shared structure. We
show that the sliding sum convolution kernels are more efficient than the
commonly used GEMM kernels on the CPU, and could even outperform their GPU
counterparts.
Roman Snytsar
# Sliding Window Sum Algorithms for Deep Neural Networks

###### Abstract

Sliding window sums are widely used for string indexing, hashing and time series analysis. We have developed a family of the generic vectorized sliding sum algorithms that provide speedup of \(O(P/w)\) for window size \(w\) and number of processors \(P\). For a sum with a commutative operator the speedup is improved to \(O(P/log(w))\). Even more important, our algorithms exhibit efficient memory access patterns.

In this paper we study the application of the sliding sum algorithms to the training and inference of the Deep Neural Networks. We demonstrate how both pooling and convolution primitives could be expressed as sliding sums and evaluated by the compute kernels with the shared structure.

We show that the sliding sum convolution kernels are more efficient than the commonly used GEMM kernels on the CPU, and could even outperform their GPU counterparts.

## 1 Introduction

A Deep Neural Network (DNN) is one of the most significant tools in the arsenal of the machine learning (ML) researcher [9]. DNNs are constructed from multiple layers that transform the data sequentially via operations such as pooling, convolution, and activation. In most successful DNNs the great portion of computational resources is consumed performing convolution.

A common approach to implementing convolutional layers is to expand the input into a column matrix (_im2col_) and then call a highly tuned General Matrix Multiplication (GEMM) procedure from the existing linear algebra library such as BLIS [13] or MKL [15]. Since the hardware optimized GEMM implementations exist for every standard CPU, graphics processing unit (GPU), or digital signal processor (DSP), the _im2col_ approach has been highly successful in DNN frameworks such as Caffe [8], ONNX [2] and Torch [5]. The popularity of the _im2col_ tactics also influences the design of the custom hardware marketed as ML accelerators that are in fact the GEMM accelerators.

The major disadvantages of the _im2col_ conversion are the greatly increased memory footprint of the input matrix and the reduced data locality. For a convolution with a filter size k, the column matrix is k times larger than the original input. A lot of effort has been put into remediating this problem. A mainstream approach is converting input to the low precision floating point or even integer representation [6]. The quantization techniques reduce the memory footprint and latency by an order of magnitude and even influence the design of the GEMM accelerators. It is important to note though that the quantization is not entangled with GEMM and could be equally successfully applied to the original convolution problem. Yet another research trend is applying the GEMM routines to the smaller intermediate data structures [1][14]or even to the original input data [16].

We propose an approach that replaces GEMM with a new kind of the computation kernel that operates on the unmodified input.
## 2 Methods

### Prefix Sum

At the foundation of our method is the concept of a _prefix sum_, and the accompanying _reduce_ and _scan_ algorithms. A prefix sum is a transformation that takes an operator \(\oplus\), and a sequence of elements

\[x_{0},x_{1},\ldots,x_{k},\ldots\]

and returns the sequence

\[y_{i}=\sum_{j=0}^{i}x_{j}=x_{0}\oplus x_{1}\oplus\ldots\oplus x_{i} \tag{1}\]

or in recurrent form

\[y_{i+1}=y_{i}\oplus x_{i+1} \tag{2}\]

Despite the data carry dependency, the \(Nth\) element of the prefix sum with an associative operator could be computed in \(O(log(N))\) parallel steps using _reduce_ algorithm. Even stronger statement is that all \(N\) elements of the prefix sum could be computed in the same \(O(log(N))\) parallel steps using _scan_ algorithm, as shown by [3].

### Sliding Window Sum

Sliding window sum (sliding sum) takes a window size \(w\) in addition to an operator \(\oplus\), and a sequence of elements, and returns the sequence

\[y_{i}=\sum_{j=i}^{i+w-1}x_{j}=x_{i}\oplus x_{i+1}\oplus\ldots\oplus x_{i+w-1} \tag{3}\]

where each sum is defined in terms of the operator \(\oplus\) and contains exactly \(w\) addends. The asymptotic complexity of a naive sliding sum algorithm is \(O(wN)\) where \(N\) is the length of the source sequence.

Every sum defined by Equation 3 is a prefix sum with operator \(\oplus\) and input sequence \(x_{i}\ldots\oplus x_{i+w-1}\). Many useful operators are associative, so the prefix scan algorithm is applicable here, reducing complexity of every sum in Equation 3 to \(O(log(w))\) and, trivially, the overall sliding sum complexity to \(O(Nlog(w))\) parallel steps.

While working on the bioinformatics applications, we have used sliding window sums to represent the minimizer seeds and have developed a family of algorithms to achieve better parallel speedups [11]. Now we will apply the same approach to the DNN operators.

### Pooling

The average pooling operator is trivially the sliding window sum with the associative operator \(+\). By analogy, the max pooling operator is a sliding window sum with the associative operator \(max\). Implementing the sliding pooling operator could be a warm-up before concentrating on the convolution.

### Dot Product

As a corollary we will show that a dot product is a prefix sum. Dot product of the two vectors of length \(M\) is defined as:

\[c=\sum_{i=0}^{M-1}a_{i}b_{i} \tag{4}\]

First, we replace vectors \(a\) and \(b\) with vectors \(\alpha\) and \(\beta\) so that:

\[\alpha_{i}=\begin{cases}1,&a_{i}=0\\ a_{i},&otherwise\end{cases},\beta_{i}=\begin{cases}0,&a_{i}=0\\ b_{i},&otherwise\end{cases} \tag{5}\]

It holds that

\[\sum_{i=0}^{M-1}\alpha_{i}\beta_{i}=\sum_{i=0}^{M-1}a_{i}b_{i} \tag{6}\]
Next, we define a sequence of \((M+1)\) pairs,

\[\gamma_{i}=\begin{pmatrix}u_{i}\\ v_{i}\end{pmatrix},\quad where\quad u_{i}=\begin{cases}1,&i=0\\ \frac{\alpha_{i-1}}{\alpha_{i}},&0<i<M\;,\quad and\quad v_{i}=\begin{cases}\beta _{i},&i<M\\ 0,&i=M\end{cases} \tag{7}\]

and the operator \(\oplus\) such that

\[\gamma_{i}\oplus\gamma_{j}=\begin{pmatrix}u_{i}\\ v_{i}\end{pmatrix}\oplus\begin{pmatrix}u_{j}\\ v_{j}\end{pmatrix}=\begin{pmatrix}u_{i}\bullet u_{j}\\ u_{j}\bullet v_{i}+v_{j}\end{pmatrix} \tag{8}\]

Operator \(\oplus\) is associative [3], and the sequence

\[\delta_{i}=\begin{cases}\gamma_{0},&i=0\\ \delta_{i-1}\oplus\gamma_{i},&0<i\leq M\end{cases} \tag{9}\]

is a prefix sum. The bottom element of the last sum \(\delta_{M}\) is the original dot product, and \(\delta_{M}\) could be evaluated using reduce algorithm in \(log(M)\) parallel steps of fused multiply-add (FMA) operations. The total work is still \(M\). The important result here is representing dot product as a prefix sum in preparations for the sliding window implementation of the convolution.

### Convolution

Convolution is a sliding window sum (dot product) with the associative operator defined by equation 8. Consequently our family of the sliding window algorithms is applicable to the convolution operator.

## 3 Algorithms

### Vector Algorithms

Our first algorithm is a vector-friendly way of calculating sliding sum assuming the input sequence elements become available one by one and are processed using the vector instructions of width \(P>w\):

```
procedureScalarInput(\(x_{0}\dots x_{n-1}\)) \(Y\leftarrow\Big{(}\underbrace{w^{-2}}_{j=0}x_{j},\sum_{j=1}^{ w-2}x_{j},\dots,x_{w-3}\oplus x_{w-2},x_{w-2},\dots,0\Big{)}\) for\(i=w-1\) to N do \(X\leftarrow\Big{(}\underbrace{x_{i},x_{i},\dots,x_{j}}_{w},0,\dots,0\Big{)}\) \(Y\gets Y\oplus X\) \(y_{i-w+1}\gets Y[0]\) \(Y\gets Y\lll 1\) endprocedure
```

**Algorithm 1** Scalar Input

Vector Y is initialized to the suffix sums with the number of elements decreasing from \(w-1\) to \(0\). Then in a loop every incoming element \(x_{k}\) is broadcast to the first \(w\) elements of vector X. After vector addition the zeroth element of Y contains the next sliding sum. Next, the vector Y is shifted left by one element, as denoted by operator \(\lll\), and the state is ready for the next iteration.

Asymptotic complexity of the scalar input algorithm is \(O(N)\) with no additional requirements on the operator \(\oplus\).

This result could be improved if we assume that the input sequence arrives packed in vectors of width \(P>w\).

At every iteration \(P\) input elements are placed into vector \(X\). X1 is filled with the prefix sums of up to \(w\) addends, and Y1 is filled with the suffix sums constructed from the elements of \(X\). Then the vector sum of \(Y\) and \(X1\) yields the
next P output elements. Finally, the suffix sums from \(Y1\) are shifted into proper positions in vector \(Y\), and it is ready for the next iteration.

The asymptotic complexity thus is \(O(N\cdot w/P)\) with the parallel speedup \(O(P/w)\) for any operator \(\oplus\). If \(\oplus\) is associative, the prefix/suffix sums could be computed in parallel using the algorithm in [3], and the complexity is reduced to \(O(N\cdot log(w)/P)\) with the speedup improving to \(O(P/log(w))\).

For example, since \(min\) is an associative operator, the sliding window minimum can be computed using the faster version of the vector input algorithm.

One might notice that the suffix sum computation utilises only \(w-1\) out of \(p\) elements of vector \(Y1\). Eliminating this inefficiency leads to the Ping Pong algorithm where both suffix and prefix sums yield some elements of the output sequence.

The Ping Pong algorithm does not offer any asymptotic improvements but is 30-50% faster in practice. It however accesses memory in strides unaligned to \(P\), and the two memory loads per iteration present a challenge while implementing boundary conditions like padding, mirroring, or periodicity.

A simpler algorithm could be formulated in terms of the Slide operation that extracts a vector from the concatenation of the two inputs. It directly maps to the **EXT** instruction from the ARM SVE instruction set [12], and is easily implemented via appropriately named **vsilideup/down** instructions of the RISC-V vector extension [10], or the **vperm*2ps** Intel AVX512 instructions [7]. Similar to the previous algorithms, if \(\oplus\) is associative, the inner loop could be replaced by the parallel reduction for maximum parallel speedup.

## 4 Experiments

We have implemented the sliding window convolution as an alternate execution path in the ONNX framework [2]. Figure 1 shows the achieved speedup when compared to the baseline --emphMlasConv procedure applied to the large 1-D input and the filters of various sizes. The speedup is approximately proportional to the logarithm of the kernel size.

Additionally, we have recreated the scenario for the very large dilated kernel described in [4]. We have achieved up to 6.8x speedup on the small data set and around 4x speedup across the board approaching the GPU performance numbers.
```
procedurePingPong(\(x_{0}\dots x_{N-1}\)) for\(i=0\) to N step 2P-w+1 do \(Y\leftarrow\left(x_{i},x_{i+1},\dots,x_{i+p-1}\right)\) \(X\leftarrow\left(x_{i+P},x_{i+P+1},\dots,x_{i+2P-1}\right)\) \(Y1\leftarrow\Big{(}\underbrace{w-1}_{j=0}^{w-1}Y[j],\sum_{j=1}^{w}Y[j],\dots, \sum_{j=P-w}^{P-1}Y[j],\sum_{j=P-w}^{P-1}Y[j],\dots,Y_{P-1}\Big{)}\) \(y_{i}\dots y_{i+P-w}\gets Y[0]\dots Y[P-w]\) \(Y1\gets Y1\lll(P-w)\) \(X1\leftarrow\Big{(}\underbrace{X[0],X[0]\oplus X[1],\dots,\sum_{j=0}^{w-2}X[j],\sum_{j=0}^{w-1}X[j],\dots,\sum_{j=p-w}^{P-1}X[j]\Big{)}\) \(Y1\gets Y1\oplus X1\) \(y_{i+P-w+1}\dots y_{i+2P-w}\gets Y[0]\dots Y[P-1]\) endfor endprocedure
```

**Algorithm 3** Ping Pong

```
procedureSlide(\(Y1,Y2,offset\)) \(Y\gets CONCAT(Y1,Y2)\) return\(\big{(}\underbrace{Y[offset],Y[offset+1],\dots,Y[offset+P-1]}_{P}\big{)}\) endprocedureprocedureVectorSlide(\(x_{0}\dots x_{N-1}\)) \(Y\leftarrow(x_{0},x_{1},\dots,x_{P-1})\) for\(i=P\) to \(N\) step \(P\)do \(Y1\leftarrow\Big{(}x_{i},x_{i+1},\dots,x_{i+p-1}\Big{)}\) \(X\gets Y1\) for\(k=1\) to \(w-1\)do \(X\gets X\oplus\)Slide(\(Y,Y1,P-k\)) endfor \(y_{i}\dots y_{i+p-1}\gets X[0]\dots X[p-1]\) \(Y\gets Y1\) endfor endprocedure
```

**Algorithm 4** Vector Slide
## 5 Conclusion

We have demonstrated excellent performance of the Sliding Sum algorithms using the commodity hardware without the need for the specialized accelerators. Despite the promising results, there is plenty of work ahead.

The most obvious next step is extending the sliding convolution approach to more than one dimension covering the majority of the DNN applications.

The most common filter sizes in the DNN applications are 3 an 5 in every dimension. With the filter this small the current sliding convolution algorithms demonstrate very modest speedup since the number of the arithmetic instructions per the memory load is low. The situation improves in the multiple dimensions but still could require custom compute kernels for the small filter sizes.

Lastly, since the accelerators for the matrix multiplication are already present in the current generation of the hardware, it would be wise re-using them. Thus, it is important to re-formulate our algorithms in terms of the small matrix multiplication completing the circle.
