\section{Introduction}
\label{sec:introduction}
\noindent
Deep neural networks (DNNs) have revitalized the field of machine learning by achieving accuracy levels beyond human perception in a variety of image, video, text and speech processing tasks~\cite{resNet,vgg16,googleNet,Venugopalan:2015,Hannun:2014,Zhang:2015,Hinton:2012,Li:2016}. Today, they have emerged as the \emph{de facto} solution approach for challenging artificial intelligence problems and are deployed in several products and services, including Google's image and voice search~\cite{Google:2013}, Facebook's DeepFace~\cite{deepface} and DeepText~\cite{Facebook:2016}, Apple's Siri~\cite{Apple:2016}, to name a few. However, one of the key challenges to the ubiquitous adoption of DNNs is their computational demand, which far outstrips the capabilities of modern computing platforms. For example, ResNet-152~\cite{resNet}, a state-of-the-art DNN for image recognition, requires $\sim$11.3 giga operations to classify a single 224$\times$224 image. As DNNs begin to pervade the spectrum of computing devices from high performance servers to low-power edge/IoT devices, the demands on compute efficiency are expected grow even more stringent. Consequently, there is an urgent need to explore new approaches to improve their implementation efficiency.

Realizing this need, prior research efforts have explored several key directions. A majority of research efforts exploit the compute and communication characteristics of DNNs to efficiently parallelize them on distributed platforms~\cite{Krizhevsky:2014,Iandola:2015,Dean:2012,Rhu:2016,Das:2016,Zlateski:2016,Nervana:2016} or design specialized accelerator architectures~\cite{Farabet:2011,Chakradhar:2010,Chen:2014,Chen:2014:1,Venkataramani:2017,Jouppi:2017,Reagen:2016,Venkataramani:2013,Eldridge:2015,Majumdar:2012,Gokhale:2014,eyeriss}. Other efforts leverage the intrinsic error resilience of DNNs to approximate selected computations~\cite{Gupta:2015,Venkataramani:2014,Zhu:2016,Courbariaux:2015,deepComp,Jaderberg:2014,Seide:2014,Zhang:2015:1} or utilize non-CMOS technologies whose device characteristics match the compute kernels of DNNs~\cite{Liu:2015,Ramasubramanian:2014,Shafiee:2016}.

Complementary to the above efforts, we focus on a different design scenario and target deeply embedded systems such as wearables and IoT edge devices where the additional area/ cost imposed by custom-accelerators is prohibitive. For example, even a low power DNN accelerator such as Eyeriss~\cite{eyeriss} occupies $12.25mm^2$ area which is 30$\times$ larger than the $0.4mm^2$ occupied by an ultra low power Cortex A35 core~\cite{arm}. Accelerating DNNs on these low-power systems then comes down to accelerating them on the indispensable general-purpose processor (GPP) cores already present in these systems. We focus on improving the execution time of DNNs on GPPs by leveraging \emph{sparsity} in different DNN data-structures, \emph{viz.}, features, weights and backpropagated errors. Sparsity in DNNs can be both static or dynamic depending on whether the zero values remain constant or vary across different inputs to the network. Sparsity in weights, introduced by pruning connections in the network after training, is static in nature. In contrast, feature and error sparsities, caused by the thresholding nature of the ReLU (Rectified Linear Unit) activation functions, are dynamic in nature.

Prior efforts that exploit sparsity to accelerate DNNs can be grouped into two classes based on whether they exploit static or dynamic sparsity, as shown in Figure~\ref{fig:relwork}. Specialized accelerators~\cite{cnvlutin, eyeriss, eie, scnn, cambriconx, Park:2016} have been proposed to exploit both forms of sparsity. However, these techniques are closely tied to a given accelerator architecture thereby preventing their applicability to GPPs. In contrast, static sparsity in weights can be exploited on GPPs through software-only approaches~\cite{BLiu:2015,Wen:2016,Yu:2017} that involve sparse encodings and sparse matrix multiplication routines. However, extending them to exploit the intermediate levels of dynamic sparsity (40\%-70\%) can be counter-productive because of the encoding overhead involved at runtime. For example, the sparse encoding overhead involved in performing sparseBLAS based matrix multiplication on the sparse conv3 layer features of AlexNet~\cite{alexNet} causes a slowdown of 2.78$\times$.
We observe across 6 image-recognition DNNs that dynamic sparsity in features results in 45.1\% of the computations being rendered redundant during inference, presenting a significant opportunity for improving performance. We believe our effort, \sagpp, is the first to successfully exploit dynamic sparsity in DNNs on GPPs. To that end, we propose micro-architectural and ISA extensions that are minimally intrusive and low-overhead. Moreover, since static sparsity is a special case of dynamic sparsity where the location of zeros remains constant across inputs, the extensions also allow \sagpp~to exploit static sparsity in weights.

\begin{figure}[htb]
	\centering
	\includegraphics[clip, width=1.0\columnwidth]{relatedWorkFig.pdf}
	\caption{Related work: Exploiting sparsity in DNNs}
	\label{fig:relwork}
\end{figure}
To exploit dynamic sparsity, GPPs need to be equipped to dynamically detect if the result of an instruction ({\em e.g.}, a load from a sparse data structure) is zero and if so, skip a set of future instructions that use its result. However, there are two key challenges: (i) the instructions to be skipped often do not immediately follow the instruction that determines whether they may be; moreover, the instructions to be skipped may be non-contiguous in the program, and (ii) to maximize performance benefits, the instructions to be skipped should be \emph{prevented from even being fetched}, as squashing the instruction in-flight would diminish the benefits by introducing a pipeline stall.

 \vspace{5pt}
{\bf \noindent Sparsity aware Core Extensions (\sagpp).} To overcome the aforementioned challenges, we propose \underline{Spar}sity aware \underline{C}ore \underline{E}xtensions (\sagpp), comprised of two key micro-architectural enhancements. First, \sagpp~contains a \emph{Sparsity Register File} (SpRF) to track which general purpose registers that contain a zero. We achieve this by augmenting the writeback stage of the processor to check if the update to a register is zero and appropriately modify the corresponding SpRF entry. Next, a \emph{Sparsity Aware Skip Address} (SASA) table is used to store instruction sequences and the conditions under which they can be skipped \emph{i.e.,} the registers in the SpRF that need to be zero for the instructions to become redundant. Whenever \sagpp~fetches an instruction, it uses the SASA table and the SpRF to \emph{pre-identify} whether the following instruction(s) can be skipped. If so, the program counter is modified to directly fetch the next irredundant instruction.

In summary, the key contributions of this work are:
\begin{itemize}
\item We propose Sparsity aware Core Extensions (\sagpp) to accelerate DNNs on general purpose processors by skipping redundant computations borne out of sparsity in the different DNN data-structures.
\item \sagpp~comprises of micro-architectural enhancements to dynamically track when the result of an instruction is zero, pre-identify future instructions that are rendered redundant, and prevent them from being fetched and executed, thereby improving performance.
\item We evaluate \sagpp~on a suite of 6 state-of-the-art image-recognition DNN benchmarks using the Caffe deep learning framework. We achieve application-level speedups of 19\%-31\%
 on low-power embedded scalar microprocessors. We also achieve speedups of 8\%-15\%
 over highly optimized baseline implementations that use OpenBLAS on an ARMv8 processor with 4-way SIMD and prefetching support.
\end{itemize}

The rest of the paper is organized as follows. Section~\ref{sec:motivation} explains the sources of sparsity in DNNs. Section~\ref{sec:design} details the key design principles of \sagpp~and demonstrates them in the context of an in-order pipelined processor. Section~\ref{sec:codeGen} shows \sagpp~in action using the ARM-BLAS GEMM routine as a case study. Section~\ref{sec:exptMeth} describes the experimental methodology and the results are presented in Section~\ref{sec:results}. Section~\ref{sec:relatedwork} presents related research efforts and Section~\ref{sec:conclusion} concludes the paper.

\vspace*{-0pt}
\section{Sparsity in DNNs: Sources and Opportunity}
\label{sec:motivation}
\noindent
In this section, we first provide a brief background on DNNs. We then explain the various static and dynamic sources of sparsity in the different DNN data-structures, and quantify the opportunity for performance improvement afforded by sparsity.

\begin{figure*}[htb]
  \centering

  \includegraphics[clip,width=\textwidth]{sparsity.pdf}
  \vspace*{-20pt}
  \caption{Different forms of sparsity in DNNs}
  \label{fig:src_sparsity}

