Durrieu

Spatial quantification of vegetation density from terrestrial laser scanner data for characterization of 3D forest structure at plot level

Keywords: Terrestrial lidar, forest canopy, 3D model, architecture, stand structure

Precise description of forest 3D structure at plot level is required for sustainable ecosystem management. However a detailed structure description from traditional field measurements is tedious. We propose an innovative method to quantify the spatial distribution in 3D of forest structure from terrestrial lidar data. The method rests on the hypothesis that the normalized number of laser returns within a given volume element is proportional to the density of vegetation material inside this volume. The developed model is based on analysis made inside Svoxels (spherical voxels) and allows calculating a spatialized vegetation density index. The model was first tested on two different scans of the same forest plot. The resulting vegetation density index well represents the vegetation structure as observed within the lidar point cloud. Quantitative analyses confirmed a global consistency of the results within and between scans. A slight decrease in the density index with distance and the distribution of between-scan differences raise the possibility of a bias that could be explained by a distance-dependant attenuation of the signal. Future work will focus on improving our algorithm and correcting biases. These results are promising for the development of quantitative measures of the 3D forest structure.

Introduction

Precise description of 3D structure of forests is useful for timber resource monitoring, ecosystem management and preservation, or improved understanding on ecosystem functioning. However the spatial complexity of forests makes structure measurement very difficult, particularly since structure is not a satisfyingly defined feature [START_REF] Fleck | Terrestrial lidar measurements for analysing canopy structure in an old-growth forest[END_REF]. Measurements from ground plots are fastidious from traditional methods to the point where a complete 3D description is not practical. Terrestrial lidars provide 3D data offering much more details compared with traditional field inventories opening up new opportunities to derive metrics closely linked to forest structure and to reduce time and costly field measurements [START_REF] Hopkinson | Assessing forest metrics with a ground-based scanning lidar[END_REF]).

Terrestrial lidars were originally developed for civil engineering (see [START_REF] Lichti | Ground-based laser scanners: operation, systems and applications[END_REF] for examples of systems and applications). Recent studies expanded their use on tree or stand structure measurements. Most of them focused on estimating traditional field-based forest parameters. [START_REF] Hopkinson | Assessing forest metrics with a ground-based scanning lidar[END_REF] first demonstrated that it is possible to locate and identify individual trees with high precision and to measure total tree height and diameter at breath height (dbh). Tree heights were however underestimated of about 1.5 m when compared with field validation data. This was mostly due to low sampling density at the upper canopy level caused by occlusion effects of the signal and a suboptimal survey protocol. Results for mean dbh differed by only 1 cm from tape measurements. Similar results were obtained by other authors for both height and dbh measurements using semi-automatic data extraction methods [START_REF] Watt | Measuring forest structure with terrestrial laser scanning[END_REF][START_REF] Fleck | Terrestrial lidar measurements for analysing canopy structure in an old-growth forest[END_REF][START_REF] Wezyk | Terrestrial laser scanning versus traditional forest inventory. First results from the Polish forests[END_REF]. Other forest parameters such as stem density, total basal area, gross and merchantable timber volume were also estimated from terrestrial lidar data with a good agreement when compared with traditional field measurements (e.g. volume estimations within 7% of the traditional field estimations [START_REF] Hopkinson | Assessing forest metrics with a ground-based scanning lidar[END_REF]). Other efforts dealt with automatic tree location and height, dbh, stand basal area or timber volume estimations [START_REF] Aschoff | Describing forest stands using terrestrial laser-scanning[END_REF][START_REF] Bienert | Tree detection and diameter estimations by analysis of forest terrestrial laser scanner point clouds[END_REF][START_REF] Király | Tree height estimation methods for terrestrial laser scanning in a forest reserve[END_REF][START_REF] Wezyk | Terrestrial laser scanning versus traditional forest inventory. First results from the Polish forests[END_REF].

The very high sampling rate of terrestrial laser systems allows generating detailed 3D models of canopy therefore opening up the possibility to analyze fine scale stand structure, foliage distribution, canopy light transfer or leaf area indices that are important to understand and model forest function and dynamic. However few studies have demonstrated the interest of such systems for ascertaining parameters beyond those from the traditional inventories. As an exception, [START_REF] Fleck | Terrestrial lidar measurements for analysing canopy structure in an old-growth forest[END_REF] proposed a method to quantify canopy projection far much precise than the 8-point canopy projection from a ground operator used in traditional inventories. As other non-traditional measures, [START_REF] Danson | Forest canopy gap fraction from terrestrial laser scanning[END_REF] proposed a method to estimate canopy directional gap fraction and [START_REF] Van Der Zande | Influence of measurement set-up of ground-based LiDAR for derivation of tree structure[END_REF] an approach for vegetation profile reconstruction. Studies using terrestrial lidar show much opportunity for developing new methods for forest canopy metrics that will take full advantage of terrestrial lidar datasets. One of the main issues will be to solve the problem of the distance-dependent varying point density from the lidar returns. This paper introduces an innovative approach to analyze the vegetation structure from 3D point clouds acquired with terrestrial lidar. The method quantifies the 3D spatial distribution of forest canopy material in volume elements (~dm level). It makes available operational calculations linking the 3D point cloud recorded by a terrestrial lidar with the spatial distribution of the vegetation. This study was also realized considering the link between airborne lidar data and field data with the aim of improving information extraction from airborne lidar data on forested areas. Indeed airborne lidars proved capable to estimate the spatial distribution of forest parameters such as height, crown area, timber volume or biomass at both tree or stand level [START_REF] Lim | Lidar remote sensing of forest structure[END_REF]. However these airborne estimates require local calibration through acquisition of field data.

Method

Study area and field data

The main study site is part of a National Environmental Observatory (ORE Draix) located in the southern part of the French Alps. It is part of the Haute-Bléone state forest, mainly composed of black pine (Pinus nigra) planted in the 1880's to protect against soil erosion. Most of the stands are even-aged and mature. Elevations range from 802 to 1263 m. Traditional field inventory was conducted during December 2007 within circular plots of 15 m and 9 m radius. Within the plots the following characteristics were measured for all the trees with dbh > 7 cm: dbh, total and timber heights, crown base height, crown diameter and tree position.

Data acquisition with the terrestrial lidar

Terrestrial lidar surveys were made on March 2008 using an ILRIS-3D system (Optech Inc, Toronto, Canada). The system measures the laser returns within a window 40° wide in both horizontal and vertical directions. The laser emits and measures light at 1,500 nm. Point density of each scan is controlled by the operator. The system can register the intensity and distance for either the first or the last backscattered signal. In our study, we selected primarily the last returns considering that they would provide a better statistical representation of the vegetation distribution compared with first returns. However, first and last returns were recorded at some particular system base stations (i.e. system location) for comparison and quality assessment. The ILRIS-3D base stations were selected outside the plot at varying distance from the plot centre and separated by an angle of about 120° relating to the plot center. Artificial targets (polystyrene spheres with 8 cm diameter) were distributed within the plot and measured using differential GPS and total station to improve the alignment (co-registration) and the georegistration of the scans acquired from different base stations.

Method developed for quantifying the spatial distribution of vegetative elements

This study's objective was to develop an algorithm capable to calculate vegetation density by relating lidar returns visible in the form of point clouds. The point density needs to be locally transformed into density of vegetation components. Vegetation density was related to a density index extracted from the lidar measurement. In order to estimate this index, we divided the plot-space into constant volume elements (voxels). For each voxel, we calculated (1) the number of lidar returns within the voxel and (2) the number of laser beams passing through the voxel. The density index of each voxel is given by the ratio (1) / (2). Our method has two spatial characteristics: a regularly spaced grid of voxel centers and the use of spherical voxels.

Regular 3D grid and spherical voxels

Voxel centers were arranged on a 3D grid regularly distributed along x, y and z axes. The grid was georeferenced in the Lambert III conformal conic coordinate system and was used to process each scan of a same plot. Computations from all scans of the plot could therefore be compared and integrated. Before processing each scan, the Lambert III grid is changed into the Cartesian system of the scan. The transformation model is computed using (1) The Lambert III coordinates of target centers, measured on the field (total station + DGPS), and (2) the Cartesian coordinates of the targets, measured on each scan by fitting a spherical shape on its corresponding point clouds. The 3D Reshaper ® software was used for that purpose. A minimum of 4 spheres was required for computing the transformation model.

Data acquisition with the terrestrial lidar follows a spherical geometry. We therefore adopted a spherical geometry to simplify computations on the resulting point cloud from lidar measurements. Lidar position was taken as the origin of the spherical system. The space illuminated by the lidar was already divided into voxels. Therefore each voxel center was associated with a spherical coordinate (r, θ, φ) and bounded with the following conditions:

1. 4 angles: θ min = θ -dθ, θ max = θ + dθ, φ min = φ -dφ and φ max = φ + dφ, 2. 2 distances: r min = r -dr, r max = r + dr, with dr set to half the grid resolution. This new volume is referred to as the spherical voxel or Svoxel (Figure 1). The geometry of the Svoxel differs only slightly from the one of its corresponding voxel if dθ and dφ are set to ensure a constant volume of Svoxels (V = R 3 , with R the 3d resolution of the grid). The resulting Svoxels have the following properties:

1. Distortion of a Svoxel compared to the reference voxel is proportional to r (cf. figure 1), 2. Distortion of a Svoxel increases when angles θ and φ increase, 3. Svoxels are not strictly contiguous. Small overlaps or gaps can occur which are more important for larger values of θ and φ, 4. For a given center point, the Svoxels generated from different base station locations will not strictly overlap due to slight changes in shape and orientation. The highest differences will occur when comparing Svoxels from scan with a 45° (modulo 90°) difference between viewing angles. Even with these properties, differences between voxels and Svoxels remain small and it is thus assumed that they are not detrimental to precise density index computation.

A b

Figure 1: Shape of a Svoxel at 1 m (a) and at 3 m (b) for 50 cm grid resolution.

Algorithm to calculate density index of the grid points

The following algorithm was implemented to calculate the density index of each Svoxel in the lidar scanning field of view:

1. Generation of a 3D regularly spaced grid in Lambert III at a resolution R, 2. Projection of the grid in the sensor Cartesian system, 3. Switch scan point cloud and grid into spherical system, 4. For each point of the grid : 1. Computation of the theoretical number of laser beams (N theorical ) passing through the Svoxel using the point density selected for the scan, 2. Identification of the number of returns inside the targeted Svoxel (N inside : points satisfying the 4 angles and 2 distances equations), 3. Identification of the number of laser beams intercepted before the targeted Svoxel (N before : points satisfying the 4 angles equation with a distance lower than r min ), 4. Computation of the vegetation density index (I density = N inside / (N theorical -N before )*100).

If N theorical -N before = 0, a nodata value is set. If N theorical -N before is lower than a given threshold Ts, results are considered as non-significant. 5. Output of results in Lambert III. 6. Steps 2 to 5 can be reapplied to other scans of a same plot acquired from other base stations.

Data analysis and validation

