Distributed Deep Learning Using Synchronous Stochastic Gradient Descent

Dipankar Das, Sasikanth Avancha, Dheevatsa Mudigere, Karthikeyan Vaidynathan, Srinivas Sridharan, Dhiraj Kalamkar, Bharat Kaul, Pradeep Dubey

Introduction

With the efficacy of large deep neural networks well established (Schmidhuber, 2015), the key challenge is to train them in a reasonable amount of time, preferably hours or even minutes. The largest networks require several Exaflops of computationVGG-A needs 33.6 GFlops per image, 43 PFlops per epoch on ImageNet-1k, and 1 ExaFlop per 25 epochs; clearly a single node or single card implementation fails to meet this challenge. Distributed (i.e., multinode) training methods are required to address this gap. Also, to stay relevant, deep learning must ride the scaling curve like it took advantage of increasing compute density and Moore’s law in recent years. In order to address this challenge, researchers have developed multinode/multi-card frameworks such as TensorFlow (Abadi et al., 2015), FireCaffe (Iandola et al., 2015), and DeepImage (Wu et al., 2015).

However, scaling synchronous Stochastic Gradient Descent (SGD) is challenging, as it is a strong-scaling problem on a small amount of work bound by compute needs for processing a minibatch of data-points (between 64-5120). Thus, many of these frameworks suffer from poor scalability, and typically do not scale beyond tens of nodes. Many variants of synchronous SGD such as 1-bit SGD (Seide et al., 2014a), elastic-SGD (Zhang et al., 2014), as well as asynchronous SGD (Chilimbi et al., 2014) have been proposed for better scaling. However, unlike these methods we do not alter hyperparameters (like minibatch or learning rate) or the algorithm, or use any compression methods, but focus on deeply understanding the vanilla SGD algorithm and scaling it.

Our approach is to systematically develop detailed system balance equations, and solve them to obtain limits for performance; we also find optimal design points for Xeon-based multinode training frameworks. We examine data-parallelism and model-parallelism, and propose a new algorithm for hybrid data- and model-parallelism. Based on this analysis we identify which strategy is best suited for different layers of a neural network.

However, before we build a distributed deep-learning training framework, we must first optimize a single node implementation to the highest possible efficiency. Again, we take an analytic approach to study the balance between computation and memory bandwidth in order to obtain optimal cache blocking strategies. Specific to x86-based architectures, we present details of threading, register blocking, as well as instruction sequence of innermost loops, and analyze how these techniques enable a deep learning framework to achieve high single-node efficiency.

We note that the analysis in this work is generic and applicable, to the best of our understanding, to other non-x86-based deep learning systems as well (including those using GPUs and accelerators).

Combining a highly efficient single node implementation and an optimized multi-node scaling strategy, we present the best time-to-train for several CNNs used for image recognition, and break all published training throughput records using Xeon based systems. We achieve efficiency of 9̃0% for several convolutional layer operations and 7̃0% for fully-connected layers on Intel Xeon E5-269Xv3 machines for a wide variety of neural networks. Our multinode algorithms and tuned implementations help us significantly surpass best published scaling efficiency for a wide gamut of neural networks.

The rest of the paper is organized as follows: in section 2 we discuss balance equations for the convolutional and fully-connected compute layers and in section 3 for the communication patterns. Section 4 describes the components of our software framework and their design. We present detailed single-node and multi-node performance results for two well-known deep learning topologies – OverFeat (Sermanet et al., 2013) and VGG-A (Simonyan & Zisserman, 2014) – on the Cori system http://www.nersc.gov/users/computational-systems/cori/ and AWS https://aws.amazon.com/.

Optimizing Computation in Neural Network Training

We can view the computation of neural network training as a task graph where each node represents a block of computation (typically one layer), and edges represent data-flow (usually a multi-dimensional tensor), which define data-dependencies between nodes. It is critical to understand the compute and memory-bandwidth needs of nodes, and identify optimal threading, cache-blocking, vectorization, and register-blocking strategies.

The compute heavy convolution and fully-connected layers take two k+2k+2-dimensional tensors as input and produce a k+2k+2-dimensional output tensor. Here kk is the number of dimensions of a feature map, or kernel. The additional two dimensions depend on the type of data: for tensors containing inputs and outputs or gradients of inputs and outputs, they represent the minibatch and feature identifier; for tensors containing weights and gradients of weights, they represent a pair of input-output features.

