Pablo Mesejo
email: pmesejo@ce.unipr.it

Stefano Cagnoni
email: cagnoni@ce.unipr.it

An experimental study on the automatic segmentation of in situ hybridization-derived images

Keywords: Biomedical Image Segmentation, In Situ Hybridization-derived Images, Histological Images, Hippocampus, Deformable Models, Soft Computing, Comparative Study

Histology is a well established branch of biology focused on the study of the microscopic anatomy of cells and tissues. Despite its importance, it is surprising to check that most literature devoted to histological image processing and analysis is focused on registration and 3D reconstruction of the whole brain. Therefore, there is a lack of research about automatic segmentation of anatomical brain structures in histological images. To bridge this gap, this paper introduces a comparative study of different segmentation techniques applied to this kind of images.

A wide and representative image set has been collected to run experiments on the hippocampus, due to the importance of this anatomical district in learning, memory and spatial navigation. Seven approaches, from deterministic to non-deterministic ones, and from recent trends to classical computer vision techniques, have been compared using different standard metrics (Dice Similariry Coefficient, Jaccard Index, Hausdorff Distance, True Positive Rate and False Positive Rate). Proper statistical tests have been performed to draw accurate conclusions about the results. The best performance on this particular problem was obtained by a combination of Active Shape Models (optimized using Differential Evolution) with a refinement step based on Random Forests. This approach achieved an average Dice Similarity Coefficient of 0.89 with a standard deviation of 0.03.

I. INTRODUCTION

Classically, image segmentation is defined as the partitioning of an image into non overlapping regions that are homogeneous with respect to some visual feature such as intensity or texture [START_REF] Pham | Current Methods in Medical Image Segmentation[END_REF]. In particular, image segmentation plays a crucial role in many biomedical imaging applications (detection of lesions, measurement of organ size/volume, quantitative tissue analysis, computer-integrated surgery,...), by automating and facilitating the delineation of anatomical structures and different regions of interest (ROIs). As one can easily check taking a look at the literature [START_REF] Maulik | Medical image segmentation using genetic algorithms[END_REF]- [START_REF] Mcinerney | Topology adaptive deformable surfaces for medical image volume segmentation[END_REF], most of the segmentation algorithms developed so far have been applied to MRI and CT images, not paying too much attention to other important imaging modalities like ultrasound or microscopy.

Among the different anatomical structures which make up the mammalian brain, the hippocampal formation (HPF) (see Figure 1) is particularly interesting. From an anatomical viewpoint, the HPF, composed by the Hippocampus and the Subiculum (SUB), is located within the medial temporal lobe. In turn, the Hippocampus is composed by the Dentate Gyrus (DG) and Ammon's Horn (CA), which is further composed by three different regions (CA1, CA2, and CA3).

The hippocampus has long been known for its crucial role in learning and memory processes [START_REF] Norman | How hippocampus and cortex contribute to recognition memory: revisiting the complementary learning systems model[END_REF]- [START_REF] Thompson | Genomic Anatomy of the Hippocampus[END_REF]. Moreover, it has recently been demonstrated that the volume of the hippocampus is an early biomarker for Alzheimer's disease [START_REF] Barnes | A meta-analysis of hippocampal atrophy rates in Alzheimer's disease[END_REF]. Therefore, there is a great interest in understanding the cellular and molecular events that take place in this structure, under both normal and abnormal conditions. From this point of view, a precise gene expression map at the cellular and subcellular level within this region can provide crucial information for understanding such biological mechanisms.

A very promising data source to derive this map has recently been provided by the Allen Brain Atlas (ABA), a huge, publicly available database that contains high-resolution images mapping the expression patterns of most genes contained in the genomes of the analyzed organisms. The first release of the ABA was focused on mouse, and contained the expression patterns of about 20,000 genes obtained by In Situ Hybridization (ISH) of mouse brain sections [START_REF]Allen Reference Atlases[END_REF]. In a more recent release, a similar resource, based on spatially mapped microarray data, has also been provided for the human brain [START_REF] Jones | The Allen Brain Atlas: 5 years and beyond[END_REF], [START_REF] Hawrylycz | An anatomically comprehensive atlas of the adult human brain transcriptome[END_REF]. The availability of such a large number of images, containing important morphological and functional information, has made it extremely important to develop completely automatic methods that can accurately and robustly extract information from the hippocampus. In the ABA, for each gene under consideration, several images are provided, corresponding to different sections of the brain, and each image is labelled according to the corresponding image of a reference atlas (comprising 132 coronal and 21 sagittal sections spaced at 100 µm and 200 µm intervals, respectively).

ISH is a technology for gene expression analysis that uses a labeled complementary DNA or RNA strand to localize a specific DNA or RNA sequence in a portion or section of tissue (in situ) [START_REF] Jin | In situ hybridization: Methods and applications[END_REF]. As its name suggests, it is a method that localizes and detects specific mRNA sequences in morphologically preserved tissue sections or cell preparations, by hybridizing the complementary strand of a nucleotide probe to the sequence of interest.

Since the task under study in this work is the automatic segmentation of the hippocampus in histological images obtained by ISH, it is important to describe the main characteristics of such images in general, and of the hippocampus region in particular. These features can be summarized as follows:

• natural variability of brain structure shapes in different subjects; • fuzziness of the hippocampus boundary;

• limited relevance of color for detecting anatomical structures: regions with similar colors may represent different structures and vice versa, depending on the dye used as well as on local image acquisition settings; • contrast variability between structures: different genes are not expressed equally in the same anatomical region, making it difficult to construct a consistent model for each landmark in all images. Moreover, grained patterns with many irregularities hamper the classification of individual pixels as belonging to the anatomical structures under consideration; • orientation issues: the different structures may be rotated or displaced on the slice with respect to a "standard" alignment; • lighting issues: within the same set of images, some are much brighter than others. • variable resolution even within the same image: highresolution regions coexist with low-resolution ones; • presence of artefacts: tears, scraps, bubbles, streaks in tissues, partial cut-off of regions; • large image size (the typical resolution of complete ABA images is about 15,000 × 7,000 pixels, and the ROIs which we work with about 2,500 × 2,000 pixels). More specifically, and from a low level point of view, these ROIs have the following average histogram (Figure 4), where the darker parts represent, among other things, the hippocampus. These ROIs are automatically extracted from the original images, as explained in [START_REF] Mesejo | Automatic Hippocampus Localization in Histological Images using Differential Evolution-Based Deformable Models[END_REF] (Section 3.1). The main contribution of this paper to the state-of-theart is threefold: first, the application of the level set method to histological images, which is not common in literature; second, the hybridization of geometric and parametric DMs (two of the methods use a parametric deformable model to initialise a level set); and third, it is one of the few comparative studies devoted to the segmentation of anatomical structures in histological images (where classical pixel-and intensitybased methods, deformable models (DM), and hybridizations with soft computing techniques, are tested and compared) 1 .

II. THEORETICAL BACKGROUND AND RELATED WORK

In this section, we describe the techniques used and their essential working, as well as the main contributions to the field found in literature.

The methods here compared are the following:

• Active Shape Models (and Iterative Otsu Thresholding Method) refined using Random Forests (ASM+RF) [START_REF] Mesejo | Automatic Hippocampus Localization in Histological Images using Differential Evolution-Based Deformable Models[END_REF], [START_REF] Mesejo | Automatic Segmentation of Hippocampus in Histological Images of Mouse Brains using Deformable Models and Random Forest[END_REF]. This method, published in 2012, uses a training set of shapes to create a model with its limits for the deformation. The approach is inspired by Active Shape Models (ASM) [START_REF] Cootes | The Use of Active Shape Models For Locating Structures in Medical Images[END_REF], [START_REF] Cootes | Active shape models-their training and application[END_REF], but takes into account only grey level intensities (not color nor texture), and uses a medial-based shape representation in polar coordinates, with the objective of creating simple models that can be managed in an easy and fast manner (the complexity of the models has been reduced as much as possible, minimizing the number of control points). Such a parametric model is moved and deformed by Differential Evolution (DE) [START_REF] Das | Differential Evolution: A Survey of the Stateof-the-Art[END_REF] according to an intensitybased similarity function between the model and the object itself. After that, Otsu's thresholding method [START_REF] Otsu | A threshold selection method from gray-level histograms[END_REF] is iteratively applied on every region identified by the located control points and, in every iteration, only the largest segmented component is kept. Therefore, the whole segmentation process can be seen as an intelligent technique for localizing areas where a fast and simple segmentation method is to be applied. Finally, in order to extend the segmentation towards the parts that have not been segmented yet (mainly because the localization did not cover completely the hippocampus), an Ensemble Classifier (in this case Random Forests [START_REF] Breiman | Random forests[END_REF]), trained using textural features, is applied to detect/classify which points belong to the hippocampus and expand, like a region growing method, the segmented area to the regions that were not properly localized. • Soft Thresholding (ST) [START_REF] Aja-Fernandez | Soft thresholding for medical image segmentation[END_REF]. This method, presented in 2010, is based on relating each pixel in the image to the different regions via a membership function, rather than through hard decisions, and such a membership function is derived from the image histogram. As a consequence, each pixel will belong to different regions with a different level of membership.

In a first stage, the normalized histogram of the image under study is calculated (using the maximum grey level to bound the histogram to the interval [0, 1]). Then, a sum of weighted known distributions is fit to the histogram, and each probability distribution represents the probability for a pixel with a certain value to belong to the corresponding region. This segmentation technique has the following advantages: (1) it is totally automatic, and does not require human intervention, which makes it suitable for automatic processes;

(2) the hard decision is postponed to the final stage. So, all the spatial operations are performed before taking into account the different memberships and obtaining the thresholded image; and (3) the spatial operations performed make the thresholding more robust to noise and artifacts.