For this preliminary study the algorithm was applied to 2 of the 8 available scans on a circular 15 m radius plot with a medium tree density (66 stems/ha) and located on a flat area. Scan density was set to 6.24 mm (resp. 7.02 mm) at 15 m for scan 1 (resp. scan 2) and the last returns were recorded. Three Svoxel resolutions were selected: 0.25, 0.5 and 1 m. Results were first evaluated from a preliminary visual assessment where Svoxels with a positive and significant density index were visualized on the lidar 3D point clouds of selected trees. Preliminary tests allowed us to adopt a value of 50 for the threshold defined for non-significant values (Ts).

Then two sets of procedures were realized: 1. In order to evaluate the result consistency inside a given scan, several stand crowns located at various distances from the base station 1 were extracted and the distribution of positive and significant density indices were analysed. Results on four black pines and one Spanish fir (Abies pinsapo) were compared (cf. Fig 2). 2. Density index values obtained from two different base stations were also compared to evaluate the consistency of the results between different scans. This preliminary analysis defined if results from multiple scans can be compared and merged.

Results and discussion

Visual analysis

The method gave visually consistent results. Figure 3 compares the point cloud from the original scan and the values of the density index for Svoxels on a vertical slice of a black pine crown. The tree shape is well described by Svoxels with density index values apparently reasonable. However tree outline description quality is getting coarser when Svoxel size increases. Highest density values are logically located along the trunk and close to large branches. At a Svoxel density of 1m, the low density index values are located at the crown periphery. The tree back part is not as well described in the scan due to occlusion effects (shaded points). However no dissymmetry in density index can be noticed between the crown part facing the scanner system and the back part of the tree even though 3D laser point density largely differs between these two tree parts (Fig 3b,c,d). This confirms that the algorithm evaluates adequately the intercepted laser beams.

Comparative analysis of different tree crowns in a same scan

Results for the five selected trees are summarized in table 1. The number of Svoxels with a positive and significant density index decreases with the tree distance from the laser scanning system (related to the mean theoretical entering beam number). Similarly occlusion effects are more predominant with distance from the scanner depending on the obstacles in the path of the light beam. For example black pines 3 and 5 are located at a similar distance from the base station but occlusions are more numerous for black pine 5 located behind trees 1 and 2 (see figure 2). The mean density index is comparable for the five trees. But it tends to decrease with distance to the base station. Since only few trees were analyzed this could be due to natural tree heterogeneity. However, our initial results raise the possibility of a small bias. This bias could originate from partial occlusion of the incident beams. These effects, which increase with distance, are likely to increase the number of last returns with a weak intensity. These returns are not recorded by the system. The anomalies related to the distance can hence result from (1) an underestimation of lidar returns from the Svoxels and (2) a lower chance to record returns from trunks or large branches, which are characterized by high vegetation density but are partially occluded by crowns. Histograms of the density index values allow to compare the distribution for the 5 selected crowns for the three grid resolutions. Figure 4 presents the histogram for a voxel grid resolution of 50 cm. The histograms are comparable for all the pines. For the Spanish fir a slight difference can be noticed on figure 4 and was observed at all the 3 resolutions: density index frequencies are higher than those from the pines for densities ranging from 20 to 50. Consequently standard deviations were similar for all the black pine crowns and were higher for the Spanish fir (table 1). A higher foliage density for this species could explain this result. This open up the possibility to classify species using density index distribution.

Table 1: For each tree crown density index mean and standard deviation were computed for 3 grid resolutions: 25, 50 and 100 cm. The theoretical entering beam number gives an indication of the crown distance from the lidar system.

Svoxel resolution

Abies

Glauca 4s

Black

Comparison of density index for two scans

Table 2 recaps the results of the comparison of the 2 studied scans for two grid resolutions (0.5 and 1 m). The total number of Svoxels was calculated for a grid including the circular plot. After merging two scans from different locations we noticed that the no-data values represented only about 12 % of the total number of Svoxels in the plot for all grid resolutions. The voxel centers, for which a significant density index value was computed from both scans, are only about 55 % of the total number of Svoxels of the grid. This low value is explained by the fact that only the bottom part of the plot was scanned in the second scan. The significant differences in the magnitude for the "Mean density index difference" and the "Mean difference for positive and significant density index values" are explained by a high number of Svoxels located in vegetation gaps. These Svoxels, with a null index value, are consistent between scans. Large differences in density index values are observed inside the vegetation elements. For the 50 cm grid resolution about 15 % of the density index values differ from less than 1 % and 45 % from less than 5% but 20 % of the Svoxels have index values with a difference higher than 20%. Part of these differences can be explained by (1) the difference between the Svoxel shapes observed from two points of view and (2) by the type of vegetation material hit by the laser beams. For example trunks or large branches can be sources of differences since they are not seen at the same place according to the base station location (back part of them, relative to view point, is occluded). Some differences may also be related to the potential bias we previously mentioned linking scan distance with density index value. All these hypotheses will have to be verified and possible bias need to be corrected before proposing a way to merge results from various scans. Lastly, we observed from the results that mean differences decrease with resolution while standard deviation increase. This tends to confirm the influence of large wooded elements present in the Svoxel on density index value differences. Actually, when grid resolution is getting coarser the proportion of large wooded elements inside the Svoxel decreases thus reducing the mean difference.

Conclusion

We proposed an innovative method to quantify spatial distribution in 3D of forest structure from terrestrial laser scanner data. The method rests on the hypothesis that the amount of laser beam returns inside a Svoxel (volume element defined in the lidar spherical coordinate system) is proportional to the density of vegetation material included inside this Svoxel. First results confirmed that our hypotheses are valid and that some adjustments can improve further the interrelationship between the lidar returns and the amount of forest components in the Svoxels.

Visual and numerical analysis allowed verifying consistency of the derived density index values. However a slight bias could be noticed with a tendency of decrease in density index with distance from the base station position. This bias could originate from partial occlusions of the incident beams leading to a relative under sampling of trunks and large branches inside the 3D lidar point cloud when distance from the base station increases. Future work will focus on improving our algorithm, refining calculations, and correcting potential biases. In-depth analysis of scans acquired in both first and last pulse modes and multi-scan comparisons at different grid resolutions also need to be tested out. Our analysis was an essential prerequisite for developing a method aiming at merging the different scans acquired on a same plot. This study was realized considering the prospect of establishing a link between airborne lidar data and field data with the aim of improving information extraction from airborne lidar data on forested areas. These results are very promising for the development of quantitative measures of the 3D forest structure that will meet the actual information needs in the fields related to forest ecology and management.

Figure 2 :

 2 Figure 2: The tree crowns selected for analysis are shown on the 3D point cloud obtained from the system position 1 and viewed from the top. Crowns represent black pine (1, 2, 3 and 5) and Spanish fir (4).

Figure 3 :

 3 Figure 3: Density index were computed for the three grid dimension (0.25, 0.5, 1m). Density index are superimposed on their corresponding Svoxel centre on the lidar 3D point cloud. The results are given for a slice cut through a tree in the scan direction (a). Density index were separated into 4 classes using quartiles.

Figure 4 :

 4 Figure 4: Histograms of density index values (positive and significant) for 5 tree crowns and for a Svoxel resolution of 50 cm.

Table 2 :

 2 Results of comparison between density index values computed for two scans.

	Svoxel	Total number	Mean difference of density indices	Mean difference for positive and
	resolution	of Svoxel in	(SD) [Number of Svoxels]	significant density index values
		the grid				(SD) [Number of Svoxels]
	1 m	58499	-1,4	(6,4)	[432175]	-4,1	(11,5)	[698]
	50 cm	467999	-0,7	(5,4)	[208531]	-3,1	(16,1)	[2266]

SilviLaser 2008, Sept. 17-19, 2008 -Edinburgh, UK

Acknowledgements:

This research was made possible though funding provided by the Centre National d'Etudes Spatiales (CNES). We are also thankful to the GIS Draix for its logistical support during the fieldwork.


Siegfried Benkner

Jesper Sabri Pllana

Larsson Träff
email: traff@par.univie.ac.at

Philippas Tsigas

Uwe Dolinsky

Cédric Augonnet

Beverly Bachmayer

Christoph Kessler

David Moloney

Vitaly Osipov

PEPPHER: Efficient and Productive Usage of Hybrid Computing Systems

The 3-year European FP7 project PEPPHER addresses efficient utilization and usage of hybrid (heterogeneous) computer systems consisting of multi-core CPUs with GPU-type accelerators. PEPPHER is concerned with two major aspects: programmability and efficiency on given heterogeneous systems, and code and performance portability between different heterogeneous systems. The PEPPHER approach is pluralistic and parallelization agnostic, aiming to support different parallel languages and frameworks at different levels of parallelism. The central idea of PEPPHER is to maintain multiple, tailored implementation variants of performance-critical components of the application and schedule these efficiently either dynamically or statically across the available CPU and GPU resources. Implementation variants are supplied incrementally by hand, by compilation, by component composition, by auto-tuning, or taken from expert-written, adaptive libraries. This paper outlines the PEPPHER performance aware component model, its means for performance prediction, the PEPPHER run-time system, and other major aspects of the project concerned with algorithm and data structure support, compilation, and hardware feedback. A larger example demonstrating performance portability with the PEPPHER approach across hybrid systems with one to four GPUs is discussed: each GPU that is added to the system brings a linear performance increase, and performance aware scheduling provides for more efficient utilization of the combined CPU-GPU resources.

I. INTRODUCTION

With the proliferation of radically different computer architectures (many-core CPUs, embedded CPUs, SIMD instruction sets, Cell Broadband Engine Architecture, Intel MIC and SCC architectures, etc.), the fusion of different architectures into hybrid systems, as well as the rapid succession of architecture generations (e.g., NVidia GPUs) ensuring both a reasonable level of performance and a sufficient degree of functional and performance portability between different hybrid systems are vital and pressing challenges to current computer science research and engineering. Due to the large architectural parameter space it is also clear that ensuring programmability and portability for such heterogeneous systems cannot be tackled manually, but must be handled or assisted by automatic means [START_REF] Hall | Self-configuring applications for heterogeneous systems: Program composition and optimization using cognitive techniques[END_REF]. A large number of current research projects are addressing these problems, and the European FP7 project PEPPHER 1 is one of them. The PEPPHER project attacks the problems of performance portability and efficient use of heterogeneous systems at several levels at the same time. PEPPHER proposes solutions that involve a combination of static and dynamic scheduling, automatic adaptation of (library) components, compilation and transformation techniques, and a resource aware run-time that is aided by performance information with feedback monitoring for gathering empirical performance information. Specific PEPPHER contributions that will be discussed in this paper include:

• making (legacy) code written in existing programming languages and parallel APIs performance aware and portable through stepwise componentization,

• a component model with meta-data for performance awareness and adaptivity,

• algorithms and data structures for hybrid, parallel architectures, library based performance portability through algorithmic auto-tuning and compositional techniques,

• performance aware, heterogeneous run-time scheduling,

• hardware feedback mechanisms through a tailored simulator, and finally

• source-to-source transformation and OpenCL compilation.

The paper explains the vision and basic premises of PEPPHER, and highlights some concrete results to substantiate our approach. In particular, we show how a numerical kernel algorithm built from available component variants can achieve a linear performance improvement with each 1 An acronym for "PErformance Portability and Programmability for HEterogeneous many-core aRchitectures". See also www. peppher.eu extra GPU added to a hybrid CPU-GPU system, and at the same time more efficiently exploit the CPU. The example demonstrates performance portability across hybrid systems with different numbers of attached GPUs. Best overall performance is achieved, not by fixed, static offloading of tasks to the GPUs, but by dynamically scheduling of component tasks on either CPU or GPUs based on availability of resources and relative efficiency of the component variants.