\end{figure*}
\subsection{DNN: Background}
DNNs are networks of primitive compute units called \emph{neurons}, organized into layers. Each layer is associated with a set of parameters called weights. DNNs operate in two phases \emph{viz.} training and inference. During the training phase, a labelled training dataset is used to iteratively refine the weights of the DNN. In the inference phase, the trained DNN is used to classify new inputs.

Computationally, DNN executions iteratively perform 3 key steps \emph{viz.} Forward Propagation (FP), Backpropagation (BP), and Weight Gradient and Update (WG), that operate on 4 data-structures, \emph{viz.}, features, weights, errors and gradients. All three steps (FP, BP, WG) are performed during the training phase, while inference involves only the FP step. In FP, inputs to the DNN are propagated through its layers to produce the DNN outputs. In each layer, the input features are operated on with its weights to produce its output features, which are then fed to the next layer and so on. In BP, errors observed at the output of the DNN are propagated backwards through each layer of the DNN. In this case, the error at the output of a layer is operated on with its weights to compute the error at its inputs. In WG, the input features and the output error of each layer are used to refine its weights.

\subsection{Sources of Sparsity in DNNs}
In practice, all major DNN data-structures - features, weights, errors and gradients - exhibit significant levels of sparsity, which can be exploited for computational savings. Three of the four data-structures (all except weight gradients) are used as inputs to multiply-and-accumulate operations in the different steps (FP/BP/WG), which become redundant when one of the input operands is zero. Among these three sparse multiply-and-accumulate operands, weights exhibit static sparsity that remains constant across different inputs while the remaining two data-structures (features and errors) exhibit dynamic sparsity that varies dynamically across different inputs. Figure~\ref{fig:src_sparsity} summarizes the sparsity in the different DNN data-structures, which we describe in the remainder of this subsection.

\subsubsection{Static Sparsity}

\textbf{Weight Sparsity.} Sparsity in weights occurs during the inference phase of the DNN. As shown in Figure~\ref{fig:src_sparsity}, after training, connections whose weights are close to zero are pruned to compress the model size~\cite{deepComp,Reagen:2016,Jaderberg:2014}. The last column in Figure~\ref{fig:src_sparsity} shows the fraction of zero weights in the different layers of the AlexNet model trained using deep compression~\cite{deepComp}. We find the sparsity to vary between 18\%-85\% across the different layers. Weight sparsity is static  in nature because of the fact that zero weights are identified before the inference phase.

\subsubsection{Dynamic Sparsity}

\textbf{Feature sparsity.} Sparsity in features stems from the thresholding nature of the activation function present at the output of each layer. As shown in Figure~\ref{fig:src_sparsity}, the predominantly used ReLU (Rectified Linear Unit) activation function clips negative inputs to zero. Figure~\ref{fig:src_sparsity} also shows the average feature sparsity exhibited by various DNN benchmarks, which ranges between $\sim$25\%-60\%. Feature sparsity has a \emph{dynamic} nature \emph{i.e.,} because neurons whose outputs are zero vary considerably across inputs. Figure~\ref{fig:ft_sparsity2} illustrates this property using feature maps obtained at the conv3 of AlexNet, produced by two different inputs from the ImageNet dataset. We find the spatial variation in white pixels, which indicate zero output, to be substantial across inputs.

\begin{figure}[htb]
  \centering
  \vspace*{-5pt}
  \includegraphics[clip,width=0.9\columnwidth]{featureSparsity2.pdf}

  \caption{Variation in feature sparsity of AlexNet CONV3 layer across input images}
  \label{fig:ft_sparsity2}

\end{figure}

\textbf{Error Sparsity.} Sparsity in the error data-structure originates from two sources. First, the derivative of the activation function, such as ReLU, is zero when the error at the output of the layer is negative. Next, when errors are propagated back through a max-pooling layer, as shown in Equation~\ref{eq:maxpool}, only one input of each pooling window is set a non-zero error value.
\begin{equation}
\label{eq:maxpool}
\frac{\partial E}{\partial y_{l}} (x+p, y+q)= \begin{cases}
0, & \text{if } y_{l+1}(x,y)\neq y_{l}(x+p,y+q)\\
\frac{\partial E}{\partial y_{l+1}}, & \text{otherwise}
\end{cases}
\end{equation}

For example, if a pooling window of size 2$\times$2 is used, at least three quarters of the error values are sparse.
Similar to feature sparsity, the error sparsity is also dynamic. The average error sparsity in different layers of a CIFAR-10 DNN is shown in Figure~\ref{fig:src_sparsity}.

\subsection{Opportunity for Computational Savings}
\begin{wrapfigure}{r}{0.55\columnwidth}
  \centering
  \vspace*{-5pt}
  \includegraphics[clip,width=0.55\columnwidth]{savingsVsBenchmarks.pdf}
  \vspace*{-17pt}
  \caption{Redundant ops}
  \label{fig:bm_opportunity}
  \vspace*{-10pt}
\end{wrapfigure}
Figure~\ref{fig:bm_opportunity} shows the fraction of multiply-accumulate (MAC) computations that are rendered redundant for each DNN benchmark due to dynamic sparsity in features during inference. We find that between 25\%-60\% (average: 45\%) of the computations can be skipped, underscoring the substantial opportunity for performance improvement. Figure~\ref{fig:inp_var_sparsity} shows how the fraction of redundant operations varies across 1000 different inputs for the AlexNet DNN. We observe $\sim$14\% variation across inputs, although each input shows considerable opportunity for reduction in execution time (minimum: 28\%).

\begin{figure}[htb]
  \centering
  \vspace*{-10pt}
  \includegraphics[clip,width=0.9\columnwidth]{savingsVsInput}

  \caption{Fraction of redundant ops across different inputs of AlexNet}
  \label{fig:inp_var_sparsity}

\end{figure}

In summary, the dynamic sparsity present in the feature and error data-structures offers a substantial opportunity to accelerate DNNs. However, the levels of sparsity are not extreme enough to completely exploit them in software, and this coupled with their dynamic nature necessitates hardware solutions to realize benefits in the context of general purpose processors.

\section{SparCE: Sparsity Aware General \\ Purpose Core Extensions}
\label{sec:design}
\noindent
To exploit the different forms of sparsity and improve DNN performance on GPPs, we propose
Sparsity aware Core Extensions (\sagpp), a set of micro-architectural and ISA extensions that are general-purpose, minimally intrusive and low-overhead. In this section, we present the key ideas behind \sagpp~and describe how they can be integrated within an in-order processor pipeline.

\subsection{Challenges}
The key challenge in exploiting sparsity is to equip the processor with the ability to dynamically detect if the result of an instruction is zero and if so, skip a list of future instructions that are rendered redundant. We illustrate this challenge using the assembly code snippet shown in Figure~\ref{fig:loop}, which computes the dot-product of two vectors, $INP$ and $KER$, each of size $N$, to produce a scalar $OUT$. Registers $r0$, $r1$ and $r2$ hold the data operands, while $p0$, $p1$ and $p2$ are pointers that hold their respective memory locations. For each instruction in the program, Figure~\ref{fig:loop} shows the instructions that can be skipped when its result is zero. For example, when the $INP$ load returns a zero (Inst. 2), the subsequent $KER$ load (Inst. 4), and the multiply and add instructions can be skipped (Insts. 6-7). It is noteworthy that the computational savings is a weighted sum of the number of instructions skipped and the cycles taken by each instruction. For instance, floating point multiply and add instructions may take 3-5 cycles to execute, while a load incurs variable cycles depending on the level of cache hierarchy accessed.

\begin{figure}[htb]
  \centering
  \vspace*{-5pt}
  \includegraphics[clip,width=1.0\columnwidth]{simpleLoop.pdf}

  \caption{Redundant instructions due to sparsity in vector dot-product evaluation}
  \label{fig:loop}

\end{figure}

The following observations highlight the challenges in detecting and benefiting from sparsity.
\begin{itemize}
\item \textbf{Location of redundant instructions.} In a program, the instructions that can be skipped may not immediately follow the instruction producing the zero result. Worse, redundant instruction sequences may be scattered non-contiguously through the program. For instance, in the program shown in Figure~\ref{fig:loop}, when inst. 2 returns a zero, 2 non-contiguous instruction sequences (Inst. 4 and Insts. 6-7) need to skipped. Hence, efficient ways to \emph{capture which instructions can be skipped} is key to leveraging sparsity.
\item \textbf{Avoiding redundant instruction fetches.} To maximize performance, the \emph{instructions that can be skipped need to be prevented from even being fetched}. This is especially critical in the context of in-order pipelines, where even if the instruction is squashed after being fetched, it would introduce a bubble in the pipeline. For example, if the condition for $r0$ or $r1$ being zero is checked after the $FMUL$ instruction is fetched (Inst. 6), it would result in a bubble flowing through the pipeline in place of Inst. 6. It is worth noting that, in the context of multi-cycle instructions, squashing the instruction after it is fetched can still improve performance.
\item \textbf{Use of SIMD instructions.} GPPs use vector units with SIMD (Single Instruction, Multiple Data) execution engines to exploit fine-grained data parallelism in workloads. SIMD instructions can be skipped only if computations performed on all the vector lanes are redundant. This constrains the sparsity to be relatively coarse grained, as irregularly scattered zero values are not beneficial.
\end{itemize}