The compute operations for forward propagation, backpropagation and determining weight gradient are identical 2k+32k+3-dimensional loops. As an example, consider a 2-D convolution forward propagation operation, where the input is a 4-D tensor over minibatch, input feature map (ifm) identifier, output feature map height (outhout_{h}) and width (outwout_{w}) respectively. The weight tensor is another 4-D tensor over input feature map, output feature map (ofm), kernel width (kwk_{w}) and kernel height (khk_{h}). In the rest of the paper, we refer to kernelkernel and weightsweights interchangeably. Algorithm 1 describes the forward propagation operation. The variable ss represents the stride for the convolution.

A fully-connected layer can similarly be written as a special case of this 7-nested loop when khk_{h}, kwk_{w}, outhout_{h}, outwout_{w} are all 1. Moreover other operations such as backpropagation and weight gradient computation have identical loops with different multiply and accumulate operations.

Backpropagation: grad_input[i0,i1,i3∗s+i5−1,i4∗s+i6−1]+=grad\_input[i_{0},i_{1},i_{3}*s+i_{5}-1,i_{4}*s+i_{6}-1]+= grad_output[i0,i2,i3,i4]∗wts[i1,i2,i5,i6]grad\_output[i_{0},i_{2},i_{3},i_{4}]*wts[i_{1},i_{2},i_{5},i_{6}]

Weight-gradient Update: gradwts[i1,i2,i5,i6]+=grad_{w}ts[i_{1},i_{2},i_{5},i_{6}]+= input[i0,i1,i3∗s+i5−1,i4∗s+i6−1]∗input[i_{0},i_{1},i_{3}*s+i_{5}-1,i_{4}*s+i_{6}-1]* grad_output[i0,i2,i3,i4]grad\_output[i_{0},i_{2},i_{3},i_{4}]

This similarity between all the operations implies that the memory access pattern for all three operations is identical, and therefore, a cache blocking strategy for one of them should apply to others as well.

2 Cache Blocking

Unless the activations and weights completely fit within the CPU cache hierarchy (which is often not the case), the loop over i3i_{3} reads outh∗outwout_{h}*out_{w} output activations, (outh∗s+kh−1)∗(outw∗s+kw−1)(out_{h}*s+k_{h}-1)*(out_{w}*s+k_{w}-1) input activations (which we denote inh∗inwin_{h}*in_{w}), and kh∗kwk_{h}*k_{w} weights from external memory (i.e., DRAM). The loop performs kh∗kw∗outw∗outhk_{h}*k_{w}*out_{w}*out_{h} multiply-and-accumulate operations. Therefore, the Bytes to FLOPs (floating-point operations) ratio is: B/F=sizedata∗(outw∗outh+inw∗inh+kw∗kh)B/F=size_{data}*(out_{w}*out_{h}+in_{w}*in_{h}+k_{w}*k_{h}) /(2∗kw∗kh∗outw∗outh)/(2*k_{w}*k_{h}*out_{w}*out_{h}).

For a convolutional layer with 12∗1212*12 output, 3∗33*3 kernel, 512 input feature maps and 1024 output feature maps (such as C5 in OverFeat-FAST), the B/F ratio is 0.54; typically the system B/F ratio is less than 0.08. On the other hand, if all data (inputs, outputs, weights) fit into the cache hierarchy, the B/F ratio reduces because with a one-time read from DRAM, all 7 loops can be computed at one go: B/F=sizedata∗(minibatch∗ofm∗outw∗outh+minibatch∗ifm∗inw∗inh+ifm∗ofm∗kw∗kh)B/F=size_{data}*(minibatch*ofm*out_{w}*out_{h}+minibatch*ifm*in_{w}*in_{h}+ifm*ofm*k_{w}*k_{h}) /(2∗minibatch∗ofm∗ifm∗kw∗kh∗outw∗outh)/(2*minibatch*ofm*ifm*k_{w}*k_{h}*out_{w}*out_{h}).

Now, the best achievable B/F ratio for C5 in OverFeat-FAST is 0.003. Clearly, the capacity (i.e., size in bytes) of the cache hierarchy determines the range of B/F ratios. Given the capacity, we formulate the cache blocking problem as a constrained minimization problem (s1,s2,…sk+2s_{1},s_{2},\dots s_{k+2} are strides):