Having been succesfully applied to CT, MRI and ultrasound, it seemed interesting to apply it also to histological imaging, to check its performance with this image modality. • Otsu Thresholding Method (OTSU) [START_REF] Otsu | A threshold selection method from gray-level histograms[END_REF]. This method, introduced in 1979, is proposed from the viewpoint of discriminant analysis, and it automatically selects an optimal threshold. Otsu's thresholding method involves iterating through all possible threshold values and calculating a measure of spread for the pixel levels on each side of the threshold. The aim is to find the threshold value that minimizes the sum of foreground and background spreads, i.e. it tries to minimize the within-class-variance by maximizing the between-class-variance (which is less computationally intensive to calculate).

The algorithm assumes that the image to be thresholded contains two classes of pixels, or that its histogram is bimodal (e.g. foreground and background), and then calculates the optimum threshold separating those two classes so that their combined spread (intra-class variance) is minimal. Since ASM+RF uses Otsu in an iterative and local way, it also seemed interesting to test the proper Otsu method directly over the whole ROI, in order to check the improvement obtained by the ASM+RF method when including prior shape knowledge, as well as texture information.

• Geodesic Active Contours (GAC) [START_REF] Caselles | Geodesic active contours[END_REF]. This method is a geodesic approach for object segmentation which connects classical 'snakes' based on energy minimization and geometric active contours based on the theory of curve evolution. The technique, introduced in 1997, is based on active contours that evolve in time according to intrinsic geometric measures of the image: the evolving contours naturally split and merge, allowing the simultaneous detection of several objects and both interior and exterior boundaries.

The Partial Differential Equation of the GAC is the following:

u t = α • div(g∇u/|∇u|)|∇u| + β • g|∇u| (1)
where g is a positive and strictly decreasing function, ∇ is the gradient operator computed on the image I, div is the divergence operator (that measures the magnitude of a vector field's source or sink at a given point) , and α and β are the contour (internal force) and expansion (external force) weights, respectively. The higher α, the more regularized/smoothed the contour. If β is positive it expands outwards, if negative it expands inwards (moving faster and crossing gradients, the smaller its value).

The problem is formalized as the minimization of the energy:

E(C) = |C|ε 0 g(|∇I(C(s))|)ds, (2)
where |C| ε is the Euclidean length of a contour C, s is the arc length of the contour, and ds is equivalent to |C (q)|dq, being C(q) a parameterized planar curve. This method is equivalent to the minimization of the length of the curve C according to a Riemannian metric, and such a metric depends here on the local gradient of the image I.

"Geodesic active contours" minimize Equation 2 via a gradient descent scheme and a level-sets representation of the curve. Unfortunately, the method is sensitive to initialization and the global minimum of Equation 2 is not always found. In this work, two implementations of GAC have been tested. The first one uses as initial contour the whole image, while the second one, called GAC+ASM, employs the localization method of ASM+RF (i.e. an ASM using DE) to create the initial contour of the geometric DM. Except the initial contour, all the remaining parameters were exactly the same in both approaches. • Chan&Vese Level Set Model (C&V) [START_REF] Chan | Active contours without edges[END_REF] 3.

∂φ ∂t = δ ε (φ)[µ div( φ | φ| )-ν-λ 1 (u 0 -c 1 ) 2 +λ 2 (u 0 -c 2 ) 2 ]
(3) The parameter µ is a scaling parameter for the length of the curve represented in terms of curvature. The smaller µ, the more the length of the curve can increase without penalizing the minimization. This allows the model to detect small objects and holes. The larger µ, the less freedom is there for the curve to increase in length, thus, it will only be able to detect larger objects. The parameter ν is also a scaling term for the area of the curve, but it can be set to 0 in the Euler-Lagrange derivation, since µ is sufficient to scale the curve according to the objects that need to be detected. Finally, λ 1 and λ 2 are weighting parameters for the forces inside the curve and outside the curve respectively. Since we want to give both forces equal weight, we have set λ 1 = λ 2 = 1 in all experiments. Therefore, the evolution of the curve is influenced by two terms: the "curvature" that regularizes the curve and makes it smooth during evolution, and the "region term" that affects its motion. The main theoretical advantages of this approach are that it automatically detects interior contours and both type of boundaries (with and without gradient), and the initial curve can be placed anywhere in the image. Also in this case, like in GAC, two implementations have been tested. The first one uses the whole image as initial contour, and the second one employs the localization results obtained by ASM to create the initial contour. Regarding the segmentation of structures in histological images, our work is related to the investigations of Senyukova et al. [START_REF] Senyukova | Automated atlas-based segmentation of Nissl-stained mouse brain sections using supervised learning[END_REF]. There, the segmentation is treated as a classification problem using RF and Markov Random Fields, which refine the results at the pixel level, but the method requires previous knowledge about the reference slice associated to that image. In [START_REF] Riklin-Raviv | Propagating distributions for segmentation of brain atlas[END_REF], the segmentation using three-dimensional Gaussian mixtures and Level Sets is carried out slice-by-slice, where the successful segmentation of one section provides a prior for the subsequent one, assuming that the segmentation of few sparsely sampled slices is done manually (so it is not a completely automatic method).

On the other hand, when CV, GAC or ASM are hybridized with DE, this implies a combination between DM and evolutionary computation techniques. There are also some examples in this regard, like [START_REF] Ballerini | Genetic snakes: active contour models by genetic algorithms[END_REF], where "Genetic snakes" are active contour models with an optimization procedure based on GA, or [START_REF] Mcintosh | Medial-based deformable models in non-convex shape-spaces for medical image segmentation using genetic algorithms[END_REF], where a GA evolves a population of medialbased shapes, using prior shape knowledge to produce feasible deformations while also controlling the scale and localization of these deformations. Also in [START_REF] Ghosh | A genetic algorithmbased level set curve evolution for prostate segmentation on pelvic ct and mri images[END_REF] a GA is used to perform level set curve evolution using texture and shape information to automatically segment the prostate in CT and MRI pelvic images. Other example can be found in [START_REF] Novo | Topological active models optimization with differential evolution[END_REF], where DE, in combination with a greedy search algorithm, is used to evolve Topological Active Nets for CT image segmentation.

III. EXPERIMENTAL RESULTS

It is important to notice that all methods used in the experimentation, apart from ASM+RF, and GAC and CV when using ASM to obtain the initial contour, are deterministic techniques, i.e. given a particular input, they will always produce the same output. This is not the case of the other three approaches, since stochastic methods, like DE, are embedded in these algorithms, and can exhibit different behaviors on different runs. Due to this, it is essential to run such algorithms several times to estimate and compare their performances.

In order to evaluate the accuracy of all methods, tests on real images were run on an Intel R Core TM i5-2410M CPU @ 2.30GHz with 4.00 GB of RAM, using MATLAB as programming language. Table II reports the results obtained using 26 real images, employing 25 runs per image (for nondeterministic methods), for a total of 2054 experiments. The results are organized per blocks of five rows each. Every block corresponds to a different metric [START_REF] Babalola | Comparison and evaluation of segmentation techniques for subcortical structures in brain mri[END_REF]: Dice Similarity Coefficient, Hausdorff Distance, Jaccard Index, True Positive Rate, and False Positive Rate. Inside every block, the average, median, standard deviation, best value, and worst values are shown. These calculations are computed over 26 or 650 values depending on the nature of the method under consideration (deterministic or non-deterministic, respectively).

A ground truth was created by manually segmenting the hippocampus in these 26 images. In order to avoid erroneous or incomplete manual segmentations, these were supervised by an expert in molecular biology. Every image was manually segmented 5 times and the intersection and union images were calculated. Nevertheless, the difference between the maximum (union) and minimum (intersection) area segmented in these (equally valid) segmentations showed an average difference of 28.44% with a standard deviation of 5.61%. Simultaneous truth and performance level estimation (STAPLE) [START_REF] Warfield | Simultaneous truth and performance level estimation (staple): An algorithm for the validation of image segmentation[END_REF] was used to create a consensus ground truth, since it is considered one of the state-of-the-art methods to create a ground truth from several manual segmentations [START_REF] Biancardi | A comparison of ground truth estimation methods[END_REF].

A post-processing phase was applied to remove isolated pixels: a median filtering (15x15), only in the case of ST and Otsu, and the removal of connected components with less than 5000 pixels in all methods, except ASM+RF. The manually tuned parameters used on each method are shown in Table I. For GAC, C&V, GAC+ASM and C&V+ASM, a median filtering (10x10) was applied during pre-processing in order to regularize the surface of the image.

The test set is composed of 26 representative genes downloaded from the ABA database (see Figure 2). In particular, the images downloaded, and used for the comparative study, were the ones associated with the following gene symbols: Wbp7, 1300018I05Rik 86, 0610010D24Rik, 1300018I05Rik 94, 5430437P03Rik, Atrx, B230215L15Rik, Nmt1, Wars, Tubb3, Trem1, 3830406C13Rik, A030009H04Rik, Atp1b2, Azin1, Camk2a, Camk2b, Camk2g, Cd9, Cutl2, Gad1, Gapdh, Gfap, Mbp, Slc17a7, and Zim1.

IV. DISCUSSION AND CONCLUSIONS

In order to check the statistical significance of the results obtained, a Kruskal-Wallis test was performed with a level of confidence of 0.01. Since the normality and homoscedasticity assumptions were not accomplished, as checked through the application of Kolmogorov-Smirnov and Bartlett's tests, nonparametric tests were used. 2 The p-value is near zero, suggesting that at least one sample median is significantly different from the others. Some conclusions can be extracted from the numerical data in Table II and from some of the representative results included in Figure 5. The first conclusion could be the difficulty of tackling satisfactorily these images, since established and wellknown segmentation methods, like C&V and GAC, did not obtain as good results as one could expect. In particular, GAC presents problems with images that have small variations of their gradient values, i.e. when the difference, in terms of intensity, between the hippocampus and the background is low.

According to Table II, the best results were obtained by ASM+RF, being the best in all metrics but in TP rate, where C&V, C&V+ASM and Otsu achieved better results. This is explained by the tendency of these three methods to oversegment the images and, therefore, to segment the hippocampus and many other regions in the image. In general terms, the best methods to segment these images have been ASM+RF, C&V+ASM, GAC+ASM, and ST.

From the values of HD interesting conclusions can also be drawn. This metric represents the largest of all the distances from a point in one set (the ground truth) to the closest point in the other set (the result of the automatic segmentation). It is interesting to check how ST has a higher (worst) HD with respect to C&V, even if it obtained better DSC results than the latter. This could be justified by the absence in the former of any kind of shape restrictions or regularization terms, that makes that the pixels are segmented only taken into account their grey level values.

On the other hand, Otsu and ST are the fastest methods since they work only, and directly, on the histogram. The difference in performance between Otsu and ASM+RF demonstrate the great improvement obtained with the introduction of prior knowledge.

Geometric DM have shown a good ability of changing topology and evolving naturally the contour, but on the other hand, they are very sensitive to initialization, pre-processing stages 3 , and parameter selection. GAC and C&V suffer with bad initializations (the use of a better initial contour favors a fast convergence to the contour and better results), but this affects specially to GAC (as can be cheked in the improvement obtained when using the localization method based on ASM). If the initial contour is located partially inside and outside of the structure of interest (as sometimes happens due to a bad localization), only the boundaries will be segmented and not the internal part of the hippocampus, because the high gradient values are located on the boundary and not in the internal part of the structure.

Regarding both geometric approaches (C&V and GAC, with and without localization stage), the results would improve in great manner if statistical shape priors for level sets were used (thickness, length, statistical information about the image to segment,...). The only piece of prior knowledge that GAC and C&V are using now is the existence of high gradient values in the boundary of the objects to segment, and the prior knowledge included in the parameter selection (weighting the different components of the energy). Another aspect to take into account is that, sometimes, the inner part of sg is completely segmented without separating the two extremes of this region. This owes to the fact that the regularization factor does not allow the evolving contour to "pass" through small narrow regions and, therefore, it fills completely the empty space of sg. Furthermore, the problem with GAC and C&V in noisy images is that, if the images are smoothed too much, it will be impossible for the method to stop at the gradients, but if they are not smoothed, the noise would provoke the segmentation of many noisy parts in the image.

As previously said, ASM+RF obtained the best results, but it is fair to underline its ad-hoc nature. It needs a training set of shapes to create the template and its possible deformations, and it also needs a training set of textural patterns for the expansion phase. Also it is not able to manage topological changes in a natural way, as geometric DMs can do. So, it is a very accurate method but it sacrifices its general applicability to segment other shapes.

Possible future works include the use of a local C&V instead of a global one [START_REF] Wang | An efficient local chan-vese model for image segmentation[END_REF]; introducing segmentation through registration algorithms in the comparison; and trying to create a new segmentation algorithm combining the good features of different methods (like the combination of the ability of geometric DM to manage, naturally, topological changes, with the global search capabilities of Evolutionary Algorithms, or initializing them using existing segmentation methods).

Fig. 1 .

 1 Fig. 1. Regions in Hippocampal Formation (HPF).

Fig. 3 .

 3 Fig. 3. Average histogram of the images under study.

Fig. 2 .

 2 Fig. 2. All images used in the experiments. From up to down, and from left to right: Gad1, Camk2a, 0610010D24Rik, Gfap, Mbp, B230215L15Rik, 3830406C13Rik, A030009H04Rik, Gapdh, 5430437P03Rik, Atp1b2, Tubb3, Azin1, 1300018I05Rik 94, Trem1, Cd9, Wars, Slc17a7, Zim1, Atrx, 1300018I05Rik 86, Cutl2, Camk2g, Camk2b, Wbp7, and Nmt1.

Fig. 4 .

 4 Fig. 4. Box-plot representing the DSC for all methods. The best results are obtained by ASM+RF.

Fig. 5 .

 5 Fig. 5. Results for some images (from up to down: Zim1, Gfap, Camk2a, and Atp1b2). The order of methods is (from left to right): OTSU, ST, C&V, C&V+ASM, GAC, GAC+ASM and ASM+RF. True Positives are represented in white, False Positives in green and False Negatives in red.

  . While many segmentation methods heavily rely in some way on edge detection, the "Active Contours Without Edges" method by Chan and Vese (2001) ignores them completely. They have proposed an active contour model that can detect objects whose boundaries are not necessarily defined by gray level gradients. Instead, the segmentation boundary is represented implicitly with a level set function, which allows the segmentation to handle topological changes more easily than explicit snake approaches, and the method is based on trying to separate the image into regions based on intensities. The Chan&Vese model is a region-based model for detecting objects in an image, based on a restriction Mumford-Shah model which divides an image into regions and represents each region by a piecewise constant (the minimal partition problem). More formally, the Euler-Lagrange equation (as derived in the original paper) representing the time motion of the curve C is shown in Equation

TABLE I PARAMETERS

 I USED IN TESTING DIFFERENT SEGMENTATION TECHNIQUES.

	ST	OTSU	ASM+RF	C&V	GAC
	L = 2 regions	No parameters to tune	Cr = 0.9	number of iterations = 500	number of iterations = 500
	Relative max normalization		F = 0.7	µ = 0.1 (length term)	β = -1 (expansion weight)
			Uniform Crossover	ν = 0 (area term)	α = 3 (contour weight)
			DE/target-to-best/1	λ1 = λ2
			Population Size = 80
			Iterations = 250
			Median Filter [25×25]
			RF with 500 trees

TABLE II SEGMENTATION

 II RESULTS USING 5 DIFFERENT METRICS: DICE SIMILARITY COEFFICIENT (DSC), JACCARD INDEX (JI), HAUSDORFF DISTANCE (HD), TRUE POSITIVES (TP), AND FALSE POSITIVES (FP). THE BEST RESULTS FOR EVERY METRIC ARE SHOWN IN BOLD LETTERS.

						Methods
			OTSU	ST	C&V	C&V+ASM	GAC	GAC+ASM ASM+RF
		Average	0.5755	0.6798	0.6164	0.6488	0.6373	0.7469	0.8882
		Median	0.5788	0.7505	0.6906	0.7577	0.6985	0.7971	0.8975
	DSC	Std	0.1820	0.2276	0.2432	0.2052	0.2723	0.1993	0.0347
		Best	0.8341	0.9084	0.9120	0.8909	0.9189	0.9198	0.9350
		Worst	0.2868	0.0687	0.2183	0.2358	0.0000	0.0214	0.7422
		Average	0.4263	0.5531	0.4869	0.5121	0.5205	0.6283	0.8006
		Median	0.4090	0.6007	0.5274	0.6099	0.5367	0.6626	0.8141
	JI	Std	0.1823	0.2330	0.2465	0.2117	0.2782	0.2080	0.0541
		Best	0.7154	0.8321	0.8382	0.8033	0.8500	0.8514	0.8779
		Worst	0.1674	0.0356	0.1225	0.1336	0.0000	0.0108	0.5901
		Average	936.2499	660.4281	566.7963	397.8478	543.2649	247.6652	110.8424
		Median	988.7435	723.4520	602.9961	378.8250	565.7718	212.2499	93.9415
	HD	Std	220.5734	354.8354	340.3507	253.4867	392.5249	174.1109	47.6210
		Best	506.7790	57.0088	33.5410	40	67.0820	30	44.7214
		Worst	1.3837e+003	1.3911e+003	1.1942e+003	1.0330e+003	1.3067e+003	721.9591	296.0152
		Average	89.4691	80.0974	95.6464	97.2764	85.1330	82.4358	87.3440
		Median	95.4543	90.2218	98.5369	99.2808	97.4342	95.2620	87.7858
	TP	Std	15.5791	23.5970	9.0898	6.8649	25.9801	25.5602	4.6220
		Best	99.9584	99.9659	100	100	99.9364	99.9364	94.9942
		Worst	35.8155	5.7835	64.1182	64.1182	0.0000	1.3607	60.5404
		Average	54.7366	34.2408	50.3108	47.9808	41.3426	28.1579	9.2846
		Median	58.7041	26.9258	46.5393	38.9488	27.4716	22.1582	7.8321
	FP	Std	19.3433	26.0209	25.0742	21.5470	30.7668	19.4360	5.5229
		Best	25.6864	1.5257	11.6259	16.0377	7.8152	3.6460	3.0877
		Worst	81.2674	91.5339	87.7290	86.5843	100	95.0183	32.4049

The methods here compared are all automatic segmentation techniques (manual/semi-automatic methods have not been included).

Note that the Kruskal-Wallis test is a nonparametric version of the classical one-way ANOVA, and an extension of the Wilcoxon rank sum test to more than two groups. It returns the p-value for the null hypothesis that all samples are drawn from the same population.

A brief example to show the influence of the type of pre-processing applied. If the size of the median filter were 20x20 instead of 10x10 the results would get worse for CV (DSC Mean: 0.586159 -DSC Median: 0.639210 -DSC Std: 0.256346) and they would improve for GAC (DSC Mean: 0.668008 -DSC Median: 0.841995 -DSC Std: 0.344408)

ACKNOWLEDGMENT Pablo Mesejo is funded by the European Commission (MIBISOC Marie Curie Initial Training Network, FP7 PEOPLE-ITN-2008, GA n. 238819). All the images were downloaded from the Allen Mouse Brain Atlas [http://mouse.brain-map.org]. Seattle (WA): Allen Institute for Brain Science. c 2009.

The authors want to thank Ferdinando Di Cunto, from the University of Torino, for facilitating the ground truth for these images.


Achal Agrawal

Samson Lasaulce

Olivier Beaude

Raphael Visoz

A Framework for Optimal Decentralized Power Control with Partial CSI

published or not. The documents may come   L'archive ouverte pluridisciplinaire

A Framework for Optimal Decentralized Power Control with Partial CSI Achal Agrawal 1 , Samson Lasaulce 1 , Olivier Beaude 1 and Raphael Visoz 2 1 L2S (CNRS -CentraleSupélec-Univ. Paris Sud 11), France 2 Orange Labs, Issy-les-Moulineux, France 1 {agrawal, lasaulce, beaude}@lss.supelec.fr, 2 raphael.visoz@orange-ftgroup.com Abstract-In this paper, we use a recent information theoretical result to develop a general framework for finding optimal power control policies in the case of interference channels. The aforementioned result characterizes the achievable payoffs for an N -Agent (transmitters in our application) coordination problem with a certain information structure. We then provide an algorithm which exploits the characterization of achievable payoffs by conditional probability distributions to find optimal decision functions for the transmitters. Due to its general nature, the developed framework is conducive for applications to diverse scenarios in wireless communications. In this article, we restrict our attention to the case of decentralized power control in interference channels for different utility functions namely sum-rate, sum-energy and sum-goodput. Our approach has the following salient features: 1) The method proposes optimal power control functions for any given utility function as opposed to ad-hoc solutions for different utilities proposed in the literature, and 2) Noise in the channel estimation is taken into account, thus providing robust optimal solutions.