II. THE PEPPHER APPROACH

The fundamental premise of PEPPHER for enabling performance portability is to provide performance-critical parts of the applications in multiple variants that are suitable for different types of cores, usage contexts and performance criteria. For this PEPPHER provides a flexible and powerful component model. Preselection and specialization of variants for a given heterogeneous architecture are performed statically as far as possible by component composition techniques, while the final selection of the most appropriate variants is delegated to a resource aware run-time system. Run-time selection is based on optimization objective, resource availability, data availability and placement, and available performance information for the variants, while respecting data dependencies between components. Component variants can be generated in part by compilation and auto-tuning mechanisms provided or enabled by the PEPPHER framework, or supplied directly by the more skilled ("expert") programmer, for instance as part of a performance portable library of algorithms and data structures. Component variants can themselves be parallel, and in that case make explicit requirements for specific, parallel resources. Preselection, composition, and run-time selection is guided by performance information that can be incrementally provided by the application developer, and that is evaluated in the context of an explicit platform model. The PEPPHER framework and methodology in this way makes it possible to gradually make an existing application more efficient for a given, heterogeneous parallel system, as well as more performance portable across different types of heterogeneous systems, namely by progressively supplying more suitable and efficient component variants, and by componentizing more and more parts of the application.

Concretely, PEPPHER introduces a flexible and extensible component model for encapsulating and annotating performance critical parts of the application. In particular, components are made performance aware by association of performance models or regressions based on performance history for predicting a desired aspect of performance (execution time, power consumption, or other). Performance aspects are parameterized and evaluated relative to an abstract platform description [START_REF] Sandrieser | Explicit platform descriptions for heterogeneous many-core architectures[END_REF]. The component model also provides for specification of resource constraints and requirements, as well as other non-functional properties of components that are nevertheless essential for the efficient execution in given contexts. The capability to maintain suitable implementation variants of components for different platforms under different circumstances and for different optimization objectives is essential. Variants can be generated automatically by compilation to different platforms and by auto-tuning techniques. For the latter, the component model makes it possible to expose tunable component parameters for which good ("optimal") values can be found by suitable auto-tuning tools. Such auto-tuning tools will, however, not be specifically developed in PEPPHER. Implementation variants can also be supplied manually by the expert programmer, e.g., targeted to different platforms. PEPPHER components may already have been parallelized using conventional parallel models and languages (OpenMP, OpenCL, Pthreads, etc.).

At run-time the components form a directed acyclic graph of component tasks. A component task variant can be scheduled when all data dependencies are resolved. Performance information, input information, optimization criteria, resource requirements and availability, data placement in the system, e.g., in main CPU or in GPU memory, are all used to determine which of the ready component task variants are scheduled on which part of the system. The execution model is parallel at multiple levels: ready component tasks can be executed in parallel on different parts of the system, and component tasks can themselves be parallel, e.g., OpenCL or CUDA variants for the GPU and multi-core parallel variants for the CPU.

Figure 1 gives an overview of the PEPPHER architecture and software stack, indicating how the elements of the approach fit together. The figure illustrates the PEPPHER approach to assist in development and generation of efficient, performance portable applications for heterogeneous many-core systems.

A. Performance Guidelines and Portability

Performance portability is an elusive notion. An application could be said to be performance portable if it executes with the same efficiency (fraction of theoretical peak performance) across different heterogeneous multi-core systems. This is a strong, absolute requirement. A different conception of performance portability is that no application restructuring be necessary when porting code from one architecture to another. In a library-based approach to performance portability this would be ensured by having implementations of the library routines targeted to different types of architectures, input configurations, and performance criteria. At each call the library would be responsible for invoking the best implementation for the given situation using an efficient lookup and selection procedure. In PEPPHER the components of the application, instead of only fixed library functionalities, are the units at which performance portability is enabled. Performance portability is supported and enforced by guidelines and requirements that application components and library routines must fulfill. An example guideline would prescribe that no library functionality can be improved just by expressing the same functionality by means of other, related or specialized library components. Such a guideline would ensure that no performance gain on a different system would be possible by the application programmer by simply reimplementing a library or component functionality in terms of other, related library functionality. For the PEPPHER framework this is an obligation to the library implementation that efficient implementation variants for different architectures are in place, and that general library functionalities select the most efficient special case implementations as appropriate. The guideline would further imply that the PEPPHER framework does the best possible selection among different available implementations of a library component. This is a non-trivial requirement in a heterogeneous and dynamic setting. By enforcing such requirements the user would be relieved from the temptation to try to do better and write the selection code himself. Other rules govern the use of component annotations: the more information is provided, the better the PEPPHER framework can do at selecting the most suited component variant. If no information is provided for the components, the PEPPHER framework produces default code. In this way PEPPHER provides an incremental approach to making applications performance portable.

B. Related work

A large number of projects are currently concerned with aspects of multi-core programmability as mentioned above. In contrast to many other European projects, e.g., HyVM, SARC, AppleCore, PEPPHER is not focusing on providing a common programming model or virtual machine type portability layer. In PEPPHER the application programmer provides performance information by annotating components and describing characteristics of the actual environment/architecture, using the most convenient API for implementation variants that are tailored to different types of CPU, GPU and other cores. Likewise PEPPHER is not concerned with automatic parallelization per se. PEPPHER is not an auto-tuning project, but enables auto-tuning techniques to be used by exposing tunable parameters of both components and parameterized, adaptive library algorithms. PEPPHER is taking a general-purpose approach in contrast to implicit parallelization and performance portability via domain specific languages, as in, e.g., [START_REF] Chafi | A domain-specific approach to heterogeneous parallelism[END_REF].

Many other projects also take the provision of implementation variants of functions, methods, or components tailored to different architectures as basic premise for addressing performance and performance portability issues. Three recent such projects are PetaBricks [START_REF] Ansel | PetaBricks: A language and compiler for algorithmic choice[END_REF], Merge [START_REF] Linderman | Merge: A programming model for heterogeneous multi-core systems[END_REF],

and Elastic computing [START_REF] Wernsing | Elastic computing: A framework for transparent, portable, and adaptive multi-core heterogeneous computing[END_REF]. PetaBricks [START_REF] Ansel | PetaBricks: A language and compiler for algorithmic choice[END_REF] is an auto-tuning project that addresses performance portability mostly across homogeneous multi-core architectures by focusing on auto-tuning methods for different types of optimization criteria. Parallelism is implicit. Merge [START_REF] Linderman | Merge: A programming model for heterogeneous multi-core systems[END_REF] also provides variants, but focuses on MapReduce [START_REF] Dean | MapReduce: Simplified data processing on large clusters[END_REF] as a unified, high-level programming model. Elastic computing [START_REF] Wernsing | Elastic computing: A framework for transparent, portable, and adaptive multi-core heterogeneous computing[END_REF] focuses on provision of large number of variants, so called elastic functions, among which the best combination is composed mostly by static means. Selection is guided by performance profiles and models. Use of algorithmic variants has previously been studied extensively in the STAPL project [START_REF] Thomas | A framework for adaptive algorithm selection in STAPL[END_REF]. Other work that has influenced PEPPHER is the Sequoia programming model for heterogeneous architectures based on a tree abstraction of the memory system [START_REF] Fatahalian | Sequoia: Programming the memory hierarchy[END_REF].

In contrast to some of these works PEPPHER is distinguished by a holistic approach, which attacks performance portability at multiple layers from high-level component based programming, compilation, library and run-time support, to hardware mechanisms for performance monitoring and feedback. Ultimately, PEPPHER aims to become language and parallelization agnostic and support different implementation languages and parallelization interfaces. Thus, the programmer can choose the most convenient language and API for implementing the components.

III. TECHNICAL ASPECTS OF PEPPHER

In this section we discuss in more detail the technical ingredients of the PEPPHER framework.

A. The PEPPHER component model

As a starting point PEPPHER has developed a comprehensive and flexible component model for specification of performance aware components, implementation variants, and tunable parameters of components. Applications for PEPPHER are currently assumed to be written in C++. PEPPHER components are identified by language external means (e.g., pragmas) rather than language extensions and employ XML schemata for their specification. In this sense the PEPPHER framework is non-intrusive, requiring little direct application modification, for instance in legacy source code.

A PEPPHER component is an annotated software module that implements a specific functionality declared in a PEPPHER interface. A PEPPHER interface is defined by an interface descriptor, an XML document that specifies the name, parameter types and the access types (read, write or both) of a function to be implemented, and in addition which performance aspects about the following properties:

• The implemented PEPPHER interface.

• The required PEPPHER interfaces, i.e., other component functionality called from this component, if any.

• The source file(s) of this component implementation.

• Deployment information such as compilation commands and options and required software modules.

• A reference to the platform consisting of the programming model/language used for the component implementation and of the target architecture.

• Type and amount of resources required for execution at the given platform.

• A reference to a performance prediction function. whose parameters are given by a context descriptor data structure.

• The structure of the context descriptor consisting of call parameter context and current resource availability information.

• Tunable parameters of the component implementation, such as buffer sizes, loop blocking factors, cut-off values and so on.

• Additional constraints for component selectability such as parameter ranges.

Performance prediction functions are usually supplied by the component developer. These can be either purely analytical, use performance data tables determined by micro-benchmarking for the target platform, or be based on historic performance data. The latter are maintained in a performance data repository. The PEPPHER run-time is able to automatically generate performance models from historic data.

The actual platform properties are defined separately in another XML document and aim at describing hardware-and software characteristics of the heterogeneous target environment. Such platform descriptions can be used at multiple levels of the PEPPHER framework [START_REF] Sandrieser | Explicit platform descriptions for heterogeneous many-core architectures[END_REF]. Lookup of specific platform properties is done by the composition tool, the run-time, or by component developers themselves.

At run-time component invocations result in tasks that are managed by the PEPPHER runtime system and executed non-preemptively. PEPPHER components and tasks are stateless, and by default so are input/output parameters which determine the dependencies among component tasks. Parameter data placement in CPU or GPU memory is tracked by the run-time system and taken into account during variant scheduling. However, parameters can have state themselves, implemented by wrapping into STL-like container data structures, and keep track of, e.g., in which memory modules of the target system which parts of the data are currently located or mirrored. The container state then becomes part of the call context information since it is relevant for performance prediction.

The PEPPHER framework keeps track of the different component implementation variants by storing descriptors in repositories that can be explored by the composition tool and the run-time.

Composition performs preselection of a subset of appropriate implementation variants based on statically available parameter and resource information. Parameter information can be either concrete or abstract values given by constraints. Composition can in the extreme case prune the number of available variants to only one, but in general several component variants will be available for selection by the run-time system. Part of the variant selection code may also be generated by the composition tool. Executable component variants are finally generated by compilation for the platforms described in the component variant descriptions.

The final selection among the remaining component variants is done by the run-time system (see Section III-B). This is the default mechanism in PEPPHER. In the special case where sufficient meta-data for performance prediction is available for all selectable component variants, composition can be done completely statically and co-optimized with the required resource allocation [START_REF] Kessler | A framework for performance-aware composition of explicitly parallel components[END_REF].