\subsection{SparCE: Overview}
Figure~\ref{fig:sagpp_overview} shows an overview of the micro-architectural and ISA enhancements proposed in \sagpp. We describe these extensions in detail, and demonstrate how they address the aforementioned challenges to leverage sparsity in DNNs.

\begin{figure}[htb]
  \centering
  \vspace{-5pt}
  \includegraphics[clip,width=\columnwidth]{overview1}
  \vspace{-20pt}
  \caption{\sagpp: Design Overview}
  \label{fig:sagpp_overview}

\end{figure}

\subsubsection{Micro-architectural states and ISA Extension}
In \sagpp, we augment the processor with two new micro-architectural states \emph{viz.} Sparsity Register File (SpRF) and Sparsity Aware Skip Address (SASA) table. The SpRF is used to dynamically track which registers in the processor's register file contain zero values. The SpRF contains one entry (few bits) corresponding to each register in the register file. When an instruction that writes to a register retires, the SpRF is updated if the result is zero. The SASA table stores information about which instructions can be skipped and under what conditions. Specifically, each entry stores the program counter ($PC$) of the instruction preceding a redundant instruction sequence, the index of the register that determines redundancy and the length of the sequence. Storing the $PC$ of the instruction preceding a redundant instruction sequence allows \sagpp~to \emph{pre-identify whether the next set of instructions can be skipped and if so, skip them before the instructions are fetched}.

\vspace{3pt}
{\bf \noindent SASA-LD Instruction.} \sagpp~enables software to explicitly identify potentially redundant instruction regions by pre-loading the SASA table at program startup, or before the program execution enters a given code region. To this end, we extend the ISA with a new instruction \emph{viz.} $SASA$-$LD$, which loads a given region of memory into the SASA table. As shown in Equation~\ref{eq:newISA},  the $SASA$-$LD$ instruction takes a register operand ($Rn$) that points to the SASA table's location in memory and an immediate operand ($size$) that denotes the size of the SASA table.
\begin{equation}
  \centering

  \label{eq:newISA}
  SASA\text{-}LD \quad [Rn], \#size

\end{equation}
At a given point in the program execution, the number of entries in the SASA table limits the number of redundant instruction regions that can be skipped by \sagpp. However, the SASA table can be periodically refreshed from memory as the program execution progresses. In the context of DNNs, we found that 20 entries in the SASA table suffice to capture all redundant instruction sequences, since the computational kernels are captured by a small number of library ({\em e.g.}, BLAS) functions.

\subsubsection{Tracking, Pre-identifying, and Skipping Redundant Instructions}
\sagpp~utilizes the SpRF and the SASA table to dynamically skip redundant instructions borne out of sparsity in the input data-structures. As shown in Figure~\ref{fig:sagpp_overview}, the micro-architecture of \sagpp~is extended to support the following functions.

\vspace{3pt}
{\bf \noindent Track Sparse Registers.} \sagpp~contains a Sparse Value Checker (\emph{SVC}), which in the processor's writeback stage compares the result of each instruction to zero and if so, updates the entry corresponding to the instruction's destination register in the SpRF.

\vspace{3pt}
{\bf \noindent Pre-identify \& Skip Redundant Instructions.}  \sagpp~is equipped with a Pre-identify and Skip Redundancy Unit (PSRU) that utilizes the SASA table to identify and skip redundant instruction regions. For each instruction, we check if its $PC$ contains an entry in the SASA table. An entry in the SASA table indicates that the instruction following the current instruction is the start of a potentially redundant instruction sequence. In this case, the PSRU checks the SpRF to identify if the registers indicated in the SASA table entry are currently zero. If so, it increments the $PC$ to the end of the redundant instruction sequence, thereby skipping instructions to benefit performance. If not, \sagpp~proceeds to execute instructions in program order.

In summary, \sagpp~uses the SpRF and the SASA table to seamlessly track sparse registers, pre-identify instruction sequences that are redundant and dynamically skip them before they are even fetched to improve performance.

\begin{figure*}[htb]
  \centering

  \includegraphics[clip,width=\textwidth]{pipeline1.pdf}
  \vspace*{-20pt}
  \caption{Block diagram of \sagpp~in-order processor architecture}
  \label{fig:pipeline}
  \vspace*{-15pt}
\end{figure*}

\subsection{In-order SparCE Processor Pipeline} \label{subsec:sapp}
We now explain how \sagpp~can be integrated into an in-order processor. Figure~\ref{fig:pipeline} shows the block diagram of the overall \sagpp~processor architecture. We start with a conventional 4-stage (fetch, decode, execute/memory, and writeback) pipelined processor architecture implementing a RISC-style instruction set with at most 2 source register operands and one destination register operand. Although the \sagpp~architecture is described in this section with a scalar execution unit for ease of illustration, it is directly applicable to vector processors with any number of SIMD execution lanes. We augment the processor with the following structures.

\vspace{5pt}
\noindent \textbf{Sparsity Register File (SpRF).} The SpRF is located at the fetch stage of the \sagpp~processor. It is a multi-ported register file, with one entry corresponding to each register in the processor's register file. Each entry in the SpRF contains only two single-bit fields - \emph{isSparse} and \emph{regUpdInFlight}. The \emph{isSparse} bit is set to 1 for registers containing zeros, and reset otherwise. The \emph{regUpdInFlight} bit indicates whether an instruction modifying the register is in flight in any stage of the pipeline. For instance, an instruction modifying register $Rd$, sets the $SpRF[Rd][regUpdInFlight]$ field when it enters the decode stage and resets it after committing its result in the writeback stage. In determining whether a future instruction is redundant, the \emph{regUpdInFlight} field ensures that we do not use a stale \emph{isSparse} value when a more recent instruction updating this register is under execution.

\vspace{5pt}
\noindent \textbf{Sparse Value Checker (SVC).} The SVC is added to the writeback stage of the \sagpp~processor. It contains a comparator that checks if the output of each instruction that updates a register is zero. It then correspondingly updates the SpRF. For example, when the output of an instruction with destination register $Rd$ is zero, then $SpRF[Rd][isSparse]$ is set and $SpRF[Rd][regUpdInFlight]$ is reset.

\vspace{5pt}
\noindent \textbf{Sparsity Aware Skip Address (SASA) table.} The SASA table is also present in the fetch stage of the \sagpp~processor. As shown in Figure~\ref{fig:pipeline}, the SASA table is an associative memory structure with three fields: (i) $preceedingPC$, which stores the $PC$ of the instruction prior to the redundant instruction sequence, (ii) $SpRFCondition$ field which stores a Boolean combination of 2 register indices in the SpRF that should be satisfied for the instruction region to be skipped, and (iii) $instsToSkip$, which contains the length of the redundant instruction sequence. As an example, for the code in Figure~\ref{fig:loop}, the SASA table entry to skip instructions 6-7 would be \{$preceedingPC$=5, $SpRFCondition$=$r0|r1$, $instsToSkip$=2\}

\begin{figure}[htb]
  \centering

  \includegraphics[clip,width=0.8\columnwidth]{flowchart1.pdf}

  \caption{Flowchart for pre-identifying and skipping redundancy}
  \label{fig:flowchart}

\end{figure}

\vspace{5pt}
\noindent \textbf{Pre-identify and Skip Redundancy Unit (PSRU)} The PSRU is also added to the fetch stage of \sagpp~processor. The PSRU determines whether an instruction region specified in the SASA table can be skipped, and appropriately modulates the $PC$. Figure~\ref{fig:flowchart} shows the flowchart depicting the operation of PSRU. Given a $PC$, an associative lookup is performed on the $preceedingPC$ field of the SASA table. If the $PC$ is a hit in the SASA table, then the next instruction marks the beginning of a potentially redundant instruction sequence. To ascertain if the instruction sequence can be skipped, the PSRU reads the register indices specified in the $SpRFCondition$ field of the SASA table (say $R_{s1}$ and $R_{s2}$) from the SpRF. If neither $SpRF[R_{s1}]$ nor $SpRF[R_{s2}]$ have their $regUpdInFlight$ field set, then PSRU computes the Boolean condition (specified in the $SpRFCondition$ field) on their respective $isSparse$ fields. If the condition is satisfied, then the instruction region is deemed redundant and the $PC$ is incremented by $instsToSkip$ to point to the instruction immediately following the end of the redundant region. If not, $PC$ is incremented by 1 and the instructions are executed in program order.
However, if the $regUpdInFlight$ field is set for either $SpRF[R_{s1}]$ or $SpRF[R_{s2}]$ and the Boolean condition in $SpRFCondition$ cannot be definitively evaluated, then the instruction region is temporarily marked as a \emph{skippable} region within the PSRU. The $PC$ is incremented by 1 and the program execution is continued. It is worth noting that continuing or aborting execution of a skippable region will \emph{not} affect program functionality.