BS=sizedata∗(b10∗…b1k+2+b20∗…b2k+2+b10∗b20∗(b12∗s1+b22−1)…(b1k+2∗sk+2+b2k+2−1))BS=size_{data}*(b1_{0}*\ldots b1_{k+2}+b2_{0}*\ldots b2_{k+2}+b1_{0}*b2_{0}*(b1_{2}*s_{1}+b2_{2}-1)\ldots(b1_{k+2}*s_{k+2}+b2_{k+2}-1)) CPB=2∗b10∗b12∗…b1k+2∗b21∗b22∗…b2k+2CPB=2*b1_{0}*b1_{2}*\ldots b1_{k+2}*b2_{1}*b2_{2}*\ldots b2_{k+2} B/F=BS/CPBB/F=BS/CPB ∀i find b1i, b2i\forall_{i}~{}find~{}b1_{i},~{}b2_{i} that minimize B/FB/F, s.t. BS<SizecacheBS<Size_{cache}

Here, BSBS is the size of the block residing in on-chip memory, CPBCPB is the amount of compute to perform on the block, and SizecacheSize_{cache} is the size of the on-chip memory/cache (with due consideration for double buffering).

We observe that traversing along consecutive blocks in any dimension, results in memory reuse and therefore, better B/F ratios. For example, consider the scenario where block size along the loop i3i_{3} is 1, and we traverse along the height-dimension in a 2-D convolution with stride=1. Here we need to read in only one row of ifmifm for each row of ofmofm, instead of khk_{h} rows of ifmifm. Similarly, traversing along the ifmifm dimension precludes reading the output-block.

We write a multithreaded program to perform a brute-force state space search over all values of loop iterators in order to find the minimum B/F ratio for different 2-D convolutional layers, given a limit on the cache size. Additionally, to compute BSBS, we constrain one of the dimensions to be a multiple of SIMD-width or warp-size since all modern high-performance processors use some form of SIMD operation. Hence the sub-problem to be solved must have one dimension (preferably the output dimension) set to a multiple of SIMD-width.

We find that with 128 KB of cache per thread in modern Xeon CPUs, a B/F ratio of ≤\leq 0.04 can be maintained for most convolutional layers even for a minibatch size of 1.

3 Data Layout and Vectorization

In addition to cache blocking, we need to vectorize the operations and perform register blocking as well. A fully optimized vectorized forward propagation operation is presented in Algorithm 2. It contains a 10-nested loop with cache blocking blocking along ifmifm, and ofmofm, and register blocking along outputhoutput_{h}, and outputwoutput_{w} dimensions. For this blocked loop structure we lay out data so that access in the innermost loops is as contiguous as possible. This results in better utilization of cache lines (and hence bandwidth) and improves prefetcher performance.

In this work we lay out all data, including activations and weights with the innermost dimension over groups of SIMD-width (SWSW) output feature maps. That is we lay out the different data structures as: and gradient of activations: N×C×H×W→N×(C/SW)×H×W×SWN\times C\times H\times W\rightarrow N\times(C/SW)\times H\times W\times SW Weights and gradients of weights: IFM×OFM×KH×KW→IFM×(OFM/SW)×KH×KW×SWIFM\times OFM\times KH\times KW\rightarrow IFM\times(OFM/SW)\times KH\times KW\times SW Transpose-weights: IFM×OFM×KH×KW→OFM×(IFM/SW)×KH×KW×SWIFM\times OFM\times KH\times KW\rightarrow OFM\times(IFM/SW)\times KH\times KW\times SW

Here, N stands for minibatch, C stands for feature-maps, H for feature map height, W for feature map width, IFM for input feature maps, OFM for output feature maps, KH for kernel-height and KW for kernel-width. This layout also enables vectorization of operations, such that the multiply-and-accumulate in Algorithm 1 can now be replaced by a broadcast and vector fused-multiply and add (Algorithm 2).

4 Register Blocking

The aim of register blocking is two-fold: firstly it improves the ratio of vector fused multiply and add (VFMA) operations to that of load/store operations. Secondly, a sequence of regreg consecutive VFMA instructions is needed to hide the latency of these instructions. The latency for a VFMA operation on the Xeon CPU core is 5 cycles, and a Xeon CPU core can execute 2 VFMA instructions per cycle. Hence in order to completely hide the latency of VFMA instructions we should have a register block size of at least 10. Hence 15≥RBh∗RBw≥1015\geq RB_{h}*RB_{w}\geq 10, as we need one register to store the weights.

In Algorithm 2 we illustrate a 2-D register block for forward propagation. Since a Xeon core can perform 2 Loads per cycle, 2 VFMAs per cycle, and 1 store per cycle, the cycles spent on load/store instructions (LSLS) and VFMA instructions (FMAFMA) for the inner loop (line 5-29 in Algorithm 2) can be computed as:

LS=(RBh∗RBw+SW∗(khend−khstart)∗(kwend−kwstart))/2+RBh∗RBwLS=(RB_{h}*RB_{w}+SW*(kh_{end}-kh_{start})*(kw_{end}-kw_{start}))/2+RB_{h}*RB_{w}