To enhance programmability and provide for more variant selection and auto-tuning possibilities, portable higher-level coordination strategies for expressing computations in terms of structured execution of PEPPHER components are being added to the framework, e.g., pipeline pattern, wave-front pattern, task farming, standard skeletons, and others. First steps have been taken in SkePU [START_REF] Dastgeer | Auto-tuning SkePU: A multi-backend skeleton programming framework for multi-GPU systems[END_REF], an auto-tunable C++ template library of data-parallel generic skeleton components such as map, reduce, or scan, each with multiple implementation variants including CUDA, OpenCL, and OpenMP, for multi-GPU based systems. By micro-benchmarking SkePU can be automatically tuned to select, depending on the call context, the expected fastest back-end and values for tunable parameters such as number of GPU threads and thread block size.

B. The PEPPHER run-time system

Execution of a PEPPHER application that consists of compiled component variants together with the parts of the application that have not been componentized is delegated to a flexible, performance and resource aware, heterogeneous run-time system. The executable forms a directed acyclic graph of component tasks with data dependencies, where each task is a set of one or more variants together with performance and other information to trigger selection. Based on system availability, resource requirements, estimated performance, execution history, and input availability the run-time system schedules the most promising component variant on the best available resource. The further development of such a performance, resource and architecture/memory aware run-time system, that can also handle scheduling of parallelized component variants over different parts of the heterogeneous system is an important aspect of the project [START_REF] Augonnet | StarPU: A unified platform for task scheduling on heterogeneous multicore architectures[END_REF].

For CPU-GPU based systems with separate memory spaces, the run-time system implements a directory based virtual shared memory system. The application registers its input and output data with the system. Scheduling of the ready component tasks is done in a centralized fashion from a CPU. At run-time the actual placement of data in either main CPU memory or GPU device memory is used together with the component performance models to decide where to launch the next component variant that uses these data. To this end, a data transfer cost model is used together with the cost estimation for the component execution. Actual data transfers between different types of memory are handled automatically by the run-time. Through so called filters structured non-consecutive as well as block distributed data can be handled.

The actual scheduling policy used by the run-time system can be modified by the PEPPHER user. The current default strategy is HEFT [START_REF] Topcuoglu | Performance-effective and low-complexity task scheduling for heterogeneous computing[END_REF], but alternate, less centralized scheduling policies that would use more hints on processor resource requirements in order to execute parallel tasks on the multi-core CPU are being developed [START_REF] Träff | Work-stealing for mixed-mode parallelism by deterministic team-building[END_REF].

The multi-level parallel framework providing for concurrency between component variants and intra-component-variant parallelism, the static, resource and architecture aware compositional techniques, and the dynamic, flexible run-time scheduling together enable the PEPPHER framework to both efficiently utilize given, heterogeneous resources, as well as to provide performance portability (by recompilation, recomposition, re-tuning relative to a different platform description)

to entirely different architectures.

C. Tunable algorithms and data structures for parallel architectures

The static compositional and dynamic run-time supported approach to performance portability is complemented by expert-written auto-tuned, architecture and context adaptive algorithmic components for further enhancing performance portability. Algorithms written by specialists with a detailed understanding of the underlying architecture and its range of tunable architecturedependent parameters can more flexibly and efficiently adapt to architectural changes, and provide for more detailed control over performance. Through libraries such highly performance portable algorithms are made available as components to the application programmer. The highly nontrivial, auto-tunable GPU sorting algorithm developed in [START_REF] Leischner | GPU sample sort[END_REF] is an example of the level of adaptable, portable performance that can be achieved by the algorithm engineering expert.

Some of the basic, performance determining parameters of the NVidia GPU Tesla and Fermi architectures are summarized in Table I. Based on these, tunable, algorithmic parameters related to these architectural features are inferred by the algorithm developer as shown in Table II. These parameters are partially interrelated through an analytical performance model of the algorithm, and must partly be determined experimentally or through auto-tuning. By finding the right settings for the free parameters, superior performance of the algorithm is achieved on both architectures as shown in Figure 3. The figure compares a previous (original) [START_REF] Leischner | GPU sample sort[END_REF] to a current, improved version version of the algorithm. The "best possible" values of the tunable parameters are taken from Table II, and the impact of changes to the four basic parameters k, t, p, M can be seen in Figure 4. Some of these parameters are determined analytically, others have been found by auto-tuning experiments. In PEPPHER concrete values for architectural parameters are looked up in the PEPPHER platform descriptors. At application deployment time Tesla and Fermi sample sort variants can thus be generated, partly through auto-tuning for determining best values for the non-determinate tuning parameters, and added to the PEPPHER component repository. The algorithm has furthermore been implemented in both CUDA and OpenCL, but only the CUDA results are shown here. We note that the CUDA and OpenCL interfaces by themselves only ensure code portability, but do not provide performance portability: this is achieved by the parameterized algorithm.

A larger selection of similarly tunable algorithms for many-core CPU architectures is incorporated in the Multi-Core STL library (MCSTL) [START_REF] Singler | MCSTL: The multi-core standard template library[END_REF] and some contributed to the Thrust GPU library. Also here, algorithms are parameterized with typical architectural features, and can thus be tuned to achieve a high degree of performance portability for the class of target architectures. Not mentioned here, algorithms and data structures for lock-free programming on CPUs and GPUs are needed and being developed in PEPPHER, both as application programmer components and for supporting the implementation of the run-time system, see, e.g., [START_REF] Gidenstam | Cache-aware lock-free queues for multiple producers/consumers and weak memory consistency[END_REF].

Choice between CPU and GPU variants for such library components is done either by specialized glue-code that takes the resource availability into account, that is by the library framework itself, or is delegated to the PEPPHER run-time. Thus, variant selection in PEPPHER is always external to the implementation of the variant.

D. Compilation techniques in PEPPHER

OpenCL is a possible portability layer for current heterogeneous systems, especially such that employ GPUs, but of a low abstraction level and not in itself solving the performance portability problem. Therefore, efficient compilation from C++ to OpenCL can support both the component model and the run-time system. Specifically, PEPPHER develops a C++ extension termed OffloadCL that allows for explicit compilation and offloading to GPUs but also to Cell SPEs with the compiler taking care of the necessary call graph duplication, functional duplication, and replication of host data. The Offload compiler interfaces directly with the PEPPHER run-time system by encapsulating offloadable code as PEPPHER component tasks [START_REF] Cooper | Offload -automating code migration to heterogeneous multicore systems[END_REF].

E. Hardware support and feedback

In order to provide a larger spectrum of heterogeneous target architectures and to be able to investigate possibilities for hardware support for performance monitoring and portability, PEPPHER develops the highly configurable hardware simulator PeppherSim. This simulator makes it possible to run benchmarks/kernels on architecture configurations that are not physically available, and experiment in a controlled way with the performance portability that can be achieved through the PEPPHER framework. The simulator enables "what if" types of experiments to determine performance bottlenecks, as well as to investigate new mechanisms for more detailed performance monitoring and feedback, especially energy consumption feedback. Also, new synchronization primitives, e.g., [START_REF] Ha | NB-FEB: A universal scalable easy-to-use synchronization primitive for manycore architectures[END_REF], and other architecture support for algorithms and run-time can be investigated with the simulator.

F. Application benchmarks

A small set of larger application benchmarks has been selected to experiment with and validate the performance portability that can be achieved with the PEPPHER framework. The benchmarks are intended to cover relevant application areas ranging from embedded, server/enterprise/generalpurpose, to high-performance computing domains. In addition, important numerical kernels have been included as manageable test cases, that will also be useful as library components. Table III summarizes the benchmarks and kernels and the target architectures for which component variants already exist.

IV. A LARGER EXAMPLE

A final example illustrates how PEPPHER can provide efficient utilization of heterogeneous compute resources and performance portability across systems with varying numbers of GPU type accelerators.

The application is a standard Tile-QR factorization algorithm based on BLAS components that has been implemented on top of the PEPPHER run-time [START_REF] Agullo | QR Factorization on a Multicore Node Enhanced with Multiple GPU Accelerators[END_REF]. The application is constructed from BLAS kernel functions as shown in Figure 5, and these kernels are turned into PEPPHER components by giving suitable XML interface specifications. Expert implementation variants of the BLAS kernels for CPUs and GPUs are already available through existing GPU and CPU libraries, namely PLASMA [START_REF] Buttari | A class of parallel tiled linear algebra algorithms for multicore architectures[END_REF] and MAGMA [START_REF] Tomov | Dense linear algebra solvers for multicore with gpu accelerators[END_REF], and are added as component variants Fig. 6. QR factorization on a 16-core AMD system with up to 4 NVidia GPUs (C1060).

in the PEPPHER implementation. Figure 6 shows the results of running the QR factorization application on CPU-GPU system with one to four GPUs. As can be seen, the SOMQR and SSSMQR kernels are particularly efficient on the GPU, and the GPU variants are therefore likely to be selected by the PEPPHER run-time.

The plot shows that the PEPPHER run-time is able to exploit the GPUs efficiently in that each additional GPU brings a linear increase in performance of about 200GFlops. In addition, the run-time is able to exploit the remaining 12 CPUs for execution of the most suitable BLAS components. The latter is interesting; by automatically avoiding to waste GPU processing power with GPU-unfriendly BLAS kernels (that are scheduled on CPUs instead), more performance can be achieved from the 12 CPUs (about 200 GFlops) than would have been possible by running the whole QR application on these 12 CPUs (only 150 GFlops) alone. In the experiment 20% of the SGQRT tasks are scheduled on the GPUs, whereas for the SSSMQR tasks over 90% are executed on the GPUs, see again the speed-up ratios in Table IV.

The experiment demonstrates efficiency and performance portability of this particular application across systems with one to four GPUs. It is important to note that it relied on existing expert written components, so no extra implementation effort was needed for the kernels.

V. CONCLUSION

This paper outlined the PEPPHER approach to achieving performance portability and programmability for hybrid (heterogeneous) many-core architectures, particularly CPU-GPU type systems, and described the current state the project. PEPPHER combines a compositional, performance aware software framework, auto-tunable, adaptive algorithmic libraries, specific compilation techniques and an efficient run-time, and is thereby independent of specific programming models, virtual machines and architectures. The combination of component based adaptation with a performance and resource aware run-time system provides for the necessary degrees of freedom, possibly not found in many auto-tuning projects, where scheduling decisions may be hard-coded at an early stage. We contend that neither best performance nor performance portability can be achieved by fixed, static offloading of promising tasks of the component based application onto preselected cores. Instead, component tasks must be scheduled dynamically on available cores for which a variant giving the best performance exists, and we aim to substantiate this with the PEPPHER project.

Fig. 1 .

 1 Fig. 1. The PEPPHER architecture and software stack. Applications written in C/C++ are annotated by the application programmer to enable generation of component variants. The more skilled, "expert" programmer or library writer can provide additional, highly performant component variant implementations for specific hardware, optimization criteria, or input configurations. Static composition performs component pruning and preselection and generates necessary glue code, while the final selection is done dynamically by the performance aware PEPPHER run-time system.

