Understanding Neural Networks Through Deep Visualization

Recent years have produced great advances in training large, deep neural networks (DNNs), including notable successes in training convolutional neural networks (convnets) to recognize natural images. However, our understanding of how these models work, especially what computations they perform at intermediate layers, has lagged behind. Progress in the field will be further accelerated by the development of better tools for visualizing and interpreting neural nets. We introduce two such tools here. The first is a tool that visualizes the activations produced on each layer of a trained convnet as it processes an image or video (e.g. a live webcam stream). We have found that looking at live activations that change in response to user input helps build valuable intuitions about how convnets work. The second tool enables visualizing features at each layer of a DNN via regularized optimization in image space. Because previous versions of this idea produced less recognizable images, here we introduce several new regularization methods that combine to produce qualitatively clearer, more interpretable visualizations. Both tools are open source and work on a pre-trained convnet with minimal setup.

Introduction
The last several years have produced tremendous progress in training powerful, deep neural network models that are approaching and even surpassing human abilities on a variety of challenging machine learning tasks (Taigman et al., 2014;Schroff et al., 2015;Hannun et al., 2014). A flagship example is training deep, convolutional neural networks (CNNs) with supervised learning to classify natural images . That area has benefitted from the combined effects of faster computing (e.g. GPUs), better training techniques (e.g. dropout ), better activation units (e.g. rectified linear units (Glorot et al., 2011)), and larger labeled datasets (Deng et al., 2009;Lin et al., 2014).
While there has thus been considerable improvements in our knowledge of how to create high-performing architectures and learning algorithms, our understanding of how these large neural models operate has lagged behind. Neural networks have long been known as "black boxes" because it is difficult to understand exactly how any particular, trained neural network functions due to the large number of interacting, non-linear parts. Large modern neural networks are even harder to study because of their size; for example, understanding the widely-used AlexNet DNN involves making sense of the values taken by the 60 million trained network parameters. Understanding what is learned is interesting in its own right, but it is also one key way of further improving models: the intuitions provided by understanding the current generation of models should suggest ways to make them better. For example, the deconvolutional technique for visualizing the features learned by the hidden units of DNNs suggested an architectural change of smaller convolutional filters that led to state of the art performance on the ImageNet benchmark in 2013 (Zeiler & Fergus, 2013).
We also note that tools that enable understanding will especially benefit the vast numbers of newcomers to deep learning, who would like to take advantage of off-the-shelf software packages -like Theano (Bergstra et al., 2010), Pylearn2 , Caffe (Jia et al., 2014), and Torch (Collobert et al., 2011) -in new domains, but who may not have any intuition for why their models work (or do not). Experts can also benefit as they iterate ideas for new models or when they are searching for good hyperparameters. We thus believe that both experts and newcomers will benefit from tools that provide intuitions about the inner workings of DNNs. This paper provides two such tools, both of which are open source so that scientists and practitioners can integrate them with their own DNNs to better understand them.
The first tool is software that interactively plots the activations produced on each layer of a trained DNN for userprovided images or video. Static images afford a slow, detailed investigation of a particular input, whereas video input highlights the DNNs changing responses to dynamic input. At present, the videos are processed live from a user's computer camera, which is especially helpful because users can move different items around the field of view, occlude and combine them, and perform other manipulations to actively learn how different features in the network respond.
The second tool we introduce enables better visualization of the learned features computed by individual neurons at every layer of a DNN. Seeing what features have been learned is important both to understand how current DNNs work and to fuel intuitions for how to improve them.
Attempting to understand what computations are performed at each layer in DNNs is an increasingly popular direction of research. One approach is to study each layer as a group and investigate the type of computation performed by the set of neurons on a layer as a whole Mahendran & Vedaldi, 2014). This approach is informative because the neurons in a layer interact with each other to pass information to higher layers, and thus each neuron's contribution to the entire function performed by the DNN depends on that neuron's context in the layer.
Another approach is to try to interpret the function computed by each individual neuron. Past studies in this vein roughly divide into two different camps: dataset-centric and network-centric. The former requires both a trained DNN and running data through that network; the latter requires only the trained network itself. One dataset-centric approach is to display images from the training or test set that cause high or low activations for individual units. Another is the deconvolution method of Zeiler & Fer-gus (2013), which highlights the portions of a particular image that are responsible for the firing of each neural unit.
Network-centric approaches investigate a network directly without any data from a dataset. For example, Erhan et al. (2009) synthesized images that cause high activations for particular units. Starting with some initial input x = x 0 , the activation a i (x) caused at some unit i by this input is computed, and then steps are taken in input space along the gradient ∂a i (x)/∂x to synthesize inputs that cause higher and higher activations of unit i, eventually terminating at some x * which is deemed to be a preferred input stimulus for the unit in question. In the case where the input space is an image, x * can be displayed directly for interpretation. Others have followed suit, using the gradient to find images that cause higher activations (Simonyan et al., 2013;Nguyen et al., 2014) or lower activations (Szegedy et al., 2013) for output units.
These gradient-based approaches are attractive in their simplicity, but the optimization process tends to produce images that do not greatly resemble natural images. Instead, they are composed of a collection of "hacks" that happen to cause high (or low) activations: extreme pixel values, structured high frequency patterns, and copies of common motifs without global structure (Simonyan et al., 2013;Nguyen et al., 2014;Szegedy et al., 2013;Goodfellow et al., 2014). The fact that activations may be effected by such hacks is better understood thanks to several recent studies. Specifically, it has been shown that such hacks may be applied to correctly classified images to cause them to be misclassified even via imperceptibly small changes (Szegedy et al., 2013), that such hacks can be found even without the gradient information to produce unrecognizable "fooling examples" (Nguyen et al., 2014), and that the abundance of non-natural looking images that cause extreme activations can be explained by the locally linear behavior of neural nets (Goodfellow et al., 2014).
With such strong evidence that optimizing images to cause high activations produces unrecognizable images, is there any hope of using such methods to obtain useful visualizations? It turns out there is, if one is able to appropriately regularize the optimization. Simonyan et al. (2013) showed that slightly discernible images for the final layers of a convnet could be produced with L 2 -regularization. Mahendran and Vedaldi (2014) also showed the importance of incorporating natural-image priors in the optimization process when producing images that mimic an entire-layer's firing pattern produced by a specific input image. We build on these works and contribute three additional forms of regularization that, when combined, produce more recognizable, optimization-based samples than previous methods. Because the optimization is stochastic, by starting at different random initial images, we can produce a set of opti- Figure 1. The bottom shows a screenshot from the interactive visualization software. The webcam input is shown, along with the whole layer of conv5 activations. The selected channel pane shows an enlarged version of the 13x13 conv5151 channel activations. Below it, the deconv starting at the selected channel is shown. On the right, three selections of nine images are shown: synthetic images produced using the regularized gradient ascent methods described in Section 3, the top 9 image patches from the training set (the images from the training set that caused the highest activations for the selected channel), and the deconv of the those top 9 images. All areas highlighted with a green star relate to the particular selected channel, here conv5151; when the selection changes, these panels update. The top depicts enlarged numerical optimization results for this and other channels. conv52 is a channel that responds most strongly to dog faces (as evidenced by the top nine images, which are not shown due to space constraints), but it also responds to flowers on the blanket on the bottom and half way up the right side of the image (as seen in the inset red highlight). This response to flowers can be partially seen in the optimized images but would be missed in an analysis focusing only on the top nine images and their deconv versions, which contain no flowers. conv5151 detects different types of faces. The top nine images are all of human faces, but here we see it responds also to the cat's face (and in Figure 2 a lion's face). Finally, conv5111 activates strongly for the cat's face, the optimized images show catlike fur and ears, and the top nine images (not shown here) are also all of cats. For this image, the softmax output layer top two predictions are "Egyptian Cat" and "Computer Keyboard." All figures in this paper are best viewed digitally, in color, significantly zoomed in. mized images whose variance provides information about the invariances learned by the unit.
To summarize, this paper makes the following two contributions: 1. We describe and release a software tool that provides a live, interactive visualization of every neuron in a trained convnet as it responds to a user-provided image or video. The tool displays forward activation values, preferred stimuli via gradient ascent, top images for each unit from the training set, deconv highlighting (Zeiler & Fergus, 2013) of top images, and backward diffs computed via backprop or deconv starting from arbitrary units. The combined effect of these complementary visualizations promotes a greater understanding of what a neuron computes than any single method on its own. We also describe a few insights we have gained from using this tool. (Section 2).
2. We extend past efforts to visualize preferred activation patterns in input space by adding several new types of regularization, which produce what we believe are the most interpretable images for large convnets so far (Section 3).
Both of our tools are released as open source and are available at http://yosinski.com/deepvis. While the tools could be adapted to integrate with any DNN software framework, they work out of the box with the popular Caffe DNN software package (Jia et al., 2014). Users may run visualizations with their own Caffe DNN or our pretrained DNN, which comes with pre-computed images optimized to activate each neuron in this trained network. Our pre-trained network is nearly identical to the "AlexNet" architecture , but with local reponse normalization layers after pooling layers following (Jia et al., 2014). It was trained with the Caffe framework on the ImageNet 2012 dataset (Deng et al., 2009).

Visualizing Live Convnet Activations
Our first visualization method is straightforward: plotting the activation values for the neurons in each layer of a convnet in response to an image or video. In fully connected neural networks, the order of the units is irrelevant, so plots of these vectors are not spatially informative. However, in convolutional networks, filters are applied in a way that respects the underlying geometry of the input; in the case of 2D images, filters are applied in a 2D convolution over the two spatial dimensions of the image. This convolution produces activations on subsequent layers that are, for each channel, also arranged spatially. Figure 1 shows examples of this type of plot for the conv5 layer. The conv5 layer has size 256×13×13, which we depict as 256 separate 13×13 grayscale images. Each of the 256 small images contains activations in the same spatial x-y spatial layout as the input data, and the 256 images are simply and arbitrarily tiled into a 16×16 grid in row-major order. Figure 2 shows a zoomed in view of one particular channel, conv5 151 , that responds to human and animal faces. All layers can be viewed in the software tool, including pooling and normalization layers. Visualizing these layers provides intuitions about their effects and functions.
Although this visualization is simple to implement, we find it informative because all data flowing through the network can be visualized. There is nothing mysterious happening behind the scenes. Because this convnet contains only a single path from input to output, every layer is a bottleneck through which all information must pass en-route to a classification decision. The layer sizes are all small enough that any one layer can easily fit on a computer screen. 1 So far, we have gleaned several surprising intuitions from using the tool: • One of the most interesting conclusions so far has been that representations on some layers seem to be surprisingly local. Instead of finding distributed representations on all layers, we see, for example, detectors for text, flowers, fruit, and faces on conv4 and conv5. These conclusions can be drawn either from the live visualization or the optimized images (or, best, by using both in concert) and suggest several directions for future research (discussed in Section 4).
• When using direct file input to classify photos from Flickr or Google Images, classifications are often correct and highly confident (softmax probability for correct class near 1). On the other hand, when using input from a webcam, predictions often cannot be correct because no items from the training set are shown in the image. The training set's 1000 classes, though numerous, do not cover most common household objects. Thus, when shown a typical webcam view of a person with no ImageNet classes present, the output has no single high probability, as is expected. Surprisingly, however, this probability vector is noisy and varies significantly in response to tiny changes in the input, often changing merely in response to the noise from the webcam. We might have instead expected unchanging and low confidence predictions for a given scene when no object the network has been trained to classify is present. Plotting the fully connected layers (fc6 and fc7) also reveals a similar sensitivity to small input changes. • Although the last three layers are sensitive to small input changes, much of the lower layer computation is more robust. For example, when visualizing the conv5 layer, one can find many invariant detectors for faces, shoulders, text, etc. by moving oneself or objects in front of the camera. Even though the 1000 classes contain no explicitly labeled faces or text, the network learns to identify these concepts simply because they represent useful partial information for making a later classification decision. One face detector, denoted conv5 151 (channel number 151 on conv5), is shown in Figure 2 activating for human and lion faces and in Figure 1 activating for a cat face. Zhou et al. (2014) recently observed a similar effect where convnets trained only to recognize different scene types -playgrounds, restaurant patios, living rooms, etc. -learn object detectors (e.g. for chairs, books, and sofas) on intermediate layers.
The reader is encouraged to try this visualization tool out for him or herself. The code, together with pre-trained models and images synthesized by gradient ascent, can be downloaded at http://yosinski.com/deepvis.

Visualizing via Regularized Optimization
The second contribution of this work is introducing several regularization methods to bias images found via optimization toward more visually interpretable examples. While each of these regularization methods helps on its own, in combination they are even more effective. We found useful combinations via a random hyperparameter search, as discussed below.
Formally, consider an image x ∈ R C×H×W , where C = 3 color channels and the height (H) and width (W ) are both 227 pixels. When this image is presented to a neural network, it causes an activation a i (x) for some unit i, where for simplicity i is an index that runs over all units on all layers. We also define a parameterized regularization function R θ (x) that penalizes images in various ways.
Our network was trained on ImageNet by first subtracting the per-pixel mean of examples in ImageNet before inputting training examples to the network. Thus, the direct input to the network, x, can be thought of as a zero-centered input. We may pose the optimization problem as finding an image x * where In practice, we use a slightly different formulation. Because we search for x * by starting at some x 0 and taking gradient steps, we instead define the regularization via an operator r θ (·) that maps x to a slightly more regularized version of itself. This latter definition is strictly more expressive, allowing regularization operators r θ that are not the gradient of any R θ . This method is easy to implement within a gradient descent framework by simply alternating between taking a step toward the gradient of a i (x) and taking a step in the direction given by r θ . With a gradient descent step size of η, a single step in this process applies the update: We investigated the following four regularizations. All are designed to overcome different pathologies commonly encountered by gradient descent without regularization.
L 2 decay: A common regularization, L 2 decay penalizes large values and is implemented as r θ (x) = (1−θ decay )·x. L 2 decay tends to prevent a small number of extreme pixel values from dominating the example image. Such extreme single-pixel values neither occur naturally with great frequency nor are useful for visualization. L 2 decay was also used by Simonyan et al. (2013).
Gaussian blur: Producing images via gradient ascent tends to produce examples with high frequency information (see Supplementary Section S1 for a possible reason). While these images cause high activations, they are neither realistic nor interpretable (Nguyen et al., 2014). A useful regularization is thus to penalize high frequency information. We implement this as a Gaussian blur step r θ (x) = GaussianBlur(x, θ b width ). Convolving with a blur kernel is more computationally expensive than the other regularization methods, so we added another hyperparameter θ b every to allow, for example, blurring every several optimization steps instead of every step. Blurring an image multiple times with a small width Gaussian kernel is equivalent to blurring once with a larger width kernel, and the effect will be similar even if the image changes slightly during the optimization process. This technique thus lowers computational costs without limiting the expressiveness of the regularization. Mahendran & Vedaldi (2014) used a penalty with a similar effect to blurring, called total variation, in their work reconstructing images from layer codes.
Clipping pixels with small norm: The first two regularizations suppress high amplitude and high frequency information, so after applying both, we are left with an x * that contains somewhat small, somewhat smooth values. However, x * will still tend to contain non-zero pixel values everywhere. Even if some pixels in x * show the primary object or type of input causing the unit under consideration to activate, the gradient with respect to all other pixels in x * will still generally be non-zero, so these pixels will also shift to show some pattern as well, contributing in whatever small way they can to ultimately raise the chosen unit's activation. We wish to bias the search away from such behavior and instead show only the main object, letting other regions be exactly zero if they are not needed. We implement this bias using an r θ (x) that computes the norm of each pixel (over red, green, and blue channels) and then sets any pixels with small norm to zero. The threshold for the norm, θ n pct , is specified as a percentile of all pixel norms in x.
Clipping pixels with small contribution: Instead of clipping pixels with small norms, we can try something slightly smarter and clip pixels with small contributions to the activation. One way of computing a pixel's contribution to an activation is to measure how much the activation increases or decreases when the pixel is set to zero; that is, to compute the contribution as |a i (x) − a i (x −j )|, where x −j is x but with the j th pixel set to zero. This approach is straightforward but prohibitively slow, requiring a forward pass for every pixel. Instead, we approximate this process by linearizing a i (x) around x, in which case the contribution of each dimension of x can be estimated as the elementwise product of x and the gradient. We then sum over all three channels and take the absolute value, computing | c x • ∇ x a i (x)|. We use the absolute value to find pixels with small contribution in either direction, positive or negative. While we could choose to keep the pixel transitions where setting the pixel to zero would result in a large activation increase, these shifts are already handled by gradient ascent, and here we prefer to clip only the pixels that are deemed not to matter, not to take large gradient steps outside the region where the linear approximation is most valid. We define this r θ (x) as the operation that sets pixels with contribution under the θ c pct percentile to zero.
If the above regularization methods are applied individually, they are somewhat effective at producing more interpretable images; Figure 3 shows the effects of each individual hyperparameter. However, preliminary experiments uncovered that their combined effect produces better visualizations. To pick a reasonable set of hyperparameters for all methods at once, we ran a random hyperparameter search of 300 possible combinations and settled on four that complement each other well. The four selected combinations are listed in Table 1 and optimized images using each are shown for the "Gorilla" class output unit in Figure 4. Of the four, some show high frequency information, others low frequency; some contain dense pixel data, and others contain only sparse outlines of important regions. We found the version in the lower-left quadrant to be the best single set of hyperparameters, but often greater intuition can be gleaned by considering all four at once. Figure 5 shows the optimization results computed for a selection of units on all layers. A single image for every filter of all five convolutional layers is shown in Supplementary Figure S1. Nine images for each filter of all layers, including each of the 1000 ImageNet output classes, can be viewed at http://yosinski.com/deepvis.  (Table 1). For all other classes, we have selected four interpretable visualizations produced by our regularized optimization method. We chose the four combinations of regularization hyperparameters by performing a random hyperparameter search and selecting combinations that complement each other. For example, the lower left quadrant tends to show lower frequency patterns, the upper right shows high frequency patterns, and the upper left shows a sparse set of important regions. Often greater intuition can be gleaned by considering all four at once. In nearly every case, we have found that one can guess what class a neuron represents by viewing sets of these optimized, preferred images. Best viewed electronically, zoomed in. In each layer, we show visualizations from 4 random gradient descent runs for each channel. While these images are hand picked to showcase the diversity and interpretability of the visualizations, one image for each filter of all five convolutional layers is shown in Figure S1 in supplementary information. One can recognize important features of objects at different scales, such as edges, corners, wheels, eyes, shoulders, faces, handles, bottles, etc. The visualizations show the increase in complexity and variation on higher layers, comprised of simpler components from lower layers. The variation of patterns increases with increasing layer number, indicating that increasingly invariant representations are learned. In particular, the jump from Layer 5 (the last convolution layer) to Layer 6 (the first fully-connected layer) brings about a large increase in variation. Best viewed electronically, zoomed in. Each of the four rows shows a linear sweep in hyperparameter space from no regularization (left) to strong regularization (right). When applied too strongly, some regularizations cause the optimization to fail (e.g. L2 decay, top row) or the images to be less interpretable (small norm and small contribution clipping, bottom two rows). For this reason, a random hyperparameter search was useful for finding joint hyperparameter settings that worked well together (see Figure 4). Best viewed electronically, zoomed in.

Discussion and Conclusion
We have introduced two visual tools for aiding in the interpretation of trained neural nets. Intuition gained from these tools may prompt ideas for improved methods and future research. Here we discuss several such ideas.
The interactive tool reveals that representations on later convolutional layers tend to be somewhat local, where channels correspond to specific, natural parts (e.g. wheels, faces) instead of being dimensions in a completely distributed code. That said, not all features correspond to natural parts, raising the possibility of a different decomposition of the world than humans might expect. These visualizations suggest that further study into the exact nature of learned representations -whether they are local to a single channel or distributed across several -is likely to be interesting (see Zhou et al. (2014) for work in this direction). The locality of the representation also suggests that during transfer learning, when new models are trained atop the conv4 or conv5 representations, a bias toward sparse connectivity could be helpful because it may be necessary to combine only a few features from these layers to create important features at higher layers.
The second tool -new regularizations that enable improved, interpretable, optimized visualizations of learned features -will help researchers and practitioners understand, debug, and improve their models. The visualizations also reveal a new twist in an ongoing story. Previous stud-  Figure 4. The third row hyperparameters produced most of the visualizations for the other classes in Figure 4, and all of those in Figure 5. ies have shown that discriminative networks can easily be fooled or hacked by the addition of certain structured noise in image space (Szegedy et al., 2013;Nguyen et al., 2014). An oft-cited reason for this property is that discriminative training leads networks to ignore non-discriminative information in their input, e.g. learning to detect jaguars by matching the unique spots on their fur while ignoring the fact that they have four legs. For this reason it has been seen as a hopeless endeavor to create a generative model in which one randomly samples an x from a broad distribution on the space of all possible images and then iteratively transforms x into a recognizable image by moving it to a region that satisfies both a prior p(x) and posterior p(y|x) for some class label y. Past attempts have largely supported this view by producing unrealistic images using this method (Nguyen et al., 2014;Simonyan et al., 2013).
However, the results presented here suggest an alternate possibility: the previously used priors may simply have been too weak (see Section S1 for one hypothesis of why a strong p(x) model is needed). With the careful design or learning of a p(x) model that biases toward realism, one may be able to harness the large number of parameters present in a discriminately learned p(y|x) model to generate realistic images by enforcing probability under both models simultaneously. Even with the simple, hand-coded p(x) models we use in this paper as regularizers, complex dependencies between distant pixels already arise (cf. the beetles with structure spanning over 100 pixels in Figure 4). This implies that the discriminative parameters also contain significant "generative" structure from the training dataset; that is, the parameters encode not only the jaguar's spots, but to some extent also its four legs. With better, learned probabilistic models over the input and activations of higher layers, much more structure may be apparent. Work by Dai et al. (2015) shows some interesting results in this direction. While the images generated in this paper are far from being photo-realistic, they do suggest that transferring discriminatively trained parameters to generative models -opposite the direction of the usual unsupervised pretraining approach -may be a fruitful area for further investigation.

Hod Lipson
HOD.LIPSON@CORNELL.EDU S1. Why are gradient optimized images dominated by high frequencies?
In the main text we mentioned that images produced by gradient ascent to maximize the activations of neurons in convolutional networks tend to be dominated by high frequency information (cf. the left column of Figure 3). One hypothesis for why this occurs centers around the differing statistics of the activations of channels in a convnet. The conv1 layer consists of blobs of color and oriented Gabor edge filters of varying frequencies. The average activation values (after the rectifier) of the edge filters vary across filters, with low frequency filters generally having much higher average activation values than high frequency filters.
In one experiment we observed that the average activation values of the five lowest frequency edge filters was 90 versus an average for the five highest frequency filters of 5.4, a difference of a factor of 17 (manuscript in preparation) 2,3 .
The activation values for blobs of color generally fall in the middle of the range. This phenomenon likely arises for reasons related to the 1/f power spectrum of natural images in which low spatial frequencies tend to contain higher energy than high spatial frequencies (Torralba & Oliva, 2003). Now consider the connections from the conv1 filters to a single unit on conv2. In order to merge information from both low frequency and high frequency conv1 filters, the connection weights from high frequency conv1 units may generally have to be larger than connections from low frequency conv1 units in order to allow both signals to affect the conv2 unit's activation similarly. If this is the case, then due to the larger multipliers, the activation of this particular conv2 unit is affected more by small changes in the activations of high frequency filters than low frequency filters.
Seen in the other direction: when gradient information is passed from higher layers to lower layers during backprop, the partial derivative arriving at this conv2 unit (a scalar) will be passed backward and multiplied by larger values when destined for high frequency conv1 filters than low frequency filters. Thus, following the gradient in pixel space may tend to produce an overabundance of high frequency changes instead of low frequency changes.
The above discussion focuses on the differing statistics of edge filters in conv1, but note that activation statistics on subsequent layers also vary across each layer. 4 This may produce a similar (though more subtle to observe) effect in which rare higher layer features are also overrepresented compared to more common higher layer features.
Of course, this hypothesis is only one tentative explanation for why high frequency information dominates the gradient. It relies on the assumption that the average activation of a unit is a representative statistic of the whole distribution of activations for that unit. In our observation this has been the case, with most units having similar, albeit scaled, distributions. However, more study is needed before a definitive conclusion can be reached.

S2. Conv Layer Montages
One example optimized image using the hyperparameter settings from the third row of Table 1 for every filter of all five convolutional layers is shown in Figure S1. 4 We have observed that statistics vary on higher layers, but in a different manner: most channels on these layers have similar average activations, with most of the variance across channels being dominated by a small number of channels with unusually small or unusually large averages (Li, Yosinski, Clune, Song, Hopcroft, Lipson. 2015. How similar are features learned by different deep neural networks? In preparation.) conv5 conv3 conv4 conv2 conv1 Figure S1. One optimized, preferred image for every channel of all five convolutional layers. These images were produced with the hyperparameter combinations from the third row of Table 1. Best viewed electronically, zoomed in.


Fast Blocked Clause Decomposition with High Quality

Any CNF formula can be decomposed two blocked subsets such that both can be solved by BCE (Blocked Clause Elimination). To make the decomposition more useful, one hopes to have the decomposition as unbalanced as possible. It is often time consuming to achieve this goal, since it has been proven to be NP-hard. So far there have been several decomposition and post-processing algorithms such as PureDecompose, QuickDecompose, EagerMover etc. We found that these existing algorithms are often either inefficient or low-quality decomposition. To build an efficient and high-quality algorithm, we improve the existing BCE, and present two new variants of PureDecompose, a new heuristic decomposition called LessInterfereDecompose, and a new post-processing algorithm called RsetGuidedDecompose. Combining these new techniques results in a new algorithm called MixDecompose. In our experiments, there is no application formula where the quality of PureDecompose+EagerMover is better than MixDecompose. In terms of speed, MixDecompose is also very fast. Our average runtime is a little longer, but the worst-case runtime is shorter. In theory, our two variants of PureDecompose requires linear time in the number of clauses. By limiting the size of the touch list used by BCE, we can guarantee always that MixDecompose runs in linear time.

Introduction
Recently, one found that blocked clause decomposition can not only efficiently find backbone variables [1] and implied binary equivalences through SAT sweeping, but also improve the performance of the state-of-the-art SAT solvers such as Lingeling [2] on hard application benchmarks [3,4]. Due to its importance, many researchers have been attracted to pay attention to this subject.
A set of clauses is said to a blocked set if it can be removed it completely by Blocked Clause Elimination (BCE) [5,6]. Any CNF formula can be decomposed into two blocked subsets. To make a blocked clause decomposition more useful, one wants always to have two blocked subsets as unbalanced as possible. The problem is that it is not easy to find the most unbalanced subsets. In theory, one has proven that finding a maximal blocked subset of a CNF formula with the largest cardinality (MaxBS for short) is NP-hard [3]. In other words, it is impossible to find the best decomposition in polynomial time unless P = N P . So far a few decomposition algorithms were proposed. However, no algorithm achieves optimization in all terms. PureDecompose [3] is the fastest, but its quality is poor. To improve the quality, Heule et al [3] presented QuickDecompose. However, QuickDecompose is time-consuming. Soon after, to improve the speed, Balyo et al [4] developed a post-processing algorithm called EagerMover. Through a exhaustive series of experiments, we noted that although the decomposition quality of PureDecompose+EagerMover (PureEager for short) and QuickDecompose can outperform PureDecompose, their quality is not high yet.
This paper aims to present an algorithm for finding efficiently a high-quality decomposition. To achieve this goal, we present two new variants of PureDecompose, a new decomposition algorithm based on clause correlation degree, and a new post-processing algorithm. In addition, we improve the existing BCE to speed up the decomposition. The algorithm resulting from integrating these new techniques is called MixDecompose, which can improve significantly the quality of decomposition. On application instances, the decomposition quality of MixDecompose is better than that of PureEager. there is no application formula where the quality of PureEager is better than MixDecompose. In terms of speed, MixDecompose is still fast. On average, it took 8.97 seconds on our machine, which is a little slower than PureEager took 7.41 seconds. However, in the worst case, MixDecompose was faster than PureEager. The latter exceeded 300 seconds in some cases, whereas the former took at most 110 seconds.

Preliminaries
In this section, we present basic concepts that will be used in subsequent algorithms for blocked clause decomposition.
CNF. It is short for conjunctive normal form. A formula in CNF is formulated as a conjunction of clauses, where each clause is a disjunction of literals, each literal being either a Boolean variable or its negation. The negation of a variable x is denoted byx or ¬x. In general, a clause C is written as is a clause. The symbols var(F ) and lit(F ) denote the sets of variables and literals occurring in a formula F , respectively.
Resolution. Given two clauses C 1 = (l ∨ a 1 ∨ · · · ∨ a m ) and C 2 = (l ∨ b 1 ∨ · · · b n ), the clause C = (l ∨ a 1 ∨ · · · ∨ a m ∨ b 1 ∨ · · · b n ) is called the resolvent of C 1 and C 2 on the literal l, which is denoted by Blocked Clauses.Given a CNF formula F , a clause C, a literal l ∈ C is said to block C w.r.t. F if (i) C is a tautology w.r.t. l, or (ii) for each clause C ′ ∈ F withl ∈ C ′ , C ′ ⊗ l C is a tautology. A clause is a tautology if it contains both x andx for some variable x. When l blocks C w.r.t. F , the literal l and the clause C are called a blocking literal and a blocked clause, respectively. BCE. It is short for blocked clause elimination, which removes blocked clauses from CNF formulas. By BCE(F ) we mean the CNF formula resulting from repeating the following operation until fixpoint: If there is a blocked clause C ∈ F w.r.t. F , let F := F − {C}. It is said that BCE can solve a formula F if and only if BCE(F ) = ∅. The seminal work in BCE is due to Kullmmann [5].

Blocked Clause Decomposition
In theory, any CNF formula can be decomposed into two blocked subsets. However, not all the decompositions are effective. In general, The larger one of the blocked sets is, the better the decomposition quality is, since the larger it is, the more it resembles the original formula. Therefore, the size difference of the two sets is considered as a measure of the decomposition quality. Nevertheless, computing the largest blocked set from a CNF formula is NP-hard. Hence, we here aim to find a fast decomposition with higher quality, rather than the highest quality.
The simplest decomposition way is called pure decomposition, which is shown in Fig. 1. Let the symbols L and R denote the left (large) subset and the right (remainder ) subset, respectively. For each variable x, this algorithm adds always the larger of F x and Fx to L and the smaller to R, where F x (Fx) is the set of clauses of F where x occurs positively (negatively). At the termination of this algorithm, we have F = L ∪ R with |L| ≥ |R|. In Fig. 1, max{F x , Fx} means the set with the larger cardinality between F x and Fx.  The advantage of this algorithm is that it can be easily implemented to run in linear time in the size of F , using a standard structure of occurrence lists. Therefore it is very fast. The drawback is that its decomposition quality is not high on many formulas. For this reason, we will improve it by combining some variants of it and the other algorithms.
By a few empirical observations, we found that the performance of PureDecompose rely significantly on the order in which variables are eliminated. Here, we present two variants of PureDecompose with different variable elimination ordering heuristics. The first variant is called min pure decomposition, which is shown in Fig. 2. One fifth of variable eliminations are to be eliminated in the same order as PureDecompose. The remaining variables are to be eliminated in order from the lowest occurrence of literals to the highest. If there are multiple literals with the lowest occurrence, the literal with the minimum total number of clauses containing it is eliminated first. The total clause size of a literal x can be formulated as C∈Fx |C|. Compared to PureDecompose, this algorithm adds only the search of variables to be eliminated. This search can be done in O(n log n) time, using an order heap, where n is the number of variables. In the actual implementation, the number of variables in min x∈lit(F ) |F x | is limited to 30000. This can guarantee that MinPureDecompose is still very fast even if n is very large. In terms of decomposition quality, this algorithm is superior to the other algorithms on some application instances such as ctl 4291 567 5 unsat pre. Now we consider the second variant of PureDecompose. The order of its variable eliminations is opposite to that of the first variant. We call this variant max pure decomposition, which is shown in Fig. 3. It always eliminate first a literal with the highest occurrence. When multiple literals have the same highest occurrence, we select a variable with the lowest different of its two literal occurrences. This can be done by computing min |Fx|=m ||F x | − |Fx||, where m is defined as max x∈lit(F ) |F x |. Actually, the first variant can introduce also this tie-break method. Unlike the first variant, MaxPureDecompose needn't compute the total clause size of each literal. So it should run faster than the first variant. In order to ensure that is still very fast even if the number of variables is very large, the range for finding a literal with the highest occurrence is limited to 5000 variables. The decomposition quality of this algorithm is superior to that of the other algorithms on some application instances such as complete-500. The above two algorithms both are based on the order of variable elimination. However, in some cases, the algorithm obtained by optimizing the order of variable elimination is not necessarily optimal. Therefore, one must find more efficient algorithms to replace them. Blow we present a new algorithm called less interfere decomposition, which is based on the order of clause elimination. Its pseudo-code is shown in Fig. 5. This algorithm may sketch the basic outline as follows: move blocked clauses in F to L by BCE, compute the candidate set S, move each clause C ∈ S ∩ F to R. These steps are repeated until F is empty. The computation of the candidate set S is based on the notion of interfering degree. The interfering degree of a clause C can be defined as ogy clause, and one otherwise. The probability that C ′ ⊗ l C is not a tautology clause is very high. To save the computing cost, we may approximate the interfering degree as Fig. 5 uses this approximation version to compute the interfering degree, and call this measure score, i.e., score[C] = C ′ ∈F |var(C ′ ) ∩ var(C)|. To get clauses with the maximum score, all the clauses in F are traversed. If only one clause with the maximum score is moved to R from F each time F is traversed, it is time consuming. So we decide to move p clauses to R one time, where p = |F | θ , where θ is a constant. In our experiment, when |F | < 8 × 10 4 , θ is set to 2300. Otherwise, θ is set to 400 for random instances, and 100 for non-random instances. When |F | θ < 18, p is set to 18. As shown in Fig. 5, the clauses with the first p highest scores are stored in S as the candidate clauses to be moved to R. In order to save time further, we compute the interfering degree produced by only literals with the lowest occurrence, not all literals.
LessInterfereDecompose calls to BCE shown in Fig. 4 to move blocked clauses to L. BCE here is different from that presented in [6]. BCE in [6] is based on a literal-based priority queue, while our BCE is based on a clause-based linear linked list. Another important difference from the usual BCE is that we do not try to test whether each literal l in C is a blocking literal when |F | ≥ 200000. We test only literal l with |Fl| < 2. That is, we replace the statement " for l ∈ C do " in the usual BCE with the statement " for l ∈ C with (|Fl| < 2 or |F | < 200000) do". Due to this simplification, our BCE is much faster than the usual BCE. Surprisingly, the decomposition quality keep unchanged in most cases. Even if it is changed, its change is still very small. In addition, our touch function is different from that in [6]. Here is our definition about it.
Fx otherwise When |F | ≥ 600000, we consider only the clauses touched by the negation of literals with the number of occurrences < 2. This can speed up the decomposition of large instances. For example, using the above touch, the runtime of decomposing on q query 3 L90 can be reduced from 600 seconds to 9 seconds. The decomposition quality keep unchanged still. The runtime of LessInterfereDecompose consists of three parts: BCE, computing scores and determining S's. The runtime of computing scores is O( |F | 2 p ) = O(θ|F |). If scores are given, determining a S can be done in a linear time in |F |, since there exist linear time algorithms for finding the p-th highest score [7,8]. The total runtime of computing scores plus determining S's does not exceed O(θ|F |). If the number of clauses touched by each clause does not exceed a constant δ, where δ is certainly smaller than the maximal number of literal occurrences times the maximal size of clauses, i.e., max x∈lit(F ) |F x | × max C∈F |C|. The total time required by all BCE 's is at most O(δ|F |). Thus, the total runtime of LessInterfereDecompose is at most O((δ + θ)|F |). In practice, δ is generally very small. Should δ is very large, we can remove a part of touched clauses to reduce the time required by BCE to test whether a clause in the touch list is blocked, or limit the size of the touch list a small constant, say 2000. Using such a policy can guarantee that the time complexity of LessInterfereDecompose is linear in |F |. Compared with EagerMover in [4], the runtime of BCE in LessIn-terfereDecompose is smaller than that in EagerMover. EagerMover calls at least four times BCE on a subset with the size of 0.75|F |. All calls to BCE on each clause C in F can be viewed as a call to BCE on the whole F . Thus, the total runtime of BCE in LessInterfereDecompose corresponds to double the runtime of BCE on a F . As long as the runtime of computing scores and determining S's is smaller than the runtime of BCE on a F , LessInterfereDecompose should be faster than EagerMover. In fact, that is true. On some instances, the former are indeed faster than the latter.
In general, the above algorithms do not achieve maximal blocked set decomposition. However, they can be improved further by post-processing. The postprocessing often used is MoveBlockedClasue algorithm shown in Fig. 7, which is to move blocked (with respect to the current L) clauses from R to L. We noted that even if this post-processing algorithm is applied, the decomposition quality can be improved still. For this reason, we present a new post-processing algorithm, called Right set guided decomposition, which is shown in Fig. 6. It is a simplified version of LessInterfereDecompose. Replacing S with R results in this algorithm. This algorithm requires that the right blocked set R must be given in advance. Hence, it is used generally as post-processing. It is faster than LessIn-terfereDecompose, since it need not computing R. Its time complexity depends mainly on that of BCE. For some benchmarks, this algorithm can improve significantly their decomposition quality. For example, to decompose SAT dat.k75-24 1 rule 3 using MinPureDecompose, the fractions of the large subset (i.e., |L| |F | ) with RsetGuidedDecompose and without it are 83.9% and 69.9%, respectively. If replacing MinPureDecompose with LessInterfereDecompose, their fractions are 87.8% and 87.3%, respectively. RsetGuidedDecompose raises still the quality by 0.5%. However, the speed difference among the three algorithms is big. On this instance, LessInterfereDecompose, RsetGuidedDecompose and MinPureDecompose spent 25, 4 and 1 seconds, respectively. The slowest LessInterfereDecompose is not suitable for huge instances with ten millions of clauses. In general, we do not know which algorithm is the best. Because each of the previous three algorithms is very fast, running the three algorithms one after another does not lose much time. Thus, we can construct an algorithm with high speed and high performance by combining the algorithms above-mentioned. The detailed implementation is shown in Fig. 7. We call this algorithm MixDecom-pose. Its basic idea is to take the maximum from three left sets outputted by three algorithms as the initial L first. If the formula to be decomposition is not large, say the number of clauses is less than 5 × 10 6 , we invoke LessInterfereDecompose to get a larger L. Finally, we enlarge the size of L by calling two post-processing algorithms: RsetGuidedDecompose and MoveBlockedClasue. Like the usual postprocessing, the task of MoveBlockedClasue is to move blocked clauses from R to L. Notice, if F is large, say |F | > 10 7 , the last post-processing can be canceled to save the running time. According to whether both subsets can be solved by BCE, a blocked clause decomposition can be classified as symmetric or asymmetric. If yes, it is symmetric. If only one of the subsets can be solved by BCE, it is asymmetric. Clearly, our two variants of PureDecompose are symmetric. Like EagerMover [4], if blockable clauses (whose definition is given below) are allowed to move to L, RsetGuided-Decompose are asymmetric, since it cannot guarantee that the left set can be solved by BCE. However, by our observation, on almost all application instances, they are symmetric.

Empirical evaluation
We evaluated the performance of decomposition algorithms under the following experimental platform: Intel Core 2 Quad Q6600 CPU with speed of 2.40GHz and 2GB memory.
Because QuickDecompose [3] requires more time than EagerMover [4] for many instances, we decided to select more competitive PureDecompose+EagerMover [9] (PureEager for short) as our comparison object. PureEager and MixDecompose both are written in C.
The large set L obtained by PureEager contains blockable clauses in addition to blocked clauses. A clause C is said to be blockable w.r.t. a blocked set L if each literal l ∈ C is not a blocking literal of any clause in L. The reason why blockable clauses are add to the blocked set is that they do not destroy the blocked property. That is, blocked sets containing blockable clauses are still satisfiable. To keep identical with the performance evaluation of PureEager, the large set L of our MixDecompose contains also blockable clauses. In addition, before calling MixDecompose, we use the same unit decomposition policy given in [4] as EagerMover to preprocess the input formula. We evaluated PureEager and MixDecompose on the 297 instances from the application track of the SAT competition 2014, except for three huge instances. Due to limited space, Table 1 lists only a part of results on the application track and a random instance from the SAT competition 2014 shown in the last row. Column |F | denotes the number of the clauses in formula F in ten thousands of clauses, where F is simplified by removing satisfied clauses, but contains unit clauses. Column |L| |F | indicates the fraction of the large set. Column Time shows the runtime in seconds. As seen in Table 1, in terms of decomposition quality, MixDecompose outperforms completely PureEager. In terms of speed, the former is sometimes faster than the latter, and vice versa. MixDecompose was able to finish the decomposition on all SAT 2014 application benchmarks excluding three huge instances within 110 seconds. However, PureEager was not able to finish on some benchmarks such as 9vliw m 9stage iq3 C1 b7 within 300 seconds.  Table 2 shows the empirical result on 297 benchmarks from SAT Competition 2014 application track. Three huge instances zfcp-2.8-u2-nh, esawn uw3.debugged and post-cbmc-zfcp-2.8-u2 were removed. The second column shows the average fraction of the large set. Column '# of best' indicates the number of the best results obtained by an algorithm. Column '# of eq' is the number of results equivalent to another algorithm. On 226 out of 297 benchmarks, the size of the large set obtained by MixDecompose is larger than that obtained by PureEager. On 71 remaining benchmarks, the quality of the two algorithms is identical. there is no application formula where the quality of PureEager is better than MixDecompose. In addition, we conducted also experiments on random benchmarks. We observed that on all random instances, the quality of MixDecompose is strictly better than that of PureEager. As seen from the last row of Table 2, MixDecompose can solve huge random instances with millions of clauses in a reasonable time. For 3-SAT random instances, it can increase the fraction of the large set by 5%.
The fifth column in Table 2 shows the average runtime taken by each algorithm in seconds. Here, computing the average runtime counts only solved instances, excluding timed-out instances. The last column in Table 2, lists the number of times the time-out was hit. The timeout for each algorithm was set to 300 seconds. MixDecompose did not time out on the tested benchmarks, while PureEager did on 7 benchmarks. MixDecompose took at most 110 seconds. Although on average, PureEager run faster than MixDecompose, in this experiment, the worst-case runtime of the former was significantly lager than the latter.
According to results of four algorithms, MixDecompose outputs the best result among them. For 297 application instances, the number of the best results obtained by LessInterfereDecompose, PureDecompose, MaxPureDecompose and MinPureDecompose are 215, 18, 30 and 34, respectively. This empirical result shows that in most cases, LessInterfereDecompose is superior to the other algorithms.

Conclusions and Future work
In this paper, we developed a new blocked clause decomposition algorithm by combining several decomposition strategies. The new algorithm not only achieves high quality decomposition, but also is fast. Even for large instances, it can ensure that the decomposition is done within 110 seconds on our machine. Because our machine is slower than the platform of SAT competition 2014, If running on the latter, the speed will be more fast.
So far we know only that we can get a higher quality decomposition than the existing algorithms such as PureEager. However, this does not mean that MixDecompose is the best. How to develop a better and more efficient than MixDecompose will be a future research topic.


Neural Networks Based Equalizer For Signal Restoration In Digital Communication Channels

One of the main obstacles to reliable communications is the inter symbol interference (ISI). An equalizer is required at the receiver to mitigate the effects of non-ideal channel characteristics and restore the transmitted signal. This paper presents the equalization of digital communication channels using artificial neural network structures. The performances of a nonlinear equalizer using multilayer perceptron (MLP) trained by

INTRODUCTION
In modern digital communications systems, there is an ever increasing need for high speed data transmission rate over the variety of limited bandwidth channels which, by their nature, distort the digital signal. The effect of this distortion is to cause the transmitted symbol to persist beyond the time interval allocated for its transmission; thus causing subsequent symbols to interfere. This distortion is called inter symbol interference (ISI) which can also arise in multipath propagation. The transmitted signal is also subject of other impairment such as nonlinear distortion and noise. At the receiver, an equalizer is used to mitigate these effects and restore the transmitted symbols [1]. Due to the variety and time varying nature of some channels, adaptive equalization techniques need to be implemented. Adaptive equalization has been an attractive area of research; many of the techniques are firmly based on linear adaptive filters, which limit the performance of the system. To improve system performances, non-linear structures have been proposed for the channel equalization task. Nonlinear equalizers are superior to linear ones in situation where the channel distortion is too severe. In particular, a linear equalizer does not perform well on non-minimum phase channels, which are difficult to equalize [2], [3]. Neural networks provide powerful nonlinear processing architectures, which have been proved to outperform traditional linear techniques in many common signal-processing applications. In particular, neural networks have been proposed for digital equalization of communication channels [4], [5], [6], and [11]. Neural networks have several properties, which make them very attractive for channel equalization; these properties are adaptive processing, universal approximation, learning and generalisation capacities etc...
The key issue in neural network is to find an appropriate structure that gives the best result. Different artificial neural network structures such as multilayer perceptron [2], [3], [5], [7] radial basis function [4] and recurrent neural network [8], [9] have been investigated in channel equalization. Among all these structures, the most widely used is the MLP structure due to its computational simplicity, finite parameterization and its stability [7].  LTD, 2015 In this paper, nonlinear channel equalization using multilayer perceptron is considered. We present a study of the MLP based equalizer and we exhibit its capacities to equalize minimum phase channel and non-minimum phase one as well. We compare and evaluate the performance of the MLP equalizer against the conventional LTE in term of the steady state MSE reached and the equalized signals for the two channels considered.

A.
Problem statement Many channels can be characterized by discrete-time equivalent model described by a finite impulse response digital filter (FIR) and an additive noise [1]. Such appropriate channel model is depicted in Figure 1.
The digital data sequence x k is passed through a linear dispersive channel of finite impulse response with a transfer function given by: Where M is length of the impulse response, the observed sequence r k is formed by adding Gaussian random noise n k to the output of the FIR filter.
When dealing with the equalization, the main goal is to recover the transmitted unobserved data symbols x k , with the highest fidelity, using the information given by the noisy observation from a set of received symbols r k . Thus, the task of the equalizer is to construct a causal approximation to the inverse of the channel. The output of the equalizer is than applied to a memory less decision device that selects the element of the transmitted alphabet. For binary data, it is the sign function.

B.
Digital channel models Many channels can be characterized by discrete-time equivalent model described by a Finite Impulse Response (FIR) digital filter [1], with a transfer function given by (1). Such appropriate channel model is depicted in Figure 2.

C. Linear Transversal Equalizer (LTE)
The structure of a linear transversal equalizer (LTE) is depicted in Figure 3.The digital data sequence observed at the input of the equaliser can be given by a convolution as following: Where x k is the binary random sequence of the transmitted symbols drawn from {-1, 1}, with equal probability and are assumed to be independent of one another. The additive noise samples are chosen independently from a Gaussian distribution with mean zero and variance σ n 2 .
Equivalently the equation above can be rewritten in a matrix vector notation as: International Letters of Chemistry, Physics and Astronomy Vol. 55 189 The coefficients of the LTE can be trained to minimise the mean square error (MSE) cost function between the sequence equalised y k and the wanted sequence x k-d as followed: The iterative method by means of LMS algorithm is used to determine the optimum solution. The coefficients of the equaliser are updating as follows: μ is the step size; its choice must be a compromise between convergence speed and fine equalizer settings.

MULTILAYER PERCEPTRON BASED EQUALIZER
Before Neural network based equalizers have been privileged in the nonlinear equalization area [2], [6] because of their nonlinear processing capacities. Especially the MLP based equalizer is the most popular equalizer because of its attractive properties [2], [3], [7]. Its basic element, the artificial neuron illustrated in Figure 4 (a), is composed of a x N linear combiner and an activation function. The neuron receives inputs from other neurons processors. The linear combiner output is the weighted sum of the inputs plus a bias term. The activation function, which can be linear or nonlinear, gives then the neuron output: Where x j is the j th input value of the neuron, W j the corresponding synaptic weight, and b the bias term. f (.) is the activation function ( Figure 4 (b)), the most commonly used is of the sigmoid type [10], defined by: Where δ define the nonlinearity degree of the sigmoid function it is generally set to one.
A multilayer perceptron (Figure 4 (c)) is composed of neurons connected to each other. The input information is processed from the input layer to the output layer. The network inputs are the inputs of the first layer. The outputs of the neurons in one layer form the inputs to the next layer. The network outputs are the outputs of the output layer. In general, all neurons in a layer are fully interconnected to neurons in adjacent layers, but there is no connection within a layer, and normally no connections bridging layers. A multilayer perceptron (MLP) consists of several hidden layers of neurons, which are capable of performing complex, non-linear mapping between the input and the out put layer. The hidden layers have the ability of performing complex, nonlinear mappings between the input and output layer. The hidden layers provide the capability to use the nonlinear sigmoid ability (Figure 4 (b)); the MLP can have several of these hidden layers [10]. The output X ik of neuron (i, k) of the MLP is given by: Where i, is the layer index, X ik is the output of neuron k of layer i, W ijk is the weight that links the output X i-1j to neuron k of layer i. N i is the number of neurons in layer i. It has been demonstrated that a two layer perceptron with a sigmoid activation function and a scalar output can approximate arbitrarily well continuous functions, provided that an arbitrarily large number of neurons is available [10]. This property is called the universal approximation property.

A. Learning Algorithm
The MLP can be trained to achieve a particular task using the back propagation (BP) algorithm. A set of input/output pairs (x(n), d(n)) trains the network to implement the desired mapping. The error back propagation learning consists of two passes through the different layers of the network: a forward and a backward passes. In the forward pass, an activity pattern is applied to the input neurons of the network and its effects propagate through the network, layer by layer, until an output is produced. Weights between neurons of successive layers are initially assigned in random. In the backward pass, the error observed between the network and the desired responses is computed and used to amend the weights. Each neuron starting from the hidden layer is muddled with a non-linear activation function. The back propagation adjusts the MLP weights to minimize any differentiable cost function, e.g. the squared error between the network output and the desired output, as below: En  dnX L n 2 . Where X L (n) is the network output at time n and d(n) is the desired output. The BP algorithm performs a gradient descent on the cost function in order to reach a minimum. The weights and thresholds levels are updated, according to the following relations:

ILCPA Volume 55
Where μ is the learning gain, α is the momentum term. δ i k is the error term calculated starting from the output layer. The output error term is given by: Where f' denotes the derivative of the activation function f ' y df y / dy . The term error is then recursively back propagated to lower layers. For the hidden neuron (i, k), the error term δ ik is given by the following: A momentum term α is used to allow for a rapid learning and to filter out high frequency variations of the weight vector. Thus, the convergence rate is much faster and the fast weight changes are smoothed out.

SIMULATION RESULTS
Simulation results were all realized using MATLAB. The digital message applied to the channel is made of uniformly distributed bipolar random numbers {-1, 1}. The channel output is corrupted by zero mean white Gaussian noise with variance σ n 2 . Two digital channel models are used for the simulation test; a minimum phase channel and a non-minimum phase channel denoted as CH 1 and CH 2 respectively. Both the conventional the LTE and the MLPE use five input samples. The MLP has 3 layers, 9 neurons in the first hidden layer, 3 neurons in the second hidden layer, 1 neuron in the output layer. The parameters of the back propagation and the LMS algorithms have been chosen as in [3] and [7].

A.
Channel caracteristics Simulations are carried out for two digital channel models as following:   The minimum phase channel, having the following transfer function n in Z domain: CH 1 z  0.87 0.44z 1  0.23z 2  And the non minimum phase channel with the transfer function in Z domain given by:

ILCPA Volume 55
These channel models are widely used to evaluate the performance of equalizers in communication systems [2], [3], [5], [7], [9], [12]. Channel characteristics are depicted in Figure 6 for the minimum phase channel CH 1 and the non-minimum phase channel CH 2 respectively. The magnitude responses of the channels (part (a)) exhibit deep notches, about -5dB for CH 1 and about -15dB for CH 2 , which are difficult to equalise specially for CH 2 .
The phase response (part(c)) is non-linear for CH 1 and linear for CH 2 . The channel impulse responses coefficients are about three as depicted in part (b). Channels roots in Z plan (part (d)) are about two zeros. Both are inside the unit circle for CH1, one zero is outside the unit circle for CH2. The out side zero may causes instability of the equalizers.

B.
Restored signal performance In this part of simulation, the performance measure is obtained by expecting the restored signal behaviour quality. Figure 7 and 8 show the restored signal features by the equalizers for the channels CH 1 and CH 2 respectively. In both figures, Part (a) shows the transmitted signal, (b) illustrates the received signal witch is completely distorted by channel ISI and additive noise; thus the detection without any equalization process will produce severe errors. Where parts (c) and (d) represent the restored signal by the LTE and the MLPE equalizers respectively, the two equalizers perform well, whoever the improvement in the restoration quality for the MLP equalizer is pretty obvious. The impairments of the channel CH2 on the International Letters of Chemistry, Physics and Astronomy Vol. 55 transmitted signal became worse. The amplitude levels of some symbols are greatly attenuated. The LTE restored signal presents some residual distortion; however the MLPE exhibit greater improvement. The MLPE restores the transmitted signal with successful manner thus allowing detection without any errors. The MLPE shows the best performance for the two channels used, this results from the nonlinear processing of the MLP structure.

C. Learning curves performance
In this part, the performance measure is obtained by learning curves that describe the evolution of the mean square error (MSE) at the output of the equalizers. Figure 9 and 10 depict the convergence behaviour of the two equalizers for the channels CH 1 and CH 2 respectively. The MLPE converge more slowly than the LTE, but it shows a clear improvement in the steady state value of averaged square error produced (-28 dB), which is lower than the noise level (-20dB). The LTE produce a steady state MSE level of about -16 dB, which is greater than the noise level. A similar improvement of the MLPE is also obtained for the non-minimum channel CH 2 . The steady state MSE reached by the LTE became worse under the great distortion effect engendered by the large ISI of the channel CH 2 . Te LTE MSE level increase to achieve -11 dB. Whereas the MSE level of the MLPE is around -25dB which steel under the noise level.

ILCPA Volume 55
The increase in the steady state MSE is about -5dB for the LTE and only -3dB for the MLPE influenced by the sever ISI introduced by the non-minimum channel CH 2 . It should be clear from these figures that the steady state MSE of MLP based equalizer is below the noise level for both the channels considered.

D.
BER performance To more investigate the consistency in performance of the MLPE, bit error rate (BER) performances are studied for the equalizer configurations and are shown in figures 12 and 13.

International Letters of Chemistry, Physics and Astronomy Vol. 55 197
The MLPE shows a very good performance, it performs better than the LTE. Figure 11, which corresponds to CH shows that a gain of about 2dB is made by the MLPE at a BER of 10 −3 when compared to the LTE. On the other hand in the for the Figure 12 which corresponds to CH 2, the MLPE still consistently behaving better than the LTE a gain of about 2,5dB is made by the MLPE over the LTE at BER of 10 −3 . It can be seen that the MLPE achieves almost the best performance even in all the scenarios considered.


An Interrupt-Driven Work-Sharing For-Loop Scheduler

In this paper we present a parallel for-loop scheduler which is based on work-stealing principles but runs under a completely cooperative scheme. POSIX signals are used by idle threads to interrupt left-behind workers, which in turn decide what portion of their workload can be given to the requester. We call this scheme Interrupt-Driven Work-Sharing (IDWS). This article describes how IDWS works, how it can be integrated into any POSIX-compliant OpenMP implementation and how a user can manually replace OpenMP parallel for-loops with IDWS in existing POSIX-compliant C++ applications. Additionally, we measure its performance using both a synthetic benchmark with varying distributions of workload across the iteration space and a real-life application on Sandy Bridge and Xeon Phi systems. Regardless the workload distribution and the underlying hardware, IDWS is always the best or among the best-performing strategies, providing a good all-around solution to the scheduling-choice dilemma.

INTRODUCTION
Most parallelism in shared-memory parallel programming comes from loops of independent iterations, i.e. iterations which can be safely executed in parallel. However, distributing the iteration space over the available computational resources of a system is not always a simple thing. Finegrained control of distribution is often associated with high overhead whereas static partitioning of the iteration space can lead to significant load imbalance. In both cases, the impact on performance is serious.
Research on an advanced for-loop scheduler was motivated by our work on PRAgMaTIc [12], a hybrid OpenMP/MPI mesh adaptivity framework. Profiling data revealed that many of PRAgMaTIc's parallel loops are highly diverse, involving irregular computations which introduce high levels of iteration-to-iteration load imbalance. Existing scheduling strategies provided by OpenMP fail to achieve good balance with low scheduling overhead, whereas adaptive mesh algorithms which constantly modify mesh topology make it impossible to balance workload a priori.
We wanted the new scheduler to be portable and easily plug-able into the widely-adopted OpenMP API, so that it can target an as wide as possible range of systems, like Fujitsu's FX-10, a SPARC64-based supercomputer [6]. Those portability requirements prohibit the use of platform-or vendor-specific threading mechanisms and parallel libraries, like Intel R Cilk TM Plus [11,5,3]. On the contrary, they call for a POSIX-compliant implementation, based on the fact that most operating systems used in scientific computing are POSIX-compliant and most compilers (e.g. Linux versions of gcc, icc, xlc, etc.) implement OpenMP threads as POSIX threads (we have found it out by experimenting with those compilers). Of course, since every OS has threading and signalling mechanisms, the new scheduler can be implemented into any compiler on any OS.
The main contributions of this article are the following: • Present an new interrupt-driven work-sharing scheduler (IDWS) which can easily be used with existing POSIX-compliant OpenMP applications • Demonstrate how OpenMP loops can be converted to IDWS loops • Describe how a compiler vendor can incorporate the new scheduler into their product • Show using a variety of benchmarks that IDWS is a good all around solution to the scheduling-choice dilemma, always being among the best-performing strategies in all benchmarks The rest of this paper is organized as follows: Section 2 provides an overview of loop scheduling options currently available, listing their advantages and weaknesses. In Section 3 we describe the new scheduler, the way it works and explain why it offers a better tradeoff between load balance and scheduling overhead compared with other alternatives. Section 4 goes into details about the current implementation in C++ and explains how it can be used to replace OpenMP for-loops in existing codes. We demonstrate the scheduler's performance in Section 5 using both a synthetic benchmark and real for-loops from PRAgMaTIc. Finally, we conclude the paper in Section 6.

BACKGROUND
OpenMP offers three different scheduling strategies for parallel loops: static, dynamic and guided [7]. There is also a more advanced scheduling technique, known as "workstealing", which is implemented by libraries such as Intel Cilk Plus, though it is not part of the OpenMP specification, nor is it supported (to the best of our knowledge) by any OpenMP implementation. In this section we will present these four options and compare them in terms of load balance, scheduling overhead and overall efficiency.

OpenMP static
Under the static scheduling scheme, the iteration space is divided into equally large chunks which are then assigned to threads. This can be seen in the first example in Figure 1.
Partitioning of iteration space is done statically at the beginning of the for-loop, so there is zero scheduling overhead.
On the other hand, this scheme can lead to significant load imbalance, especially in a highly diverse loop.

OpenMP dynamic
Dynamic scheduling is a first approach to the problem of load imbalance. Instead of a static partitioning of the iteration space, chunks of work are assigned to threads dynamically. Once a thread finishes a chunk, it takes the next available from the iteration space. This is shown in the middle example in Figure 1. Access to chunks is done via atomic updates of the loop counter; a thread acquiring a chunk reads the current value of the loop counter and increments it atomically by the chunk size.
Dynamic scheduling solves imbalance problems as threads proceed to the next iterations of the for-loop in a fine-grained way. As an immediate consequence, good load balance comes at a cost. The loop counter is updated atomically and this constitutes a 2-way source of overhead. The two components of overhead are related to instruction latency and thread competition. The time it takes to execute an atomic instruction can vary anywhere between a standard update in L1 (if the thread performing the update is running on the same physical core as the thread which last updated the shared variable) and an update in RAM (if the last thread to update the shared variable is running on another socket in case of NUMA systems). This may not be a problem in short for-loops, but becomes easily a hotspot in loops with millions of iterations and little work per iteration (i.e. when atomic instruction latency is comparable to the loop body itself). Secondly, as the number of threads increases, so does the competition for the shared variable, leading to either (depending on the architecture) increased locking or increased number of failed update attempts, thus making atomic instruction latency even longer.
It could be argued that this overhead can be mitigated by increasing the chunk size, therefore lowering the number of times a thread will need to access the loop counter. On the other hand, increasing the chunk size can introduce load imbalance once again. Additionally, it is usually impossible to know the optimal chunk size at compile time and/or it can vary greatly between successive executions of an algorithm. Besides, relying on the chunk size for performance optimization puts an extra burden on the programmer.
We have found that using dynamic scheduling over guided in PRAgMaTIc can increase the execution time of specific algorithms by up to three times, as will be shown in Section 5. Following that, dynamic scheduling was rejected as an option for that framework.

OpenMP guided
The guided scheme is an attempt to reduce dynamic scheduling overhead while retaining good load balance. The key difference between the two strategies is how chunks of work are allocated. Whereas in dynamic scheduling the chunk size is constant, the guided scheme starts with large chunks and the size is reduced exponentially as threads proceed to subsequent iterations of the for-loop. This can be seen in the last example of Figure 1. Initial large chunks account for reduced atomic accesses to the loop counter while the more fine-grained control towards the end of the loop tries to maintain good load balance.
For the most part, guided scheduling works well in PRAg-MaTIc, yet there are cases where we have observed significant load imbalance. This can happen if, for instance, most of the work in an irregular loop is accumulated in a few of the initial big chunks. In a case like that, threads processing the "loaded" chunks are busy for long while others go through the remaining "light" iterations relatively quickly and reach the end of the for-loop early, waiting for the busy workers to finish as well.

Work-stealing
Work-stealing ( [4,13]) is a more sophisticated technique aiming at balancing workload among threads while keeping scheduling overhead as low as possible. The generic workstealing algorithm for a set of tasks [4] can be summarized as follows. Each thread keeps a deque (double-ended queue) of tasks to execute. While the deque is full, the thread pops workitems from the front. Once the deque is empty, the thread becomes a thief, i.e. it chooses a victim thread randomly and steals a chunk of workitems from the back of the victim's deque.
For a parallel for-loop with a predefined number of iterations N the deques can simply be replaced with pairs of indices < istart, i end > corresponding to the range in the iteration space [istart, i end ) which has been assigned to each thread, 0 ≤ istart, i end < N . In this case, every thread executes iterations by using istart as the loop counter whereas thieves steal work by decrementing a victim's i end .
Accesses to those pairs of indices can lead to race conditions, so they need to be accessed with atomics. Following that, work-stealing for for-loops comes close to OpenMP's dynamic scheduling with some chunk size > 1, with a major difference being that in work-stealing threads do not compete all together for atomic access to the same shared vari-static dynamic guided Thread 0 Thread 1 Thread 2 Thread 3 Figure 1: Example with four threads of the three scheduling strategies offered by OpenMP: static, dynamic (chunk=2) and guided. Note that under the guided scheme chuck size is reduced exponentially.
able (the common loop counter); instead, congestion is rare and happens only if two thieves try to steal from the same victim.
Performance can still suffer from load imbalance and scheduling overhead when using work-stealing. The main drawback of the classic work-stealing algorithm is that thieves choose victims randomly. There is no heuristic method to indicate which threads are more suitable victims (i.e. have more remaining workload) than others. Stealing comes at a cost and picking victims with too little or no remaining work is inefficient, as it leads to the need for frequent stealing which induces some overhead. Additionally, failed attempts do not help balance the workload. As an example of an extreme case, a single thread becomes the sole remaining worker while the rest waste time trying to steal from each other in vain.
Mitigating the effects of random choice was our main concern when designing the new for-loop scheduler. We devised a low-overhead heuristic method for finding appropriate victims. At the same time, we tried to reduce scheduling overhead by eliminating the need to use atomics when accessing each thread's < istart, i end > pair. The following section describes in detail how the scheduler is implemented.

INTERRUPT-DRIVEN WORK SHARING
Our new scheduler differs from existing work-stealing approaches in two major ways. First of all, as was mentioned in Section 2.4, every worker constantly "advertises" its progress so that thieves can find suitable victims which have been left behind. Secondly, a thief does not actually steal work from the victim in the classic sense; instead, it interrupts the chosen victim by sending a POSIX signal. The signal handler executed by the victim encapsulates the code with which the victim decides what portion of its remaining workload can be given away. As it becomes apparent, the new scheduling algorithm is much closer to work-sharing than work-stealing, therefore we call it Interrupt-Driven Work-Sharing (IDWS). Nonetheless, we will use work-stealing terminology throughout this article.
The abstract description of this scheme can be split into three parts: Algorithm 1 Parallel loop executed by each thread for i = istart; i < i end ; i ← i + 1 do flush i ⊲ from register file to memory, so that ⊲ thieves can see this thread's progress execute i th iteration flush i end ⊲ from memory to register file, as it may ⊲ have been modified by the signal handler end for Algorithm 2 Work-stealing for all threads tn do remainingt n ← i end,tn − istart,t n − it n end for let T ← tn for which remainingt n = max send signal to victim T wait for answer update own istart, i end execute loop chunk Loop execution (Algorithm 1). Every thread executes the iterations of the chunk it has acquired in the same way as it would using OpenMP's static scheduling scheme. Initially, the iteration space is divided statically into chunks of equal size and every thread tn is assigned one chunk. The chunk's boundaries for thread tn, referred to as istart,t n and i end,tn , are globally visible variables accessible by all threads. Compared to static scheduling, the important addition here is some necessary flushing of the loop counter it n and the loop boundary i end,tn . More precisely, the value of it n has to be written back to memory (instead of being cached in some register) at the beginning of every iteration so that potential thieves can monitor tn's progress, calculate how much work is left for tn and decide whether it is worth stealing from it. Similarly, the end boundary i end,tn has to read by tn from memory (instead of caching it in some register) before proceeding to the next iteration because i end,tn might have been modified by the signal handler if a thief interrupted tn while the latter was executing an iteration of the for-loop.
Choosing suitable victims (Algorithm 2). By flushing their loop counters, threads advertise their progress so potential thieves can find where to steal from. When a thread becomes a thief, it calculates the remaining workload for all other threads by reading the associated values it n , istart,t n and i end,tn . This way, we have a heuristic method for finding which thread has the most remaining work, thus being a more suitable victim than others. This heuristic may not be optimal, but is an improvement over random choice. Once the thief has spotted its victim, it sends a signal and waits for an answer. The victim executes the signal handler and replies with the boundaries (a pair of < istart, i end >) of the chunk it wants to give away. Finally, the thief becomes a worker once again and moves on to process the newly acquired chunk.

Signal handler (Algorithm 3). When a victim is inter-
rupted by the signal, control is transferred by the operating system to the associated handler. Inside that code, the thread calculates how much work it can give away (if any), replies to the thief with the boundaries of the donated chunk, re-adjusts the boundaries of its own chunk and finally returns from the signal handler.
It is clear that there are no races and no need for atomic access to any loop variables during the stealing process, as the donor is the one who decides what will be donated. Of course, switching from user to kernel mode to execute the signalling system call and busy-waiting for a reply from the victim involves some overhead; however, as will be shown in the results section, this method seems to be more efficient than classic work stealing.

C++ IMPLEMENTATION AND USAGE
This section describes how the IDWS scheduler is implemented and how it can replace existing OpenMP for-loops.

IDWS namespace
IDWS is a namespace encapsulating all necessary data structures and functions used by the new scheduler. Its declaration can be seen in Code Snippet 1. struct IDWS::thread_state_t. The heart of IDWS is a struct named thread_state_t, which encapsulates all variables involved in parallel loop execution and work-stealing. Each participating thread has its own instance of this struct, which is accessible by all other threads. The struct can be seen in Code Snippet 2. • processed: Is used by a thread to advertise its progress through the loop • current_ctx: IDWS loops are nowait loops, which means that a thread can proceed to the rest of the program without synchronising with other threads. In order to know whether two threads work inside the same loop, so stealing work from one another is valid, a counter current_ctx is used, which is incremented each time a thread finishes a loop. Here we assume that all threads will go through all loops of the program.
• active: Indicates whether the thread is inside the loop; this variable is used by thieves to skip immediately threads which have also become thieves.
• signal_arg: POSIX signals can only have two arguments, what signal is to sent and to whom. The victim needs to know, however, who the thief is, so sig-nal_arg is used by the thief to send its ID to the victim.
• comm_lock: In order to avoid needless busy-waiting by other thieves while one thief has already sent a signal to its victim, we use a lock (in form of a mutex); while this lock is held by a thief, other thieves will choose other victims to steal from.
• ptid: POSIX ID of the thread; it is used by the thief to raise the signal.
• request_mutex and wait_for_answer: POSIX mutex and condition variables which assist the process of sending the signal and waiting for a reply.
Note that we need two separate mutexes and cannot use re-quest_mutex in place of comm_lock. The former is implicitly released by the thief in order to enable the victim to signal the condition variable; in the meantime, before the victim locks request_mutex, another thief might acquire the lock and destroy the process.

Signal handler. A victim decided what portion of its chunk
can be donated by executing the signal handler. The way it is done is described in Code Snippet 3. The victim first checks that the thief works in the same context. Then, it calculates how much work it can give away using start, end and processed, also leaving a safety margin due to an uncertainty regarding the true value of processed. In case of success, the victim updates both the thief's and its own start and end and sets sig_arg=1 to indicate successful donation (otherwise, sig_arg is set to another value). Finally, the victim signals the condition variable to let the thief know that the signal handler is over.

Prologue and epilogue macros
The new scheduler is defined in two parts, using macros IDWS_prologue and IDWS_epilogue. These macros must surround the loop body.
IDWS_prologue macro. Before entering a loop, the iteration space is split into equal chunks which are assigned to threads. After that, each thread begins the execution of its chunk. Compared to a standard for-loop, a IDWS for-loop is defined slightly differently. Apart from checking for the end of the loop and incrementing the counter after every iteration, in IDWS we must also enforce the compiler to flush the counter back to memory and load the updated value of i end from memory (which might have been modified by the signal handler), as indicated by Algorithm 1. Flushing is done selectively for those two variables by casting them to volatile datatypes. Using #pragma omp flush would flush the entire shared program state, which is not efficient. A pseudo-code of how the macro expands is given in Code Snippet 4. Parameters TYPE, NAME and SIZE correspond to the datatype of the loop counter, its name and the size of the iteration space, respectively. In the current implementation of the new scheduler we assume that loops run from 0 to SIZE with increments of 1 and that the loop counter is of an unsigned integral datatype.
IDWS_epilogue macro. After a thread finished its chunk, it becomes a thief. That means it has to enter the stealing process, as described in Algorithm 2. The IDWS_epilogue macro serves this purpose. The way the macro expands can be seen in Code Snippet 4. The thief calculated for all active workers the amount of remaining work. Then, starting from the worker with the highest remaining workload, it tries to acquire the worker's comm_lock. If no suitable worker is found, then the thief exits the IDWS loop and proceeds to the rest of the code. Otherwise, the thief locks the victim's mutex, sends the signal and waits on the victim's condition variable for an answer. The answer comes via sig_arg.
If sig_arg==-1, then the victim has set the thief's start and end variables, so the thief becomes a worker again.
If any other answer has been sent back, then the thief exits the IDWS loop. It is important to note that a memory fence is necessary on the thief's side between setting the victim's signal argument sig_arg and raising the signal, so that we make sure that the victim will see the correct value of sig_arg. Locking the victim's mutex before sending the signal works as an implicit memory fence. The new scheduler can be used directly with virtually any C++ OpenMP application written for any POSIX-compliant operating system (provided that the compiler used implements OpenMP upon pthreads). A prerequisite for converting an OpenMP loop to a IDWS one is that the former is written as shown in Code Snippet 5, i.e. the loop must be inside an omp parallel region. Conversion to IDWS loops is shown in Code Snippet 6. The user needs to include header file "IDWS.h" which can be downloaded from PRAg-MaTIc's page on Launchpad 1 . This header file defines the IDWS namespace and the prologue and epilogue macros. Compared to the initial version, we need to define:

OpenMP to IDWS
• int nthreads=omp_get_max_threads(): shared variable outside the parallel region • int tid=omp_get_thread_num(): thread-private variable inside the parallel region , remove the #pragma omp for directive and the for-loop declaration and, finally, surround the loop-body with the IDWS_prologue and IDWS_epilogue macros.  In order to measure the performance of our new scheduler, we ran a series of tests using both synthetic benchmarks and real kernels from the PRAgMaTIc framework. We used three systems: a dual-socket Intel Xeon E5-2650 (Sandy Bridge, 2.00GHz, 8 physical cores per socket, 2 hyperthreads per core, 32 threads in total), a dual-socket Intel Xeon E5-2643 (Sandy Bridge, 3.30GHz, 4 physical cores per socket, 2 hyperthreads per core, 16 threads in total) and an Intel Xeon Phi 0b/01 board (1.2GHz, 61 physical cores, 4 hyperthreads per core, 244 threads in total). The two Xeon systems run Red Hat Enterprise Linux Server release 6.4 (Santiago). Both versions of the code (intel64 and mic) were compiled with Intel Composer XE 2013 SP1 using the -O3 optimisation flag. The benchmarks were run using Intel's thread-core affinity support with the maximum number of available threads on each platform. Additionally, we ran a second series of benchmarks on Xeon Phi using half the available number of threads and more specifically using all 61 cores with 2 hyperthreads per core. We did so because we have observed that for some codes Xeon Phi performs best when using this threading configuration.

EXPERIMENTAL RESULTS
The synthetic benchmark was designed to be compute-bound with minimal memory traffic and no thread synchronization. Our purpose is to show how the different scheduling strategies compare to each other in terms of achievable load balance and incurred scheduling overhead without being affected by other factors (such as memory bandwidth, data locality etc.  butions of workload and total amount of work. Each population has been given a name: • Regular: All elements of states are set equal to 2. This is a distribution corresponding to a regular loop which does the exact same thing in every iteration.
• Random: states is populated with random values following a uniform distribution. This sub-benchmark corresponds to real-life distributions in problems like graph coloring or the swap and smooth kernels in PRAg-MaTIc.
• Dense End: Most of the workload is accumulated towards the end of the iteration space, where states[i]=3, while the beginning is populated with a uniform mixture of values [0..3). The rest of the iteration space is set to 0, i.e. no work. This is a distribution closely related to the refinement kernel in PRAgMaTIc.
• Dense Start: Mirrored distribution of Dense End. Closely related to the PRAgMaTIc's coarsening kernel. This is an example of workload distribution for which OpenMP guided scheduling is a bad choice.
• Periodic: There is a repeating pattern of states throughout the iteration space. It is particularly bad for static scheduling with interleaved allocations of iterations (i.e. with some chunk size).
Apart from the synthetic benchmark, we also ran PRAg-MaTIc using the various scheduling options in order to see how each strategy performs in a real-life scenario, where compute capacity is not the only performance-limiting factor. It should be noted that PRAgMaTIc is build upon OpenMP, so there are no results for Cilk+ in this case. Table 1, Table 2 and Tables 3 & 4 show the execution time on the three platforms, respectively, using six scheduling strategies for each distribution of the synthetic benchmark and the four PRAgMaTIc kernels. The strategy named "OMP static,1" is static scheduling with chunk size equal to 1. As can be seen, IDWS is either the fastest scheduling option or very close to the fastest for each benchmarkplatform combination. Additionally, it clearly outperforms Cilk Plus, with the performance gap becoming wider as the number of threads increases and Cilk's design to pick victims in a random fashion becomes inefficient. Those results are a good indication that IDWS is future-proof and ready for the thousand-core era.
Regarding PRAgMaTIc, IDWS's major competitor seems to be OpenMP's guided scheduling. Despite not being very suitable for certain kernels (coarsening) theoretically, in practice it performs just as well as IDWS. A notable exception is the 244-thread case on Xeon Phi, where guided scheduling is the worst choice among the available options.
A comparison of relative performance between the three major competitors (IDWS, OpenMP guided and Cilk Plus) is shown in Figure 2 (Intel Xeon E5-2650 system), Figure 3 (Intel Xeon E5-2643 system), Figure 4 (Intel Xeon Phi with 122 threads) and Figure 5 (Intel Xeon Phi with 244 threads).
For each benchmark, we compare the relative execution time between IDWS, OpenMP guided and Cilk Plus (for PRAg-MaTIc kernels there is only IDWS vs OpenMP guided comparison). Reference execution time per benchmark, i.e. the one which corresponds to 1.0 on the y-axis, is execution time of the fastest scheduler.

CONCLUSIONS AND FUTURE WORK
We have presented an Interrupt-Driven Work-Sharing forloop scheduler which is based on work-stealing principles and tries to address major problem of the original work-stealing algorithm: random choice of victims. The first implementation of IDWS works very efficiently, outperforming Intel Cilk Plus, while being from slightly slower to considerably faster than the best (per kernel) OpenMP scheduling strategy. These results indicate that IDWS could become the universal default scheduler for OpenMP for-loops, freeing the programmer from tricky and disruptive management of load balance.
Two main points of focus for further work should be data locality and efficiency of work-sharing. Work on locality issues has been published by several groups ([9, 8, 10, 1]), whereas Adnan and Sato have presented interesting ideas on efficient work-stealing strategies [2], some of which could be applicable to our work-sharing scheduler.