Reformulating the direct convolution for high-performance deep learning inference on ARM processors

Journal of Systems Architecture 135 (2023) 102806 Available online 15 December 2022 1383-7621/© 2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/). Contents lists available at ScienceDirect Journal of Systems Architecture journal homepage: www.elsevier.com/locate/sysarc Reformulating the direct convolution for high-performance deep learning inference on ARM processors Sergio Barrachina a, Adrián Castelló b, Manuel F. Dolz a,∗, Tze Meng Lowc, Héctor Martínez d, Enrique S. Quintana-Ortí b, Upasana Sridhar c, Andrés E. Tomása aUniversitat Jaume I, Castellón de la Plana, Spain bUniversitat Politècnica de València, Valencia, Spain cCarnegie Mellon University, Pittsburgh, USA dUniversidad de Córdoba, Córdoba, Spain ARTICLE INFO Keywords: Convolution Direct algorithm Deep learning High performance ARMv8 architecture ABSTRACT We present two high-performance implementations of the convolution operator via the direct algorithm that outperform the so-called lowering approach based on the im2col transform plus the gemm kernel on an ARMv8based processor. One of our methods presents the additional advantage of zero-memory overhead while the other employs an additional yet rather moderate workspace, substantially smaller than that required by the im2col+gemm solution. In contrast with a previous implementation of a similar zero-memory overhead direct convolution, this work exhibits the key advantage of preserving the conventional NHWC data layout for the input/output activations of the convolution layers. 1. Introduction The convolution (operator) is a key computational kernel that concentrates a significant part of the arithmetic cost for the type of deep neural networks (DNNs) that are frequently leveraged in signal processing and computer vision tasks [1,2]. Therefore, it is natural that a considerable effort has been dedicated to the efficient implementation of this operator. The most common implementation is the lowering (or im2col-based) approach which transforms the input tensor into an augmented matrix in order to cast the operator in terms of a cache-friendly general matrix–matrix (gemm) multiplication [3–5]. This approach has the advantage of exploiting the excellent performance of optimized matrix product implementations, but it requires a very large amount of memory. In contrast, the direct algorithm, consisting of 6–7 nested loops around a multiply-and-accumulate instruction, provides arguably the simplest implementation while requiring much less extra memory [2]. In this paper, we initially follow the ideas in [6] to design a cachefriendly, gemm-like formulation of the direct convolution, making the following new contributions1with respect to that work: •In contrast with [7], our solution preserves the standard NHWC layout for the input tensor and only requires specialized packing for the filter tensor. To achieve this, (1) we decouple the dimension of the microkernel from the cache configuration parameters, ∗Corresponding author. E-mail address: [email protected] (M.F. Dolz). 1The source code of the direct convolution variants presented in this paper is available at https://github.com/hpca-uji/convDirect. which has the positive side effect of improving performance; and (2) we integrate a row-wise oriented microkernel, similar to that introduced in [8]. •In addition, we propose an alternative blocked variant which, during the execution of the convolution operator, packs a small piece of the input tensor into a buffer, in order to accesses with ensure unit stride memory to both the filter and input tensor during the microkernel execution. This approach permits the use of the architecture-specific microkernels for gemm in the BLIS library. •Third, we conduct a complete performance analysis of the two new algorithms on an ARMv8 core for two relevant convolutional DNNs (GoogleLeNet and ResNet-50 v1.5), combined with the ImageNet dataset, showing their advantages in comparison with high-performance implementations of the convolution operator based on the im2col approach and the Winograd transform. Concretely, our experiments with GoogleNet and ResNet-50 show that new algorithms are 13% to 85% faster than the traditional im2col-based approach. The rest of the paper is structured as follows. Section 2provides a brief summary of different convolution algorithms and optimization approaches for the direct convolution variant. In Section 3we review https://doi.org/10.1016/j.sysarc.2022.102806 Received 13 May 2022; Received in revised form 25 November 2022; Accepted 12 December 2022 Journal of Systems Architecture 135 (2023) 102806 2 S. Barrachina et al. the basics of the convolution operator and several scalar and blocked algorithms for its computation, including the high-performance implementation proposed in [7]. Next, in Section 4, we present our new two NHWC-preserving algorithms for the direct convolution, emphasizing the differences with respect to previous solutions. Finally, in Section 5 we report the results from a complete evaluation, and in Section 6we close the paper with a few concluding remarks. 2. Related work Among the different methods proposed in the literature for the convolution operator, we can list (1) the direct algorithm, usually implemented as six nested loops around a multiply-and-add instruction [7]; (2) the lowering approach, which applies either im2col or im2row to the activation input matrix to then obtain the result from a gemm [3,4]; (3) the FFT-based algorithm, which shifts the computation into the frequency domain in order to reduce the arithmetic requirements [9–11]; and (4) the Winograd-based convolution, which leverages the Winograd minimal filtering algorithm to decrease the arithmetic cost of the convolution [9,12]. With regards to the direct algorithm for the convolution, the implementation by Zhang et al. [7] proposes a reorganization of the algorithm loops around a microkernel, mimicking the conventional implementation of high-performance instances of gemm in libraries such as GotoBLAS, OpenBLAS, and BLIS [13–15]. That work also includes a specialized, cache-friendly data layout for both the input and filter tensors to operate with standard data layouts (e.g., NCHW or NHWC). For that, however, the input/output tensors have to be transformed to/from that format before the first DNN layer and after the last one, respectively. This need for a specialized layout may offer a reason why the direct convolution algorithm is not a widely-adopted. To deal with that, in this work we refine this algorithm to preserve the NHWC input/output tensor format. In this line, we can also find optimization approaches for the direct convolution algorithm based on exploiting data reuse and the design of near I/O-optimal data layout strategies. For instance, the work by Zhang et al. [16], applies an aggressive design of auto-tuning based on I/O lower bounds which outperform the optimal solutions by TVM. Similarly, the work by Li et al. [17] introduces a block hardwarefriendly direct convolution implementation that can completely avoid the off-chip memory transfers of intermediate feature maps on memorylimited devices, e.g., FPGAs (field programmable gate arrays) or MCUs (micro-controller units). Other trends for deriving optimized direct convolution implementations use meta-programming techniques, such as those proposed by Zlateski et al. in [18], or auto-tuning compilation stacks, such as Apache TVM [19]. For instance, the work by Liao et al. [20] uses TVM to analyse the sparsification opportunities for the acceleration of the direct convolution by pruning certain some CNNs weights. In a similar line, the work by Georganas et al. [21] proposes JIT-based optimization approaches for the direct convolution kernels on x86 architectures. 3. Algorithms for the direct convolution 3.1. Simple implementation of the direct convolution Let us consider the convolution operator 𝑂=Conv(𝐹 , 𝐼),(1) where 𝐼is a 4D input tensor consisting of 𝑁input images of size 𝐻𝑖×𝑊𝑖×𝐶𝑖each, 𝐻𝑖×𝑊𝑖denote the image height ×width, and 𝐶𝑖 stands for its number of input channels (or features). Furthermore, 𝐹is a 4D input tensor consisting of 𝐶𝑜filters of height ×width ×channels =𝐻𝑓×𝑊𝑓×𝐶𝑖each. The operator in (1) then produces a 4D output tensor 𝑂, with 𝑁outputs of size 𝐻𝑜×𝑊𝑜×𝐶𝑜each, where 𝐻𝑜×𝑊𝑜stand for the output height ×width, and 𝐶𝑜for the number of output channels 1void ConvDirect_Naive( 𝐼[𝑁][𝐻𝑖][𝑊𝑖][𝐶𝑖],, 2𝐹[𝐻𝑓][𝑊𝑓][𝐶𝑖][𝐶𝑜], 3𝑂[𝑁][𝐻𝑜][𝑊𝑜][𝐶𝑜]) 4{ 5for (ℎ= 0;ℎ<𝑁;ℎ++ ) 6for (𝑖= 0;𝑖<𝐶𝑖;𝑖++ ) 7for (𝑗= 0;𝑗<𝐶𝑜;𝑗++ ) 8for (𝑘= 0;𝑘<𝑊𝑜;𝑘++ ) 9for (𝑙= 0;𝑙<𝐻𝑜;𝑙++ ) 10 for (𝑚= 0;𝑚<𝑊𝑓;𝑚++ ) 11 for (𝑛= 0;𝑛<𝐻𝑓;𝑛++ ) 12 𝑂[ℎ][𝑙][𝑘][𝑗] + = 𝐼[ℎ][𝑙+𝑛][𝑘+𝑚][𝑖] ∗ 𝐹[𝑛][𝑚][𝑖][𝑗]; 13 } Listing 1: Naive algorithm for the direct convolution. (which, as stated previously, is equal to the number of filters). For this purpose, each of the 𝐶𝑜individual filters combines (or convolves) a subtensor of the inputs, with the same dimension as the filter, to render a single scalar value (entry) for one of the 𝐶𝑜outputs. By repeatedly applying the filter to the whole input, in a sliding window manner, the convolution operator produces the complete entries of this single output; see [2] and Fig. 1. For simplicity, hereafter we will consider that the filter is applied with unit vertical/horizontal strides; and the output is not padded so that 𝐻𝑜=𝐻𝑖−𝐻𝑓+ 1, 𝑊𝑜=𝑊𝑖−𝑊𝑓+ 1. A naive algorithm that computes the convolution is given in Listing 1. Overall, the algorithm traverses the dimensions of the convolution operator in the order 𝑁→𝐶𝑖→𝐶𝑜→𝑊𝑜(𝑊𝑖)→𝐻𝑜(𝐻𝑖)→𝑊𝑓→𝐻𝑓, from the outermost to the innermost loop, which results in a certain pattern of memory accesses depending on the layout of the tensors in memory. Here we will assume that the operands are laid out following a generalization of the row-major matrix format for the multidimensional tensors, which implies that the dimensions are closer in memory from right to left (closest to farthest; see Fig. 2). Furthermore, we will adopt the notation 𝐼mem[𝑁][𝐻𝑖][𝑊𝑖][𝐶𝑖], 𝑂mem[𝑁][𝐻𝑜][𝑊𝑜][𝐶𝑜]. to specify this layout in memory for the 4D input/output tensors. Note that, this corresponds to the standard NHWC format for 4D tensors in deep learning (DL). Also, the 4D filter tensor is stored following the standard RSCK format; that is, 𝐹mem[𝐻𝑓][𝑊𝑓][𝐶𝑖][𝐶𝑜]. We close this short review of the direct convolution by highlighting that the 7 loops in the algorithm are independent of each other as well as from the memory layout of the input, output and filter tensors. In consequence, the algorithm loops can be reorganized in any order while still producing the correct result. 3.2. Blocking the direct convolution In [7], the authors propose a complete reorganization of the naive algorithm for the direct convolution in Listing 1with the enhancements described in the following paragraphs. Reordering the loops to saturate computations. The two innermost loops of the algorithms are chosen to traverse the 𝐶𝑜and 𝑊𝑜dimensions of the convolution operator (from innermost to outermost) in order to: 1. Favor that the processor floating-point arithmetic units (FPUs) experience no stalls due to data dependencies between consecutive writes to the same entry of the output tensor; and 2. Accommodate SIMD (single instruction, multiple data) vectorization in the innermost loop by exploiting the number of channel outputs, 𝐶𝑜. (As 𝐶𝑜is a user-specified parameter, it can be selected to be a multiple of the size of the SIMD vector.) Journal of Systems Architecture 135 (2023) 102806 3 S. Barrachina et al. Fig. 1. Sliding window of 𝐶𝑜filters over an input image producing the output of the convolution operator for each output channel. Fig. 2. Memory layout for the 4D output tensor 𝑂assuming 𝑁= 1 and a generalization of the conventional row-major layout for matrices and the NHWC format for the tensor. The same layout applies to the 3D input tensor 𝐼. Fig. 3. Convolution-friendly memory layout for the 4D output tensor proposed in [7], assuming 𝑁= 1. 1void ConvDirect_Reordered( 𝐼[𝑁][𝐻𝑖][𝑊𝑖][𝐶𝑖],, 2𝐹[𝐻𝑓][𝑊𝑓][𝐶𝑖][𝐶𝑜], 3𝑂[𝑁][𝐻𝑜][𝑊𝑜][𝐶𝑜]) 4{ 5for (ℎ= 0;ℎ<𝑁;ℎ++ ) 6for (𝑙= 0;𝑙<𝐻𝑜;𝑙++ ) 7for (𝑛= 0;𝑛<𝐻𝑓;𝑛++ ) 8for (𝑚= 0;𝑚<𝑊𝑓;𝑚++ ) 9for (𝑖= 0;𝑖<𝐶𝑖;𝑖++ ) 10 for (𝑘= 0;𝑘<𝑊𝑜;𝑘++ ) 11 for (𝑗= 0;𝑗<𝐶𝑜;𝑗++ ) 12 𝑂[ℎ][𝑙][𝑘][𝑗] + = 𝐼[ℎ][𝑙+𝑛][𝑘+𝑚][𝑖] ∗ 𝐹[𝑛][𝑚][𝑖][𝑗]; 13 } Listing 2: Algorithm for the direct convolution with re-ordered loops. Ordering the loops to optimize data reuse. After having reordered the two innermost loops, the next three loops in the sequence from innermost to outermost are set to traverse the 𝐶𝑖, 𝑊𝑓, 𝐻𝑓dimensions of the operator so as to ensure that the data from the input and output tensors are accessed in the same order. For convolutional DNNs, this is beneficial because the output of one convolution layer is often the input to the next one. The initial algorithm with the reordered loops in [7] is reproduced in Listing 2. Blocking for the memory hierarchy. In this stage: 1. Due to the limited number of vector registers in any processor architecture, tiling [22] with block sizes 𝐶𝑜,𝑏, 𝑊𝑜,𝑏 is applied to the two innermost loops of the reordered algorithm. This is performed to avoid register spilling that can result in degraded performance. 2. Furthermore, the loop traversing the 𝑊𝑜∕𝑊𝑜,𝑏 dimension is placed between the 𝐻𝑜and 𝐻𝑓loops in order to mimic the ordering of the loops traversing the filter weights. 3. Also, the loop traversing the 𝐶𝑜∕𝐶𝑜,𝑏 dimension is set as the outermost loop to facilitate parallelization. 4. In addition, the 𝐶𝑖dimension of the operator is partitioned with block size 𝐶𝑖,𝑏 so as to split the input dataset into smaller blocks that fit into the memory hierarchy. The resulting blocked algorithm in [7] is illustrated in Listing 3. Convolution-friendly data layout. To ensure that the output of the convolution operator is accessed with unit stride: 1. The output tensor is rearranged in memory into sequential blocks of 𝑁×𝐻𝑜×𝑊𝑜×𝐶𝑜,𝑏 elements, where the entries of each block are first arranged in the channel dimension, and then following the row-major order into 𝐻𝑜×𝑊𝑜blocks of length 𝐶𝑜,𝑏. 2. For compatibility between the outputs and inputs of consecutive layers in convolutional DNNs, the same layout applies to the input tensor; see Fig. 3. In our notation, this corresponds to: 𝐼mem[𝑁][𝐶𝑖∕𝐶𝑖,𝑏][𝐻𝑖][𝑊𝑖][𝐶𝑖,𝑏], 𝑂mem[𝑁][𝐶𝑜∕𝐶𝑜,𝑏][𝐻𝑜][𝑊𝑜][𝐶𝑜,𝑏]. 3. Finally, the filter tensor is rearranged in memory with its dimensions organized as [7]: 𝐹mem[𝐶𝑜∕𝐶𝑜,𝑏][𝐶𝑖∕𝐶𝑖,𝑏][𝐻𝑓][𝑊𝑓][𝐶𝑖,𝑏][𝐶𝑜,𝑏]. Journal of Systems Architecture 135 (2023) 102806 4 S. Barrachina et al. 1void ConvDirect_Blocked( 𝐼[𝑁][𝐻𝑖][𝑊𝑖][𝐶𝑖], 2𝐹[𝐻𝑓][𝑊𝑓][𝐶𝑖][𝐶𝑜], 3𝑂[𝑁][𝐻𝑜][𝑊𝑜][𝐶𝑜]) 4{ 5for (ℎ= 0;ℎ<𝑁;ℎ++ ) 6for (𝑗′= 0;𝑗′<𝐶𝑜∕𝐶𝑜,𝑏;𝑗′++ ) 7for (𝑖′= 0;𝑖′<𝐶𝑖∕𝐶𝑖,𝑏;𝑖′++ ) 8for (𝑙= 0;𝑙<𝐻𝑜;𝑙++ ) 9for (𝑘′= 0;𝑘′<𝑊𝑜∕𝑊𝑜,𝑏;𝑘′++ ) 10 for (𝑛= 0;𝑛<𝐻𝑓;𝑛++ ) 11 for (𝑚= 0;𝑚<𝑊𝑓;𝑚++ ) 12 // Micro -kernel for C+=A*B 13 for (𝑖𝑖 = 0;𝑖𝑖 <𝐶𝑖,𝑏;𝑖𝑖++ ) 14 for (𝑘𝑘 = 0;𝑘𝑘 <𝑊𝑜,𝑏;𝑘𝑘++ ) 15 for (𝑗𝑗 = 0;𝑗𝑗 <𝐶𝑜,𝑏;𝑗𝑗 ++ ) 16 𝑂[ℎ][𝑙][𝑘′⋅𝑊𝑜,𝑏 +𝑘𝑘][𝑗′⋅𝐶𝑜,𝑏 +𝑗𝑗] 17 + = 𝐼[ℎ][𝑙+𝑛][𝑘′⋅𝑊𝑜,𝑏 +𝑘𝑘 +𝑚][𝑖′⋅𝐶𝑖,𝑏 +𝑖𝑖] 18 ∗𝐹[𝑛][𝑚][𝑗′⋅𝐶𝑜,𝑏 +𝑗𝑗][𝑖′⋅𝐶𝑖,𝑏 +𝑖𝑖] 19 } Listing 3: Blocked algorithm for the direct convolution. Fig. 4. Microkernel for 𝐶+ = 𝐴⋅𝐵. The microtile 𝐶(in brown) is initially retrieved into the processor registers and updated at each the 𝐶𝑖,𝑏 iterations of a loop via an outer product between a column of the micropanel 𝐴(in red) and a column of the micropanel 𝐵(in blue). Upon completion, 𝐶is copied back into the memory. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Gemm-based microkernel. The three innermost loops in the blocked algorithm for the convolution, traversing the 𝐶𝑖,𝑏, 𝑊𝑜,𝑏, 𝐶𝑜,𝑏 dimensions of the operator (see Listing 3), perform a small gemm, of the form 𝐶+ = 𝐴⋅𝐵, where 𝐶is a 𝑊𝑜,𝑏 ×𝐶𝑜,𝑏 microtile of 𝑂;𝐴is a 𝑊𝑜,𝑏 ×𝐶𝑖,𝑏 micropanel of 𝐼; and 𝐵is a 𝐶𝑖,𝑏 ×𝐶𝑜,𝑏 micropanel of 𝐹. Following the usual implementation of this type of microkernels in high-performance instances of the BLAS (basic linear algebra subprograms [23]), this is done by first loading the contents of the microtile 𝐶into the processor (vector) registers, to then update these elements via a sequence of 𝐶𝑖,𝑏 outer products, each involving a single column/row of the micropanels 𝐴∕𝐵consisting of 𝑊𝑜,𝑏∕𝐶𝑜,𝑏 elements; see [15,24] and Fig. 4. For high performance, vector fused-multiply-add (FMAs) instructions are used to perform these updates using the processor SIMD FPUs. The size of the microkernel integrated in the direct convolution is determined based on the latency and throughput of these type of instructions, as well as on the available number of architecture vector registers. Specifically, given a fixed number of SIMD registers available on the target processor, the design of the microkernel seeks to maximize the number of output elements computed while staying within the fixed number of available registers. The layout chosen for the output (𝑂) and filter tensors (𝐹) in [7] ensures that the entries of the microtile 𝐶and the rows of the micropanel 𝐵can be accessed from the microkernel with unit stride. In contrast, the layout of the input (𝐼) implies that the micropanel 𝐴is stored with consecutive columns but with a row stride (or row leading dimension) ldA=𝐶𝑖,𝑏. For the microkernel, this results in non-unit stride accesses 1void ConvDirect_New( 𝐼[𝑁][𝐻𝑖][𝑊𝑖][𝐶𝑖], 2𝐹[𝐻𝑓][𝑊𝑓][𝐶𝑖][𝐶𝑜], 3𝑂[𝑁][𝐻𝑜][𝑊𝑜][𝐶𝑜]) 4{ 5// Loops for h, j’, i’, l, k’, n 6// remain the same as in the blocked 7// algorithm in Listing 3and therefore 8// they are omitted for brevity 9for (𝑚= 0;𝑚<𝑊𝑓;𝑚++ ) 10 for (𝑘𝑘 = 0;𝑘𝑘 <𝑊𝑜,𝑏;𝑘𝑘 + = 𝑀𝑟) 11 for (𝑗𝑗 = 0;𝑗𝑗 <𝐶𝑜,𝑏;𝑗𝑗 + = 𝑁𝑟) 12 // Micro -kernel 13 for (𝑖𝑖 = 0;𝑖𝑖 <𝐶𝑖,𝑏;𝑖𝑖++ ) 14 for (𝑗𝑟=𝑘𝑘;𝑗𝑟<𝑘𝑘 +𝑀𝑟;𝑗𝑟++ ) 15 for (𝑖𝑟=𝐶𝑜,𝑏;𝑖𝑟<𝐶𝑜,𝑏 +𝑁𝑟;𝑖𝑟++ ) 16 𝑂[ℎ][𝑙][𝑘′⋅𝑊𝑜,𝑏 +𝑗𝑟][𝑗′⋅𝐶𝑜,𝑏 +𝑖𝑟] 17 + = 𝐼[ℎ][𝑙+𝑛][𝑘′⋅𝑊𝑜,𝑏 +𝑗𝑟+𝑚][𝑖′⋅𝐶𝑖,𝑏 +𝑖𝑖] 18 ∗𝐹[𝑛][𝑚][𝑗′⋅𝐶𝑜,𝑏 +𝑖𝑟][𝑖′⋅𝐶𝑖,𝑏 +𝑖𝑖] 19 } Listing 4: Blocked algorithm for the direct convolution with decoupled micro-kernel loops. to the columns of the micropanel 𝐴. However, 𝐶𝑖,𝑏 can be chosen to be small enough to ensure that the consecutive elements reside in the same cache line. 4. New blocked NHWC-preserving algorithms for the direct convolution In this paper, we build upon the work in [7] by maintaining part of the blocking strategy proposed there, but making the major changes described in the following paragraphs. Decouple the microkernel dimension from the cache blocking parameters. The blocked algorithm in [7] ties the blocking parameters 𝐶𝑖,𝑏,𝑊𝑜,𝑏,𝐶𝑜,𝑏 to the dimensions of the small gemm performed inside the microkernel. In our work, we avoid this by introducing two additional loops that traverse the 𝑊𝑜,𝑏, 𝐶𝑖,𝑏 dimensions of the operator in steps of size 𝑀𝑟, 𝑁𝑟 respectively, as shown in Listing 4. As a result, the microkernel now performs a small gemm of the form 𝐶+ = 𝐴⋅𝐵, where 𝐶is a microtile of size 𝑀𝑟×𝑁𝑟, and 𝐴, 𝐵 are micropanels of size 𝑀𝑟×𝐶𝑖,𝑏 and 𝐶𝑖,𝑏 ×𝑁𝑟, respectively. Preserving the input/output NHWC layout. In contrast with [7], we preserve the original NHWC layout for the input and output tensors. Accesses with non-unit stride to 𝐶from the microkernel can be expected to be well amortized due to the 𝐶𝑖,𝑏 ratio between the number of memory accesses to load/store the microtile 𝐶(2𝑀𝑟𝑁𝑟memory accesses) and the number of floating-point arithmetic operations (2𝑀𝑟𝑁𝑟𝐶𝑖,𝑏 flops) to update it. The costs of accessing 𝐵and 𝐴from the microkernel are respectively tackled by arranging the filter tensor following a particular layout in memory, and preparing the data for the microkernel, as detailed next. Ensuring unit-stride accesses to 𝐵: data layout for filter tensor. We aim at maintaining a micropanel 𝐵, of dimension 𝐶𝑖,𝑏 ×𝑁𝑟, in the L1 cache during the execution of the microkernel by arranging the entries of the filter tensor by blocks (tiling), as: 𝐹mem[𝐻𝑓][𝑊𝑓][𝐶𝑜∕𝐶𝑜,𝑏][𝐶𝑖][𝐶𝑜,𝑏], combined with an appropriate selection of the cache configuration parameters 𝐶𝑖,𝑏, 𝑁𝑟, which match the size of the L1 cache and its associativity degree [24]. This has the additional advantage of ensuring that the micropanel of 𝐵is properly arranged in memory so as to accommodate unit-stride accesses from the microkernel. Also, the cost Journal of Systems Architecture 135 (2023) 102806 5 S. Barrachina et al. 1#define VFMA(Creg , Breg, Creg , offset)\ 2Creg = vfmaq_laneq_f32(Creg , Breg, Areg , offset) 3#define VLD vld1q_f32 4#define VST vst1q_f32 5 6#define ACCESS_C(VOP) {\ 7VOP(&C[0], C00); VOP(&C[4], C01);\ 8VOP(&C[8], C02);\ 9// ... Omitted for brevity\ 10 VOP(&C[ldC*6], C60); VOP(&C[4+ldC*6], C61);\ 11 VOP(&C[8+ldC*6], C62);} 12 13 #define LOAD_B(j) {\ 14 B0 = VLD(&Bptr[j]); B1 = VLD(&Bptr[4+j]);\ 15 B2 = VLD(&Bptr[8+j]);} 16 17 #define UPDATE_C(j) {\ 18 VFMA(C00 ,B0,A0,j); VFMA(C01,B1,A0,j); VFMA(C02,B2,A0,j);\ 19 // ... Omitted for brevity \ 20 VFMA(C60 ,B0,A6,j); VFMA(C61,B1,A6,j); VFMA(C62,B2,A6,j);} 21 22 // Load micro -tile of C 23 ACCESS_C(VLD) 24 25 // Iterate inside loop from kr = 0 to Cib -3 26 baseB = 0; 27 for ( kr = 0; kr < Cib -3; kr += 4 ) { 28 // Load 7 x 4 block of A 29 A0 = VLD(&A[kr]); 30 // ... Omitted for brevity 31 A6 = VLD(&A[kr+ldA*6]); 32 // Load 1x12 row of B; Update micro -tile w.r.t column of A 33 LOAD_B(baseB); UPDATE_C(0);LOAD_B(baseB+12);UPDATE_C(1); 34 LOAD_B(baseB+24);UPDATE_C(2);LOAD_B(baseB+36);UPDATE_C(3); 35 baseB += 4*Nr; 36 } 37 38 // Last iterations , from kr+1 to Cib-1, omitted for brevity 39 // Store micro -tile of C 40 ACCESS_C(VST) Listing 5: Structure of the 7 × 12 microkernel for algorithm A, implemented using ARM NEON intrinsics. of this reorganization is low and only needs to be done once, as the filters of a DNN model do not vary during the inference process. We next describe two distinct strategies to deal with the non-unit accesses to the elements of 𝐴from the microkernel, leading to two distinct algorithms, described in the next two subsections. 4.1. New Algorithm A Amortizing the cost of non-unit stride accesses to 𝐴with a specialized microkernel. To illustrate this technique in detail, via a concrete example, let us consider the IEEE single-precision floating-point (FP32) format as the baseline datatype for the tensor data and arithmetic, and consider a target ARMv8 architecture with 32 (SIMD) vector registers of width 128 bits (i.e., with capacity for 4 FP32 numbers). We can follow [8] to implement a microkernel that operates with a microtile of dimension 𝑀𝑟×𝑁𝑟= 7 × 12, while unrolling the loop that traverses the 𝐶𝑖,𝑏 dimension with a factor 4, with the following implications: •21 (7 ×3) vector registers are employed to store the entries of the microtile 𝐶. As each vector register can store 4 FP32 numbers, this results in the required 7 × (3 ⋅4) = 7 × 12 microtile. •7 vector registers are dedicated to storing 4 entries each of 7 consecutive rows of 𝐴. •3 vector registers are employed to store (3 ⋅4) = 12 entries of 𝐵. •Each iteration of the microkernel uploads 7 rows of 𝐴, to then repeat four times the following: load a row of 𝐵(12 elements) and then update the microtile (7 ×12) 𝐶with respect to one of the block columns uploaded from 𝐴and the row of 𝐵. •Vector FMAs are used for the updates. The general structure of the microkernel is illustrated in Listing 5, showing the use of ARM NEON intrinsics for the vector loads, stores and FMAs. Considering the data layout for the data tensors, we remark the following details: •As corresponds to the NHWC format, the entries of the microtile 𝐶are stored in row-major order with row stride ldC=𝐶𝑜. The complete microtile is loaded before the microkernel loop, using 21 vector instructions, and stored after it, also with 21 vector instructions. The costs of these memory accesses are amortized over the 𝐶𝑖,𝑏 iterations of the loop. •For compatibility, the entries of the micropanel 𝐴are also stored in row-major order, in this case with row stride LDA=𝐶𝑖. By unrolling the loop with a factor 4, at each iteration of the loop in the microkernel and for each row of 𝐴we upload 4 consecutive entries of the micropanel using a vector instruction. This Journal of Systems Architecture 135 (2023) 102806 6 S. Barrachina et al. Fig. 5. Memory layout for the 𝑊𝑜,𝑏 ×𝐶𝑖,𝑏 buffer 𝐴𝑐. compensates the cost of accessing non-consecutive rows of 𝐴in memory. •The entries of the micropanel 𝐵are stored in row-major order, with row stride 𝑁𝑟= 12. Each iteration of the loop in the microkernel retrieves (4⋅12) = 48 consecutive entries of the panel, yielding an access with unit stride. Other microkernel dimensions. This technique is indeed not exclusive for the 𝑀𝑟×𝑁𝑟= 7× 12 microkernel. For example, the same principles can be applied to an 𝑀𝑟×𝑁𝑟= 6 × 16 microkernel where: •24 (6 ×4) vector registers are employed to store the microtile of 𝐶. •6 vector registers are dedicated to 4 entries each of 6 consecutive rows of 𝐴. •2 vector registers are employed to store (2⋅4) = 8 entries of 𝐵, but these vector registers are reutilized four times during the unrolled loop iteration. This yields a total utilization of 32 vector registers in the variant with a 6 ×16 microkernel versus the 31 vector registers employed in the 7×12 microkernel. In general, the dimensions of the micro-kernel are constrained by the number of vector registers in the processor architecture. Here we follow the principles in [24] to select the dimensions that maximize the usage of vector registers without incurring in register spilling during the execution of the micro-kernel. This led us to choose several configurations among which, we will only report results for those that delivered the best general performance for each algorithm. Amortizing the cost of loading/storing the microtile of 𝐶.We next note that the microtile of 𝐶updated in the new Algorithm A is independent of the loops traversing the 𝐻𝑓, 𝑊𝑓filter dimensions (indexed by 𝑛, 𝑚 in our algorithms). Therefore, it is possible to move these two loops inside the microkernel, so as to amortize the costs of these accesses to the 𝐶operand, which is expected to reside in the main memory and, therefore, can be expensive to access. 4.2. New Algorithm B Ensuring unit-stride accesses to 𝐴with via packing. The previous algorithm requires a specialized microkernel, different from that available in the BLIS library for matrix multiplication. Developing this microkernel for each particular processor architecture can be a complex and costly process, as this is usually done using assembly instructions to carefully prefetch the data into the processor vector registers in order to completely overlap communication with computation [15,25,26]. The alternative strategy adopted in the algorithm proposed next follows the BLIS strategy to pack a 𝑊𝑜,𝑏 ×𝐶𝑖,𝑏 block of 𝐷=𝐴into a temporary buffer 𝐴𝑐, in blocks of 𝑀𝑟rows as illustrated in Fig. 5. This packing is done inside the loop that traverses the 𝑊𝑜dimension of the problem in Listing 4. For illustrative purposes, Listing 6shows a BLIS-like implementation of an 8 ×12 microkernel for algorithm B that assumes that 𝐴is packed into a contiguous buffer 𝐴𝑐. There: •24 (8 ×3) vector registers are employed to store the entries of the 8 ×12 microtile 𝐶. •2 vector registers are dedicated to store 8 consecutive entries of a column of 𝐴𝑐. •3 vector registers are employed to store 12 entries of a row of 𝐵. •Each iteration of the microkernel uploads a column of the buffer 𝐴𝑐and 12 entries in a row of 𝐵to then update the 8 ×12 microtile of 𝐶with respect to them. •Vector FMAs are used for the updates. Compared with all the previous algorithms for the direct convolution, including our new algorithm A, on the negative side: 1. This variant introduces a memory overhead due to the additional workspace for the 𝐴𝑐buffer. 2. There is some cost due to the copying of the data into 𝐴𝑐. On the positive side: 1. This new algorithm B can leverage the natural BLIS microkernel to perform the small gemm. 2. The cost of data copying is likely to be well amortized over enough computations in the inner loops. 3. The rearrangement ensures unit-stride accesses to the micropanel of 𝐴from the microkernel; see the algorithm in Listing 6. Algorithm B presents an additional advantage for convolution operators with horizontal and/or vertical strides larger than 1. Concretely, this type of access to the input tensor can be accommodated when packing the entries of the 𝑊𝑜,𝑏 ×𝐶𝑖,𝑏 of 𝐴into the buffer 𝐴𝑐so that the high-performance implementation of the microkernel does not need to be modified. Amortizing the cost of loading/storing the microtile of 𝐶.A similar idea to that described at the end of the previous subsection can be applied to Algorithm B. However, that would imply that the use of the BLIS microkernel is no longer an option and, therefore, we do not consider it further for this algorithm. 4.3. Vector intrinsics versus assembly The implementations of the 7 ×12 microkernel for Algorithm A and the 8 ×12 microkernel for Algorithm B, respectively shown in Listing 5and 6, employed ARM NEON intrinsics. This leaves in the hands of the compiler the decision on (1) how to translate each intrinsic into a specific collection of assembly instructions; and (2) which architecture vector registers to use in order to store each variable, which often results in large variability in the attained performance. This is because the characteristics of individual instructions (e.g., latency and throughput) are significant factors in the design of the microkernel. As such, we chose to implement the microkernel using inline assembly to specifically use certain types of instructions. This is illustrated by redefining the VFMA macro utilized in the two microkernels as in Listing 7. For the second issue, we noticed that the compiler tends to produce code that resulted in a significant amount of register spilling, resulting in superfluous load and store instructions being generated. We overcame this issue by implementing the entire microkernel in assembly. We have done so for some of our simpler cases, including the most used microkernels for Algorithms A and B. Here it is worth mentioning that the development of microkernels is a systematic process, which can be automated to a certain extent with good performance results [27]. We leave the application of this research line to the convolution operator as part of future work. Journal of Systems Architecture 135 (2023) 102806 7 S. Barrachina et al. 1#define VFMA(Creg , Breg , Creg , offset)\ 2Creg = vfmaq_laneq_f32(Creg , Breg , Areg, offset) 3#define VLD vld1q_f32 4#define VST vst1q_f32 5 6#define ACCESS_C(VOP) {\ 7VOP(&C[0], C00); VOP(&C[4], C01);\ 8VOP(&C[8], C02);\ 9VOP(&C[ldC], C10); VOP(&C[4+ldC], C11);\ 10 VOP(&C[8+ldC], C12);\ 11 // ... Omitted for brevity\ 12 VOP(&C[ldC*7], C70); VOP(&C[4+ldC*7], C71);\ 13 VOP(&C[8+ldC*7], C72);} 14 15 // Load micro -tile of C 16 ACCESS_C(VLD) 17 18 // Iterate inside loop from kr = 0 to Cib -1 19 baseA = 0; baseB = 0; 20 for ( kr = 0; kr < Cib; kr++ ) { 21 // Load column of A 22 A0 = VLD(&A[baseA]); 23 A1 = VLD(&A[baseA+4]); 24 // Load row of B 25 B0 = VLD(&B[baseB]); 26 B1 = VLD(&B[baseB+4); 27 B2 = VLD(&B[baseB+8]); 28 // Update micro -tile w.r.t column of A and row of B 29 VFMA(C00 ,B0,A0 ,0); VFMA(C01,B1,A0,0); VFMA(C02 ,B2,A0 ,0);\ 30 VFMA(C10 ,B0,A0 ,1); VFMA(C11,B1,A0,1); VFMA(C12 ,B2,A0 ,1);\ 31 // ... Omitted for brevity\ 32 VFMA(C70 ,B0,A1 ,3); VFMA(C71,B1,A1,3); VFMA(C72 ,B2,A1 ,3);} 33 baseA += Mr; baseB += Nr; 34 } 35 36 // Store micro -tile of C 37 ACCESS_C(VST) Listing 6: Structure of the 8 × 12 microkernel for algorithm B, implemented using ARM NEON intrinsics. 1#define VFMA(Creg , Breg , Areg , offset)\ 2__asm__ volatile\ 3(\ 4" fmla %[c_reg].4s, " \ 5" %[b_reg].4s, " \ 6" %[a_reg].s[ " #offset " ] \n\t " \ 7: [c_reg] " +w " (Creg)\ 8[b_reg] " w " (Breg),\ 9: [a_reg] " w " (Areg)\ 10 );\ Listing 7: Use of assembly inlining. 5. Experimental results 5.1. Hardware setup All the experiments in this section were performed using IEEE single-precision (FP32) on a single core of an NVIDIA Carmel (ARMv8.2) processor. This architecture is equipped with a private 4way associative (data) L1 cache of 64KiB per core, a 16-way associative L2 cache of 2 MiB shared between each pair of cores, a 16-way associative L2 cache of 4 MiB shared by all 8 cores, and 4GiB of RAM. To avoid the performance distortions caused by the utilization of distinct power modes (and associated processor core frequencies), the operating core frequency was set to 2.3GHz. The theoretical peak performance for a single core of this architecture, operating at such frequency, is 36.8 (FP32) GFLOPS (billions of floating-point operations, or flops, per second). Moreover, the single worker thread was pinned to a specific core of the NVIDIA Carmel processor. There are several approaches to parallelize the convolution algorithms targeted in this work by exploiting multi-threading on a multicore processor. In [28], we evaluate different parallel approaches and conclude that the best strategy is to run one instance of the model per core in the system, and then distribute the batch among the different instances. From that point of view, we believe that an evaluation using a single thread is illustrative of the behavior that could be attained with this parallelization strategy. 5.2. Data setup For the experiments, we selected the convolution layers present in two popular DNNs: GoogleLeNet [29] and ResNet-50 (v1.5) [2,30]. From these, we discarded those layers with horizontal/vertical strides greater than 1, because they are more difficult to tackle efficiently by the blocked algorithm in [7] and the new algorithm A described in Section 4. Overall, the convolution layers included in the following experimental study represent around 70% of the total execution time of the CNN on the target architecture (though we note that our codes do not exploit other orthogonal optimization techniques, such as fusion, Journal of Systems Architecture 135 (2023) 102806 8 S. Barrachina et al. which could increase the contribution of the convolution operators to the total cost of the inference). The dataset that was selected in both cases was ImageNet [31] and the batch size was set to 𝑁= 1 in order to reflect the single stream scenario of the ML Commons benchmark for inference on the edge.2 5.3. Algorithms In the comparison, we include five options for the calculation of the convolution operators appearing in the corresponding layers of the two selected DNNs: •lowering: The lowering approach with the im2col transform implemented as an optimized C routine; and the subsequent gemm performed by invoking BLIS (v0.8.1), with the configuration of this library adjusted to the NVIDIA Carmel processor. (We also ran the same option linked with OpenBLAS version v0.3.19 and ARMPL version v.21.1, in general obtaining lower performance. Therefore, we do not consider their gemm instances for the lowering alternative in the following discussion.) •blocked: The blocked algorithm for the direct convolution described in [7]. In this case, we developed a microkernel with 𝑂 resident in registers of dimension 𝑊𝑜,𝑏 ×𝐶𝑜,𝑏 = 7 × 12, following the ideas in [8], and optimized it using ARM NEON intrinsics (see Listing 5). We remind that, because of the connection between consecutive convolution layers, this algorithm enforces the constraint that 𝐶𝑖,𝑏 =𝐶𝑜,𝑏. •new-A: The blocked NHWC-preserving Algorithm A for the direct convolution introduced in Section 4of this paper, with 𝑀𝑟×𝑁𝑟= 7 × 12 and 6 ×16 microkernels optimized using ARM NEON intrinsics (see, e.g., Listing 5). For this particular case, we also developed our own assembly versions of these two microkernels. For this algorithm, we select a value for 𝐶𝑖,𝑏 which favors that a𝐶𝑖,𝑏 ×𝑁𝑟micropanel of 𝐵(corresponding to a block of the filter tensor) remains in the L1 cache during the execution of the microkernel [24]. •new-B: The blocked NHWC-preserving Algorithm B for the direct convolution in Section 4, with the cache configuration parameters selected via the analytical model in [24]. In this case, we either (1) employ our own manually encoded 𝑀𝑟×𝑁𝑟= 8×12 and 4 ×20 microkernels using ARM NEON intrinsics; (2) utilize our own equivalent microkernels fully developed in assembly following the approach in Listing 6; or (3) rely on this specific component in BLIS for the ARM processor. •Winograd: The implementation of this algorithm in the NNPACK library,3with a 3 ×3 filter. In most of the following plots, the results are reported in terms of GFLOPS, using a cost of 2𝐶𝑖𝐶𝑜𝐻𝑜𝑊𝑜𝐻𝑓𝑊𝑓flops for all the previous algorithmic options. The GFLOPS rate is inversely proportional to the execution time of the algorithm and presents the advantage of being bounded by the peak performance of the target platform (36.8 GFLOPS). The utilization of this count also for the Winograd-based algorithm (even though it performs a considerably smaller number of flops,) allows a direct comparison with the other approaches. The execution of each implementation was repeated for at least 2 minutes, and the results were then averaged. We do not include results for the (Apache) AutoTVM tool TVM because optimizing the convolution operator using it requires a considerable execution time. Concretely, AutoTVM explores about 1,000 options per convolution layer, resulting in between 40 min and 1 h to optimize each layer. In contrast, our direct convolution algorithms are straight-forward to tune for a particular architecture, for almost 2See https://mlcommons.org/. 3https://github.com/Maratyszcza/NNPACK any convolutional layer, independently of its specific parameters. Given these reasons, we prefer to omit the comparison against TVM as we would like to preserve the focus of the paper on those optimizations which can be obtained by mimicking the techniques usually applied in linear algebra libraries for the matrix multiplication. 5.4. Performance evaluation Microkernels for the new variants. The performance of new-A and new-B strongly depends on the implementation of the underlying microkernel. For this particular work, due to the dimensions of the convolution operators appearing in the two target models, we ‘‘manually’’ developed microkernels that employ ARM NEON intrinsics (or, alternatively, assembly instructions) for the two variants, with 𝑀𝑟×𝑁𝑟= 7 × 12 (see Listing 5) and 6 ×16 for new-A; and 𝑀𝑟×𝑁𝑟= 8 × 12 (see Listing 6) and 4 ×20 for new-B. In addition, for new-B, we consider the native microkernel for ARM from BLIS (which also sets 𝑀𝑟×𝑁𝑟= 8 × 12). As an initial point for our evaluation, we therefore assess the impact of the selected microkernel on the performance of the variants. The results in Figs. 6 and 7for new-A show that the 6 ×16 microkernel is consistently the best option for GoogleLeNet-ImageNet and dominates the performance in 24 out of 46 layers for ResNet-50-ImageNet. (The layers of the GoogleLeNet-ImageNet and ResNet-50-ImageNet models were respectively split into 3 and 2 plots in these figures to improve visibility.) For new-B,Figs. 8 and 9report the clear performance gain when using the 8 ×12 microkernel for GoogleLeNet-ImageNet model and, in general, superior performance for the same microkernel for ResNet-50-ImageNet (in 31 out of 46 layers). In addition, for variant new-B, the manual microkernels offer superior performance to that observed when using the native implementation in BLIS. This may be surprising but in an independent experiment we could deduce that the reason for this behavior is that the BLIS microkernel is highly optimized when the actual dimensions of the microtile that is to be processed, say 𝑚𝑟×𝑛𝑟, match the dimensions of the native microkernel. However, the BLIS implementation suffers a serious drop in performance for the ‘‘border’’ cases, when 𝑚𝑟< 𝑀𝑟= 8,𝑛𝑟< 𝑁𝑟= 12, or both. To increase the readability of the following plots, hereafter we will refer with the labels new-A or new-B to the implementations that respectively integrate the optimal microkernel into Algorithm A or B, which depends on the particular model and layer. Winograd. While we recognize that the Winograd-based convolution may offer competitive performance, depending on the layer properties, this alternative approach reduces the number of arithmetic operations at the cost of sacrificing numerical accuracy and flexibility. Therefore, a comparison based only on raw performance is incomplete. Nonetheless, we evaluated the Winograd-based algorithm in NNPACK for ResNet50 model. The first aspect to note is that the package could be used with very few convolution layers of the model: 17 out of 46 (around 37% only). The second issue is that the performance in all these cases varied between 7.5 and 9.5 GFLOPS, therefore being much lower than that observed for the direct convolution implementations applied to the ResNet-50 model with ImageNet presented earlier. Therefore, we exclude the Winograd-based algorithm from the global evaluation in the following experiments. Global comparison. Figs. 10 and 11 report the performance rates for the distinct implementations of the convolution operator and the two DNN models considered in this work. These results show fair benefits for the new-A and new-B implementations, which outperform the other options for most layers of both models. In particular, for the GoogleLeNet-ImageNet pair, the new variants new-B and new-A outperform the (best) lowering approach by an average speed-up of 1.22 and 1.09 respectively. For new-B, the highest speed-up is 1.88 (layer #84) and the lowest one is 0.49 (layer #1). In comparison, for new-A the highest speed-up is 1.85 (layer #14) and the lowest one is 0.5 (layer #1). Journal of Systems Architecture 135 (2023) 102806 9 S. Barrachina et al. Fig. 6. Performance evaluation of the new-A variant for GoogleLeNet. Fig. 7. Performance evaluation of the new-A variant for ResNet-50. For RestNet-50, the average speed-up of new-A/new-B with respect to the (best) lowering approach is 1.08 and 1.23 respectively, with the highest speed-up being 1.87 (layer #167) and the lowest one 0.90 (layer #92) for new-B; and 1.62 (layer #160 and #170) and 0.79 (layer #80) for new-A, respectively. While the speed ups attained for the individual layers for GoogleLeNet-ImageNet or ResNet-50-ImageNet can be regarded as a useful piece of information, to put them into perspective we need to consider their impact on the cost of each layer. In order to do this, we integrated the distinct convolution algorithms into our PyDTNN framework for deep learning [32,33], and performed a complete evaluation of the inference process. Fig. 12 displays the absolute (aggregated) time of the two options: lowering and the new-B algorithm. (The newAalgorithm could not be considered in this comparison as it neither