In the case when a $PC$ is not present in the SASA table, the PSRU checks if the $PC$ is part of an active skippable region. For the registers present in the $SpRFCondition$, the $regUpdInFlight$ fields are re-checked from the SpRF. If they are reset, the Boolean condition in $SpRFCondition$ is evaluated. If the skippable region is determined to be redundant, then the remaining instructions in the region are skipped by appropriately modifying the $PC$. Further, instructions belonging to the skippable region prior to the current $PC$ are squashed if they are still in flight. If the Boolean condition evaluates to a \emph{false}, then the remaining instructions in the skippable region are executed. Finally, for a $PC$ that misses the SASA table and is not part of any active skippable region, the $PC$ is incremented by 1 to fetch the next instruction.

We note that the logic introduced in \sagpp~to pre-identify redundant instructions executes in parallel with the instruction cache (ICache) access. The additional logic does not impact the latency of the fetch stage, as both the SASA table and the SpRF are significantly smaller structures compared to the ICache.

In summary, by using the SpRF and the SASA table, \sagpp~dynamically tracks sparse registers, pre-identifies if an instruction region is redundant and skips instructions before they are even fetched to improve performance. Thus \sagpp~enables DNN acceleration on GPPs by exploiting the sparsity resident in their data-structures.

\section{Software for SparCE \\ Processsors}
\label{sec:codeGen}
\begin{figure*}[htb]
  \centering
  \vspace*{-10pt}
  \includegraphics[clip,width=\textwidth]{BlasCode.pdf}
  \vspace*{-20pt}
  \caption{Zero skipping for \emph{sgemm\_kernel} subroutine in BLAS}
  \label{fig:blas_code}
  \vspace*{-15pt}
\end{figure*}

To extract maximum performance from \sagpp, the software needs to suitably leverage the sparsity-aware micro-architectural extensions. In this section, we outline the key principles behind software design for \sagpp, and demonstrate them in the context of a highly optimized implementation of matrix multiplication (GEMM) from the OpenBLAS library.

\subsection{SparCE Software Design Steps}
The following steps need to be performed in software to leverage sparsity on \sagpp~processors.

\vspace{3pt}
\noindent \textbf{Identifying sparse data-structures and redundant instructions.} One of the key requirements on software is to indicate (through the SASA table) which instruction regions are potentially redundant due to sparsity. To this end, sparse data-structures in the workload need to be first identified by the programmer. Next, when the application is compiled, the registers that hold the sparse data-structures and the instructions that load them from memory are identified. Then, a static dependency analysis of the instruction stream reveals the instructions that are affected by the sparse data-structures, which are then marked as potentially redundant instruction sequences. Following this, the condition for skipping each instruction sequence is derived from: which sparse data-structures affect the region, the registers containing them and the type of operation performed. Finally, the application assembly is instrumented with appropriate $SASA$-$LD$ instructions.

\vspace{3pt}
\noindent \textbf{Separate redundant instructions from instructions causing redundancy.} To identify if an instruction region can be skipped, the instruction whose result makes them redundant should have completed execution. Therefore, in the context of in-order processors, they need to be spaced at least few instructions (3 in our as case as SVC is located in writeback stage) apart in the program, so that redundant instructions are not fetched and no pipeline bubble is introduced due to squashing redundant instructions. To this end, the instruction stream is re-ordered by inserting independent non-redundant instructions where needed to enable sufficient separation.

\vspace{3pt}
\noindent \textbf{Mapping sparse data-structures on vector processors.} In vector processors, typically one of the input operands is shared by all SIMD lanes, while the other is different across lanes. A vector instruction can be skipped only if computations performed on all SIMD lanes are redundant. Hence, it is better to map a sparse data-structure as the shared input operand, since the likelihood of all non-shared inputs being zero is low. If both data-structures are sparse, then the data-structure that exhibits the most block-wise sparsity is mapped as the non-shared input operand.

\subsection{Case Study: SparCE GEMM Routine} \label{sec:sparce_gemm}
Popular deep learning frameworks, such as Caffe~\cite{caffe}, tensor flow~\cite{tensorflow}, \emph{etc.,} execute DNNs as a series of matrix multiply operations, and leverage highly optimized software libraries such as BLAS (Basic Linear Algebra Subprograms) to realize them. Therefore, we develop a \sagpp~version of the GEMM (Generalized Matrix Multiply) routine from the OpenBLAS library~\cite{openblas} targeting an ARM processor. We utilize the \sagpp-GEMM routine to quantify the reduction in DNN execution time achieved on a \sagpp~processor.

Figure~\ref{fig:blas_code} shows the assembly program for the \emph{sgemm\_kernel\_l1\_M16\_22} subroutine, which performs single-precision floating point matrix multiplication ($B$$\times$$A$$=$$C$). The \emph{sgemm} routine executes two smaller subroutines \emph{viz.} \emph{kernel16x4\_M1} (which we abbreviate as $M1$) and \emph{kernel16x4\_M2} ($M2$) sequentially in a loop. The subroutines utilize vector registers $v8$ and $v12$ to hold the first input operand ($B$) and registers $v0$-$v7$ to hold the second operand ($A$). The intermediate results are computed in registers $v16$-$v31$. To optimize performance, each subroutine prefetches data for the other. For example, $M1$ fetches data for operand $A$ into the \emph{odd} registers, which are then used by $M2$ in the subsequent iteration. The memory addresses are provided by scalar registers $pB$ and $pA\_0$-$pA\_3$. The subroutines also prefetch data in the L1-cache using the \emph{prfm} instruction.

We identify redundant instruction sequences in the \emph{sgemm\_kernel\_l1\_M16\_22} subroutine assuming one of the input operands (say $B$) is sparse. The analysis can be easily extended to cover the scenario when either $A$ or both $A$ and $B$ are sparse. Since $B$ is sparse, it is beneficial to map it as the operand shared across the SIMD lanes, which is already the case in Figure~\ref{fig:blas_code}. Even when one of the words in a vector register containing $B$ ($v8$ or $v12$) is zero, 4 $fmla$ instructions can be skipped. For instance, when $v12.s[1]$ equals $0$, instructions 9,10,13,14 in $M2$ can be skipped. This forms 2 redundant instruction sequences (9-10 and 13-14), each of size 2. This amounts to 16 $fmla$s being skipped when the entire vector register ($v8$ or $v12$) is zero. It is worth noting that, had the sparse data-structure ($B$) been mapped as the non-shared SIMD operand in the program, only 4 $fmla$s could be skipped even when the entire vector register is zero.

In addition to the $fmla$ instructions being skipped, when the vector register is fully zero, the load instructions for the second operand can also be skipped. For example, when $v12$ is zero, $ld1$ instructions (7, 11, 15 and 19) for operand $A$ can be skipped in $M1$. Also, note that we did not re-order any instruction in the \emph{sgemm} routine, as the redundant instruction sequences were sufficiently spaced apart from the instruction that triggers their skipping. In the context of $fmla$ instructions, the vector loads happen in a different subroutine owing to pre-fetching, and in the case of the $ld1$ instructions, they were naturally spaced $>$3 instructions apart.

\begin{figure}[htb]
  \centering

  \includegraphics[clip,width=0.8\columnwidth]{SASAtable.pdf}

  \caption{SASA table entries for \emph{kernel16x4\_M1} subroutine}
  \label{fig:sasa_m1}
  \vspace*{-15pt}
\end{figure}

Based on the above analysis, Figure~\ref{fig:sasa_m1} shows the SASA table corresponding to the $M1$ subroutine. We find a total of 12 entries in the SASA table, 8 entries of size 2 for the 16 $fmla$s and 4 entries of size 1 for the $ld1$ instructions.