(

  such as mean execution time or energy consumption) the prediction functions of component implementations must provide. The component model is illustrated in Figure2. Interfaces can be generic in static entities such as parameter types or code; genericity is resolved statically by expansion as with C++ templates. From the interface descriptor, header files in the various component implementations' source languages can be generated.Several component variants may implement the functionality defined in a PEPPHER interface by different algorithms or for different execution platforms; also, further component implementation variants may be generated automatically from a common source module, e.g., by special compiler transformations or by instantiation of tunable parameters. The latter would be the task of auto-tuning tools. These variants differ by their resource requirements and performance behavior, and thereby become alternative choices for composition or run-time selection whenever the interface function is called. In order to prepare and guide variant selection, component implementations need to expose their relevant properties explicitly to composition tool and runtime as will be briefly described later. Each PEPPHER component implementation variant thus provides its own component descriptor, an XML document that contains information (meta-data)

Fig. 3 .Fig. 4 .

 34 Fig. 3. Comparison of original and improved GPU sorting algorithm on Tesla and Fermi GPUs, number of sorted elements per time unit as function of number of elements to be sorted.

PEPPHER Run-time Drivers (CUDA, OpenCL, OpenMP) CPU& GPU APU Peppher Sim PEPPHER Task graph Transformation & Composition PePU

		Applications		C/C++ source code with annotated
	Embedded	General Purpose	HPC	components
		Components
	C/C++, OpenMP, CUDA, OpenCL,		Component implementation.
					Variants for different cores,
					architectures, algorithms,..
					Component glue code generation.
					Static variant selection, composition
				Performance Models	Component task graph with explicit data dependencies
					Performance-aware, data-aware
	Scheduling Strategy 1	Scheduling Strategy N	dynamic scheduling of "best" component variants onto available
					execution units
					Single-node heterogeneous manycore
					systems

Offload Autotuned Algorithms Data Structures

  Fig. 2. A PEPPHER component consists of an abstract interface describing the functionality of the component, and a number of implementation variants targeted for different architectures and with different performance characteristics. Interface and variant meta-data descriptions are described in XML documents.

		Abstract Interface -Functional Definition
	Component Variant	Component Interface	Component Variant
	Variant Meta-Data	Interface Meta-Data	Variant Meta-Data
		abstract function interface
	source codes references	supported performance aspects	source codes references
	compilation information	...	compilation information
	execution target platform		execution target platform
	performance aspect models, prediction functions		performance aspect models, prediction functions
	deployment information		resource requirements
	tunable parameters		...
	Concrete Implementation for Interface		Concrete Implementation for Interface

TABLE II THE

 II TUNING PARAMETERS FOR THE TWO ALGORITHM VARIANTS OF THE GPU SAMPLE SORT ON NVIDIA TESLA

AND

FERMI GPUS.

TABLE IV KERNEL

 IV COMPONENTS OF THE QR APPLICATION. THE TABLE SHOWS THE PERFORMANCE OF THE KERNELS ON CPU AND GPU, RESPECTIVELY.

		1100	1GPU
		1000	2GPUs 3GPUs
			4GPUs
		900	All
		800
	Performance [GFLOPS]	400 500 600 700
		300
		200
		100
		0
			0	5000	10000	15000	20000	25000	30000	35000	40000
							Matrix Size

  The test platform is composed of four quad-core AMD Opteron 8358 SE CPU cores (16 cores total) running at 2.4 GHz with 32 GB

	of memory divided into four NUMA nodes. It is accelerated with four NVIDIA Tesla C1060
	GPUs of 240 cores each (960 GPU cores total) running at 1.3 GHz with 4 GB of GDDR3
	memory (102 GB/s) per GPU. One CPU is reserved for each GPU, leaving 12 CPUs available

for computational work. The performance of the kernels of the QR application is shown in Table

IV

, which also shows the speed-up achievable on the GPU relative to a single CPU core.

  Jesper Larsson Träff received an M.Sc. in computer science in 1989, and, after two years at the industrial research center ECRC in Munich, a Ph.D. in 1995, both from the University of Copenhagen. He spent four years as a Research Associate in the Algorithms Group of the Max-Planck Institute for Computer Science in Saarbrücken, and the Efficient Algorithms Group at the Technical University of Munich. From 1998 until late 2009 he was working at the NEC Laboratories Europe in Sankt Augustin, Germany on efficient implementations of MPI for NEC vector supercomputers. This work led to a doctorate (Dr. Scient.) from the University of Copenhagen in 2009. Since 2010 he is Professor for Scientific Computing at the University of Vienna. His research interests are broadly in parallel processing and include interfaces, algorithms, and architectures. He is currently scientific coordinator for the European FP7 project PEPPHER. With Martti Forsell he organizes the annual Euro-Par Workshop on Highly Parallel Processing on a Chip (HPPC).

	Nordbergstrasse 15/C3
	1090 Vienna
	Austria
	Phone: +43 1 4277 39411
	Fax: +43 1 4277 9394
	Email: pllana@par.univie.ac.at
	Contact Information:
	University of Vienna
	Faculty of Computer Science
	Research Group of Scientific Computing
	Nordbergstrasse 15/3C
	1090 Vienna
	Austria
	Phone: 0043 1 4277 39432
	Fax: 43 1 4277 9394

Email: traff@par.univie.ac.at

Acknowledgment

This work was supported by the EU as part of the FP7 Project PEPPHER under grant 248481.

We gratefully acknowledge the contributions of all people involved in the PEPPHER project to the work presented in this paper. The authors also thank the anonymous reviewers for constructive feedback.


J L Fouquet

H Thuillier

J M Vanherpe

A P Wojda

ON MINIMUM (K q , k) STABLE GRAPHS

Keywords:

A graph G is a (Kq, k) vertex stable graph (q ≥ 3) if it contains a Kq after deleting any subset of k vertices (k ≥ 0). We are interested by the (Kq, κ(q)) vertex stable graphs of minimum size where κ(q) is the maximum value for which for every nonnegative integer k < κ(q) the only (Kq, k) vertex stable graph of minimum size is K q+k .

Introduction

In [START_REF] Katona | Extremal stable graphs[END_REF] Horwárth and Katona consider the notion of (H, k) stable graph: given a simple graph H, an integer k and a graph G containing H as subgraph, G is a a (H, k) stable graph whenever the deletion of any set of k edges does not lead to a H-free graph. These authors consider (P n , k) stable graphs and prove a conjecture stated in [START_REF] Frankl | Extremal k-edge hamiltonian hypergraphs[END_REF] on the minimum size of a (P 4 , k) stable graph. In [START_REF] Dudek | stable graphs with minimum size(H, k)[END_REF], Dudek, Szymański and Zwonek are interested in a vertex version of this notion and introduce the (H, k) vertex stable graphs. In [START_REF] Fouquet | On (Kq, k) vertex stable graphs with minimum size[END_REF] we have characterized (K q , k) vertex stable graphs with minimum size for q = 3, 4, 5 (where K q denotes the clique on q vertices). Definition 1.1. [START_REF] Dudek | stable graphs with minimum size(H, k)[END_REF] Given an integer k ≥ 0 and a graph H, a graph G containing a subgraph isomorphic to H is said to be a (H, k) vertex stable graph if, for every subset S of k vertices, G -S contains (a subgraph isomorphic to) H.

In this paper, we are only interested by (H, k) vertex stable graphs and, since no confusion will be possible, a (H, k) vertex stable graph shall be simply called a (H, k) stable graph. Definition 1.2. A (H, k) stable graph with minimum size (i.e. with minimum number of edges) is called minimum (H, k) stable graph. The minimum size of a (H, k) stable graph shall be denoted by stab (H, k).

It is clear that if G is a (H, k) stable graph with minimum size then the graph obtained from G by addition or deletion of some isolated vertices is also minimum (H, k) stable. Hence we shall asume that all the graphs considered in the paper have no isolated vertices. Here we consider (K q , k) stable graphs. We have proved in [START_REF] Fouquet | On (Kq, k) vertex stable graphs with minimum size[END_REF]

that stab(K q , k) = q+k 2
for q ≥ 3 and k = 1, 2. Moreover K q is the only minimum (K q , 0) stable graph for q ≥ 2 and K q+1 is the unique (K q , 1) stable graph for q ≥ 4. Dudek, Szymański and Zwonek have proved the following result. Theorem 1.3. [START_REF] Dudek | stable graphs with minimum size(H, k)[END_REF] For every q ≥ 4, there exists an integer k(q) such that stab

(K q , k) ≤ (2q -3)(k + 1) for k ≥ k(q).
As a consequence of this last result, for every k ≥ k(q) the graph K q+k is not minimum (K q , k) stable. Definition 1.4. For every integer q ≥ 4, we denote by κ(q) the greatest integer such that for 1 ≤ k < κ(q) the only minimum (K q , k) stable graph is K q+k .

In two previous papers we have proved the following. Theorem 1.5. [START_REF] Fouquet | On (Kq, k) stable graphs with small k[END_REF][START_REF] Fouquet | On (Kq, k) vertex stable graphs with minimum size[END_REF] κ(3) = 1, κ(4) = 3, κ(5) = 4 and for q ≥ 6 κ(q) > q 2 + 1. In this paper we give an upper bound for the value of κ(q). Theorem 1.6. For every q ≥ 4, if κ(q) is even then κ(q) < 2(q -1)(q -2) and if κ(q) is odd then κ(q) < 1 + 2(q -1)(q -2)

We prove that these upper bounds are reached for values of q such that there exists a minimum (K q , κ(q)) stable disconnected graph (note that it is the case for q = 4 and q = 5).

Theorem 1.7. Let q ≥ 4 and suppose that there exists a disconnected minimum

(K q , κ(q)) stable graph. Set ρ(q) = 1 2 (q -1)(q -2) -1. If 1 2 (q -1)(q -2) > ρ(q) 2 + ρ(q) then κ(q) = 2ρ(q) + 1. If 1
2 (q -1)(q -2) ≤ ρ(q) 2 + ρ(q) then κ(q) = 2ρ(q). Proofs of Theorems 1.6 and 1.7 shall be given in subsection 3.2.

If there is no minimum disconnected (K q , κ(q)) stable graph then, by definition of κ(q), there exists a connected minimum (K q , κ(q)) stable graph G q which is not complete. We think that it never happens, so we propose the following conjecture.

Conjecture 1.8. If G is a minimum (K q , k) stable graph then every component of G is complete.

If this last conjecture is true then Theorem 1.7 gives the exact value of κ(q) for every q ≥ 4.

Notations and general results

For terms not defined here we refer to [START_REF] Bondy | Graph Theory[END_REF]. As usually, the order of a graph G is the number of its vertices and the size of G is the number of its edges (it is denoted by e(G)). The disjoint union of two graphs G 1 and G 2 is denoted by G 1 + G 2 . The union of p mutually disjoint copies of a graph G is denoted by pG. For any set A of vertices, we denote by G[A] the subgraph induced by A and by G -A the subgraph induced by V (G) -A. If A = {v} we write G -v for G -{v}. For any set F of edges, we denote by G -F the spanning subgraph (V (G), E(G) -F ). If F = {e} we write G -e instead of G -{e}. A complete subgraph of order q of G is called a q-clique of G. The complete graph of order q is denoted by K q . When a graph G contains a q-clique as subgraph, we say "G contains a K q ".

Lemma 2.1. [2] Let G be a (H, k) stable graph with k ≥ 1. Then, for every vertex v, G -v is (H, k -1) stable.
A set of vertices of G that intersects every subgraph of G isomorphic to H is called a transversal of all the subgraphs isomorphic to H or simply a H-transversal of G. A H-transversal of G having the minimum number of vertices is said to be a minimum H-transversal of G. The number of vertices of a minimum H-transversal