FMA=(SW∗(khend−khstart)∗(kwend−kwstart)∗RBh∗RBw)/2FMA=(SW*(kh_{end}-kh_{start})*(kw_{end}-kw_{start})*RB_{h}*RB_{w})/2

In practice RBhRB_{h} is often 1 for forward propagation, since most feature map width are ≥12\geq 12 for CNNs. For RBw=12RB_{w}=12 and khstart=khend=0kh_{start}=kh_{end}=0 and kwstart=kwend=3kw_{start}=kw_{end}=3 (as in case of OverFeat-FAST C5 layer), we compute the efficiency to be: 88%. The loop for backpropagation is similar with blocking along inpwinp_{w}, inphinp_{h} and ifmifm, instead of outwout_{w}, outhout_{h}, and ofmofm.

Unlike forward and backpropagation, weight gradient computation has the weight (-gradient) kernel as the output. These kernels are often small (3x3, 5x5, 7x7, 11x11), and even two dimensional blocking will only yield a theoretical peak efficiency of 75% for a 3x3 kernel. Hence there is often a need to perform blocking along the ifmifm dimension as well. Indeed we propose using specific tailored strategies for different kernel sizes.

3x3 kernel: register block with one row (3 SIMD-elements) of 4 consecutive kernels along the input-feature-map dimension.

5x5 and 7x7: register block with one row of 2 consecutive kernels along the input-feature map dimension.

11x11: one dimensional register block along the kernel-width dimension.

5 Threading and Work Partitioning

We perform fine grained partitioning of work across threads. For the forward and backpropagate operations, we partition the work across multiple minibatches into jobs, each for one row of the output/input across SWSW output/input features. These jobs are then equally distributed across the different threads (iterations in lines 1-4). For weight update, we treat weight kernels for SWSW input- and output-feature pairs to be the basic unit of work, and these jobs are subsequently distributed across multiple threads. In case the number of jobs created in this way is low (like for C1 layers), we additionally partition the problem along the minibatch dimension, and then privtize and reduce weight gradient computation.

Optimizing Communication

In this work we perform strong scaling of the synchronous minibatch stochastic gradient descent algorithm. We scale the computation for one iteration across multiple nodes, such that the multi-threaded, and multi-node parallel implementation is equivalent to a single-node single-threaded serial implementation. We present a detailed theoretical analysis of computation and communication balance equations, and determine strategies for work partitioning between nodes.

We analyze the algorithmic computation-to-communication balance of data parallelism, wherein the work in an iteration is partitioned across the minibatch. For our analysis we consider a butterfly-reduce operation.

Consider a convolutional layer with ofmofm output feature maps each of size: outw∗outhout_{w}*out_{h} (width and height), ifmifm input feature maps, ss stride, and kernel of size kw∗khk_{w}*k_{h}. The amount of computation (for MBnodeMB_{node} data points assigned to a node) in the number of FLOPS in this layer for forward, backward, and weight gradient computation steps is: Comp=3∗2∗MBnode∗ifm∗ofm∗kw∗kh∗outw∗outhComp=3*2*MB_{node}*ifm*ofm*k_{w}*k_{h}*out_{w}*out_{h}

Similarly, we estimate the total communication per iteration for a data-parallel approach. In each iteration, a node sends partial weight gradients to other nodes, and receives updated weights from them. Therefore, the total communication volume is: Comm=sizedata∗ifm∗ofm∗kw∗kh∗(2−overlap))Comm=size_{data}*ifm*ofm*k_{w}*k_{h}*(2-overlap))

overlapoverlap is the amount of overlap that the software achieves between send and receive operations. Assuming floating point data representation, and overlap=1overlap=1 of sends and receives, the algorithmic communication-to-computation ratio for data parallel implementation of a single layer is: comp_comm=1.5∗outw∗outh∗MBnodecomp\_comm=1.5*out_{w}*out_{h}*MB_{node}

We observe that the algorithmic computation-to-communication ratio does not depend either on the kernel size or number of input and output feature maps or stride, but solely on the size of the output feature-map and the number of data-points assigned per node.

An advantage of data parallelism is that it can be overlapped with computation of current and previous layers. Moreover, the weight update function needs weight gradient values which are available immediately after the backpropagation of a given layer kk and the updated weights are not needed until before the forward propagation step for layer kk in the next iteration. Therefore, we estimate the scalability of layers L0,L1,L2,...,Lk−1L_{0},L_{1},L_{2},...,L_{k-1}, as follows:

ocompi=∑j<icompj+compi/3ocomp_{i}=\sum_{j<i}comp_{j}+comp_{i}/3 ocommsi=∑j<=icommsjocomms_{i}=\sum_{j<=i}comms_{j} bubblei=ocommsi/commsys−ocompsi/compssysbubble_{i}=ocomms_{i}/comm_{sys}-ocomps_{i}/comps_{sys}

Scaling efficiency is then estimated as the ratio of two terms: (∑i<kbubblei+(∑i<kcompi)(\sum_{i<k}bubble_{i}+(\sum_{i<k}comp_{i}) and (compsys)/((∑i<kcompi)/compsys)(comp_{sys})/((\sum_{i<k}comp_{i})/comp_{sys})

The best compute-communication overlap implies that ∀i>0bubblei=0\forall_{i>0}bubble_{i}=0, as we cannot avoid the communication bubble between the weight-gradient update and forward propagation steps of the first layer L0L_{0}. It is notable that the first layer need not perform backpropagation and can only perform weight gradient computation. Moreover, we perform weight gradient computation before backpropagation in order to allow for some more computation to overlap communication (which explains the compi/3comp_{i}/3 term).

We note that in typical convolutional networks (and otherwise), the size of feature maps keeps on reducing monotonically with increase in layer-id. Also, we noted earlier that the algorithmic compute-to-communication ratio depends on only the feature map size; thus, if the communication of layer ll cannot be completely overlapped with compute, then the communication of layer l+1l+1 cannot be overlapped as well. Therefore, to estimate the best possible compute-communication overlap we check if bubblek<0bubble_{k}<0, where layer LkL_{k} is the last layer in the data-parallel regime. For convolutional neural networks, this is usually the last convolutional layer. The number of nodes (N) to which the algorithm scales is: N≤minibatch∗(commssys/compsys)∗(ocompk/ocommsk)N\leq minibatch*(comms_{sys}/comp_{sys})*(ocomp_{k}/ocomms_{k}) The algorithmic computation-to-communication ratio convolutional layers of OverFeat-FAST and VGG-A are 208, and 1456 respectively.

We examine the smallest number of data-points per node for a training run of OverFeat-FAST and VGG-A on two platforms: dual-socket 16-core Xeon E5-2698v3 + FDR Infiniband, and 2-socket 9-core 2.9 GHz Xeon E5-2666 v3 + 10GigE Ethernet (Table 1). Based on this, we estimate the scaling for data parallel parts of a 256 minibatch training run. We estimate that the convolutional layers can be scaled to 128 nodes for OverFeat-FAST and 256 nodes for VGG-A. Note that in a CNN there are several fully connected layers which do not scale much in practice as compared to convolutional layers.

2 Analyzing Model Parallelism

We first consider a simple model parallel approach where each node operates on a part of the model of size: ifmb∗ofmbifm_{b}*ofm_{b} input- and output-feature maps. In this case, the computation for the forward pass, backward pass, or weight-gradient update is given as:

Computation=2∗ifmb∗ofmb∗kernelw∗kernelh∗outputw∗outputh∗minibatchComputation=2*ifm_{b}*ofm_{b}*kernel_{w}*kernel_{h}*output_{w}*output_{h}*minibatch

For the forward pass the amount of data received is:

commsrecv=sizedata∗ifmb∗inputw∗inputh∗minibatch∗(ifm/ifmb−1)comms_{recv}=size_{data}*ifm_{b}*input_{w}*input_{h}*minibatch*(ifm/ifm_{b}-1)

The amount of data sent out by the previous layer is:

commssend=sizedata∗ifmb∗inputw∗inputh∗minibatchcomms_{send}=size_{data}*ifm_{b}*input_{w}*input_{h}*minibatch

Hence the total volume of communication data is: sizedata∗ifm∗inputw∗inputh∗minibatchsize_{data}*ifm*input_{w}*input_{h}*minibatch

Hence the time taken for a forward pass with no compute and communication overlap for a given layer is:

Computation/SysFlops+(commsrecv+commssent)/CommBW+SWlatComputation/SysFlops+(comms_{recv}+comms_{sent})/CommBW+SWlat

It is notable that the communication-bandwidth is dependent on the message size, and along with the impact of software overheads, the performance of model parallel communication pattern falls sharply with decrease in size of the feature map and the number of features.

A question of interest is to determine when is model parallelism preferred to data parallelism. We do a simplified analysis to compare the amount of communicated data in each method. Model parallelism is better if:

sizedata∗ifm∗ofm∗kernelw∗kernelh∗(1+(1−overlap))>sizedata∗ifm∗inputw∗inputh∗minibatchsize_{data}*ifm*ofm*kernel_{w}*kernel_{h}*(1+(1-overlap))>size_{data}*ifm*input_{w}*input_{h}*minibatch

We can simplify this further to: ofm∗kernelw∗kernelh∗(2−overlap)>inputw∗inputh∗minibatchofm*kernel_{w}*kernel_{h}*(2-overlap)>input_{w}*input_{h}*minibatch

In convolutional layers ofmofm is typically less than 1024 and kernelhkernel_{h} and kernelwkernel_{w} are 3, or 5, while inputhinput_{h} and inputwinput_{w} are greater than 14, and minibatch size >> 64. For this case, only for a large kernel size and small minibatch does model parallelism become better. For fully connected layers, where kernelw,kernelh,inputwandinputhkernel_{w},kernel_{h},input_{w}andinput_{h} are 1 whenever ofm>minibatchofm>minibatch model parallelism is better than data parallelism. This is typically the case for most fully connected layers, unless we have large minibatches (>5000>5000) as in case of ASR networks.

3 Analyzing Hybrid Parallelism

Beyond vanilla data- and model-parallelism we explore a hybrid scheme. One may view data-parallelism as partitioning work along the ”minibatch” dimension and model parallelism as partitioning along the ”feature map” dimension. Clearly we can partition work along both ”minibatch” as well as ”feature map” dimensions.

For this we partition nodes into node groups, such that nodes within a group follow a model-parallelism regime while corresponding nodes across node groups follow a data-parallelism regime.

In this scheme the minibatch is partitioned into GG groups each containing N/GN/G nodes and responsible for mbgroup=minibatch/Gmb_{group}=minibatch/G data-points. Model parallelism on this subgroup, for forward as well as back-propagation, leads to an exchange of commsmodelcomms_{m}odel amount of data: commsmodel=2∗sizedata∗ifm∗inputw∗inputh∗mbgroupcomms_{model}=2*size_{data}*ifm*input_{w}*input_{h}*mb_{group}

Communication due to data-parallelism (send and receive weights for a 1/Gth1/G^{th} fraction of the weights): commsdata=sizedata∗ofm∗ifm∗kernelw∗kernelh∗(2−overlap)/Gcomms_{data}=size_{data}*ofm*ifm*kernel_{w}*kernel_{h}*(2-overlap)/G

By changing GG and hence mbgroupmb_{group} we can balance the volume of communication for data and model parallelism. Recall that for data parallelism we have the option to overlap the same across all previous layers, thereby hiding time taken for communication. Trading off some model-parallelism for data-parallelism also helps in increasing message sizes, thereby improving network performance. None the less it is of interest to ask if hybrid data- and model-parallelism can yield any benefits in terms of overall communication volume:

commshybrid=2∗sizedata∗ifm∗inputw∗inputh∗mbgroup+sizedata∗ofm∗ifm∗kernelw∗kernelh∗(2−overlap)∗(G/N)comms_{hybrid}=2*size_{data}*ifm*input_{w}*input_{h}*mb_{group}+size_{data}*ofm*ifm*kernel_{w}*kernel_{h}*(2-overlap)*(G/N) for G>1G>1 2∗sizedata∗ifm∗inputw∗inputh∗minibatch2*size_{data}*ifm*input_{w}*input_{h}*minibatch for G=1G=1

We can find the minimum value for the total communication volume by differentiating the expression for overall communication volume over G, and then solving for d(commshybrid)/d(G)=0d(comms_{hybrid})/d(G)=0. For a fully connected layer, and assuming FP32 data-type (and no overlap), we can find the optimal point by solving:

d(8∗ifm∗(minibatch/G+ofm∗G/N))/d(G)=0d(8*ifm*(minibatch/G+ofm*G/N))/d(G)=0 Or, −minibatch/G2+ofm/N=0-minibatch/G^{2}+ofm/N=0, hence G=(N∗minibatch/ofm)G=\sqrt{(}N*minibatch/ofm) or G=1G=1. For a layer with ofm=4096, minibatch=256, N=64, we have G = 3, and the communication volume is: 8∗ifm∗2138*ifm*213. For the pure model parallelism case of G=1, the communication volume is 8∗ifm∗2568*ifm*256. Clearly hybrid parallelism offers better overall communication volume than data- or model-parallelism.

As far as parallelism is concerned we believe that hybrid parallelism (with data- and model-parallelism as special cases) is often sufficient approach to parallelizing neural network computation. Indeed one can argue that while hybrid parallelism partitions work along both the minibatch and features (possibly 2) dimensions, partitioning work across other tensor dimensions will always be sub-optimal and should be used only if enough parallelism is not afforded by partitioning across these dimensions.

4 Deep Learning Communication Primitives

The hybrid parallel approach can be implemented using two simple multi-node data transfer operations which we call: part-reduce and part-broadcast. Given a group of nodes NgN_{g}, and a tensor τ\tau, the part-reduce operation performs reduction over partial τ\tau computed locally on each node of NgN_{g} and then scatters the reduced τ\tau to all the nodes in NgN_{g}. The two steps of reduction and scatter can be fused and best represented by MPI_Reduce_scatter() as shown in Figure 1. In data parallelism, this primitive is required in between local weight gradient computation and performing SGD to calculate updated weights and in model parallelism, this primitive is required in forward propagation before partial activation of one layer is used as input for next layer. Similarly, in part-broadcast operation each node of NgN_{g} broadcasts its locally owned strip of τ\tau to all other nodes of NgN_{g} and can be best represented by MPI_Allgather() operation as shown in Figure 2. This primitive is used to populate updated weights to all the nodes of NgN_{g} after performing distributed SGD in the back propagation phase. It is also used to construct full gradient with respect to inputs in the back pass of model parallelism.

PCL-DNN Software Framework

The PCL-DNN software framework consists of three primary modules: data handling, optimized compute library of core CNN/DNN functions optimized for x86 and an optimized MPI-based communications library to enable PCL-DNN execute on a large-scale distributed system. The data handling module functions as the data layer in our framework, ensuring a continuous stream of input data (e.g., images, speech) to the optimized compute library. This library drives the training or scoring process for a given application, as required. It executes the various computations – forward propagation, backpropagation and weight updates – on the underlying hardware.

The main role of the data handling module is to pre-process input data that the compute library uses for its functions. An important requirement we place on this module is that it must not become the bottleneck to the overall throughput of the framework, either for training or classification. That is, it must ensure continuous availability of pre-processed data to the compute library. Further, it must not compete with the latter for hardware resources (i.e., threads and cores) to perform its job. To meet these requirements, we make two design choices: one, the data handling module executes on a dedicated hardware thread and two, we access the disk via the Linux File I/O interfaces, relying on Lustre File System (LFS) to provide high-performance disk access.

The compute library processes the input data in a layer-by-layer manner according to the specified network topology. It consists of and executes primitives for convolutional and fully-connected operations with high-efficiency on the x86 architecture using AVX2 vector instructions. The library supports convolutions with filters of various sizes (e.g., 3x3, 5x5, 11x11) with different strides. To perform the compute associated with fully-connected layers, the optimized library implements highly efficient block-SGEMM functions as well as data layout transformations.

The optimized communications library, similar to the data handling module, executes on a dedicated thread, thereby ensuring clear separation between compute and communication resources on the underlying hardware. Similar to the data handling module, it is critical that the communications library not be a bottleneck to the compute library. To this end, it performs two critical tasks: one, overlap communications during backpropagation due to hybrid parallelism with compute in the forward propagation, including message reordering and ensuring contiguous access to message buffers; two, provide a lock-free command queue (Vaidyanathan et al., 2015) that enables the compute library to submit communication commands in a non-blocking manner (i.e., submit-and-forget).

Experimental Results

In this section we present performance results. We use the Cori phase I system for the for both single and multi-node experiments. This is a Cray XC machine with 1630 Intel Xeon E5-2698v3 HSW dual socket CPUs with 16 cores (supporting upto 32 threads) per socket and 128 GB of memory per node. This has Cray Aries high speed ”dragonfly” topology interconnect.

We report the single-node performance of the PCL-DNN framework for the three topologies listed above. We show results for both scoring (labeled FP in the figure) and training (FP+BP in the figure) across five minibatch sizes, 16, 32, 64, 128 and 256. From Figure 3, we observe that our framework delivers approximately 315 and 95 images/s for Overfeat-FAST and VGG-A, respectively, for scoring; for training, the throughputs are approximately 90 and 30, respectively. We also observe that the throughput PCL-DNN delivers across minibatch sizes remains nearly the same for the largest topology VGG-A for both scoring and training. For OverFeat, which is approximately 3x smaller than VGG-A, we observe that the throughput for the smaller minibatch sizes (i.e., 16 and 32) is lower, than that of the largest minibatch (256) for training, whereas there is no significant variation in the throughput for scoring. We attribute lower training throughput for smaller minibatch sizes to load imbalance on the system. These results clearly indicate that PCL-DNN executes with high-efficiency on the x86 architecture even for small minibatch sizes – a fact that is critical to achieve high scaling efficiency across a large cluster.

2 Scaling Results

PCL-DNN delivers best till date scaling performance for deep learning applications, without using specialized HW (accelerators, NW elements, fabrics …). We present scaling performance for the VGG-A on the NERSC Cori phase I cluster running upto 128 nodes. Fig 4 show VGG-A performance as we scale from 1 to 128 CPU nodes on the NERSC Cori machine. For minibatch sizes of 256 and 512, PCL-DNN scales almost linearly with the number of node. For 512 minibatch, on 128 nodes PCL-DNN scales by 90X90X, with a throughput of 2510 images per second this corresponds to scaling efficiency of 70%70\% at 128 nodes. This significantly brings down the training time, to under 10 minutes per epoch for the Imagenet-1K dataset(Deng et al., 2009). Further we show that even for a smaller minibatch of 256, PCL-DNN scaled well upto 64 nodes with a efficiency of 82%82\%. Fig 5 shows the training of the VGG-A network on 32 and 64 nodes. Since we parallelize SGD retaining its synchronous nature, and there are no hyperparameter changes, the convergence of the distributed algorithm is identical to the single node version, and the behavior for different distributed versions is identical. We see this in Fig 5, where the Top5 validation and training accuracy overlap for both 32 node and 64 node runs.

3 Scaling on Cloud

This section presents the performance and scaling of PCL-DNN on AWS EC2. The purpose of these experiments is to demonstrate the applicability of our performance optimizations in multi-tenant cloud environments that are not as high performing as dedicated HPC clusters. Our experiments were performed on a cluster of 16-node c4.x8large RHEL 7.1 instances. Each instance consists of 2-socket 9-core 2.9 GHz Xeon E5-2666 v3 with 60 GB memory. Instances are connected via 10 Gigabit Ethernet. Unlike Cori, a dedicated bare-metal HPC cluster, CPU and network resources in EC2 are virtualized incurring higher overheads.

We enabled SR-IOV network virtualization support (aka ’enhanced networking’ by AWS). Further, we dedicated a core for handling interrupt requests for the network transmit and receive queues, resulting in 30%-40% better network performance.

Figure 6 presents the performance (in images/second) and scaling of PCL-DNN on AWS EC2 for Overfeat and VGG-A topologies with mini-batch size of 256. We observe 1027 and 397 images per second on 16 nodes for Overfeat and VGG-A respectively. This translates to 11.9x and 14.2x speedup on 16 nodes respectively. We observe better speedups for VGG-A given its higher flops per network byte requirements.

4 Automatic Speech Recognition

We examine DNNs in the Automatic Speech recognition (ASR) context, using Context-Dependent Deep Neural Networks (CD-DNN) HMMs (Seide et al., 2011). CD-DNN-HMMs combine classic artificial-neural-network HMMs with traditional tied-state triphones and gives a 33% relative WER reduction over a discriminatively trained GMM-HMM on Hub5’00 switchboard dataset. This network consists of 7 fully connected hidden layers each with 2048 neurons. For this network a detailed performance study (Seide et al., 2014b), compares performance for both on single node and scaling. For the CD-DNN network, on a Xeon E5-2697v3 HSW CPU (with 14x2 cores with 1.7 TFLOPS/s SP peak) PCL-DNN delivers 4600 frames/s. This is 4X4X better than best reported CPU performance and the 2-node performance with PCL-DNN betters the performance reported for an 80-node Xeon cluster (Seide et al., 2014b). Performance scaling of the CD-DNN network is shown in fig.7. This surpasses the 3-card K20x performance (Seide et al., 2014b) on an 4 HSW nodes (delivering 13K frames/sec), and continues to scale with a performance of 29.5K frames/sec on 16-nodes.Scaling DNN is far more challenging than the CNNs presented earlier, owing to higher communication to compute ratios and these results show that PCL-DNN delivers best-in-class performance on CPUs and demonstrating the generality of this optimized framework.

Conclusions

We demonstrate that deep learning training can be performed at scale using synchronous SGD at high thoughput on CPUs. We extend the state of the art in terms of scaling and time-to-solution, and present a detailed insights and analysis of multinode training.

Acknowledgements

The authors would like to acknowledge support from Prabhat Kumar, Wahid Bhimji, and Lisa Gerhardt of the National Energy Research Scientific Computing Center for helping with access and support for the Cori supercomputer.

References