\vspace{3pt}
{\bf \noindent SparCE in action.} Figure~\ref{fig:blas_action} depicts the sequence of events that leads to instructions being skipped by \sagpp. First, when register $v12$ is loaded by the $M1$ subroutine, its SpRF entry is updated. Note that the $isSparse$ field is a bit vector, one bit for each word in the vector register. Next, when instruction $M2.12$ is fetched, it sees a hit in the SASA table and the reads $SpRF[12]$ to ascertain if the bit corresponding to $SpRF[12][1]$ is sparse. Since that is the case, the $PC$ is incremented to directly fetch $M2.15$.

\begin{figure}[htb]
  \centering

  \includegraphics[clip,width=\columnwidth]{blasAction2.pdf}

  \caption{\sagpp~in action for \emph{sgemm} routine}
  \label{fig:blas_action}
  \vspace*{-10pt}
\end{figure}

In summary, with minimal changes to software, \sagpp~processors can leverage sparsity to improve performance.

\vspace*{-0pt}
\section{Experimental Methodology}
\label{sec:exptMeth}
\noindent
In this section, we present the methodology adopted in our experiments to evaluate \sagpp.

\subsection{Performance Evaluation}
We modeled the micro-architectural extensions proposed in \sagpp~using the cycle-accurate \emph{gem5} architectural simulator~\cite{gem5}. The gem5 simulator was tightly integrated with the popular Caffe~\cite{caffe} deep learning framework, wherein the matrices corresponding to each layer and each input batch was formed in Caffe and fed into the gem5 simulator to perform the matrix computations. The results were fed back to Caffe to form the inputs for the next layer (or input batch), and so on. Figure~\ref{fig:table}(a) shows the gem5 system configuration used in our experiments. All experiments were run in full-system mode. We evaluate \sagpp~by measuring application level execution times under two scenarios. The first scenario targets embedded scalar processors that are present in ultra-low power edge/IoT devices and lack support for high performance libraries. We chose a scalar ARM v8 in-order processor architecture as the baseline. We leverage the modular nature of gem5 to cater to this scenario, wherein we disable support for advanced architectural features such as SIMD processing and pre-fetching in the ARM v8 processor. We then prototyped a direct convolution routine (which we call \emph{Dir-Conv-Scalar}) that does not utilize the disabled features and used it in our experiments. The second scenario targets a reasonably sophisticated mobile processor for which we chose the ARM v8 in-order processor architecture with 4-way SIMD as the baseline. In this case, the matrix computations were realized using the highly-optimized OpenBLAS~\cite{openblas} based GEMM routines described in Section~\ref{sec:sparce_gemm}. We refer to this implementation as \emph{OpenBLAS-SIMD4}.

\subsection{Power and Area Evaluation}
We implemented the hardware extensions of \sagpp~at the Register Transfer Level (RTL) using Verilog HDL and synthesized them to IBM 45nm technology using Synopsys Design Compiler to measure its power and area overheads. For the configuration shown in Figure~\ref{fig:table}(a), the area overhead is 0.4\% of the ARM Cortex A35 core~\cite{arm}. Thus the area overhead of \sagpp~is quite minimal allowing its deployment in the resource-constrained embedded platforms

\begin{figure}[htb]
  \centering
  \includegraphics[clip,width=1.0\columnwidth]{expt_table.pdf}

  \caption{(a) Gem5 simulation parameters (b) Application benchmarks}
  \label{fig:table}

\end{figure}

\subsection{Benchmarks}
Our benchmark suite, listed in Figure~\ref{fig:table}(b), consists of 6 state-of-the-art image-recognition DNNs \emph{viz.} CIFAR-Caffe DNN using the CIFAR-10 dataset~\cite{cifarData}, and AlexNet~\cite{alexNet}, VGG-16~\cite{vgg16}, ResNet-50~\cite{resNet}, GoogleNet~\cite{googleNet} and Deep Compression-AlexNet~\cite{deepComp} using the ImageNet dataset. These benchmarks contained 5-50 layers and took 0.01-15.4 Billion scalar operations to classify an image. We utilized pre-trained models from the Caffe Model Zoo to evaluate \sagpp~in the context of inference. For training, we utilized only the smaller CIFAR-10 benchmark, as training ImageNet models on gem5 was prohibitively time consuming. It is noteworthy that all benchmarks exhibited dynamic sparsity in features and errors, while only Deep Compression-AlexNet exhibited static parsity in weights.

\section{Results}
\label{sec:results}
\noindent
In this section, we present the results of our experiments that highlight the advantages of \sagpp.

\subsection{Performance and Energy Improvement}
Figure~\ref{fig:timeBenefits} shows the normalized execution time benefits of \sagpp~over the baseline processor for both inference and training. In the context of Dir-Conv-Scalar, the reduction in \emph{application runtime} ranges between 19\%-31\% across the benchmarks. In contrast, OpenBLAS-SIMD4 demonstrates benefits in the range of 8\%-15\% reduction in runtime. This is because $fmla$ instructions occupy a much smaller fraction of their runtime, as their execution engines are more sophisticated - multiple SIMD lanes, low floating point instruction latency \emph{etc.} Also, since they  support features such as prefetching where the data is already fetched into the higher levels of the memory subsystem, avoiding redundant data fetches has a less prominent impact on performance.

\begin{figure}[htb]
  \centering
  \includegraphics[clip,width=1.0\columnwidth]{timeBenefits.pdf}

  \caption{Improvement in execution time at the application level}
  \label{fig:timeBenefits}

\end{figure}

Among the benchmarks, the execution time benefits are largely proportional to the amount of sparsity that they exhibit (Figure~\ref{fig:src_sparsity} in Section~\ref{sec:motivation}). The Deep Compression-AlexNet benchmark achieves the most benefits because both its feature and weight data-structures are sparse, as opposed to other benchmarks whose weight data-structure is dense. In the context of training, the backpropagation step achieves more improvement compared to forward propagation. This stems from the fact that the error data-structure is more sparse compared to features.

\vspace{3pt}
{\bf \noindent Execution Time Breakdown.} To better appreciate the improvements achieved by \sagpp, Figure~\ref{fig:time_brk} hierarchically breaks down the application execution time (in the context of both Dir-Conv-Scalar and OpenBlas-SIMD4 implementations) into components that can and cannot be impacted by \sagpp. At the top level, the solid yellow and gray colors represent the execution time fraction that cannot be improved by \sagpp. This is primarily constituted by auxiliary DNN operations such as activation functions, subsampling and others, which represent 1.9\% and 12.2\% of the runtime for Dir-Conv-Scalar and OpenBlas-SIMD4 implementations, respectively. It is noteworthy that although these operations occupy $<$1\% of the total DNN FLOPs, they occupy a substantially larger fraction of the runtime for the OpenBLAS-SIMD4 implementation. This is owed to the fact that they are typically memory-bound (higher Bytes/FLOP ratio), which is further amplified as matrix multiply operations are significantly optimized by the GEMM subroutine. Also, since the inputs to the DNN are typically dense, the first DNN layer exhibits little redundancy. This occupies 14.3\% and 16.9\% of the total runtimes of AlexNet for Dir-Conv-Scalar and OpenBLAS-SIMD4 implementations, respectively. The fraction grows smaller in deeper networks such as ResNet and VGG.

In Figure~\ref{fig:time_brk} the green color bars denote the computations that can accelerated by leveraging sparsity ($\sim$71\%). For Dir-Conv-Scalar implementation, this is limited to 83.6\% of the baseline AlexNet runtime. Since AlexNet contains $\sim$36\% redundant computations (Figure~\ref{fig:bm_opportunity} in Section~\ref{sec:motivation}), \emph{the best case benefits are limited to 29.8\%}, of which \sagpp~achieves 22.3\%. For the OpenBLAS-SIMD4 implementation, the underlying GEMM involves supplementary operations like memory allocate, copy and free operations, which as marked by the dotted patterns and consumes $27\%$ of the total execution time. This constraints the opportunity for \sagpp~to $\sim$44\% of AlexNet runtime as shown by the diagonally hatched portions in Figure~\ref{fig:time_brk}. Since AlexNet contains $\sim$36\% redundant computations (Figure~\ref{fig:bm_opportunity} in Section~\ref{sec:motivation}), \emph{the best case benefits are limited to $\sim$16\%}, of which \sagpp~achieves 12\% improvement as other control operations such as pointer arithmetic, prefetching \emph{etc.} present within the loop body cannot be avoided.

\begin{figure}[htb]
  \centering
  \includegraphics[clip,width=\columnwidth]{timeBreakdown.pdf}
  \vspace{-10pt}
  \caption{Execution time breakdown for AlexNet}
  \label{fig:time_brk}
  \vspace{-5pt}