is denoted by τ H (G). Remark that G is (H, k) stable if and only if τ H (G) ≥ k + 1 Definition 2.2. Let G be a (H, k) stable graph. If G has a minimum H-transversal having exactly k + 1 vertices, G is said to be exactly (H, k) stable. Remark 2.3. H is the unique (H, 0) stable graph with minimum size. Lemma 2.4. [2] Let G be a (H, k) stable graph with k ≥ 1 and e ∈ E(G) such that G -e is not (H, k) stable. Then G is exactly (H, k) stable and G -e is exactly (H, k -1) stable. Definition 2.5. [2] Let G be a (H, k) stable graph. If G -e is not (H, k) stable for every edge e ∈ E(G), G is said to be minimal (H, k) stable.
Remark 2.6. In [START_REF] Dudek | stable graphs with minimum size(H, k)[END_REF] Theorem 2.9. Let G be an exactly (H, k) stable graph, and let G 1 , G 2 , ..., G r , with r ≥ 1, be its components. Then, there exist integers

k 1 , k 2 , ..., k r , with 0 ≤ k i ≤ k, sucht that i) for every i, with 1 ≤ i ≤ r, G i is exactly (H, k i ) stable, ii) r i=1 k i + (r -1) = k G is minimal (H, k) stable if and only if for every i, 1 ≤ i ≤ r, G i is minimal (H, k i ) stable. Moreover, if G is minimum (H, k) stable then for every i, 1 ≤ i ≤ r, G i is minimum (H, k i ) stable. Proof For each i, 1 ≤ i ≤ r,
∈ {1, • • • r} such that |S i0 | ≤ k i0 = |T i0 | -1. Then, G i0 -S i0
contains a subgraph isomorphic to H, that is, G is exactly (H, |T | -1) stable, and we have

r i=1 k i + (r -1) = k.
Let e be an edge of G and let G i be the component containing e.

Claim . G -e is (H, k) stable if and only if

G i -e is (H, k i ) stable. Proof Suppose that G i -e is (H, k i ) stable. Let U be a H-transversal of G -e. Set U i = U ∩ V (G i -e) = U ∩ V (G i ) and for every j = i, U j = U ∩ V (G j ). Since (G i -e) -U i and each G j -U j , j = i, contain no subgraphs of G -e isomorphic to H, we have for every j, 1 ≤ j ≤ r, |U j | ≥ k j + 1. Then, |U | = r j=1 |U j | ≥ k + 1. Hence, for every set S of k vertices (G -e) -S contains a subgraph isomorphic to H, that is, G -e is (H, k) stable. Conversely, suppose that G i -e is not (H, k i ) stable. Let T i be a H-transversal of (G i -e) -T i having k i vertices. For every j = i let T j be a H-transversal of G j having k j + 1 vertices. The set T = ∪ r j=1 T j has k vertices and is a H-transversal of G -e, that is, G -e is not (H, k) stable.
Thus, G is minimal (H, k) stable if and only if for every i,

1 ≤ i ≤ r, G i is minimal (H, k i ) stable.
Note that, by replacing a minimal (H, k i ) stable component G i by any minimal (H, k i ) stable graph G i (connected or not), we obtain again a minimal (H, k) stable graph. Thus, if G is minimum (H, k) stable then for every i,

1 ≤ i ≤ r, G i is minimum (H, k i ) stable. Remark 2.10. Let r be an integer ≥ 2, k 1 , • • • , k r be r non negative integers and k = r i=1 k i + (r -1). If for every i, 1 ≤ i ≤ r, G i is a minimum (H, k i ) stable graph then the disjoint union G 1 + G 2 + • • • + G r
may not be a minimum (H, k) stable graph. For example, K q is minimum (K q , 0) stable, 2K q and K q+1 are minimal (K q , 1) stable, but since e(2K q ) > e(K q+1 ), for q ≥ 4 the graph 2K q is not minimum (K q , 1) stable.

Minimum (K q , k) stable graphs

In this section we are interested by (K q , k) stable graphs with minimum size (q ≥ 3). Recall that stab(K q , k) = M in{e(G) | G is (K q , k) stable}.

3.1. Some known results. We give here some known results about this topic. By Remark 2.6 and Lemma 2.7 we have: Properties 3.1. [START_REF] Dudek | stable graphs with minimum size(H, k)[END_REF] A minimal (K q , k) stable graphs G has the following properties:

P 1 ) G is exactly (K q , k) stable. P 2 ) For every edge e, G -e is exactly (K q , k -1) stable. P 3 ) For every vertex v, G -v is exactly (K q , k -1) stable. P 4 ) Every vertex of G belongs to some q-clique of G. P 5 ) Every edge of G belongs to some q-clique of G.

Remark 3.2. For any two integers q ≥ 3 and k ≥ 1, K q+k is minimal (K q , k) stable.

Proposition 3.3. [4]

For every integer q ≥ 4, K q+1 is the unique minimum (K q , 1) stable graph.

Proposition 3.4. [START_REF] Fouquet | On (Kq, k) vertex stable graphs with minimum size[END_REF] For every integer q ≥ 4, K q+2 is the unique minimum (K q , 2) stable graph.

Proposition 3.5. [START_REF] Fouquet | On (Kq, k) vertex stable graphs with minimum size[END_REF] For every integer q ≥ 5, K q+3 is the unique minimum (K q , 3) stable graph.

Theorem 3.6. [START_REF] Dudek | stable graphs with minimum size(H, k)[END_REF] For every k ≥ 1, stab(K 3 , k) = 3k + 3 and stab(K 4 , k) = 5k + 5.

Theorem 3.7. [START_REF] Fouquet | On (Kq, k) vertex stable graphs with minimum size[END_REF] Let G be a minimum (K 3 , k) stable graph, with k ≥ 0. Then G is isomorphic to sK 4 + tK 3 , for any choice of s and t such that 2s + t = k + 1.

Theorem 3.8. [START_REF] Fouquet | On (Kq, k) vertex stable graphs with minimum size[END_REF] Let G be a minimum (K 4 , k) stable graph, with k ≥ 1. Then G is isomorphic to sK 6 + tK 5 , for any choice of s and t such that 3s + 2t = k + 1 Theorem 3.9. [START_REF] Fouquet | On (Kq, k) vertex stable graphs with minimum size[END_REF] For every k ≥ 5, stab(K 5 , k) = 7k + 7.

Theorem 3.10. [START_REF] Fouquet | On (Kq, k) vertex stable graphs with minimum size[END_REF] Let G be a minimum (K with k ≥ 0, q ≥ 3, 1 ≤ r ≤ k + 1 as the family of graphs consisting of r complete graphs

K ij with i 1 ≥ • • • ≥ i r ≥ q satisfying the condition r i=1 (i j -q) + (r -1) = k and they prove that every graph in A (Kq,k) r is minimal (K q , k) stable. We observe that if G is a (K q , k) stable graph disjoint union of r ≥ 1 cliques K ij , 1 ≤ j ≤ r, then by Theorem 2.9, G ∈ A (Kq,k) r . They defined a graph G ∈ A (Kq,k) r
as a balanced union if |i j -i l | ∈ {0, 1} for every j and l in {1, 2, • • • , r} and they proved that given q, k and r there is exactly one balanced union

B (Kq,k) r in A (Kq,k) r
, and that B (Kq,k) r has the minimum number of edges among the graphs in A (Kq,k) r . In [START_REF] Dudek | stable graphs with minimum size(H, k)[END_REF] the following lemma has been given. We give its proof for completeness. Lemma 3.11. [START_REF] Dudek | stable graphs with minimum size(H, k)[END_REF] Let G 0 be a (K q , k 0 ) stable graph (k 0 ≥ 0) which has the minimum size among all graphs beeing a disjoint union of r cliques (r ≥ 1) G i ≡ K q+kj with 1 ≤ j ≤ r, k j ≥ 0. There exists a positive integer s and a nonnegative integer

k such that k ≥ 0, s ≤ r, G 0 = sK q+k + (r -s)K q+k-1 with rk + s = k 0 + 1 and e(G 0 ) = (r(q -1) + k 0 + 1 -s)(r(q -2) + k 0 + 1 + s)) 2r .
Proof Suppose, without loss of generality, that k 1 ≥ k 2 ≥ • • • ≥ k r and that there exist two components G i and G j with i < j such that k i -k j ≥ 2. By substituting

G i ≡ K q+ki-1 for G i and G j ≡ K q+kj +1 for G j , we obtain a new (K q , k) stable graph G 0 such that e(G 0 ) = e(G 0 ) -(k i -k j -1) < e(G 0
), which is a contradiction. Thus, for any i and any j, 0 ≤ |k i -k j | ≤ 1. Hence, either for any i and any j k i and k j have the same value k and we have G 0 = rK q+k with k ≥ 0, or there exist distinct k i and k j and we have

G 0 = sK q+k + (r -s)K q+k-1 with k ≥ 1 and 1 ≤ s ≤ r -1. If G 0 = sK q+k + (r -s)K q+k-1 then a minimum K q -transversal of G 0 has k 0 + 1 = s(k + 1) + (r -s)k = s + rk vertices. Note that r divides k 0 + 1 -s. We have 2e(G 0 ) = s(q + k)(q + k -1) + (r -s)(q + k -1)(q + k -2). Since k = k0+1-s r ,
we obtain e(G 0 ) = (r(q-1)+k0+1-s)(r(q-2)+k0+1+s)) 2r .

3.2. Minimum (K q , k) stable graph for small k. Recall that for every integer q ≥ 2, K q is the unique minimum (K q , 0) stable graph. Remark 3.12. By Propositions 3.3 and 3.4, K 5 is the unique minimum (K 4 , 1) stable graph, K 6 is the unique minimum (K 4 , 2) stable graph and, by Theorem 3.8, 2K 5 is the unique minimum (K 4 , 3) stable graph. Remark 3.13. By Propositions 3.3, 3.4 and 3.5, for every k ∈ {1, 2, 3} the graph K 5+k is the unique minimum (K 5 , k) stable graph. Lemma 3.14. [START_REF] Fouquet | On (Kq, k) stable graphs with small k[END_REF] K 9 and K 6 + K 7 are the only minimum (K 5 , 4) stable graphs. Theorem 3.15. [START_REF] Fouquet | On (Kq, k) stable graphs with small k[END_REF] Let G be a minimum (K q , k) stable graph, where q ≥ 6 and k ≤ q 2 + 1. Then G is isomorphic to K q+k . Recall that for every integer q ≥ 4, κ(q) is the greatest integer such that for 1 ≤ k < κ(q) the only minimum (K q , k) stable graph is K q+k . Then, either e(K q+κ(q) ) > stab(K q , κ(q)) or e(K q+κ(q) ) = stab(K q , κ(q)) but there is a minimum (K q , κ(q)) stable graph G such that K k+κ(q) ≡ G. Lemma 3.16. κ(4) = 3 and κ(5) = 4.

Proof By Remark 3.12, κ(4) = 3 (and 2K 5 is the unique minimum (K 4 , 3) stable graph). By Remark 3.13 and by Lemma 3.14, κ(5) = 4 (and K 9 and K 6 + K 7 are the minimum (K 5 , 4) stable graphs).