I. INTRODUCTION

One of the key contributions of this paper is to exploit the recent theorem derived in [START_REF] Larrousse | Coordinating partiallyinformed agents over state-dependent networks[END_REF] to find power control functions which may exploit the available knowledge optimally (in the long run). As an example, we treat the well studied problem of optimal power control in interference channels ( [START_REF] Zander | Performance of optimum transmitter power control in cellular radio systems[END_REF], [START_REF] Yates | A framework for uplink power control in cellular radio systems[END_REF], [START_REF] Stanczak | Fundamentals of resource allocation in wireless networks: theory and algorithms[END_REF]) and show that not only does our framework attain various payoffs identical to the state-of-the-art, it also provides the optimal power control functions in diverse scenarios. The scenarios considered resemble those of [START_REF] Gohary | Robust IWFA for open-spectrum communications[END_REF], [START_REF] Anandkumar | Robust rate maximization game under bounded channel uncertainty[END_REF]. These works propose ad-hoc solutions for specific utilities but not a generic framework which directly provides a power control function that aims at exploiting the available and arbitrary information as well as possible. Indeed our approach is very general and can be used to treat more complex and interesting scenarios, e.g. by considering different kinds of channel state information available at transmitters, different common objectives, robust power control taking the noisy communications into account, vectorial optimization etc.

The generality of the approach is by virtue of exploitation of some recent information theoretical results [START_REF] Larrousse | Coordinating partiallyinformed agents over state-dependent networks[END_REF]. A similar approach, albeit with a different information structure and for only 2 transmitter-receiver pairs was proposed in [START_REF] Larrousse | Implicit coordination in 2-agent team problems. application to distributed power allocation[END_REF]. The information structure considered in this paper is more adapted for the application proposed under realistic assumptions, the important ones being 1) No knowledge of future realizations of the nature state w.r.t. which the agents are coordinating, 2) Any number of agents with only a partial and noisy information about the current nature state, and 3) Decision functions must be found in a decentralized manner leading to simpler design architectures to implement.

Indeed, for the case of power control in wireless interference channels, these assumptions have been paid some attention. The necessity of minimizing the information required at the transmitters for coordination led to approaches proposed in [START_REF] Paul De Kerret | Best-response team power control for the interference channel with local CSI[END_REF], [START_REF] Gjendemsjø | Binary power control for multi-cell capacity maximization[END_REF], [START_REF] Gjendemsj | Binary power control for sum rate maximization over multiple interfering links[END_REF]. They consider only local channel state information at the transmitters (CSIT). This is also important for reducing the complexity of the optimizations to be performed. For sum-rate maximization, the scenario considered here is comparable to a single-carrier version of the iterative water filling algorithm (IWFA) [START_REF] Scutari | Asynchronous iterative water-filling for gaussian frequency-selective interference channels[END_REF]. In the case of sum-energy utility, channel inversion strategy was proposed by [START_REF] Goodman | Power control for wireless data[END_REF], whereas the multi-carrier version of the problem was considered by [START_REF] Meshkati | A game-theoretic approach to energy-efficient power control in multicarrier CDMA systems[END_REF]. Decentralization of the power control, or for that matter for any wireless network design, is the reason why game theory has been extensively applied to such problems [START_REF] Lasaulce | Game theory and learning for wireless networks: fundamentals and applications[END_REF]. We thus compare our results to game theoretic equilibria analysis as most of the decentralized optimization literature only guarantees convergence to Nash equilibria.

The article follows the structure from more general results to a specific application of the framework proposed. To this effect, Sec. II describes the performance limits of a general coordination problem characterizing the optimal achievable payoffs. The characterization lays the foundation for an algorithm proposed in Sec. III for finding decentralized decision policies for the coordinating transmitters. The preceding sections lay the general framework. In Sec. IV-A we describe the specific problem of power control in interference channels. A numerical analysis is conducted for diverse scenarios and the results are compared to existing literature in Sec. IV-B. Concluding remarks are provided in Sec. V.

II. PERFORMANCE LIMITS OF AN N -TRANSMITTERS

INTERFERENCE COORDINATION PROBLEM Consider a coordination problem where N transmitterreceiver pairs are trying to coordinate the power emitted by them X i ∈ X i , i ∈ N = {1, ..., N }. The coordination is done with respect to a nature state represented by X 0 ∈ X 0 with the goal of optimizing a common payoff function w(x 0 , x 1 , x 2 , ...x n ) over a long time period T . The nature state X 0 for applications to wireless communications can be channel gain coefficients g ij = |h ij | 2 , with g ij being the channel gain of the link between transmitter i and receiver j. X 0 = (g 11 , ..., g N N ) is a random state which affects the common payoff function for the system and is not controlled by the coordinating transmitters. The realizations of the nature state X 0,t at each time instant t are i.i.d. and follow a probability distribution ρ. In wireless communications ρ is typically an exponential distribution for each channel gain g ij . Our aim is to characterize all the achievable expected payoffs under a certain information structure over a long time period T , T → ∞.

Firstly, we need to define the information structure under consideration. At every instant t, transmitter i is assumed to have an image or a partial observation S i,t of the nature state X 0,t with respect to which all transmitters are coordinating. In the case of our application, this could be knowledge of only local CSIT, i.e. transmitter i observes a noisy version (in general) of only the direct link channel gain g ii . One could imagine other kinds of information available at the transmitters: for example, transmitter i observes all the links g ji , ∀j. The observations S i,t are assumed to be generated by a memoryless channel whose transition probability is denoted by i (S i,t |X 0,t ). All transmitters have to choose the power emitted by them X i,t based on this information received. Formally, the sequence of decision functions for agent i, σ i,t , is defined as:

σ i,t : S t i × U -→ X i (1) (s i (1), s i (2)...s i (t), u(t)) -→ x i (t) (2)
where S t i = S i (1)×S i (2)...×S i (t) is the discrete observation alphabet till the instant t and X i is the power emitted for transmitter i with |S i |, |X i | < ∞. U is the alphabet of the auxiliary variable U which is discussed in more detail later. s i (t), u(t) and x i (t) are the realizations of the corresponding variables at instant t.

The problem is said to be decentralized as each transmitter chooses its power independently based on the information received by it. Since the channel gains are generated by a random process, and thus in general also the corresponding chosen power levels, the quantity to be optimized is the expected objective function:

E Q [w(X)] = x ∈ X w(x 0 , x 1 , ..., x N )Q(x 0 , x 1 , ..., x N ) (3) where X = (X 0 , X 1 , ..., X N ), x = (x 0 , x 1 , ..., x N ), and X = X 0 × N i=1 X i . Q(x 0 , x 1 , ..., x N )
is the joint probability distribution of the variables affecting the payoff. An important point to note is that since expectation is a linear operator, optimizing the expected payoff is equivalent to finding the optimal distribution Q(x 0 , x 1 , ..., x N ). However, the optimization problem is not so trivial as indeed there are certain restrictions on the distributions Q that are implementable given the imposed information structure. We now define the notion of an implementable distribution.

Definition 1 (Implementability). Let the information structure be as defined in [START_REF] Larrousse | Coordinating partiallyinformed agents over state-dependent networks[END_REF]. The probability distribution Q(x 0 , x 1 , ..., x N ) is implementable if there exist decision functions σ i,t such that as T → +∞, we have for all x ∈ X ,

1 T T t=1 Q X0...X N ,t (x 0 , ..., x N ) -→ Q(x 0 , ..., x N ) (4) where Q X0X1...X N ,t = Q X1,...X N |X0,i × ρ is the joint distribution induced by σ i,t at stage t.
As seen before, the expected payoff is characterized by the probability distribution Q over all the variables that intervene in the payoff function. Thus, the time averaged expected payoff w is said to be achievable, if and only if the corresponding distribution Q is implementable. The following theorem characterizes precisely the distributions that are implementable under the information structure (1).

Theorem 1.

[1] Assume the random process X 0,t to be i.i.d. following a probability distribution ρ and the available information to the transmitters S i,t to be the output of a discrete memoryless channel obtained by marginalizing the conditional probability (s 1 , ..., s N |x 0 ) . An expected payoff w is achievable in the limit T → ∞ if and only if it can be written as:

w = x 0 ,x 1 ,...x N , u,s1,...s N ρ(x 0 )P u (u) (s 1 , ...s N |x 0 )× N i=1 P Xi|Si,U (x i |s i , u) w(x 0 , x 1 , ...x N ) (5)
where U is an auxiliary variable which can be optimized.

The auxiliary variable U is an external lottery known to the transmitters beforehand, which can be used to achieve better coordination. An example for how this variable U can help coordinate better in the case of power control is shown in Sec. IV-A.

Theorem 1 helps us find the best achievable payoffs for the expected utility (4) with arbitrary information to the transmitters as long as it follows the structure (in terms of when the information is available to whom) defined in equation [START_REF] Larrousse | Coordinating partiallyinformed agents over state-dependent networks[END_REF]. However, this theorem does not provide optimal sequences of decision functions σ i,t for the time averaged expected payoff. In the following section, we shall give an algorithm which helps find a suboptimal solution for the decision functions. The suboptimality is due to the high complexity of a multilinear optimization problem. Nonetheless, In Sec.IV-B we see that the solutions obtained using this algorithm match the state-of-the-art solutions proposed for different utilities discussed.

III. ALGORITHM FOR FINDING SUB-OPTIMAL LOW-COMPLEXITY DECISION FUNCTIONS

As seen from Theorem 1, optimal performance for the coordination problem depends only on the vector of conditional probabilities (P X1|S1,U , ..., P X N |S N ,U ). Thus, it suffices to find an optimum vector of lotteries for every action X i possible to obtain the best achievable payoff. However, one can see that this task is computationally very demanding. Instead, we apply a suboptimal approach by searching for the conditional probabilities in a distributed manner, thus reducing the complexity of the search. To this effect, we apply sequential best-response dynamics [START_REF] Lasaulce | Game theory and learning for wireless networks: fundamentals and applications[END_REF]. This consists of each transmitter choosing the best corresponding power given that others keep their powers constant, and each transmitter doing so sequentially within an iteration. Note that by doing so, the optimization problem becomes linear, and hence its solution lies at a vertex, i.e. either (P Xi|S1,U ) = 0 or (P Xi|S1,U ) = 1. Thus the search for conditional probability distributions simplifies to search for decision functions

f i : S i × U -→ X i (s i (t), u(t)) -→ x i (t)
An important point to note is that we search for a stationary strategy, i.e. f i does not depend on the time instant t. This simplification is practical as it leads to simpler design structures and easier to find in terms of complexity. Also, the auxiliary variable U , which serves as a coordination key, is exchanged offline beforehand and thus does not entail any signalling cost. The procedure outlined below is performed many times till convergence is achieved, typically in a few iterations. The number of iterations required for convergence of course depends on the number of agents coordinating, but it scales up very slowly as best-response dynamics normally converges very fast [START_REF] Lasaulce | Game theory and learning for wireless networks: fundamentals and applications[END_REF]. Furthermore, since we are considering a common 'team' payoff w(x 0 , x 1 , ..., x N ), the convergence of best-response dynamics is guaranteed [START_REF] Monderer | Potential Games[END_REF]. Proposition 1. Algorithm 1 converges for a common performance criteria w(x 0 , x 1 , ..., x N ).

Proof: The result can be proved by induction or more generally by calling for an exact potential game property (the latter argument may hold in the more general case in which the transmitter have different performance criteria).

For further justification of the simplification made by finding functions f i instead of conditional probability distributions, note that the optimization problem is linear w.r.t. each component of the vector (P Xi|S1,U ) that we are searching. The problem is therefore multilinear and Algo 1 will converge to the vertex of the polytope defined by the constraints. Note however that we might lose out on finding the global optimum by doing so. Nonetheless, the power control functions obtained by the simplification are sufficient to compare with the existing literature which typically provide power control functions and not mixed randomized strategies.

To apply best-response dynamics, we rewrite the expected utility given by Theorem 1 in the following way:

w = x 0 ,x 1 ,...x N , u,s1,...s N ρ(x 0 )P U (u) (s 1 , ..., s N |x 0 )× (6)
N i=1 P Xi|Si,U (x i |s i , u) w(x 0 , x 1 , ...x N ) = ai,bi,u θ ai,bi,u P Xi|Si,U (x ai |s bi , u) (7)
where a i , b i , u are the respective indices of x i , s i , u and

θ ai,bi,u   a0 ρ(x a0 ) i (s i |x a0 ) a-i w(x a0 , x a1 , ...x a N )× b-i i =i i (s b i |x a0 ) i =i P X i |S i ,U (x a i |s b i , u)   P U (u) (8)
where a -i , b -i are the indices which represent a i , b i being constant, while all the other indices are summed over.

Written in this form, the linearity of the problem w.r.t. any one of the variables P Xi|Si,U (x i |s i , u) is even more explicit. Note that we have also assumed the independence of the observation channels of all the agents, i.e.

(s|x 0 ) = 1 (s 1 |x 0 ) × . . . × N (s N |x 0 ).
For interference channels, this assumption is reasonable as all transmitters receive independent feedback of the channel gains.

Algorithm 1: Proposed decentralized Algorithm for finding a control function inputs : An important point to note is that the transmitters can run this algorithm offline as all they need to know are the channel statistics ρ. This helps in using the power control functions for the entire timeslot during which the channel statistics remain constant. Thus we exploit the channels from the very start of a timeslot as opposed to taking some time initially to find the optimal functions online. The latter case is considered in some decentralized solutions such as the algorithms based on water-filling techniques proposed in [START_REF] Gohary | Robust IWFA for open-spectrum communications[END_REF], [START_REF] Anandkumar | Robust rate maximization game under bounded channel uncertainty[END_REF].

X i , ∀i ∈ {0, ..., N }, w(x 0 , x 1 , ...x N ), ∀x, ρ(x 0 ), S|X0 (s|x 0 ), ∀x 0 , f init i , ∀i ∈ {1...N }, P U (u), ∀u ∈ U output: f * i (s i , u i ), ∀i ∈ {1...N } Initialization: f 0 i = f init i , iter = 0, iter max = 100 while ∃i : f iter-1 i -f iter i 2 ≥

IV. APPLICATION TO POWER CONTROL IN INTERFERENCE CHANNELS

In the previous sections, we developed the theoretical framework, giving intuition as to what the defined quantities might represent in diverse optimization problems of wireless communications. In this section, we concentrate on the problem of power control in interference channels [START_REF] Lasaulce | Game theory and learning for wireless networks: fundamentals and applications[END_REF] and provide numerical analysis of some special cases for the above problem.

A. System Model

Consider N single-antenna Transmitter-Receiver pairs Transmitter i , Receiver i , i ∈ N = {1, 2..., N } with a single band interference channel with the channel gains being g ij , i, j ∈ {1, 2..., N } and g ij ∈ G. G is considered to be a discrete set and represents the alphabet of the possible channel gains. A well accepted model of statistics for the channel gains g ij = |h ij | 2 is Rayleigh fading. In this model, due to central limit theorem, the real and the imaginary components of h ij follow a normal distribution, and thus the channel gain g ij follows exponential distribution (|h ij | follows a Rayleigh distribution). We also assume that all channel gain distributions are independent of each other and that their realizations are i.i.d. The transmitters transmit at discrete power levels P i ∈ P i quantized uniformly (in dB) between [0, P max ]. The utility for the pair Transmitter i , Receiver i depends on its signal to interference plus noise ratio (SINR) γ i (P , g). Here P is the vector with each component P i being the power emitted by the i th transmitter. g is the vector containing all the channel gain coefficients g ij . More precisely,

γ i (P , g) = P i g ii σ 2 + j =i P j g ji (9)
where σ 2 is the noise variance. Without loss of generality, we will consider σ 2 = 1, and change P max to regulate the signal to noise ratio (SNR). SNR is defined to be the ratio

P max /σ 2 .
We shall consider only local CSIT, i.e. the direct channel gain g ii available at each transmitter i, based on which they have to choose the power level to emit at. For evaluation of the solutions found using our method, we shall consider three different payoff functions. These payoff functions have vastly different proposed solutions in the literature. Our method helps find the state-of-the-art solutions for all the different payoffs by simply using a different payoff function for the Algo. 1. The following common payoffs are considered:

• Sum-rate: w R (P , g) = i∈N log(1 + γ i ) ; • Sum-goodput: w G (P , g) = i∈N Ω(γ i ) ; • Sum-energy: w E (P , g) = i∈N Ω(γ i ) P i .
Typical functions for Ω are Ω(γ) = e -c γ [START_REF] Belmega | Energy-efficient precoding for multiple-antenna terminals[END_REF], where c > 0 or Ω(γ) = (1 -e -γ ) M where M ≥ 1 is the packet length [START_REF] Goodman | Power control for wireless data[END_REF]. In simulations, we chose the former function and further supposed c = 1 as it only changes the optimal solutions by a multiplicative constant. Note that Proposition 1 holds for the above payoffs as we are considering sums of individual payoffs which trivially satisfy the required condition for convergence.

The following table explicits the relation between the input variables for Algo. 1 and the corresponding quantities for the application to power control defined above. In Theorem 1,

General model Power Control Application

Nature State -X 0 Channel state -(g 11 , ..., g N N ) Decision of Transmitter i -X i Power emitted

P i ∈ P i w(x 0 , x 1 ...x N ) w R (P , g), w G (P , g), w E (P , g) ρ(x 0 ) N i=1 N j=1 e -g ij , g ij = E(g ij ) (s i |x 0 ) s i = g ii or s i = g ii f init i f i (s) = Pmax, ∀s ∈ S