\end{figure}

{\bf \noindent Layer-wise Benefits.} We now present the layer-wise breakdown of the benefits quantified in terms of the execution time, instructions and data cache (D-Cache) accesses skipped for the convolutional layers of AlexNet. Figure~\ref{fig:brk1} shows the benefits in the context of both SIMD and scalar processor implementations. As shown in Figure~\ref{fig:brk1}, we are able to achieve, on an average, 39.4\% reduction in instruction count and 35.1\% reduction in D-Cache accesses for Dir-Conv implementations on low-power embedded scalar processors. The reduction in instruction count and D-Cache accesses amount to 30.5\% and 14\% respectively for OpenBLAS-SIMD4 implementations. We also observe the benefits are typically larger for layers deeper in the DNN, as they typically exhibit more sparsity.

\begin{figure}[htb]
  \centering
  \includegraphics[clip,width=0.9\columnwidth]{benefitsBreakdown3.pdf}
  \vspace{-10pt}
  \caption{Layer-wise benefits breakdown for AlexNet}
  \label{fig:brk1}

\end{figure}
\vspace{3pt}
{\bf \noindent Energy Benefits.} Power evaluation of the \sagpp~reveals that it consumes 1.74 mW at 1 GHz, which amounts to 1.9\% of the 90 mW power consumed by even the most power-efficient baseline ARM v8 processor, the Cortex A35 processor~\cite{arm}. Accordingly, the execution time benefits translate to benefits in the range of 16.9\%-28.7\% reduction in application-level energy for a Dir-Conv-Scalar implementation. In the context of OpenBLAS-SIMD4 implementations, the reduction in energy ranges between 6.1\%-13.2\% across the benchmarks.

\subsection{Performance Scaling with Sparsity}
We now study how the performance of \sagpp~scales with increasing levels of sparsity. To this end, we considered a matrix multiplication problem ($B$$\times$$A$$=$$C$), wherein the dimensions of the input matrices $B$ and $A$ were 169$\times$3456 and 3456$\times$384 respectively.
We varied the sparsity of the $B$ matrix by constraining the number of zero entries. The location of the zeros and other entries of the matrices were chosen at random. Figure~\ref{fig:varSparsity} shows how the execution time and the fraction of instructions executed varies with sparsity in the context of both Dir-Conv-Scalar and OpenBLAS-SIMD4 implementations. We find that both implementations exhibit strong performance scaling with sparsity, outlining the ability of \sagpp~to efficiently skip computations. We find the number of instructions executed to be larger than ideal (dotted line in Figure~\ref{fig:varSparsity}) due to the presence of control instructions for pointer arithmetic, loop counts \emph{etc.,} in the program, which cannot be skipped. Also, the disparity in the fraction of instructions executed and the resultant execution time benefits is more pronounced for the OpenBLAS-SIMD4 implementation. We attribute this to the intelligent instruction ordering in the GEMM routine utilized in the implementation, wherein computations are aggressively overlapped with data-fetches. Therefore, even if computations are skipped, the improvement in performance is limited by the time taken for the data-fetches.

\begin{figure}[htb]
  \centering
  \includegraphics[clip,width=1.0\columnwidth]{varSparsityResult.pdf}
  \vspace*{-10pt}
  \caption{\sagpp~performance scaling with sparsity}
  \label{fig:varSparsity}

\end{figure}

\subsection{Operand Ordering in SparCE OpenBLAS-SIMD4 Implementations}
As described in Section~\ref{sec:codeGen}, based on the amount of sparsity exhibited by the data-structures, mapping the right data-structure as the shared-SIMD operand can have a considerable impact on performance. In the case of all benchmarks other than Deep Compression-AlexNet, only the feature data-structure is sparse. Therefore, mapping it as the shared-SIMD operand (Matrix $B$ in the \emph{sgemm} subroutine) would yield the best benefits. Figure~\ref{fig:deepComp} shows the performance improvement achieved when operands are ordered in both ways \emph{viz.} Features$\times$Weights and Weights$\times$Features. In the context of AlexNet, we find that mapping features as the shared-SIMD operand yields 1.86$\times$ better benefits ($\sim$12\% \emph{vs.} 6.5\%) compared to mapping the non-sparse weight data-structure.
\begin{wrapfigure}{r}{0.5\columnwidth}
  \centering
  \includegraphics[clip,width=0.5\columnwidth]{deepComp.pdf}

  \caption{Impact of operand ordering on performance}
  \label{fig:deepComp}

\end{wrapfigure}
 For the Deep Compression-AlexNet network, since both feature and weight data-structures are sparse, the disparity in performance due to operand ordering is relatively small ($<$2\%). Even in this case, we find that choosing features as the shared-SIMD operand is beneficial. This is attributed to the fact that some of the weight matrices have high degree of sparsity, and their zero entries are typically clustered. Therefore, using weights as the non-shared SIMD operand has less of an adverse impact on performance compared to using features.

\vspace*{-0pt}
\section{Related Work}
\label{sec:relatedwork}
\noindent
Prior research efforts that target improving the computational efficiency of DNNs can be grouped into the following broad classes.

\vspace{3pt}
\noindent \textbf{Software parallelization on multi-cores and GPUs.} A large number of previous efforts have been directed towards developing techniques for efficient parallelization of DNNs on multi-core servers and GPUs~\cite{Krizhevsky:2014, Iandola:2015, Dean:2012, Rhu:2016, Das:2016, Zlateski:2016, Nervana:2016}. However, the scalability of these techniques is often limited by synchronization and communication bottlenecks.

\vspace{3pt}
\noindent \textbf{Specialized accelerators.} Developing specialized hardware architectures has been an attractive approach to improve the computational efficiency of DNNs. These accelerators~\cite{Chen:2014,Chen:2014:1,Venkataramani:2017,Jouppi:2017,Majumdar:2012,Gokhale:2014,Farabet:2011,eyeriss,Chakradhar:2010, Reagen:2016, Eldridge:2015} utilize specialized processing cores, interconnect network and other hardware-software co-design methodologies to leverage the different forms of compute and data reuse patterns in DNNs.

\vspace{3pt}
\noindent \textbf{Approximate computing.} DNNs and the applications that use them are intrinsically resilient to errors in their underlying computations. Approximate computing techniques, such as low-precision implementations~\cite{Courbariaux:2015,Gupta:2015,Zhu:2016,Venkataramani:2014}, model compression~\cite{deepComp,Jaderberg:2014} and others~\cite{Seide:2014,Zhang:2015:1}, leverage this property to improve the computational efficiency of DNNs.

\vspace{3pt}
\noindent \textbf{Post-CMOS technology.} Post-CMOS technologies such as memristors and spintronics have succeeded in realizing the computational primitives of DNNs through their intrinsic device characteristics. Implementations of small scale DNNs using memristor crossbar arrays~\cite{Liu:2015,Shafiee:2016} and spintronic devices~\cite{Ramasubramanian:2014} have demonstrated significant promise.

All the above efforts are complementary to \sagpp, as unlike above efforts, we attempt to leverage sparsity in the different DNN data-structures to improve efficiency on cost constrained embedded platforms comprising of mainly a GPP core.

\vspace{3pt}
\noindent \textbf{Exploiting sparsity.} Prior efforts that exploit sparsity to improve DNN efficiency can be grouped to two classes based on whether they exploit static sparsity or dynamic sparsity, as shown in Figure~\ref{fig:relwork}. Specialized sparse architectures that are capable of exploiting both static and dynamic sparsity incorporate a variety of compression techniques and zero-skipping schemes to reduce storage requirement and avoid redundant multiplications in accelerators~\cite{cnvlutin, eyeriss, eie, scnn, cambriconx, Park:2016}.

On the other hand, software approaches that exploit weight sparsity on GPPs take advantage of sparse matrix libraries. These sparse libraries usually yield performance improvement only under extreme levels of sparsity ($>$95\%). Since DNNs naturally exhibit sparsity in the range of 40\%-70\%, a few research efforts force more weight sparsity into DNNs using sparse decomposition methods, \emph{etc.}~\cite{BLiu:2015,Wen:2016} or customize the pruning to match the underlying hardware organization~\cite{Yu:2017}. These invariably come at the cost of training overhead or loss of functional accuracy. In contrast, we propose micro-architectural extensions to GPPs that can exploit both dynamic and static sparsity while being effective even under intermediate levels. Thus, \sagpp~is able to exploit feature sparsity in several state-of-the-art dense DNN models as well as both feature and weight sparsity in existing pruned DNN models.