In the following, if no confusion is possible, we simply denote the integer κ(q) by κ. Lemma 3.17. Suppose that q ≥ 4. If κ is even then stab(K q , κ -1) < e(2K q+ κ 2 -1 ) and stab(K q , κ) ≤ e(K q+ κ

2 + K q+ κ 2 -1 ) . If κ is odd then stab(K q , κ-1) < e(K q+ κ-1 2 +K q+ κ-3 2
) and stab(K q , κ) ≤ e(2K q+ κ-1

2

) .

Proof Recall that, by definition of κ, K q+κ-1 is the only one minimum (K q , κ-1) stable. If κ is even then 2K q+ κ 2 -1 is exactly (K q , κ -1) stable and

K q+ κ 2 + K q+ κ 2 -1 is exactly (K q , κ) stable. If κ is odd then K q+ κ-1 2 + K q+ κ-3 2
is exactly (K q , κ -1) stable and 2K q+ κ-1 2 is exactly (K q , κ) stable. Lemma 3.18. Let q ≥ 3 and p ≥ 0 be two integers. Then, e(K q+2p ) < e(K q+p + K q+p-1 ) if and only if p 2 + p < 1 2 (q -1)(q -2) and e(K q+2p ) = e(K q+p + K q+p-1 ) if and only if p 0 = 1 2 ( 1 + 2(q -1)(q -2) -1) is an integer and p = p 0 . e(K q+2p+1 ) < e(2K q+p ) if and only if (p + 1) 2 < 1 2 (q -1)(q -2) and e(K q+2p+1 ) = e(2K q+p ) if and only if p 1 = 1 2 ( 2(q -1)(q -2) -1) is an integer and p = p 1 .

Proof It is easy to check that e(K q+2p )-e(K q+p +K q+p-1 ) = p 2 +p-1 2 (q-1)(q-2) and e(K q+2p+1 ) -e(2K q+p ) = (p + 1) 2 -1 2 (q -1)(q -2). These polynomials of degree 2 in p have respectively p 0 = 1 2 ( 1 + 2(q -1)(q -2) -1) and p 1 = 1 2 ( 2(q -1)(q -2) -1) as positive roots.

Proof of Theorem 1.6. If κ = 2p then, by Lemma 3.17, stab(K q , κ -1) < e(2K q+ κ 2 -1 ). Since κ -1 = 2(p -1) + 1, by Lemma 3.18, p 2 < 1 2 (q -1)(q -2), that is, κ < 2(q -1)(q -2).

For κ = 2p we have 1 2

(q +2p-1)(q +2p-2) < (q +p-1)(q +p-2) < (q +p-1) 2 ≤ 1 2 (q +2p)(q +2p-1) .

This implies that (B) 2p 2 < 1 2 (q -1)(q -2) ≤ p 2 + p .

Combining (A) and (B) yields

p 2 < 1 2 (q -1)(q -2) ≤ (p + 1) 2 .
This implies that 1 2 (q -1)(q -2) -1 ≤ p < 1 2 (q -1)(q -2) .

Hence, p = ρ(q) = 1 2 (q -1)(q -2) -1. By inequalities (A) and (B), position of 1 2 (q -1)(q -2) in comparison to ρ(q) 2 +ρ(q) determines the parity of κ. Hence, if 1 2 (q -1)(q -2) > ρ(q) 2 + ρ(q) then κ = 2ρ(q) + 1 = 2 1 2 (q -1)(q -2) -1 else κ = 2ρ(q) = 2 1 2 (q -1)(q -2) -2 If there is no minimum disconnected (K q , κ(q)) stable graph then, by definition of κ(q), there exists a connected minimum (K q , κ(q)) stable graph G q distinct from a clique. Note that if such a graph exists then e(G q ) ≤ M in{e(K q+κ(q) ), e(K q+ κ 2 + K q+ κ 2 -1 )} if κ(q) is even or e(G q ) ≤ M in{e(K q+κ(q) ), e(2K q+ κ-1

2

)} if κ(q) is odd .

Conjecture 1.8 states that there is no such graph G q .

  "minimal (H, k) stable graphs" are called "strong (H, k) stable graphs" by the authors. Note that a (H, k) stable graph G is minimal (H, k) stable if and only if for every e ∈ E(G) the graph G -e is exactly (H, k -1) stable. Moreover, a minimal (H, k) stable graph is exactly (H, k) stable.

	Remark 2.8. Clearly, every minimum (H, k) stable graph is minimal (H, k) stable.
	One may ask what happens for components of a (H, k) stable graph. The following
	theorem gives us an answer.

If there exists an edge e of an (H, k) stable graph G such that there are no subgraph isomorphic to H containing e then G -e is a (H, k) stable graph. Hence, we have the following. Lemma 2.7.

[START_REF] Dudek | stable graphs with minimum size(H, k)[END_REF]

Every edge of a minimal (H, k) stable graph is contained in a subgraph isomorphic to H. Consequently, every vertex of a minimal (H, k) stable graph is also contained in a subgraph isomorphic to H.

  let us consider a minimum H-transversal of G i , say T i , and set k i = |T i | -1. Clearly, for each i the graph G i is exactly (H, k i ) stable and the set T = 1≤i≤r T i is a minimum H-transversal of G. Note that the number of elements of T is |T | = r i=1 k i + l and we have |T | > k. Let S be any set of vertices of G such that |S| ≤ |T | -1 and for every i denote by S i the set S ∩ V (G i ).Clearly, there exists i 0

The research of APW was partially sponsored by polish Ministry of Science and Higher Education.

If κ = 2p + 1 then by Lemma 3.17, stab(K q , κ -1) < e(K q+ κ-1 2 + K q+ κ-3

2

).

Since κ -1 = 2p, by Lemma 3.18, p < 1 2 ( 1 + 2(q -1)(q -2) -1), that is, κ < 1 + 2(q -1)(q -2). Lemma 3.19. For every integer q ≥ 4 and κ = κ(q) we have e((κ + 1)K q ) > e(K q+κ ) Proof We have 2(e((κ + 1)K q ) -e(K q+κ )) = κ(q 2 -3q -κ + 1). By Theorem 1.6, κ < 1 + 2(q -1)(q -2) < 3q 2 , and hence κ(q 2 -3q -κ + 1) > 0.

Theorem 3.20. Let q ≥ 4 and suppose that there exists a minimum

By Theorem 2.9, there are integers

We shall prove first that r = 2. In fact, it is clear that r ≥ 2. Let us suppose that r ≥ 3. We have

Hence, e(K q+kr+kr-1+1 ) < e(K q+kr ) + e(K q+kr-1 ) and the graph

and the proof follows.

Note that Theorem 3.20 implies that there exists at most one disconnected (K q , κ) stable graph and this graph, if it exists, is

Proof of Theorem 1.7 By Lemma 3.17 and Theorem 3.20, if κ is odd then

) ≤ e(K q+κ ) (note that, by Lemma 3.18, it may be possible that e(2K q+ κ-1 2 ) = e(K q+κ ) for some values of q), if κ is even then

) ≤ e(K q+κ ) (note that, by Lemma 3.18, it may be possible that e(K q+ κ 2 + K q+ κ 2 -1 ) = e(K q+κ ) for some values of q). For κ = 2p + 1 we have 1 2

(q + 2p)(q + 2p -1) < (q + p -1) 2 < (q + p)(q + p -1) ≤ 1 2 (q + 2p + 1)(q + 2p) .

This implies that (A) p 2 + p < 1 2 (q -1)(q -2) ≤ (p + 1) 2 .


Antonia Maria Masucci
email: antonia.masucci@supelec.fr

Øyvind Ryan
email: oyvindry@ifi.uio.no

Merouane Debbah
email: merouane.debbah@supelec.fr

Merouane Debbah

Polynomial

⋆ Øyvind

Mérouane Debbah

POLYNOMIAL EXPANSION DETECTOR FOR UNIFORM LINEAR ARRAYS

Keywords: Uniform Linear Array, Polynomial Receiver, Random Vandermonde Matrix

In this paper, we design a low complexity linear MMSE decoder to recover the signal transmitted by M mobile users to a base station equipped with N receiving antennas, arranged as a uniform linear array (ULA). The angles of arrival are supposed to be uniformly distributed. As the dimension increases, it is hard to invert the system and therefore recent results on random matrix theory and polynomial expansion detectors are helpful to approximate the optimum linear MMSE receiver. Simulation results confirm the validity of this approximation.

INTRODUCTION

In the last two decades, multiuser detection has shown to be a useful design technology for detecting desired signals from interference and noise. The first optimum multiuser detector was investigated by Verdú in [?] for asynchronous Gaussian multiple access channels based on maximum-likelihood detection. The receiver is considered optimal in the sense that the performance in the absence of Gaussian noise approaches that of a single-user system. Since the computation of the inverse matrix is complex, the authors of [2] introduced a multistage linear receiver to design different linear detectors, such as decorrelating detectors and minimum mean square error detector, approximating the inverse matrix by a polynomial expansion of the correlation matrix. This expansion is formed by a linear combination of the outputs of the individual receiver stages. The polynomial approximations are helpful in practice only if the weights can be calculated more easily than performing matrix inversion. As the optimum weights depend on the eigenvalues of the correlation matrix, they are not easy to calculate either. The combination of polynomial expansion detectors and convergence results from random matrix theory gives in [START_REF] Muller | Design and analysis of low-complexity interference mitigation on vector channels[END_REF] low-complexity detectors to mitigate interference in multiple-input multiple-output (MIMO) systems. In this paper, we design a multiuser receiver using the linear minimum mean square error (MMSE) as a signal detector in a system model where M mobile users are communicating in the uplink with a base station represented by a uniform linear array. The angles of arrival are supposed to be uniformly distributed and the channel matrix can be modeled by a random Vandermonde matrix. This scenario can be met in large surveillance system using radar. Using known results on random matrix theory, in particular on the moments of random Vandermonde matrices with entries on the unit circle [START_REF] Ryan | Asymptotic behaviour of random vandermonde matrices with entries on the unit circle[END_REF], we are able to find the optimum polynomial receiver, optimum with respect to the mean square error (MSE). Note that similar approaches could be used to find optimum receiver for compact rotating MIMO antenna [START_REF] Bains | On sampling issues of a virtually rotating mimo antenna[END_REF]. In the next section, we present the model under observation. In section 3, we describe a low complexity linear MMSE receiver and in section 4 we use recent asynptotic results on Vandermonde matrices to estimate weights. In section 5, simulations show the validaty of the MMSE reciver design and in the last section give the conclusions. Notation. In the following, upper (lower) boldface symbols will be used for matrices (column vectors), whereas lower symbols will represent scalar values. (.) H will represent the hermitian transpose operator, (.) T denotes the transpose operator and || • ||F denotes the Frobenius norm operator. We denote by IM the M × M identity matrix. We let Tr be the (non-normalized) trace for square matrices, defined by Tr(A) = n i=1 aii, where aii are the diagonal elements of the M × M matrix A. We also let trM be the normalized trace, defined by trM (A) = 1 M Tr(A).

SYSTEM MODEL

We consider M mobile users, each with a single antenna, communicating with a base station equipped with N antennas, arranged as a uniform linear array (ULA), as in Figure 1. The N × 1 output signal at the base station is given by