TABLE I CORRESPONDENCE BETWEEN THE GENERAL FRAMEWORK AND ITS

APPLICATION TO POWER CONTROL we introduced an auxiliary variable U which could help in achieving better coordination. For better intuition, we show how a simply constructed external lottery helps transmitters coordinate in the case of power control.

Consider only N = 2 transmitter-receiver pairs communicating over a single band interference channel while maximizing the sum-rate w R (P , g). From [START_REF] Gjendemsjø | Binary power control for multi-cell capacity maximization[END_REF], we know that the optimal solution (P * 1 , P * 2 ) ∈ {(P max , 0), (0, P max ), (P max , P max )}. However, if transmitters try to maximize their individual utility log(1 + γ i ), the Nash equilibrium is (P max , P max ) [START_REF] Lasaulce | Game theory and learning for wireless networks: fundamentals and applications[END_REF]. Now consider the following probability distribution for a random variable U with the realizations u ∈ U = {u 1 , u 2 , u 3 }.

P U (u) =    P(u = u 1 ) = α 1 , P(u = u 2 ) = α 2 , P(u = u 3 ) = α 3 = 1 -(α 1 + α 2 ) . (10)
A large sequence of realizations for U are drawn with the probability distribution specified in [START_REF] Gjendemsj | Binary power control for sum rate maximization over multiple interfering links[END_REF] and all transmitters share this drawn sequence. Now consider a trivial function common to all transmitters ∆ : U → P which maps every realization of U , u i to one of three possible optimal power control vectors {(P max , 0), (0, P max ), (P max , P max )} which form the set of possible actions P . Under this setup, one can do Monte-Carlo simulations for a large number of channel gain realizations, and find the optimal probabilities (α * 1 , α * 2 ) by doing an exhaustive search. This is done by searching over a discrete space α i ∈ {0, 0.01, 0.02, ..., 0.99, 1}, i ∈ {1, 2}. In Fig. 1, we plot the average payoff w for the coordinated policy using the pre-drawn lottery U versus the SNR (which is equal to P max due to normalization of channel noise variance σ 2 ). We compare the average payoffs for the coordinated policy as well as for Nash equilibrium and see that at higher SNRs, it is indeed preferable to coordinate using U . This is not surprising at all, as at higher SNRs, the interference generated is too high, whereas time-sharing helps avoid this.

In the following section, we provide the power control simulation setup and discuss the power control functions given by our approach. We do not consider the pre-drawn lottery U in our analysis, although in general, U will only help improve the coordination performance. This is because the solution where transmitters ignore the coordination key is always possible.

B. Numerical Analysis

1) Simulation Setup: We implement Algo. 1 for the case of power control in a single band interference channel [START_REF] Lasaulce | Game theory and learning for wireless networks: fundamentals and applications[END_REF]. To this effect, we use the correspondence established in Sec. IV-A. We consider the case when N = 2. Note that exponential distributions for channel gain statistics g ij are characterized only by their mean g ij . We shall consider the standard scenario where ∀i ∈ N , g ii = 1, ∀i, j ∈ N , i = j, g ij ) = 0.3.

As mentioned before, the noise variance of the channel has been normalized to one. Thus, the SNR is regulated via P max = 20 dB. The possible power levels emitted by transmitters (P i ) are 50 uniformly discrete points between

P i ∈ [-20, 20] dB .
Since the setup considered is symmetric for all Transmitter -Receiver pairs, we plot the power control functions f i (s i ), only for the transmitter i. We suppose that the observation alphabet S i = G ii . Note that since we do not consider the auxiliary variable u from here on, the function f i is represented with only one argument s i . In particular, it is seen that for maximizing sum-energy, the obtained power control function exhibits a threshold under which transmission should not occur for a given transmitter.

2) Results: In Fig. 2, we plot the power control functions f i obtained using Algo. 1 for the case N = 2 against g ii the direct channel gain between the i th Transmitter-Receiver pair. We do not consider noise in the channel estimation for the moment, and thus in this case s i = g ii . Three different utilities are considered, sum-rate, sum-energy, and sum-goodput.

For the case of sum-rate, it is known that binary power control P i ∈ {P min , P max } is optimal for 2 Transmitter-Receiver pairs [START_REF] Gjendemsjø | Binary power control for multi-cell capacity maximization[END_REF]. Moreover, as shown in [START_REF] Paul De Kerret | Best-response team power control for the interference channel with local CSI[END_REF], optimal power control functions with only local CSIT (g ii ) amounts to P i (g ii ) = P min if g ii ≤ g * and P i (g ii ) = P max otherwise. It is reassuring to find that our results verify this. Thus, in the case of sum-rate with local CSIT, we find exactly the same results as state-of-the-art.

In the case of sum-energy, we see that all the available power P max is not used. Indeed, even in the case of only one Transmitter-Receiver pair, the optimal power control function is c/g 11 . We see that there is also a threshold value of g ii below which, P i = P min . Above the threshold value, the function is similar to the optimal solution obtained in the case of only one Transmitter-Receiver pair. The threshold function is also seen as solution in the case of sum-packet rate. This is because the utility function is not monotonous w.r.t. SINR. To the best of our knowledge, the power control functions found for sum-energy are new. Solutions in literature do not propose thresholding of power control functions to help reduce interference in case of bad direct channel gains. Admittedly, this requires greater analysis and comparisons with the state-of-the-art.

To investigate the function for sum-energy further, in 3 we compare the function obtained using our algorithm with that proposed by [START_REF] Goodman | Power control for wireless data[END_REF]. We see that unlike [START_REF] Goodman | Power control for wireless data[END_REF], we find a threshold below which emitting no power is more optimal. We also see that the presence of noise makes the power control function more uniform. Our approach is shown to be much more robust to noise than the one proposed by [START_REF] Goodman | Power control for wireless data[END_REF] N=2 Fig. 3. Our algorithm reveals the shape of good power control functions in the presence of interference. In contrast with related works on energyefficiency [START_REF] Goodman | Power control for wireless data[END_REF], our work shows that thresholding is required to manage interference efficiently.

function 1/g ii is very sensitive to noise at low values of g ii . Relative performance gain (w f -w N ash /w N ash in % ) of our algorithm when compared to average payoffs obtained by Nash equilibrium for sum-rate and sum-energy. The curve for sum-energy saturates as at high SNR, Pmax is not utilized as power emitted is much lower.

In Fig. 4, we compare the performance of the power control function with that of Nash equilibrium for different SNRs in the case of sum-energy and sum-rate. Thus, on the y -axis we see the relative performance gain (w fw N ash /w N ash in % ) of our algorithm when compared to average payoffs obtained by Nash equilibrium. Current stateof-the-art proposes solutions which converge to the Nash equilibrium. Thus, the solutions are generally not optimal while maximizing common utilities that are sum of the individual utilities. The evaluations for obtaining average expected payoffs was done for 10 6 channel realizations. They were drawn according to the distribution N i=1 N j=1 e -g ij . The channel gain statistics g ij were the same as the ones used for running the algorithm.

We see that for sum-energy, there is no performance gain at higher SNRs. This is because, as seen from Fig. 2, the optimal power used is very low for all g ii . However, we see that when compared to Nash equilibrium, we obtain a significant gain of around 55%. In the case of sum-rate, performance obtained matches with [START_REF] Paul De Kerret | Best-response team power control for the interference channel with local CSI[END_REF] since the power control functions obtained are identical. Not suprisingly, at higher SNRs, Algo. 1 does much better than Nash equilibrium, which is all transmitters transmitting at P max generates too much interference.

One of the advantages of our framework is that it provides optimal power control functions even for noisy channel estimates. We illustrate this in Fig. 5 where power control functions for different levels of noise in channel estimation are plotted. This noise simulates the error in estimation due to noise in feedback transmission, or just simply statistical estimation error. The noise for the simulations is gaussian, i.e. s i = g ii +z where Z ∼ N (0, σ 2 z ), with σ z ∈ {0, 1, 3}. As expected, the power control functions become more uniform at higher noise levels, as the information received is less reliable, and thus transmitters emit at a power level which maximizes the utility for an average case. Influence of channel estimation noise on the power control functions. We see that as noise increases, the optimal power emitted becomes more uniform. This is intuitive as higher the noise, higher the uncertainty in the observation leading to same power emitted for all channel gains observed in the asymptotic case.

V. CONCLUSIONS

The proposed framework was shown to be relevant in diverse scenarios of single band interference channel for finding optimal power control functions. Moreover, the power control functions depend only on local CSIT, thus having the merit of being implementable in a completely decentralized manner. Also, the solutions obtained take noise in the estimation of the channel gain into account. All the above features illustrate the generality of our approach in tackling problems of power control for maximizing sum-utility functions.

However, the framework can be exploited further for tackling other problems in wireless communications as well. For example, one could consider the problem of power allocation in a multi-band interference channel. Also, the auxiliary variable U was not exploited, which in general will only make the solution more optimal. One could also consider different information available at transmitter i: g ij , ∀j.

Fig. 1 .

 1 Fig.1. The figure shows the performance in terms of sum-rate with: the proposed coordination mechanism using a public key shared beforehand; the single-carrier counterpart of the iterative water-filling algorithm. Without additional signalling, a significant gain is obtained.

Fig. 2 .

 2 Fig.2. The figure represents the power control functions f i provided by Algo. 1 for the three performance metrics under consideration in this paper. In particular, it is seen that for maximizing sum-energy, the obtained power control function exhibits a threshold under which transmission should not occur for a given transmitter.

  as

  Fig.[START_REF] Stanczak | Fundamentals of resource allocation in wireless networks: theory and algorithms[END_REF]. Relative performance gain (w f -w N ash /w N ash in % ) of our algorithm when compared to average payoffs obtained by Nash equilibrium for sum-rate and sum-energy. The curve for sum-energy saturates as at high SNR, Pmax is not utilized as power emitted is much lower.

  Fig. 5.Influence of channel estimation noise on the power control functions. We see that as noise increases, the optimal power emitted becomes more uniform. This is intuitive as higher the noise, higher the uncertainty in the observation leading to same power emitted for all channel gains observed in the asymptotic case.