\section{Conclusion}
\label{sec:conclusion}
\noindent
As DNNs pervade the spectrum of computing devices, new approaches to improve their computational efficiency on resource-constrained IoT/ edge devices becomes critical. In this work, we accelerate DNNs on GPPs, which are an indispensable part of IoT/ edge devices, by exploiting sparsity in the different DNN data-structures. To this end, we propose sparsity aware general purpose core extensions (\sagpp) that enable GPPs to efficiently leverage sparsity, while being minimally intrusive and low-overhead. \sagpp~comprises of two key micro-architectural enhancements. First, a Sparsity Register File (SpRF) dynamically tracks zero-valued registers in the processor. Next, a Sparsity Aware Skip Address (SASA) table indicates potentially redundant instruction sequences and the conditions under which they can be skipped. A Pre-identify and Skip Redundancy Unit (PSRU) combines the information from the SpRF and the SASA table to dynamically \emph{pre-identify} if an instruction sequence can be skipped, and if so masks it from being fetched and executed. We evaluate \sagpp~on 6 image-recognition DNNs in the context of both training and inference. Our evaluations reveal that \sagpp~is a promising design that allows us to exploit all forms of static and dynamic sparsity to accelerate DNNs on GPPs.

\newpage
\bibliographystyle{unsrt}


\section{Introduction}
\label{sect:intro}