y(t) = M i=1 v(θi)p 1/2 i xi(t) + n(t) = V(θ)P 1/2 x(t) + n(t)
where x is the M × 1 input signal transmitted by the M users

x(t) = [x1(t), . . . , xM (t)] T ,
satisfying E x(t)x(t) H = IM , and n(t) is the additive Gaussian noise such that E n(t)n(t) H = σ 2 IM . We suppose that the components in x(t) and n(t) are independent. The elements of the M × M matrix represent the power with which the users send information. In the case of a line of sight between the mobile users and the base station, the N × M matrix V has the following form

P =      p1 0 • • • 0 0 p2 • • • 0 . . . . . . . . . . . . 0 • • • 0 pM     
V(θ) = 1 √ N [v(θ1), . . . , v(θM )] ,
where v(θ) = 1, e -j2π d λ sin(θ) , . . . , e -j2π d λ (N-1) sin(θ) T , and θ1, . . . , θM are the angles of the users (see Figure 1) and are supposed to be i.i.d. and uniform on ] -π 4 , π 4 [ , d is the interspacing distance between the antennas of the ULA, and λ is the wavelenght of the signal. Assuming a memoryless channel, we can write y = VP 1/2 x + n.

(1)

RECEIVER DESIGN

The linear minimum mean square error (MMSE) detector [?] for the channel in (1) is given by

x P 1/2 V H VP 1/2 + σ 2 I -1 P 1/2 V H y. (2)
Recovering the data in (2), hence, requires the inversion of the matrix between brackets which is a difficult task when the dimensions of this matrix are large. In order to compute this, we remember the Caley-Hamilton theorem, which states that any square matrix satisfies its own characteristic equation

det P 1/2 V H VP 1/2 + σ 2 IM -(µi(pi, λi) + σ 2 )IM = 0,
where λi are the eigenvalues of V H V, pi the eigenvalues of P, and µi(pi, λi) the eigenvalues of P 1/2 V H VP 1/2 (for large matrices, µi can be viewed as functions of pi, λi, due to the almost sure convergence of products of Vandermonde matrices and deterministic matrices [START_REF] Ryan | Convolution operations arising from Vandermonde matrices[END_REF], which is why µi(pi, λi) has been subsituted for the eigenvalues of P 1/2 V H VP 1/2 ). The characteristic polynomial can be written as

M i=0 αi P 1/2 V H VP 1/2 + σ 2 I i = 0
where αi = αi(λ1, . . . , λM ; σ 2 ; p1, . . . , pM ). We obtain

P 1/2 V H VP 1/2 + σ 2 I -1 = = - M i=1 αi α0 P 1/2 V H VP 1/2 + σ 2 I i-1 = - M i=1 αi α0 i-1 j=1 i -1 j σ 2(i-j-1) P 1/2 V H VP 1/2 i-1 ,
which means that we can write the inverse as a matrix polynomial of degree M -1 in P 1/2 V H VP 1/2 such that

P 1/2 V H VP 1/2 + σ 2 I -1 = M -1 i=0 βi P 1/2 V H VP 1/2 i ,
where the coefficients βi = α i α 0

i-1 j=1

i-1 j σ 2(i-j-1) and they only depend on the eigenvalues of V H V, the power P and the noise variance σ 2 . In order to reduce the complexity, this suggests the introduction of multistage receivers that approximate the inverse matrix of the LMMSE receiver by a polynomial expansion of degree

K < M in P 1/2 V H VP 1/2 : xMMSE = K i=1 βi(P 1/2 V H VP 1/2 ) i-1 P 1/2 V H y. (3)
The weight vector β = β1, . . . , βK T can be choosen to optimize some performance measure. We derive the vector of weights such that minimizes the mean square error (MSE) of the estimated vector

xMMSE β = arg min β E   x- K i=1 βi(P 1/2 V H VP 1/2 ) i-1 P 1/2 V H y 2   .
(4) Expressing equation (3) as a linear combination of the vectors zi = (P 1/2 V H VP 1/2 ) i-1 P 1/2 V H y, we can write the expectation in (4) as

E x -Z β 2 ,
where the matrix Z = [z1, . . . , zK ]. We minimize respect to β and we obtain that the optimum weight vector satisfies

β = E Z H Z -1 E Z H x = Φ -1 φ.
The elements of the K × K matrix Φ are given by

Φ(i, j) = E z H i zj = E (x H P 1/2 V H + n H )VP 1/2 (P 1/2 V H VP 1/2 ) i-1 × ×(P 1/2 V H VP 1/2 ) j-1 P 1/2 V H (VP 1/2 x + n) = 1 M Tr P 1/2 V H VP 1/2 i+j + + σ 2 M Tr P 1/2 V H VP 1/2 i+j-1 = 1 M Tr PV H V i+j + σ 2 M Tr PV H V i+j-1 (5)
and the elements of the K × 1 vector φ are

φ(i) = E z H i x = E (x H P 1/2 V H + n H )VP 1/2 (P 1/2 V H VP 1/2 ) i-1 x = 1 M Tr P 1/2 V H VP 1/2 i = 1 M Tr PV H V i . (6)

ASYMPTOTIC ANALYSIS

We need recall some results from random matrix theory. We approximate the optimum weights by using results on the asymptotic moments of random Vandermonde matrices [START_REF] Ryan | Asymptotic behaviour of random vandermonde matrices with entries on the unit circle[END_REF]. In particular, we assume that the size of V grows to infinity while the ratio of the number of columns to rows tends to a limit r > 0.

Definition 1. An N × M random Vandermonde matrix with unit entries has the form

V = 1 √ N      1 • • • 1 e -jω 1 • • • e -jω M . . . • • • . . . e -j(N-1)ω 1 • • • e -j(N-1)ω M     
, where ω1, . . . , ωM are i.i.d. random variables.

We will denote by P(n) the set of all partitions of {1, . . . , n}, and will write W1, . . . , Wr for the blocks of a given ρ ∈ P(n).

Definition 2. For ρ ∈ P(n), we define

Kρ,ω,N = 1 N n+1-|ρ| (0,2π) |ρ| n k=1 1 -e jN(ω b(k-1) -ω b(k) ) 1 -e j(ω b(k-1) -ω b(k) ) dω1 . . . dω |ρ| , ( 7
)
where ωW 1 , . . . , ωW |ρ| are i.i.d. (indexed by the blocks of ρ) with the same distribution of ω and where b(k) is the block of ρ which contains k. If the limit

Kρ,ω = lim N→∞ Kρ,ω,N (8)
exists, it is called a Vandermonde mixed moment expansion coefficient.

In [START_REF] Ryan | Asymptotic behaviour of random vandermonde matrices with entries on the unit circle[END_REF], it has been proved that the limit in (8) exists if the density of ω is continous. The calculation is based on combinatorial computation using crossing partitions since the matrices are not free. We will need the following result [START_REF] Ryan | Asymptotic behaviour of random vandermonde matrices with entries on the unit circle[END_REF]: Theorem 1. Let ω be as in Section 2, i.e. ω = 2π d λ sin(θ) with θ uniform on ] -π 4 , π 4 [ . We have that all limits Kρ,ω exist. Also,

mn = lim N→∞ E trM PV H V n = ρ∈P(n)
Kρ,ωr |ρ|-1 Pρ (9)

when M N → r > 0, Pn = limN→∞ tr(P n ), Pρ = k i=1 PW i .
The quantities (9) are also computable [START_REF] Ryan | Documentation for the random matrix library[END_REF]. Using the previous result in Theorem 1, we obtain an estimation of the optimum weights for the considered model, when N → ∞, M N → r, as follows

= Φ -1 asy φasy, (10)
where the generic element of Φasy is given by

Φasy(i, j) = ρ∈P(i+j) Kρ,ωr |ρ|-1 Pρ + σ 2 ρ∈P(i+j-1) Kρ,ωr |ρ|-1 Pρ = mi+j + σ 2 mi+j-1, and φasy(i) = ρ∈P(i)
Kρ,ωr |ρ|-1 Pρ = mi.

This means that we can write the equation (10) in the following form

     m1 m2 . . . mK      =      m2 + σ 2 m1 • • • mK+1 + σ 2 mK m3 + σ 2 m2 • • • mK+2 + σ 2 mK+1 . . . . . . . . . mK+1 + σ 2 mK • • • m2K + σ 2 m2K-1            β(1) asy β(2) asy . . . β(K) asy      
Note that we have substituted E tr PV H V n with tr PV H V n , which is valid due to almost sure convergence [START_REF] Ryan | Convolution operations arising from Vandermonde matrices[END_REF].

NUMERICAL RESULTS

In this section, we present simulations of the above result. Simulations show the validity of our approximation. Documentation of code for computing the moments of Vandermonde matrices can be found in [START_REF] Ryan | Documentation for the random matrix library[END_REF].

In the Figure 2, we have plotted, increasing M and N , the error

||G -Gasy|| 2 F where G = P 1/2 V H VP 1/2 + σ 2 I -1 P 1/2 V H and Gasy = K i=1 β(i) asy (P 1/2 V H VP 1/2 ) i-1 P 1/2 V H .
The ratio between the increasing M mobile users and N receiving antennas is assumes to be equal to 1/5, 1, 7/5. In Figure 3,||G|| 2 F and ||Gasy|| 2 F are plotted with increasing M and N with ratio and with approximation order K = 5.

In Figure 5, we have plotted the SINR per user, increasing the number of receiving and mobile users with given ratio equal to 1/5, 1, 7/5. In particular, the SINR of user ℓ (ℓ = 1, . . . , M )

reads SINR ℓ = |g ℓ h H ℓ | 2 M i=1, i =ℓ |g ℓ h H i | 2 +σ 2
where g ℓ is the ℓ-th row of the matrix Gasy and h ℓ is the ℓ-th column of the matrix product P 1/2 V. Figure 4 shows good results for order of approximation K = 4. In all simulations, we assume the power P = IM , variance noise σ = √ 0.5, wavelenght λ = 2, and distance between receiving antennas d = 1.

CONCLUSION

In this paper, we have designed a low-complexity linear MMSE receiver through random matrix theory and polynomial expansion detectors. In particular, we have used recent asymptotic results on the moments of random Vandermonde matrices with entries on the unit circle to approximate the optimum weights of the polynomial receiver. Simulation results are presented in order to confirm the validity of our approximation.

[1] S. Verdú, "Minimum probability of error for asynchronous gaussian multiple-access channels," IEEE Transactions on Information Theory, vol. 32, no. 1, pp. 85-96, 1986.

[2] S. Moshavi, E. G. Kanterakis, and D. L. Schilling, "Multistage

  This work was supported by Alcatel-Lucent within the Alcatel-Lucent Chair on flexible radio at Supélec. The work of Antonia Maria Masucci was supported by Microsoft Research through its PhD Scholarship Programme.

Fig. 1 .

 1 Fig. 1. Graphical representation of the system model.

Fig. 2 .Fig. 3 .

 23 Fig. 2. ||G -Gasy|| 2F is plotted increasing M and N with ratio 1/5, 1, 7/5.

Fig. 4 .Fig. 5 .

 45 Fig. 4. SINR per user is plotted for K = 1, 2, 3, 4 increasing M and N with a ratio 3/5.