ACKNOWLEDGMENT

Authors of this paper and particularly A. Agrawal would like to thank project LIMICOS -ANR-12-BS03-0005 for financing the project.


Youssef S G Nashed
email: nashed@ce.unipr.it

Pablo Mesejo
email: pmesejo@ce.unipr.it

Roberto Ugolotti

Jérémie Dubois-Lacoste
email: jeremie.dubois-lacoste@ulb.ac.be

Stefano Cagnoni
email: cagnoni@ce.unipr.it

A comparative study of three GPU-based metaheuristics

Keywords: Global Continuous Optimization, Particle Swarm Optimization, Differential Evolution, Scatter Search, GPGPU

In this paper we compare GPU-based implementations of three metaheuristics: Particle Swarm Optimization, Differential Evolution, and Scatter Search. A GPU-based implementation, obviously, does not change the general properties of the algorithms. As well, we give for granted that GPU-based implementation of both algorithm and fitness function produces a significant speed-up with respect to a sequential implementation. Accordingly, the main goal of this work has been to fairly assess the efficiency of the GPU-based implementations of the three metaheuristics, based on the statistical analysis of the results they obtain in optimizing a benchmark of twenty functions within a prefixed limited time.

Introduction

Modern graphics hardware has gained an important role in the area of parallel computing, since it has been used to accelerate general computations (General Purpose Graphics Processing Unit -GPGPU -programming), in addition to playing its natural role. CUDA TM (Compute Unified Distributed Architecture) is a parallel computing environment by nVIDIA TM which exploits the massively parallel computation capabilities of its massively parallel GPUs. In particular, CUDA-C [START_REF][END_REF] is an extension of the C language that allows development of GPU routines (termed kernels), that can be executed in parallel by several different CUDA TM threads, following the Single Instruction Multiple Thread model.

Among the stochastic approaches to continuous optimization, Evolutionary Algorithms (EAs) [START_REF] Eiben | Introduction to Evolutionary Computing[END_REF] and Swarm Intelligence [START_REF] Bonabeau | Swarm Intelligence: From Natural to Artificial Systems[END_REF] algorithms offer a number of attractive features: robust and reliable performance, global search capability, virtually no need of specific information about the problem to solve, easy implementation, and, above all, implicit parallelism.

In this paper, we compare the GPU implementations of three real-valued population-based optimization techniques: Particle Swarm Optimization (PSO) [START_REF] Kennedy | Particle Swarm Optimization[END_REF], Differential Evolution (DE) [START_REF] Storn | Differential Evolution -a simple and efficient adaptive scheme for global optimization over continuous spaces[END_REF], and Scatter Search (SS) [START_REF] Glover | Heuristics for integer programming using surrogate constraints[END_REF] 1 . They have been tested on a benchmark of 20 numerical problems (also implemented on GPU), comprising unimodal/multimodal and separable/non separable functions (see Table 3). The main contributions of this paper are: i. the first implementation, to the best of our knowledge, of Scatter Search in CUDA TM ; ii. a novel parallel version of DE that solves some of the problems of previous implementations; iii. the study of these three metaheuristics on a broader benchmark than those usually adopted to test GPU implementations; iv. an unbiased evaluation of the effectiveness of the GPU implementation, as each metaheuristic has been subjected to the same parameter optimization method and run for a pre-fixed limited amount of time.

2 Basics of the three metaheuristics

Particle Swarm Optimization

Particle Swarm Optimization (PSO) is a stochastic optimization algorithm which simulates the behavior of bird flocks. A set of particles (or solutions) move through a "fitness" function domain (search space) seeking the function optimum (best fitness value). Each particle's motion is described by two simple discrete-time equations which regulate the particles' velocity and position:

v n (t) = w • v n (t -1) + c 1 • rand() • (BP n -P n (t -1)) + c 2 • rand() • (BLP n -P n (t -1)) P n (t) = P n (t -1) + v n (t)
where P n (t) and v n (t) are the position and velocity of the n th particle at time t; c 1 , c 2 and w (inertia factor) are positive constants, rand() returns random values uniformly distributed in [0, 1], BP n is the best-fitness location visited so far by the particle, and BLP n is the best-fitness location visited so far by any particle in its neighborhood, which may include a limited set of particles or even coincide with the whole swarm.

Differential Evolution

Differential Evolution has recently been shown to be one of the most successful EAs for global continuous optimization [START_REF] Das | Differential Evolution: A Survey of the State-of-the-Art[END_REF][START_REF] Vesterstrom | A comparative study of differential evolution, particle swarm optimization, and evolutionary algorithms on numerical benchmark problems[END_REF]. Unlike traditional EAs, DE perturbs the current population members with the scaled differences of randomly selected and distinct individuals. In the first iterations the elements are widely scattered in the search space and have a great exploration ability. As optimization proceeds, the individuals tend to concentrate in the regions of the search space with better fitness values, so the search automatically focuses on the most promising areas. In DE, every element acts as a parent vector, for which a donor vector is created. In the original version of DE, the donor vector for the i th parent (X i ) is generated by combining three random and distinct elements X r1 , X r2 and X r3 . The donor vector V i is calculated as:

V i = X r1 + F • (X r2 -X r3 )
where F (scale factor) is a parameter that strongly influences DE's performances and typically lies in the interval [0.4, 1]. After mutation, every parent-donor pair generates one child (trial vector) by means of a crossover operation. Two kinds of crossover are typically used: binomial (uniform) and exponential. Crossover is applied with a certain probability Cr (crossover rate) that, like F , is one of the control parameters of DE. Then, the trial vector is evaluated and its fitness is compared to its parent's: the better survives while the other is discarded.

Scatter Search

Scatter Search is based on a systematic combination between solutions (instead of randomized, as is usual in EAs) taken from a subset of the population, named the "reference set", that is usually significantly smaller than a typical EA population. SS is composed of five structural "blocks" or methods (see Figure 1):

1. Diversification Generation: a population of solutions P is generated, having a certain degree of quality and diversity. The reference set R is then drawn from P , and includes the |R 1 | solutions with best fitness, and the |R 2 | solutions from the reference set that are farthest from P according to the Euclidean distance (hence, |R| = |R 1 | +|R 2 |); the evolution process acts only on R; 2. Solution Combination: in most problems a specific method to combine solutions is needed, which can be applied to all solutions or only to selected ones (e.g., the best solutions, and/or randomly selected ones). In many cases an existing crossover operator, borrowed from other EAs, is employed; 3. Subset Generation: the procedure deterministically generates subsets of R, to which the combination method is applied. 4. Improvement: an improvement method (typically a local search) is applied to the original solutions and/or to combined solutions; 5. Reference Set Update: once a new solution is obtained it replaces the worst solution in R only if it improves the quality of the reference set in terms of fitness and/or diversity;

Parallel implementation

The first parallel versions of PSO relied on multiprocessor parallel machines or cluster of computers. With the introduction of GPUs, research shifted towards GPU-based parallel PSO (GPU PSO) to alleviate multi-processor and cluster systems inefficiencies, such as network overhead, shared memory access, etc. In 2009 and 2010, respectively, the first implementations of PSO and DE using nVIDIA CUDA TM were developed [START_REF] De Veronese | Swarm's flight: Accelerating the particles using C-CUDA[END_REF][START_REF] De Veronese | Differential evolution algorithm on the GPU with C-CUDA[END_REF]. Also in 2009, a hybrid between GPU PSO and pattern search, aimed at enhancing the convergence of PSO, was presented in [START_REF] Zhou | GPU-based parallel particle swarm optimization[END_REF]. After that, other implementations of DE have been developed [START_REF] Krömer | Many-threaded implementation of differential evolution for the CUDA platform[END_REF][START_REF] Zhu | Massively parallel differential evolution-pattern search optimization with graphics hardware acceleration: an investigation on bound constrained optimization problems[END_REF], and fast versions of PSO have been implemented by relaxing the synchronicity constraints between particles' evaluation and update [START_REF] Mussi | GPU-based asynchronous particle swarm optimization[END_REF]. The early GPU PSO implementations suffered from a coarse-grained parallelization (one thread per particle), that neglected the opportunity to compute the fitness function, usually the most time-consuming process, in parallel over the problem dimensions. This aspect has been improved by the PSO implementation evaluated in this work, first introduced in [START_REF] Mussi | Evaluation of parallel particle swarm optimization algorithms within the CUDA architecture[END_REF]. In it, a thread manages a single dimension of each particle, adding a further level of parallelism. Similar inefficiencies characterized the early implementations of DE, such as a partially sequential implementation of the fitness function and random number generation. These problems were addressed by [START_REF] Krömer | Many-threaded implementation of differential evolution for the CUDA platform[END_REF] using four kernels executed sequentially. Our present implementation uses only one kernel for generating the trial vectors, and another for fitness evaluation and migration. In addition, we offer three different mutations and two kinds of crossovers, while early GPU-based DE considered only one mutation strategy (DE/rand/1) and one kind of crossover. Regarding SS, to the best of our knowledge, ours is the first parallel implementation of this metaheuristic.

The performance of CUDA TM code depends chiefly on the thread configuration, number of kernels, and memory access schemes used. All methods presented here aim at exploiting fast-access local and shared memory instead of slower global memory to the greatest possible extent, considering the number of kernels as the main criterion to assess how well an algorithm can be parallelized.