AS Cam is a main-sequence eclipsing binary star (B8V+B9.5V components), its orbital period is $\approx 3.43$ days, the orbital eccentricity is $e\approx 0.17$, the maximum visual brightness is $\approx 8.57^\mathrm{m}$ (Simbad database \footnote{http://simbad.u-strasbg.fr}). It was found in photographic plates by \citet{strohmeier1968}. \citet{hilditch1969,hilditch1972} conducted photoelectric observations of AS Cam, obtained its radial velocity curve and calculated the system's absolute parameters.

\citet{khaliullin1983} discovered the apsidal motion in AS Cam using the {\it WBVR} photometry. The obtained rate of the periastron movement $\dot\omega_\mathrm{obs}=16^{\circ}$ per century was almost 3 times less than the expected theoretical value $\dot\omega_\mathrm{th} = 44^{\circ}$ per century. This discovery was independently cofirmed by \citet{maloney1991}, \citet{wolf1996}. AS Cam became the second (after DI Her, see \citealp{martynov1980}) eclipsing system in which the apsidal motion proved to be much slowly than it was predicted by the theory. To explain the discrepancy between observational estimations and theoretical calculations of the apsidal motion rate in both binaries (DI Her and AS Cam) different authors introduced different hypothesises \citep{shakura1985,moffat1989,claret1997,claret1998}. Most of their hypothesises were already discussed by \citet{maloney1991,claret1997,claret1998}, so we did not describe them here. \citet{zakharov1988}, \citet{khaliullin1991} explained the observed anomalies in the frames of the classical and relativistic mechanics: the gravitational influence of the third companion on the central binary in case of non-coplanar orbits can slow down the apsidal motion. \citet{borkovits2007} significantly improved this idea and enriched it by numerical and analytical computations.

According to \citet{claret2010} there is still no evidence for the existence of the third companion in DI Her, its light equation cannot be found within existing errors of observational data. A leading idea that explains the system's slow apsidal motion is a non-coplanar axial rotation of the stars, \citet{albrecht2009}, \citet{albrecht2011} observed the Rossiter-McLaughlin effect in DI Her and NY Cep. This model requires very high equatorial velocities (up to 300 km s$^{-1}$). In case of AS Cam \citet{kozyreva1999} found an evidence in favour of the existence of the third body. They obtained high precision light curves of AS Cam in 1992-1996, accumulated 10 primary and 13 secondary minima and used them to find variations of times of minima. Based on these data and times of minima from the literature \citet{kozyreva1999} discovered cyclical in-phase variations of primary and secondary times of minima and explained their result by the influence of the third companion, its orbital period was found to be about 805 days, the eccentricity was $\approx 0.5$, the amplitude of the light time effect was $\approx 0.50$ astronomical units.

After the discovery of the light equation in AS Cam in 1999 a lot of new minima times were obtained. So, there is a possibility to re-test the system on the existence of the third companion with higher precision.

\section{Observations}
\label{sect:Obs}

Photometric observations of AS Cam were conducted in 2002, 2004, and 2017 in Tien Shan observatory, Fesenkov Astrophysical Institute (Kazakhstan). In 2002 and 2004 we used the 50 cm AZT-5 telescope with the photomultiplier tube (PMT) model 79 and {\it V} filter, the comparison star was HD 34463. In 2017 we obtained new CCD observations in {\it B}, {\it V}, and {\it R} filters using the Zeiss-1000 telescope equipped with the Apogee U900 CCD camera. For the latest set of observations we used TYC 4347-452-1 (comparison star) and TYC 4347-682-1 (control star) as reference stars. Usual exposure
times were about 10 seconds.

To process raw CCD data we used the MAXIM-5 program. The aperture was constant during one night, its differences from night to night was not significant. Maximum errors for a single exposure were in the range $0.003^m-0.006^m$ for the different nights. Reference stars were assumed to be constant during observations. For the values in Tables \ref{ourmoments} and \ref{parevolution} we used only the best light curves in {\it V} filter with standard deviation less than $0.007^m$. Standard dark and flat
field corrections were made. In order to obtain the maximum possible precision of times of minima we used only full light curves between their maxima, our values obtained for 2002, 2004, and 2017 observations are shown in Table \ref{ourmoments}.

\section{Light equation}
\label{sect:light}

A computer code was used to find orbital elements and system parameters. A description of the method used in the code to solve light curves can be found in Section ``Algorithms and models'' in the paper by \citet{kozyreva2001}, a similar model was described by \citet{khaliullin1984}. The code seeks for the photometric parameters and orbital elements using a simple model of two spherical stars (with a linear limb darkening law) that move around a common center of masses in elliptical orbits. The parameters are: the radii of the primary and secondary components $r_{1, 2}$, the limb darkening coefficients for the components $u_{1, 2}$, the luminosities of components in fractions of the system's total luminosity $L_{1, 2}$, the inclination of the orbit of the binary with respect to the plane of the sky $i$, the orbital eccentricity $e$, the longitude of periastron of the primary star's orbit $\omega$, the epoch of the primary minimum corresponding to the epoch of the observations analysed, the system's third light parameter $L_3$. AS Cam is a pair of definitely detached stars, therefore the model of two spherical stars with limb darkening is quite reasonable. The solution was accepted as adequate only when ``observed minus calculated'' value was the lowest and it had no systematic deviations within the minima. The values of parameters were found in a free search excluding the limb darkening coefficients. The values of limb darkening coefficients ($0.46\pm 0.03$ and $0.31\pm 0.05$ for the primary and secondary respectively) were taken from the paper by \citet{khaliullin1983}, they were computed using the high quality light curves of AS Cam obtained in 1981. The values of other parameters are slightly different from year to year being within the error bars of their values in Table III in \citet{khaliullin1983}, except $i$ and $e$ (see Table \ref{parevolution}).

We took into account all times of minima from the Table \ref{ourmoments} and from the B.R.N.O. database\footnote{http://var2.astro.cz/ocgate/index.php?lang=en}. The AS Cam apsidal motion is slow, its period is about 2400 years. During 55 years of observations this rate yields only $\approx 8^{\circ}$ of the total cycle ($360^{\circ}$), therefore it is possible to use linear approximations instead of sinusoidal changes of times of minima to investigate the possible light equation in the system. We obtained following ephemerides for the primary and secondary minima with the same orbital period of the central binary:

\begin{eqnarray}
C_1(\textrm{Min}\ \textrm{I}) & = & \textrm{HJD}2440125.60300 + 3.43096730\times E, \\
C_1(\textrm{Min}\ \textrm{II}) & = & \textrm{HJD}2440123.67395 + 3.43096730\times E.
\end{eqnarray}

\noindent where $E$ is the number of orbital cycles since the initial epoch, HJD is the Heliocentric Julian Date of the initial epoch.

The presence of the third body in an eclipsing binary system can observationally appear as the periodical variations of its times of minima in comparison with the system's linear ephemerides. Such variations arise from the motion of the center of masses of the binary star around the center of masses of the triple system. The amplitude of variations for primary minima is given by a light equation:

\begin{equation}
(O-C)_1 = \frac{a_3\sin{i_3}}{c}(1-e_3\cos{E_3})\sin(v_3+\omega_3),
\label{lighteq}
\end{equation}

\noindent where $v_3$ is the true anomaly of the third companion's orbit, $E_3$ is its eccentric anomaly, $a_3$ is the semi-major axis of the third companion's orbit, $i_3$ is the inclination of this orbit with respect to the plane of the sky, $e_3$ is its eccentricity, and $\omega_3$ is the pericentric longtitude. Figure \ref{figure} shows the reference frame for the third body's orbit. The value of $E_3$ is connected with other elements as follows:

$$
\frac{2\pi}{P_3}(t-T_3)=(E_3-e_3\sin E_3),
$$
\noindent where  $T_3$ is the time of the periastron passage by the third body, $t$ is the time, $P_3$ is the orbital period of the third body.

$(O-C)_1$ experiences periodical variations in the time scale $\approx 2$ years (see Figures \ref{figure1} and \ref{figure2}), these variations have the same phase for the primary and secondary minima. It is essential that both sets of observations (1968-1973 in Figure \ref{figure1} and 1980-2017 in Figure \ref{figure2}) can be described by the same light equation curve. The estimations of the parameters in Equation (\ref{lighteq}) that we obtained (we applied the least-squares method for the $(O-C)_1$) are listed in Table \ref{lighteqparams}, and the ephemerides for the third companion are following:

\begin{equation}
\label{eph3}
\textrm{Min}\ \textrm{III} = \textrm{HJD} 2444265.036 + 805.9\times E_1.
\end{equation}

\noindent where $E_1$ is the number of orbital cycles of the third body since the initial epoch.

The mass function

\begin{equation}
f(m) = \frac{(M_3\sin i_3)^3}{(M_1+M_2+M_3)^2} = \frac{(a_3\sin i_3)^3}{P_3^2}.
\label{mf}
\end{equation}

\noindent gives (after the subtraction of the parameters of the light equation) the lower limit of the third body's mass $M_3 \sin i_3 \approx 1.1  M_{\odot}$.

The light equation method allows to estimate only the lower limit of the third body's mass, and the real value depends also on the angle between the plane of the sky and the orbital plane of third companion. If the angle $i_3\leq 30^{\circ}$, the mass of the third body should be $\geq 2.2 M_{\odot}$ (this value is comparable with the secondary star of the central tight binary). Our photometric solutions for the light curve can be plausible only if the luminosity of the third companion is no more than 3.5\% of the total luminosity of the system. If this body is a non-degenerate main-sequence companion, its mass is less than $1.5 M_{\odot}$ and $i_3\geq 43^{\circ}$. Also the spectral lines of the third companion were not found, therefore the mass of a hypothetical main-sequence companion has an upper limit. In general, the suggested body can be a compact remnant (a white dwarf or even a neutron star) or can be a very close binary star like YY Gem.

\section{AS Cam apsidal motion rate re-estimation}

The difference between periods of primary and secondary minima is the indication of the apsidal motion rate (see Figure \ref{figure3}). Using the least square method we computed following ephemerides:

\begin{eqnarray}
C_2(\textrm{Min}\ \textrm{I}) & = & \textrm{HJD}2444939.24519 + 3.43096365\times E, \\
C_2(\textrm{Min}\ \textrm{II}) & = & \textrm{HJD}2444937.32569 + 3.43097095\times E.
\end{eqnarray}

These equations correspond to the rate of the apsidal motion $\dot{\omega}_\mathrm{obs}= 15.5^{\circ}\pm 1.5^{\circ}$ per century. The theoretical apsidal motion rate for AS Cam was found to be $\dot{\omega}_\mathrm{th} = \dot{\omega}_\mathrm{cl} + \dot{\omega}_\mathrm{rel} = 44^{\circ}$ per century \citep{maloney1991}.

The parameters obtained in this study coincide with the parameters of the light equation and of the apsidal motion rate calculated by \citet{kozyreva1999} within the error bars, despite the fact that the results in 1999 were found from much less observational data.

\section{AS Cam slow apsidal motion possible causes}

The influence of the third companion on the apsidal motion rate in general was studied by \citet{khaliullin1991}, \citet{khodykin1997}, \citet{khodykin2004}. The low apsidal motion rate in DI Her and in AS Cam probably could be explained by the perturbations cause by the third body. This idea was improved by \citet{borkovits2007}, they conducted theoretical investigations of configurations of the central binary's orbit and the third companion's orbit using AS Cam as an example. \citet{borkovits2007} considered four variants of the mutual disposition of both orbits: the angle between them equals $0.8^{\circ}$, $20^{\circ}$, $60^{\circ}$, and $90^{\circ}$. The results were shown in figures 6-9 by \citet{borkovits2007} for the fast (in a time scales of about $\leq 100$ years) and for the slow (in a time scales of about 3000 years) evolution of orbital parameters.

It is possible to compare the values of parameters obtained in past (since 1981 for our data) and modern (2017) observations for the short time scale, a large amount of light curves within primary and secondary minima was accumulated. We took light curves with the precision better than 1\% within light curve minima in our data, and the light curve published by \citet{hilditch1972}. It allowed to compute orbital elements, e. g. \citep{kozyreva2001}. We used our light curves of AS Cam obtained in 1981, 1992,1993,1994,1995,1996,2002, 2017, and light curves obtained by \citet{hilditch1972} in 1967 and 1968. The dataset covered 50 years of observations and allowed to trace the change of the eccentricity and the inclination of the central binary's orbit (see Table \ref{parevolution}) and to compare the observed changes with the theory of \citet{borkovits2007}. The eccentricity change is $\approx 0.03 \pm 0.01$ during 50 years of observations, the inclination change during the same period of time is $\approx 1 \pm 0.5^{\circ}$. Our value of the eccentricity change is higher than its theoretical value for mutual inclination of the orbits $i'=0$ and $i'=20^{\circ}$, it is almost the same as the theoretical value calculated for $i'=60^{\circ}$, see \citet{borkovits2007}, Figure 2. Their paper contains only four values of mutual inclination of the orbits, therefore it is possible that $i'$ is less than $60^{\circ}$. For the case of $i'=60^{\circ}$ \citet{borkovits2007} showed that the dependence of the difference between primary and secondary minima is more complicated than a simple difference between two sinusoidal functions shifted relative to each other by $180^{\circ}$. The dependence between times of minima is complicated and has non-uniform character, so the apsidal motion rate can change in different periods of time, see Figure 3 by \citet{borkovits2007}.  Such periods with different apsidal motion rates can have duration up to several centures, and now we can see AS Cam in a form of slow apsidal motion.

There is another explanation of the slow apsidal motion rate, it is the spin-orbit misalignment \citep{albrecht2009,albrecht2011}. \citet{pavlovski2011} obtained spectral lines of AS Cam and found that the lines are narrow in comparison with the expected line width in supposition of synchronous rotation of the stars. They concluded that the axial and orbital rotation are not aligned. \citet{shakura1985} found that the apsidal line can undegro retrograde motion if the
rotation axises of the stars are not aligned to the orbital axis, with the largest effect when the axises are perpendicular. If the stars were rotating at three times the synchronous rate around axises that lie almost in the orbital plane, and if one takes into account the projected rotational velocities, the stars' rotational axes should be tilted by $82^{\circ}$ and $87^{\circ}$  with respect to the orbital axis \citep{pavlovski2011}. If the orbit of that body is highly inclined with respect to the orbit of the inner two stars, then the Lidov-Kozai mechanism \citep{lidov1962,kozai1962} may be invoked to explain the spin-orbit misalignment.

Probably, both mechanisms (the influence of the third body and the spin-orbit misalignment) of the deceleration of the apsidal motion rate are plausible for the AS Cam binary.

\section{Conclusions}
\label{sect:conclusion}

During the 50 years time interval the orbit of the close binary star AS Cam changed its eccentricity by $\approx 0.03$. This fact indicates that the system possesses the third body in a highly inclined orbit with respect to the central binary's orbit. Such configuration can induce the precession of rotational axises of the stars in the binary, and trigger of their spin-orbit misalignment. Both mechanisms can slow down the apsidal motion rate in comparison the system without a third companion, and if the rotational axises were aligned with the orbital axis.

The unresolved question is how the third body with the mass comparable to the masses of the stars in central binary does not show itself neither in spectra of the stars, nor in photometric elements obtained as a solution for the light curves. This body can be a degenerate star or a very close binary system like YY Gem consisting of low luminosity main-sequence stars.

\begin{acknowledgements}
This study was supported by by the funding program BR05236322 of the Ministry of Education and Science RK and by a grant for Leading Scientific Schools NSh-9670.2016.2 (Russia). We are grateful to A. I. Zakharov for the software and to M. A. Krugov for his efforts to establish the collaboration with Tien Shan observatory. We thank the anonymous referee for useful comments.
\end{acknowledgements}