PSO is divided into three kernels described in [START_REF] Mussi | Evaluation of parallel particle swarm optimization algorithms within the CUDA architecture[END_REF], while DE, as mentioned earlier, can be implemented as two kernels. Each thread of the first kernel performs the following instructions:

generate two or three distinct random numbers on the GPU, according to the mutation strategy; -calculate an element of the donor vector from the population members randomly selected in the previous step; -decide whether to include the donor or the parent element in the trial vector, based on the type of crossover and the crossover rate, Cr.

The second DE kernel evaluates all trial vectors simultaneously in shared memory and, if the fitness has improved, it replaces the parent with the offspring. Clearly, SS is not as inherently parallel as the two other metaheuristics (see Figure 1). In SS a diverse population is first initialized and evaluated; diversity is simulated by generating uniform random values for each dimension over the whole search space. Then, to build the reference set R, a parallel sort operation is required to find R 1 , followed by another kernel that calculates pairwise Euclidean distances between solutions in P -R and R, sequentially adding the solutions that are farthest from the reference set for |R 2 | iterations. As for selection and crossover, a kernel selects all solution pairs in the reference set for mating, and combines them through the BLX-α crossover , generating two distinct solutions chosen with α set to (0.5 + λ) and (0.5λ), respectively. The combined solutions make up the pool, to which a parallel implementation of the Solis & Wets search method [START_REF] Wets | Minimization by random search techniques[END_REF] is then applied as improvement method. For the last step, we compared two methods for updating the reference set, one of which considers both quality and diversity as in [START_REF] Duarte | Hybrid scatter tabu search for unconstrained global optimization[END_REF], while the other updates the reference set with the best |R| solutions in (R ∪ pool). The latter yielded better results in terms of both speed and accuracy, as proven by the automatic tuning process described in the following section.

Experiments

We give for granted that the parallel metaheuristics we consider are faster than the corresponding sequential versions on sufficiently large problems, while their accuracy is the same for identical configurations, because they implement the same algorithm. Many comparisons between the accuracy of the sequential versions of the algorithms [START_REF] Duarte | Hybrid scatter tabu search for unconstrained global optimization[END_REF][START_REF] Vesterstrom | A comparative study of differential evolution, particle swarm optimization, and evolutionary algorithms on numerical benchmark problems[END_REF] have already been made, but they gave insights on their intrinsic features rather than on the computational efficiency of possible implementations. In this work, we evaluate both quality and speed of their parallel versions, analyzing the accuracy they can achieve in a limited amount of time, to assess the degree of parallelization that each of them allows to reach.

The algorithms we compared have a number of parameters that affect both accuracy and parallelism. "Manual" parameter tuning is time consuming and may introduce a bias in comparing an algorithm with a reference, due to better knowledge of the algorithm under consideration and to possible different time spent tuning each of them. Therefore, the automatic tuning of all three algorithms was performed using the irace software package [START_REF] López-Ibáñez | The irace package, iterated race for automatic algorithm configuration[END_REF], to find the configurations that yielded the best results in a given time: we set this time to one second, since it is generally short enough to avoid reaching full convergence with all three methods, allowing one to compare their short-term performances.

The tuner was run on all 20 functions with a budget of 30000 experiments, each being one run of one configuration on one function with a termination criterion of one second. Since the functions have different fitness ranges, a rank-based test is prferable to based on the solutions' mean values. Accordingly, the Friedman test was used to discard significantly worse configurations. We tuned the parameters for 30-dimensional problems, and assumed that such configurations are good also for lower-sized ones. Table 1 displays the parameters that have been tuned for each algorithm, and the best corresponding values.

We compared our results to the values that are most commonly used in literature. For instance, the authors in [START_REF] Das | Differential Evolution: A Survey of the State-of-the-Art[END_REF] suggest F ∈ (0.4, 0.95) and Cr ∈ (0.9, 1) for multimodal separable functions (the most common ones in our benchmark); we obtained similar results. Regarding PSO, in most papers, c 1 = c 2 = 2.0 [START_REF] Poli | Particle swarm optimization[END_REF], while our automatic tuning set them to slightly smaller values.

To evaluate both the effectiveness and the efficiency of the three parallel implementations, tests on 20 numerical benchmark functions (see Table 3) were run on a 64-bit Intel(R) Core TM i7 CPU running at 2.67GHz using CUDA TM v. 4.1 on a nVidia GeForce GTS450 graphics card with 1GB of DDR memory and compute capability 2.1 [START_REF][END_REF]. Table 2 reports the results obtained executing 100 runs per function (6000 independent runs) and setting 1 second as the only termination criterion. The first column is the function under consideration. The following ones are divided into two blocks according to the number of dimensions (10 and 30). Within each block, the mean best fitness and the standard deviation over all runs are reported for each method. Results reported on a grey background highlight those case in which the median over 100 runs obtained by the method is significantly better than the other methods, according to the Kruskal-Wallis test, with a confidence level of 0.01.

Discussion

The results reported in Table 2 and Figure 2 allow one to draw some conclusions about the behaviour of the three parallel metaheuristics. Conforming with previous results obtained by sequential implementations, DE obtained the best results, sometimes tied with some other method, in 35 out of the 40 experiments performed, while PSO was the best method, sometimes tied with some other method, in 20 out of 40 functions, its main drawback being its tendency to stagnate and find sub-optimal solutions more often than DE, even if a higher number of function evaluations is run. Regarding SS, whose first parallel implementation is presented here, it obtained the best result in 12 out of 40 problems; however, this metaheuristic, which is not as parallelizable as the other methods, as reflected by the number of kernels, has achieved better performance over multimodal non-separable problems and time-consuming fitness functions, like Katsuura.

All tests were run with a temporal limit of one second, a short time in which all three methods can generally obtain results close to the optima without reaching full convergence. Figure 2 shows that PSO requires almost three times as many fitness function evaluations as DE to converge on 30-dimensional problems. It is important to notice that the population size in PSO is also almost three times as large as in DE, which justifies the larger number of fitness evaluations. However, this may represent a shortcoming only if applied to larger-dimensional functions than those considered in this work. where yi = 1 +

x i -1

Fig. 1 :

 1 Fig. 1: Block diagram of the Scatter Search Algorithm.

F((yi - 1 ) 2

 12 [START_REF] Mussi | GPU-based asynchronous particle swarm optimization[END_REF][START_REF] Duarte | Hybrid scatter tabu search for unconstrained global optimization[END_REF].12] n f griewank (f rosenbrock ) (xi, xi+1), xn = x0 where F (x, y) = 0.5 + Ai -Bi(x)) 2 , x = [x0, . . . , xn-1] Ai, Bi(x) defined as in[START_REF] Suganthan | Problem definitions and evaluation criteria for the CEC 2005 special session on real-parameter optimization[END_REF]. 10 sin 2 (πyi + 1) + (yn-1 -1) 2 1 + 10 sin 2 (2πyn-1)

Table 1 :

 1 Automatically-tuned parameter values used to test different optimization techniques.

Table 2 :

 2 Results on the 20 functions.

			10 dimensions					30 dimensions
	DE		PSO		SS		DE		PSO		SS
	Avg	Std	Avg	Std	Avg	Std	Avg	Std	Avg	Std	Avg	Std

Table 3 :

 3 Benchmark functions. For every function, the table shows the name, the range of the search space, the formula, the multimodality (multimodal, unimodal) and the separability (separable, non separable). All minima are in {0} n . B} , i = 0, . . . , n -1, x = [x0, . . . , xn-1], Ai, B defined in[START_REF] Suganthan | Problem definitions and evaluation criteria for the CEC 2005 special session on real-parameter optimization[END_REF].

		Name	Range	Formula
	f0 f1 f2 Sum of Squares Sphere Elliptic f3 HyperEllipsoid f4 Schwefel 2.22 f5 Zakharov f6 Schwefel 1.2	[-100, 100] n [-100, 100] n [-1, 1] n [-1, 1] n [-10, 10] n [-10, 10] n [-100, 100] n	+	n-1 i=0 n-1 i=0 n-1 i=0 n-1 i=0 n-1 i=0 n-1 x 2 i (10 6 ) ix 2 i i D-1 x 2 i-1 i i 2 • x 2 |xi| + n-1 i=0 i=0 xi 2 + n-1 i=0 0.5 • i • xi 2 4 |xi| n-1 i=0 n-1 i=0 i j=0 xj 2	0.5 • i • xi 2 2	+	U S U S U S U S U S U S U NS
	f7	Schwefel 2.6	[-100, 100] n	max {Aix -U NS
	f8	Dixon-Price	[-10, 10] n	(x0 -1) 2 +	n-1 i=1	i • 2xi 2 -xi-1	2	U NS
	f9	Rastrigin	[-5.12, 5.12] n		n-1 i=0	x 2 i -10 • cos(2πxi) + 10	M S
	f10 Schwefel 2.26	[-500, 500] n 418.9829 • n +	n-1 i=0	xi • sin	|xi|	M S
	f11	Katsuura	[-1000, 1000] n		n-1 i=0	1 + (i + 1)

The code can be downloaded from http://sourceforge.net/p/libcudaoptimize[START_REF] Nashed | libCudaOptimize: an Open Source Library of GPU-based Metaheuristics[END_REF]

, i = 0, . . . , n -1 M NS

Acknowledgments

Youssef S. G. Nashed, Pablo Mesejo and Jérémie Dubois-Lacoste are funded by the European Commission (Marie Curie ITN MIBISOC, FP7 PEOPLE-ITN-2008, GA n. 238819). Roberto Ugolotti is funded by Compagnia di S.Paolo and Fondazione Cariparma. The authors want to thank Thomas Stützle for his comments and Enrico Viappiani for the help in implementing CUDA TM code.