Dmitry Kuzovkin
email: dmitry.kuzovkin@technicolor.com

Tania Pouli
email: tania.pouli@technicolor.com

Rémi Cozot
email: remi.cozot@irisa.fr

Olivier Le Meur
email: olivier.le_meur@irisa.fr

Jonathan Kervec
email: jonathan.kervec@technicolor.com

Kadi Bouatouch
email: kadi.bouatouch@irisa.fr

Image Selection in Photo Albums

Keywords: Photo selection, Image aesthetics, Image quality

The selection of the best photos in personal albums is a task that is often faced by photographers. This task can become laborious when the photo collection is large and it contains multiple similar photos. Recent advances on image aesthetics and photo importance evaluation has led to the creation of different metrics for automatically assessing a given image. However, these metrics are intended for the independent assessment of an image, without considering the possible context implicitly present within photo albums. In this work, we perform a user study for assessing how users select photos when provided with a complete photo album-a task that better reflects how users may review their personal photos and collections. Using the data provided by our study, we evaluate how existing state-of-the-art photo assessment methods perform relative to user selection, focusing in particular on deep learning based approaches. Finally, we explore a recent framework for adapting independent image scores to collections and evaluate in which scenarios such an adaptation can prove beneficial.

INTRODUCTION

Casual photography has in recent years become an essential part of everyday life, where people tend to document each moment of their life through a photo. Not restricted by storage limitations or camera availability, users often take numerous photos of the same life moment. This is further exaggerated by functionalities such as the burst mode available on most modern smartphones, which produce dozens of nearly identical images to ensure that the perfect moment or expression is optimally captured. As a result, users often evade the responsibility of taking one best shot during the action of photo capture itself. Eventually, they might obtain an extremely large collection of photos, where for each significant moment, the best shot has to be chosen among many similar photos.

Assessment of photographs in a photo album is a non-trivial task. The selection of the best photos is a subjective process that is also highly affected by the photo album context. For instance, a particular photo might appear of low quality when observed in an isolated manner, while inside the belonging album, the same photo may be the best candidate when compared with other similar photos. Thus, the selection of photos in the album is, to a large extent, a comparison-based process.

To facilitate and eventually automate this time-consuming task, computational modelling of these human decisions would be necessary. Recent progress in computer vision and machine learning techniques has led to a wealth of image assessment techniques [START_REF] Datta | Studying Aesthetics in Photographic Images Using a Computational Approach[END_REF][START_REF] Kong | Photo Aesthetics Ranking Network with Attributes and Content Adaptation[END_REF][START_REF] Lu | RAPID: Rating Pictorial Aesthetics Using Deep Learning[END_REF][START_REF] Luo | Photo and Video Quality Evaluation: Focusing on the Subject[END_REF][START_REF] Mai | Composition-Preserving Deep Photo Aesthetics Assessment[END_REF][START_REF] Marchesotti | Assessing the aesthetic quality of photographs using generic image descriptors[END_REF], where an image is usually assigned a ranking score, or a label of high or low aesthetics. Going a step further, certain works have attempted to analyze and understand the image features and general characteristics that affect people's decisions regarding aesthetic quality or beauty of a photograph [START_REF] Isola | What makes an image memorable[END_REF][START_REF] Marchesotti | Learning beautiful (and ugly) attributes[END_REF][START_REF] Tina | Smart Photo Selection: Interpret Gaze As Personal Interest[END_REF]. Among other applications, automatic image assessment can be used to assist different tasks, such as image retrieval [START_REF] Geng | The role of attractiveness in web image search[END_REF][START_REF] Ionescu | Result diversification in social image retrieval: a benchmarking framework[END_REF] or automatic video thumbnailing [START_REF] Song | To click or not to click: Automatic selection of beautiful thumbnails from videos[END_REF].

Although such approaches could help guide automatic photo selection or rating decisions, a few major drawbacks limit their usefulness when applied to users' photo collections. By their nature, such methods are typically trained or optimized on large, general collections accumulating photographs from multiple users. As such, they inherently represent average user preferences modeled over a large variety of content. Further, the models or features learned in such methods are often biased towards professional level photographs, as these types of images are likely to be preferred by the average user when faced with a varied selection of photographs. Finally, the evaluation of each photo is performed independently, where possible connections to other similar photos in the collection are not taken into account.

Given the above limitations, classic approaches for photo aesthetics assessment can be directly applied to more general tasks, such as image retrieval, but may be less capable to reproduce user selections within a photo album, which is the focus of our work. A framework was recently proposed for adapting a general image quality or aesthetics score to the context of a collection [START_REF] Kuzovkin | Context-aware Clustering and Assessment of Photo Collections[END_REF]. By clustering images according to their degree of similarity, individual image scores could be scaled, such that final selection preferred images that were the best in their cluster, even if they were assessed as being low quality independently. Nevertheless, this framework was demonstrated only using a sharpness metric as an independent assessment criterion, which covers only one of a multitude of image characteristics considered both by users and image aesthetics models [START_REF] Marchesotti | Learning beautiful (and ugly) attributes[END_REF][START_REF] Wolters | The Art of Deleting Snapshots[END_REF]. To better define how to assess images within the context of their surrounding collection, a deeper understanding of the behavior of users when faced with this task is necessary. To that end, we perform a user study, where, for a variety of photo albums, users have to place themselves in the role of the photographer and identify the images they would like to keep. Our study considers photo albums covering typical scenarios from vacation albums to specific events such as weddings or birthday parties, allowing us to analyze user behavior and cross-user agreement across different situations. Based on this data, we evaluate existing approaches to assess their suitability for selecting photos within collections in a manner consistent with user preferences. Our evaluation considers several state-of-the-art deep-learning methods for image assessment applied to images independently, as well as their adaption with a clustering-based framework [START_REF] Kuzovkin | Context-aware Clustering and Assessment of Photo Collections[END_REF].

RELATED WORK

Photo selection is a complex task, where a selection decision depends on objective characteristics of image quality (e.g. sharpness, dynamic range and presence of artifacts), on subjectively perceived attributes, which make an image attractive (e.g. photo composition and color style), and on semantic aspects, such as presence of important people in a photo and their face expressions [START_REF] Ceroni | Investigating human behaviors in selecting personal photos to preserve memories[END_REF][START_REF] Marchesotti | Learning beautiful (and ugly) attributes[END_REF][START_REF] Tina | Smart Photo Selection: Interpret Gaze As Personal Interest[END_REF].

To model human preferences in image evaluation, machine learning techniques are often used as they can learn from pre-evaluated data. Hand-crafted features and different types of image descriptors have been used by several approaches to guide the image analysis and evaluation [START_REF] Datta | Studying Aesthetics in Photographic Images Using a Computational Approach[END_REF][START_REF] Luo | Photo and Video Quality Evaluation: Focusing on the Subject[END_REF][START_REF] Marchesotti | Assessing the aesthetic quality of photographs using generic image descriptors[END_REF][START_REF] Mavridaki | A comprehensive aesthetic quality assessment method for natural images using basic rules of photography[END_REF][START_REF] Redi | The beauty of capturing faces: Rating the quality of digital portraits[END_REF][START_REF] Sandhaus | Employing aesthetic principles for automatic photo book layout[END_REF][START_REF] Tang | Content-Based Photo Quality Assessment[END_REF]. Similar to the factors that affect user decisions, such hand-crafted features may be inspired by photography practices or can take the form of more objective quality metrics. However, these features are not fully capable of modeling and predicting the complex and highly subjective notion of aesthetics, which tends to limit their performance in image evaluation.

To encompass higher level concepts that may influence how people evaluate images, methods for assessing image memorability and interestingness were proposed. The term of memorability is linked with the likelihood that a user will recognize the same photograph after a certain time delay [START_REF] Isola | What makes an image memorable[END_REF][START_REF] Khosla | Understanding and Predicting Image Memorability at a Large Scale[END_REF]. The term of interestingness is linked with the ability of a certain image or video to draw attention of a user to its content and keep this attention for an extent of time [START_REF] Demarty | Predicting Interestingness of Visual Content[END_REF][START_REF] Gygli | The Interestingness of Images[END_REF]. Such methods are primarily concerned with evaluating the effect a new, unseen image might have to a user, and hence, they might be not directly applicable to the image selection task.

A more recent trend in photo aesthetics assessment is represented by the wealth of deep learning based methods [START_REF] Chang | Automatic triage for a photo series[END_REF][START_REF] Jin | Image aesthetic predictors based on weighted CNNs[END_REF][START_REF] Kong | Photo Aesthetics Ranking Network with Attributes and Content Adaptation[END_REF][START_REF] Lu | RAPID: Rating Pictorial Aesthetics Using Deep Learning[END_REF][START_REF] Lu | Deep Multi-patch Aggregation Network for Image Style, Aesthetics, and Quality Estimation[END_REF][START_REF] Mai | Composition-Preserving Deep Photo Aesthetics Assessment[END_REF][START_REF] Talebi | NIMA: Neural Image Assessment[END_REF][START_REF] Wang | Event-specific image importance[END_REF], since the task of photo assessment can largely benefit from the abstract feature modeling achieved by convolutional neural networks (CNN). The first CNN-based methods were representing the entire image by a fixed-size cropped patch [START_REF] Lu | RAPID: Rating Pictorial Aesthetics Using Deep Learning[END_REF] or multi-patch aggregation models [START_REF] Lu | Deep Multi-patch Aggregation Network for Image Style, Aesthetics, and Quality Estimation[END_REF], while providing a binary label of low or high aesthetics as an output. Later proposals were able to handle the input images directly, while preserving their aspect ratio, providing a ranking score to an image [START_REF] Jin | Image aesthetic predictors based on weighted CNNs[END_REF][START_REF] Mai | Composition-Preserving Deep Photo Aesthetics Assessment[END_REF] instead of a binary label. Certain methods combine the estimation of technical image quality along with aesthetics, such as the proposal by Talebi et al. [START_REF] Talebi | NIMA: Neural Image Assessment[END_REF]. The model proposed by Kong et al. [20] can assess multiple meaningful photographic attributes from the content, and then estimate an overall ranking score from them.

Despite the number of studies conducted in the area of independent image assessment, little attention has been paid to contextdependent and personalized assessment of images. While the majority of automatic approaches assess the analyzed photo against non-related photos (or features) from the learned dataset, in a reallife scenario each photo is usually evaluated in comparison with similar photos from the same album. In addition, the actual process of photo selection in photo albums is largely affected by individual preferences, the capturing conditions and other properties associated with a particular album and a particular user.

The usage of individual preferences for image aesthetics assessment was demonstrated in a few recent works. In the work of Yeh et al. [START_REF] Yeh | Personalized Photograph Ranking and Selection System[END_REF], the influence of each extracted feature is weighted by the user's adjustments, either provided manually or learned from a photo example. Adaptation of the general photo assessment model was also proposed by Park et al. [START_REF] Park | Personalized Image Aesthetic Quality Assessment by Joint Regression and Ranking[END_REF], where a ranking model for personal preferences is learned from a subset of test images assessed by users, and is then used for adaptation. However, existing user-adaptive methods still require some amount of user interaction to be able to model specific preferences. In absence of such information we can still discover and employ useful patterns in photo albums, such as relations between similar photos taken in the same scene and characteristics of other photos from the entire album.

The characteristics of the associated photo context can be extracted by considering and comparing photos within coherent clusters [START_REF] Ceroni | To Keep or Not to Keep: An Expectation-oriented Photo Selection Method for Personal Photo Collections[END_REF][START_REF] Alexander | Automated event clustering and quality screening of consumer pictures for digital albuming[END_REF], which are typically constructed by detecting the natural boundaries in the caputured image series. Existing approaches for collection clustering are usually based on temporal information [START_REF] Cooper | Temporal Event Clustering for Digital Photo Collections[END_REF][START_REF] Platt | PhotoTOC: automatic clustering for browsing personal photographs[END_REF] or image similarity [START_REF] Ceroni | To Keep or Not to Keep: An Expectation-oriented Photo Selection Method for Personal Photo Collections[END_REF][START_REF] Chu | Automatic Selection of Representative Photo and Smart Thumbnailing Using Near-duplicate Detection[END_REF][START_REF] Kuzovkin | Context-aware Clustering and Assessment of Photo Collections[END_REF][START_REF] Alexander | Automated event clustering and quality screening of consumer pictures for digital albuming[END_REF]. A versatile technique for clustering data where the number of clusters is unknown is hierachical clustering, which was demonstrated as advantageous in photo collection based applications [START_REF] Epshtein | Hierarchical Photo Organization Using Geo-relevance[END_REF][START_REF] Krishnamachari | Image browsing using hierarchical clustering[END_REF][START_REF] Kuzovkin | Context-aware Clustering and Assessment of Photo Collections[END_REF].

Once a clustering is obtained, several approaches may be considered for analyzing the cluster contents and assessing images in their newly defined context. The method proposed by Ceroni et al. [START_REF] Ceroni | To Keep or Not to Keep: An Expectation-oriented Photo Selection Method for Personal Photo Collections[END_REF] utilizes the features collected both on intra-and inter-cluster levels, where a Support Vector Machine based prediction model is learned to predict the selection probabilities for unseen photos. Although their method takes the characteristics of each cluster into account, the learned model is not completely album-adaptive, as it is learned over numerous non-related albums and users. The nature of selection decisions within a given group of similar photos taken in the same scene was studied in the recent work by Chang et al. [START_REF] Chang | Automatic triage for a photo series[END_REF]. In their proposed method, pairwise comparisons are learned with a Siamese CNN, and a relative ranking of images in the group is produced. To handle a similar problem on the album-wise level, the method by Wang et al. [START_REF] Wang | Event-specific image importance[END_REF] complements their Siamese CNN architecture with event type information. Another approach to the context assessment was proposed by Kuzovkin et al. [START_REF] Kuzovkin | Context-aware Clustering and Assessment of Photo Collections[END_REF], where an independent score provided by an external method is adapted to the multi-level photo context extracted with hierarchical clustering.

Although the problem of automatic photo assessment has attracted a lot of attention in recent years, it has been difficult to evaluate the relative merits of the proposed methods in real-life scenarios, such as image selection in complete photo albums. In part, this is due to the absence of appropriate ground truth data for image selection within albums, and, from the other side, due to the lack of knowledge on how users perform this task themselves. In this work, we attempt to address these matters through a user study, which allows us to compare different methods and provides useful insights on the user decision patterns.

USER PHOTO SELECTION ASSESSMENT

To better understand how users select their preferred photographs within a collection, we perform a user study evaluating user selection decisions on a series of different photo collections. For each album, average user preference scores are calculated, indicating how often a particular image was chosen by the study participants, as well as user confidence scores, giving an indication of user agreement for each decision.

User Study Design

Photo Album Data: For the purpose of this user study, we selected six photo albums covering a variety of typical scenarios where amateur photographers may opt to take a large number of photos. Photo albums were selected from several different sources, including PEC dataset [START_REF] Bossard | Event Recognition in Photo Collections with a Stopwatch HMM[END_REF], YFCC100M dataset [START_REF] Thomee | YFCC100M: The New Data in Multimedia Research[END_REF], CUFED dataset [START_REF] Wang | Event-specific image importance[END_REF] and personal albums of the authors. We have limited our search to collections that were not altered by image processing software, and where no evident pre-selection was applied before, thus possibly containing multiple similar and near-duplicate photos, and reflecting a typical modern photo album taken with a digital camera or a smartphone. As the initial collections vary in the number of photos, and to limit the duration of the experiment, we have extracted 50 photos from each album, in their original consecutive order.

Each user was presented with a pair of albums, where one given album represented a typical family event, such as wedding or birthday, while the second represented a travel photo collection. Two example albums are shown in Figure 1 In this case, one pair of similar images is always scored higher than another pair, by each tested method. Thus, a global ranking of photos in an album would not provide an expected selection, if the presence of different scenes is not taken into account.

The analyzed albums demonstrate different characteristics, which allow us to study various real-life scenarios and also identify particular behavior of assessed methods:

• Family Event 1 is a wedding photo album, with a moderate number of repetitive photos (from 1 to 3 photos for each captured moment on average), where the pictures are taken with a semi-professional camera. • Family Event 2 is another wedding photo album captured with a semi-professional camera, but with a higher number of repetitive photos (in some cases, more than 5 photos for each same moment). • Family Event 3 album represents a family birthday gathering taken mostly indoors, with a point-and-shoot camera. This collection presents a large number of fuzzy shots, where the points of interest are not well defined, with a moderate-tohigh number of repetitive photos (from 3 to 6 photos on average). • Travel Album 1 represents a common scenario of vacation photos, where photos of landscapes are mixed together with photos of people posing in front of a landscape. The photos are taken with a semi-professional camera, and multiple highly similar photos of the same moment are taken (more than 5 photos on average). • Travel Album 2 consists only of landscape photos, taken with a semi-professional camera, and the number of repetitions is moderate: 2 to 3 photos are taken for each captured moment. • Travel Album 3 represents photographs taken during an amusement park visit. In this album, no people are present: it consists of multiple pictures of architecture, landscapes and objects (usually from 2 to 5 photos per same scene). The photos are taken with a point-and-shoot camera and various cases of blurred or under-exposed photos are present.

Participants: In total, 30 participants took part in our study (7F/23M), with ages ranging between 24 and 55 years. Each pair of albums was evaluated by 10 different users. All participants could be characterized as amateur or casual photographers, with varying levels of photographic experience and interest.

Task: For each shown photo album, users were presented with a browser-based interface as shown in Figure 2, and were tasked with putting themselves in the role of the photographer of that collection to select the best, more representative, or most important photos in their opinion. No limit was placed on the number of photos selected in each album. Before the start of the experiment, each user was presented with two practice collections of 12 photos each, in order to get familiar with a task and the interface. Then, once the user was ready, they could proceed to selecting photos in the two complete collections assigned to them.

All photos of each album were simultaneously visible as thumbnails, while a larger version of the examined photo was also shown. Users could navigate within the collection freely, with a possibility of viewing all photos before making any selection. Each album had to be completed before moving to the next one. Overall, users took around 20 minutes to complete the entire task.

Analysis of User Selection Results

User selections were recorded and averaged across the ten users assessing each album, obtaining a normalized user preference score for each image. Higher values in this case indicate that an image was selected more often, with a value of one identifying images that were selected by all users. The user preference score for an image i is computed as

s u i = N sel i N user s , (1)
where N sel i is the number of times image i was selected, and N user s is the total number of users.

To better visualize agreement between users, we additionally compute a user confidence score which indicates how decisive the preference score for the particular image. The confidence score is calculated using the inverse of the triangular function over the preference score

c u i = 1 tri(2s u i -1) , (2)
which produces higher confidence when the preference is closer to 1 (every user has selected the image) or to 0 (no user has selected the image). The per-image preference and confidence scores for two example albums used in our study can be seen in Figure 1.

As photo selection is a subjective process, the provided user selections can and do vary between users for each album. To assess the consistency of the selections of different users, the interagreement between observers for such data can be computed in different ways. A common choice to estimate the inter-rater agreement is the Kendall's W coefficient [START_REF] Maurice | The problem of m rankings[END_REF] or Cohen's kappa statistics [START_REF] Cohen | A coefficient of agreement for nominal scales[END_REF]. The Kendall's W coefficient is generally used for ordinal ratings, while kappa statistics are applied for nominal ratings, which is our case (selected/not selected labels). As in our case ten users rate each album, we employ the modified Fleiss' kappa measure [START_REF] Joseph | The equivalence of weighted kappa and the intraclass correlation coefficient as measures of reliability[END_REF], which is a generalization of the original kappa for more than two raters. Values for this measure can be interpreted as follows according to [START_REF] Landis | The Measurement of Observer Agreement for Categorical Data[END_REF]: κ < 0.2 indicates slight agreement, 0.2 ≤ κ < 0.4 indicates fair agreement, 0.4 ≤ κ < 0.6 indicates moderate agreement. The computed Fleiss' kappa values for each album are given in the first row of Table 2.

Most albums assessed lead to a fair to moderate agreement between users, with the exception of two albums where slight agreement was found (Travel album 2 and Travel album 3). Looking at the content of the albums, several interesting conclusions may be drawn. We observe that albums with higher agreement contain a larger number of people portraits, with repetitive similar photos of the same person or group (including Travel album 1, which contains multiple people portraits taken in front of landscapes). At the same time, Travel album 2 and Travel album 3 do not contain people's photos and consist mostly of landscapes and architecture photos. These latter albums demonstrate a larger variance in user selections: the notion of an attractive landscape appears to vary much more than the understanding of a well-captured portrait or group photo.

This suggests that for users it may be easier to perform photo selection of people's photos within an album, even when the presented people are unknown. In fact, closer observation of users' selections during the study reveals that facial expressions were a critical factor guiding their decisions when multiple photos of the same people were present. On the other hand, unique photos of people were almost always selected, irrespective of the quality or expression present.

Further, in Travel album 2 we find a particular example of a photo sequence, where the concept of best photo selection may not be directly applicable: this photo sequence present a panoramiclike capture of the surrounding landscape (seen in the second and third rows of the second album in Figure 1). In this scenario, it is unlikely that a user would want to keep a single photo, as the series is intended for a particular type of post-processing. Indeed in this case we note that users have shared their selection across the series with no particular photo showing higher selection preference. Another similar scenario could occur when capturing a bracketed series for later construction of a high dynamic range image. Such use cases can frequently occur given the general availability of advanced photo processing tools even on mobile devices.

Another challenging example where user selections become divided is when nearly identical photos are present, with no discernible differences in quality. In such cases, we found that user votes were approximately equally shared between the photos in question, meaning that no single image led to a higher preference, despite users wanting to keep at least one representative image of such scenes.

EVALUATION OF IMAGE ASSESSMENT METHODS

As discussed in the previous section, photo selection is a highly subjective task, where, depending on the type of images and collections assessed, even user agreement may be low. As such, automating this task is a daunting challenge. Although a wealth of approaches exist for rating or assessing images, it is unclear how well they perform for selecting images within a photo collection. To evaluate the applicability of different methods for this task, we compare their assigned image scores with our experimental findings. Given the promising advances in this field with recent deep learning based methods, we opt for evaluating the following, CNNbased methods, which were mentioned in Section 2:

• The approach by Kong et al. [START_REF] Kong | Photo Aesthetics Ranking Network with Attributes and Content Adaptation[END_REF], which has recently shown the state-of-the-art performance in the independent image aesthetics assessment. In their method, different photographic attributes are estimated and weighted by the image content, giving proper relevance to what should be considered important in an image. It is possible that their proposed content adaptation might be applicable in photo albums as well. • The method proposed by Jin et al. [START_REF] Jin | Image aesthetic predictors based on weighted CNNs[END_REF], which introduced another fine-tuning scheme with sample weights that should allow to assess images spanning a wide range of aesthetic quality. In addition, it was demonstrated that their method can be applied for automatic image cropping, therefore we can expect it to cope better with similar repetitive images often present in photo albums. • The NIMA method proposed by Talebi et al. [START_REF] Talebi | NIMA: Neural Image Assessment[END_REF], which was designed to estimate both technical image quality and aesthetic attractiveness of an image. In our analysis, we use their NIMA MobileNet [START_REF] Andrew G Howard | Mobilenets: Efficient convolutional neural networks for mobile vision applications[END_REF] version of the CNN architecture.

Performance of Analyzed Methods

The methods assessed provide a ranking score for each processed image, which we compare against the average user preference scores determined in our experiment. To estimate the performance of each analyzed method, we compute the Pearson correlation coefficient between the method's photo assessment scores and the user preference scores. With this approach, even if the scores from each method are not computed on the same scale, we can estimate the extent of correlation between the computed independent image scores and the user evaluations. The computed correlation values are detailed in Table 2 and visualized in Figure 3.

Several observations arise from this correlation analysis. The method by Jin et al. [START_REF] Jin | Image aesthetic predictors based on weighted CNNs[END_REF] often performs poorly in the complete photo albums, despite its potential for dealing with image crops. Nevertheless, for Travel album 1 it shows the best correlation with the user preferences, among all the methods. This album consists of a number of very similar photos of landscapes and landscapes with people, which is possibly the scenario where the approach by Jin et al. [START_REF] Jin | Image aesthetic predictors based on weighted CNNs[END_REF] demonstrates its best performance. In addition, their approach appears to respond to the presence of blur in images, and our observations show that in this album the users often selected the most sharp photos among similar ones.

The NIMA method [START_REF] Talebi | NIMA: Neural Image Assessment[END_REF] performs worse in the albums where people's photos are present, but it shows its best performance in landscape-focused Travel album 2 and Travel album 3. However, the degree of inter-observer agreement is relatively low for these collections, therefore no certain conclusions can be made.

On average, the highest performance is demonstrated by the approach of Kong et al. [START_REF] Kong | Photo Aesthetics Ranking Network with Attributes and Content Adaptation[END_REF]. Despite their primary aim of addressing general aesthetics scoring of photos, their produced score demonstrates a noticeable correlation with the user preferences. Possibly, this is due to their proposed content-dependent weighting scheme, which provides better estimation of different type of the image content. Thus, their approach could be potentially suitable for aiding the selection process as well.

CLUSTERING-BASED SCORE ADAPTATION

As it was shown in the previous section, the applicability of independent image assessment methods is limited in complete photo collections. Although in some cases this could be explained by the method's occasional failure on certain images with complex content, another important reason is that the photo context is not considered by such methods. The illustration of this point is given in Table 1: it can be observed that photos from different scenes in the same album often receive non-comparable scores, which can lead to inaccurate ranking of photos in the album, if it is performed in a global manner.

The direct approach of using independent image scores in a photo album would be to rescale and normalize the scores linearly, in accordance with other scores computed in the album. However, we have observed earlier that the original scores given by methods are often poorly correlated with user preferences. One reason for this could be the elimination of the context present in photo albums, as each photo is assessed independently, without consideration of its surrounding photos. For example, even if multiple pictures of the same scene were taken, suggesting that the user found that scene important, each picture from the scene could receive a low score and be potentially rejected. For this reason, we attempt to utilize the notion of context in the photo assessment within an album.

As discussed in Section 2, a possible approach to define the photo context is through album clustering, where the entire collection is clustered into groups of similar photos. For this purpose, we apply the framework proposed in the approach of Kuzovkin et al. [START_REF] Kuzovkin | Context-aware Clustering and Assessment of Photo Collections[END_REF], where the entire photo collection is clustered into clusters of different similarity levels to define the relevant context for each photo.

The context of a photo is modeled using three enclosed hierarchical levels: collection level (containing the whole collection), scene cluster level (containing photos depicting the same scene), and near-duplicate cluster level (reflecting very similar images). This context hierarchy is obtained by combining time-based clustering together with similarity clustering based on the SIFT descriptors.

To evaluate how the presence of context affects the performance of the state-of-the art methods tested in the previous section, each album from our user study was clustered and its scores were adapted, according to [START_REF] Kuzovkin | Context-aware Clustering and Assessment of Photo Collections[END_REF]. To adapt an image-based score to its associated context, z-scores [START_REF] Kreyszig | Advanced engineering mathematics[END_REF] are computed for each image and for each clustering level, considering both the image score and statistics of the cluster, as follows:

z I L = s I -µ L σ L , (3)
where µ L and σ L denote mean and standard deviation of the scores, computed on one of three levels L ∈ C, SC, N D, which define collection level, scene cluster level and near-duplicate cluster level, respectively. After z-scores for each level are computed, they are combined into the global score Z I , as an average of the adapted scores from three levels.

In Figure 3, we also demonstrate the result of the z-score context adaptation for the analyzed methods. It can be seen that for the methods by Jin et al. [START_REF] Jin | Image aesthetic predictors based on weighted CNNs[END_REF] and NIMA method [START_REF] Talebi | NIMA: Neural Image Assessment[END_REF] the correlation with user preferences increases for the most of the albums, after the performed adaptation. Especially noticeable increase in performance can be observed for NIMA method [START_REF] Talebi | NIMA: Neural Image Assessment[END_REF], which is also reflected in the average correlation increase, shown in Table 3. For the method of Kong et al. [START_REF] Kong | Photo Aesthetics Ranking Network with Attributes and Content Adaptation[END_REF], the performance gain is not particularly evident. In this case, the performance has increased for the albums where similar repetitive sequences of photos are largely present, such as Family event 1 and Travel album 1 (also the kappa user agreement is the highest for these albums). However, for other albums the effect of adaptation is opposite, with a decrease in the correlation. This could be also due to the nature of these albums, as some of them contain a lower number of repetitive photos.

Original scores

Clustering-adapted scores [START_REF] Kuzovkin | Context-aware Clustering and Assessment of Photo Collections[END_REF] Kong et al. [START_REF] Kong | Photo Aesthetics Ranking Network with Attributes and Content Adaptation[END_REF] 0.274 0.261 Jin et al. [START_REF] Jin | Image aesthetic predictors based on weighted CNNs[END_REF] 0.230 0.260 NIMA [START_REF] Talebi | NIMA: Neural Image Assessment[END_REF] 0.169 0.217 Table 3: Average correlation across all albums for original scores and scores after clustering-based adaptation.

It is also important to note that the utilized clustering is less robust in Travel album 2 and Travel album 3, due to occasional large viewpoint changes for the same captured scenes and presence of severe blur in some photos. Additionally, the album Family Event 2 represents a special case for all methods, where the correlation before and after adaptation is rather low, even in presence of multiple similar photos. We found that the number of low-quality photos (such as blurred ones) is smaller in this album, while, as we previously observed, the user preferences in this album were often guided by more complex factors, such as face expressions or people's poses. Due to this, the original scores may not always be accurate, which can in turn lead to an unreliable adaptation. Moreover, the employed adaptation approach sometimes may be too simplistic to model the subtleties of user employed criteria.

CONCLUSION

Our work studies the performance of image assessment methods when applied to the task of photo selection in photo albums. Despite the wealth of work available for assessing the quality or aesthetics of images, most existing methods consider images independently, without knowledge of their surrounding collection or context. To understand how users perform this task, we have collected a selection of photo albums covering different events and quality levels, and conducted a user study on them where users were asked to select which photos they would like to keep in each collection. Our findings suggest that users consider several elements in their decisions, varying from the quality of images, to the depicted scene or people. More interestingly, we find that users often show very different selection decisions between them, highlighting the difficulty of this task for an automated method. The low agreement for some albums also confirms that in certain cases a personalized user preferences modeling would be necessary.

Our comparisons of image evaluation scores from several stateof-the-art methods relative to the results of our user study show that, in most cases, independent image assessment solutions correlate to a limited degree with user selections. At the same time, most methods performed better for albums where user agreement was higher. To assess whether additional knowledge of the context of images could improve the results of automatic image evaluation approaches, we adapted independent image scores using an album clustering approach [START_REF] Kuzovkin | Context-aware Clustering and Assessment of Photo Collections[END_REF]. Although this adaptation showed some benefit in higher user-agreement albums, no clear conclusions could be drawn in lower agreement cases.

Given the large quantities of photographs captured for any key event in our lives, efficient solutions that can aid users in the cumbersome task of photo selection are likely to become increasingly necessary. In this work, we show that existing methods can work for some scenarios, but are globally far from being able to predict user decisions in the context of their own albums and photo collections. Nevertheless, our study provides some insights in user behavior that we hope may serve as a basis for future work, such as further computational analysis which image features affect user selection decisions within photo albums.

Figure 1 :

 1 Figure 1: Demonstration of two albums from our user study (Family event 1 and Travel album 2), with a visualization of user preferences shown around each image and confidence bar charts, which indicate confidence of user preferences. (a) Family event with high confidence of user selections. As it can be seen, in the albums with people's photos, users tend to agree more in their selections. (b) Travel album with lower confidence of user selections. In the albums with landscape photos, users tend to agree less.

Figure 2 :

 2 Figure 2: Interface of the conducted selection user study. User can freely browse through the entire collection and perform an image selection on key press.

Figure 3 :

 3 Figure 3: Per-album correlation between computed scores and user preferences, given for the original independent scores and the scores after clustering-based adaptation. The albums are sorted from the lowest to the highest kappa user agreement.

Table 1 :

 1 . Scoring of different photos in the same album by the analyzed image assessment methods.

		0.55	0.58	0.54	0.53
	Jin et al. [17]	3.13	4.59	4.86	5.06
	NIMA [38]	4.87	5.01	4.76	4.54

Table 2 :

 2 Travel album 2 Travel album 3 Family event 2 Family event 3 Family event 1 Travel album 1 User selections agreement and performance comparison for analyzed methods. The kappa agreement values are given in the first row (the albums are sorted from the lowest to the highest user agreement). The performance values for each method are computed as the Pearson correlation coefficient between user preferences and scores provided by analyzed methods.

	Kappa agreement	0.179	0.210	0.334	0.351	0.393	0.472
	Kong et al. [20]	0.302	0.120	0.175	0.263	0.384	0.397
	Jin et al. [17]	0.258	0.009	0.080	0.075	0.236	0.721
	NIMA [38]	0.269	0.196	0.115	0.044	0.128	0.260


Stéphane Brull
email: brull@univ-tln.fr

Jacques Schneider
email: jacques.schneider@univ-tln.fr

A new approach of the Ellipsoidal Statistical Model

.

1 Introduction.

Solving the Boltzmann equation [START_REF] Cercignani | The Boltzmann equation and its applications[END_REF] involves computing a collisional term which is a very hard task. Different approaches can be considered: either find approximating models of the interaction term that are easier to compute or/and try to improve numerical approximations. We do not want here to give a huge review of all works that have been done in the two directions but to follow a way that was initiated by Bhatnagar, Gross and Krook [START_REF] Bhatnagar-E | A model for collision processes in gases[END_REF] and their famous BGK model. Namely the BGK model consists in replacing the interaction term with the relaxation term

R(f ) = λ(M -f ) (1)
where f = f (t, x, v) is the distribution function with (t, x, v) ∈ R + × R 3 × R 3 , M = M(t, x, v) is the local equilibrium function (Maxwellian distribution) defined by

M(v) = 1 (2πT ) 3/2 exp - (v -u) 2 2T (2)
with

ρ(t, x) = R 3 f dv, u(t, x) = 1 ρ R 3 vf dv, T (t, x) = 1 3ρ R 3 |v -u| 2 f dv. (3)
Here ρ, u, T are the local mass, density and temperature of the distribution function and λ = λ(ρ, T ) is a relaxation coefficient. The interest of such a model is that it inherits most important physical properties of the true one: mass, momentum and energy are conserved, the (mathematical) entropy decays and equilibrium states are the Maxwellian distributions. Unfortunately the BGK model does not allow to recover the true Navier-Stokes equation at the hydrodynamic limit. That is a correct Prandtl number P r which defines the ratio between the viscosity µ and thermal conductivity κ coefficients P r = 5 2 µ κ .

Different models have been proposed to recover the correct Prandtl number. The most famous is the Ellipsoidal Statistical model (ES or ES-BGK model) [START_REF] Holway | Kinetic theory of shock structure using an ellipsoidal distribution function[END_REF]. Other models were proposed later on in [START_REF] Bouchut | A BGK model for small Prandtl number in the Navier-Stokes approximation[END_REF] and finally in [START_REF] Struchtrup | The BGK model with velocity dependent collision frequency[END_REF]. For numerical comparisons between the BGK model, the ES-BGK model and eventually other models one refers to [START_REF] Mieussens | Numerical comparison of BGK-models with proper Prandtl number[END_REF], [START_REF] Andries | Perthame Numerical comparison between the Boltzmann and the ES BGK models for rarefied flows[END_REF].

A new approach of relaxation models can also be found in [START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF]. Roughly speaking the idea is to add different relaxation terms toward "generalized" Maxwellians that are associated to different relaxation coefficients.

Finally one of the author [START_REF] Schneider | Entropic approximation in kinetic Theory[END_REF] has proposed in the case of Maxwellian molecules a relaxation term of the form ρ(G -f ) that fits exactly to the Boltzmann collision operator in the following weak sense

ρ(G -f )m(v)dv = Q(f, f )m(v)dv ∀m ∈ P (4)
where P can be any polynomial space. This model allows to recover the right Prandtl number. Unfortunately this model cannot be extended to other potentials. Here G is the solution to the minimization problem of the entropy functional H(f ) where

H(f ) = R 3 f (ln f -1) dv (5)
under the moment constraints [START_REF] Bhatnagar-E | A model for collision processes in gases[END_REF]. Such minimization problems are the basis of the work by D. Levermore [START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF]).

In this paper we rather look for a single relaxation term of the form λ(G-f ) that can generalize the BGK-type models. The idea is the following: G is defined through a minimization principle of the entropy under moments constraints. More precisely let P be a given polynomial space of maximum order p with vector basis m(v) = (m 1 (v), ..., m N (v)) T , λ and (λ i ) i=1,...,N be some nonnegative real number and f ∈ L 1 be a nonnegative function such that f (v) (1 + |v| p )dv < +∞. Then we introduce the following minimization problem "Find G solution to the minimization problem

G = min g∈C f H(g) (6)
where C f is the set of functions g ≥ 0 s.t the following equations hold

R 3 λ(g -f )m i (v) dv = -λ i R 3 f m i (v)dv, ∀i = 1, ..., N." (7)
λ together with (λ i ) i=1,...,N are seek so as to satisfy the physical properties of Q(f, f ). In particular (λ i ) i=1,...,N are parameters to be chosen later for fitting the hydrodynamic limit of the Boltzmann equation at different level with respect to the Knudsen number (see [START_REF] Cercignani | The Boltzmann equation and its applications[END_REF]) . Notice that the idea is different from the one developed in [START_REF] Schneider | Entropic approximation in kinetic Theory[END_REF] where the relaxation model is set through an approximation theory principle. The question whether this problem admits a solution or not is subjected to the value of λ and (λ i ) i=1,...,N and to the so-called realizability moments problem (see [START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF], [START_REF] Junk | Maximum entropy for reduced moment problems[END_REF], [START_REF] Schneider | Entropic approximation in kinetic Theory[END_REF]).

In this paper we restrict our approach to the case of

P = span[1, v, v ⊗ v].
It is divided as follows. After setting the problem we construct in section 3 a relaxation operator. This relaxation operator shows up to be the ES-BGK model henceforth giving a new approach of it. Section 4 is devoted to the study of the entropy dissipation of our model. We obtain a different proof of the entropy decay than that in [START_REF] Andries | Perthame Entropy condition for the ES BGK model of Boltzmann equation for mono and polyatomic gases[END_REF]. Finally we compare in section 5 the different relaxation models -ES-BGK, the one of Bouchut and Perthame [START_REF] Bouchut | A BGK model for small Prandtl number in the Navier-Stokes approximation[END_REF] and and the one of Struchtrup ([12], [START_REF] Mieussens | Numerical comparison of BGK-models with proper Prandtl number[END_REF])-enjoying an entropy decay and giving the right Prandtl number. We show a general variational principle that gives new proofs for the existence of those models ( [START_REF] Bouchut | A BGK model for small Prandtl number in the Navier-Stokes approximation[END_REF], [START_REF] Struchtrup | The BGK model with velocity dependent collision frequency[END_REF], [START_REF] Mieussens | Numerical comparison of BGK-models with proper Prandtl number[END_REF]) and we end up with a discussion comparing the different models.

2 Setting of the problem.

Recall that in our paper we look for a relaxation model of the form

R(f ) = λ(G -f ) ( 8
)
where G is the solution to the following minimization problem

G = min g∈C f H(g). (9)
C f is the set of functions g ≥ 0 s.t the following equations hold

R 3 (1, v, |v| 2 ) gdv = R 3 (1, v, |v| 2 ) f dv, (10)
R 3 λ(g -f )A(V ) dv = -λ 1 R 3 f A(V )dv, (11)
where A(v) are the Sonine polynomials defined by

A(v) = v ⊗ v - 1 3 |v| 2 Id, (12)
Id being the matrix identity in R 3 and

V = v -u √ T . (13)
u and T are defined by [START_REF] Bagland-P.Degond-M | Moment systems derived from relativistic kinetic equations[END_REF]. λ is the relaxation coefficient to G while λ 1 is the relaxation rate of the moment

R 3 f A(V )dv.
This means that we restrict our approach to the case where the polynomial space is

P = span[1, v, v ⊗ v].
The equations [START_REF] Junk | Maximum entropy for reduced moment problems[END_REF] just express the conservation laws of mass, momemtum and energy. The relation [START_REF] Schneider | Entropic approximation in kinetic Theory[END_REF] can be interpreted in the following sense

R 3 f A(V )dv = R 3 L( f M ) L -1 (A(V )) M(v)dv,
where L(f ) is the linearized Boltzmann operator [START_REF] Cercignani | The Boltzmann equation and its applications[END_REF]. The interest is that

L -1 (A(V )) = α(T, |V |)A(V
) is somehow an "eigenfunctions" of the linearized collision operator. Then R(f ) is chosen to mimic the behaviour of L(f ) or equivalently the behaviour of Q(f, f ) closed to equilibrium. This explains why λ and λ 1 are seek so as to match the hydrodynamic limit of the true Boltzmann equation. Now introducing the stress tensor Θ,

Θ = 1 ρ R 3 c ⊗ c f dv, c = v -u, (14)
(11) can be rewritten in the form

1 ρ R 3 c ⊗ c gdv = (1 - λ 1 λ )Θ + λ 1 λ T Id. ( 15
)
We see that this equation only depends on the fraction λ 1 λ which must be dimensionless.

In order to compare the present result with the other works on the Ellipsoidal Statistical Model, we set ν = 1 -λ 1 λ and (15) takes the form

1 ρ R 3 c ⊗ c gdv = νΘ + (1 -ν)T Id = T . ( 16
)
3 Construction of the Ellipsoidal Statistical Model.

In this section we construct the ES-BGK model performed in three steps.

1. We solve exactly the minimization problem for different values of the ratio between λ and λ 1 relaxation coefficients.

2. By using the hydrodynamic limit of our formal model through a Chapmann-Engskog expansion up to the order 1 one can define a relaxation rate λ 1 for given moments (namely the Sonine polynomials A(v) defined below) of the distribution function independently of the relaxation coefficient λ.

Then the relaxation coefficient λ is defined by equating the different terms appearing in the Chapmann-Engskog expansion at order 0.

3. λ is set so as to match with the Navier-Stokes equation at the hydrodynamic limit with the true Prandtl number P r = 2 3 . From this we recover the ES-BGK with a new interpretation.

3.1 Realisability conditions for the existence of the solution to the minimization problem [START_REF] Mieussens | Numerical comparison of BGK-models with proper Prandtl number[END_REF].

Our first purpose is to construct the relaxation function G. As we are going to see a necessary and sufficient coefficient for the problem (9) to possess a solution is ν ∈ [-1 2 , 1[. This improves the result found in [START_REF] Andries | Perthame Entropy condition for the ES BGK model of Boltzmann equation for mono and polyatomic gases[END_REF]. So let us consider the weighted space

L 1 2 = {f ∈ L 1 s.t. (1 + |v| 2 )f ∈ L 1 }.
Then we have the following result Theorem 1. For all nonnegative functions f ∈ L 1 2 and ν ∈ [-1 2 , 1[, the tensor T is symetric positive definite and the problem (9) admits a unique solution G defined by

G(v) = 1 det(2πT ) exp - 1 2 c, T -1 c . ( 17
)
Conversely, if (9) admits a unique solution for all nonnegative function f ∈ L 1 2 , then ν ∈ [-1 2 , 1[. The existence of a solution to the minimization problem (9) requires the condition C f = ∅. In [START_REF] Andries | Perthame Entropy condition for the ES BGK model of Boltzmann equation for mono and polyatomic gases[END_REF] this is done by exhibiting the function G (17) for T semi definite positive and then proving that this function is the solution to the minimization problem. Nevertheless this can be done by using more general arguments [START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF], [START_REF] Junk | Maximum entropy for reduced moment problems[END_REF], [START_REF] Schneider | Entropic approximation in kinetic Theory[END_REF]. As shown in those papers a necessary and sufficient condition is C f = ∅ for P of maximal degree two.

Lemma 1. Let χ be a symetric definite tensor of order two and consider the set C of functions g ≥ 0 s.t.

gdv = ρ, 1 ρ vgdv = u, 1 3ρ |v -u| 2 gdv = T (18) for some (ρ, u, T ) in R * + × R 3 × R * + and 1 ρ R 3 c ⊗ c gdv = χ.
then C is not empty.

The proof of lemma 1 is given in appendix.

Proof. (Theorem 1)

For the sake of clarity we recall here why C f = ∅ (see [START_REF] Andries | Perthame Entropy condition for the ES BGK model of Boltzmann equation for mono and polyatomic gases[END_REF]). By denoting with θ 1 , θ 2 , θ 3 the eigenvalues of Θ the eigenvalues of T are

(1 -ν) 3 (θ 1 + θ 2 + θ 3 ) + νθ i , i = 1, 2, 3.
These values are obviously nonnegative for ν ∈ [0, 1]. The previous expression writes for i = 1

1 + 2ν 3 θ 1 + (1 -ν) 3 (θ 2 + θ 3 )
and is therefore nonnegative for ν ≥ -1 2 . So T is symmetric definite positive. Then C f is non empty according to Lemma 1. Since the condition (16) is of power 2 in v then the moments (10, 15) are realizable ( [START_REF] Junk | Maximum entropy for reduced moment problems[END_REF], [START_REF] Schneider | Entropic approximation in kinetic Theory[END_REF]) and the minimization problem 9 has a unique solution which writes (17).

Conversely let f be a nonnegative function in L 1 2 and Θ be defined as in (14). Θ is symetric and one can diagonalize it as Θ = P t ∆P where the diagonal terms of ∆ are θ 1 , θ 2 , θ 3 with 0 < θ 1 , θ 2 , θ 3 < 3T and θ 1 +θ 2 +θ 3 = 3T . Now if ( 9) admits a solution one can compute this solution which is exactly given by (17). Both Θ and T -1 can be diagonalized in the same basis such that G writes

G(v) = 1 det(2πT ) exp - 1 2 c, P t D -1 P c .
Clearly each term of the diagonal matrix D -1 or equivalently D satisfies the same conditions as ∆. That is writing (16) in the diagonalization basis

θ i ν + (1 -ν)T ∈ ]0, 3T [ ∀i = 1, 2, 3. (19)
Since those conditions must be valid for all f then letting for example θ 1 tend to 0 gives ν ≤ 1 and letting θ 1 tend to 3T gives ν ≥ -1 2 .

Definitions of λ 1 and λ: the Ellipsoidal Statistical model

The classical way to obtain the hydrodynamic limit of (1) consists in expanding f in power of ε around the local equilibrium function M (Chapman-Engskog expansion ).

f = M(1 + εf (1) + ε 2 f (2) + . . . ), (20)
and then inserting it in [START_REF] Andries | Perthame Numerical comparison between the Boltzmann and the ES BGK models for rarefied flows[END_REF]. In our case we just have to replace the interaction term Q(f, f ) with λ(G-f ) and since G depends itself on f one has to expand G as

G = M(1 + ε ∂G ∂Θ : Θ (1) + . . . ).
Then comparing the Chapman-Engskog expansion of the kinetic equation

∂ ∂t + v • ∇ x f = λ ε (G -f ).
at order 0 gives

∂ ∂t + v • ∇ x M = λ( ∂G ∂Θ : Θ (1) -f (1) ) M.
And so this provides as usual

A(V ) : D(u) + B(V ) √ T • ∇ x T M = λ( ∂G ∂Θ : Θ (1) -f (1) )M, ( 21
)
where the stress tensor D(u) is defined as

D(u) = 1 2 (∇ x u + ∇ x u t ) - 1 3 div(u)Id and
B(V ) = V 2 (V 2 - 5 2
).

One can also expand equation ( 11) (replacing g with G) and compare the left and right hand sides at order 1.

λ

R 3 ∂G ∂Θ : Θ 1 -f (1) Mdv = -λ 1 R 3 f (1) MA(v)dv. (22)
Inserting it in equation ( 21) gives after multiplication with A(V ) and integration

R 3 (A(V ) : D(u))MA(V )dv = -λ 1 R 3 A(V )Mf (1) dv. Hence R 3 A(V )Mf (1) dv = - 1 λ 1 R 3 A(V ) : A(V )Mdv D(u).
This sets a relationship between the relaxation rate λ 1 and the viscosity µ

λ 1 = ρT µ (23)
which defines λ 1 . Remark that this definition does not depend on λ. Now

λ = λ 1 1 -ν = ρT µ(1 -ν) (24)
which entirely defines our relaxation model

R(f ) = ρT µ(1 -ν) (G -f ) (25)
with

G(v) = 1 det(2πT ) exp - 1 2 c, T -1 c ( 26
)
and T is defined by (16) (recall also the equality (24)). This turns out to be the Ellipsoidal Statistical model [START_REF] Holway | Kinetic theory of shock structure using an ellipsoidal distribution function[END_REF].

Thermal conductivity and Prandtl number.

Multiplying (21) with B(V ) and integrating with respect to v yields

R 3 B(V )f (1) Mdv = - 1 3λ R 3 B(V ) • B(V )Mdv ∇ x T √ T = - 5 2λ ∇ x T √ T = - κ ρT ∇ x T √ T .
where κ is the thermal conductivity of our model. This sets the other relationship between λ and κ as

λ = 5 2 ρT κ .
This means that the Prandtl number of this model is

P r = 5 2 µ κ = λ λ 1 = 1 1 -ν
Remark that two important cases of this model are those where 1-ν = -1 2 which gives the important physical case P r = 2 3 that holds for monatomic gas, 2-ν = 0 which gives the classical BGK model [START_REF] Bhatnagar-E | A model for collision processes in gases[END_REF] whose Prandtl number is P r = 1 and is not physical.

H-Theorem.

In this section we want to prove the entropy dissipation law for any -1 2 ≤ ν ≤ 1. Though this has already been done in [START_REF] Andries | Perthame Entropy condition for the ES BGK model of Boltzmann equation for mono and polyatomic gases[END_REF] we give here another proof. We denote with G ν the anisotropic Gaussian function defined by (26) Theorem 2. For all -1 2 ≤ ν ≤ 1 the entropy dissipation term satisfies

D(f ) = -(G ν -f ) ln f dv ≤ 0 where G is defined in (26). Moreover D(f ) < 0 for -1 2 ≤ ν < 1 with equality iff f = M.
Proof. The convexity of x ln x -x yields the inequality

-(G ν -f ) ln f dv ≤ H(G ν ) -H(f ) ≤ H(G ν ) -H(G 1 )
The last inequality is a direct consequence of the minimization problem taking ν = 1. The computation of H(G ν ) gives

H(G ν ) = ρ ln ρ √ 2πdetT - 5 2 ρ.
This function is a strictly convex function of ν since it is the minimum of the entropy functional under linear constraints in ν. Therefore it is enough to prove that G -1 2 ≤ G 1 . Remark that the minimum of H(G ν is the physical entropy

H(M) = 3 2 ρ log ρ 2/3 T -Cρ and is obtained in ν = 0. Next H(G ν ) -H(G 1 ) = 1 2 ρ log detΘ detT
is of the same sign as

p(ν) = detΘ -det(νΘ + (1 -ν)T Id)
Let us recall that θ 1 , θ 2 , θ 3 are the eigenvalues of the symmetric definite positive matrix Θ so that on one side detΘ = i θ i while on the other side

det(Θ -λId) = -i (λ -θ i ) implies detT = i (νθ i -(ν -1)T ).
Then one has

p(ν) = i θ i + i (ν(T -θ i ) -T )
Setting θi = θ i /3 then we have

p(ν) = T 3 ( i θi + i (ν(1 -θi ) -1)).
In particular

p(- 1 2 ) = T 3 3 8 (3 θ1 θ2 θ3 -θ1 θ2 -θ2 θ3 -θ1 θ3 )
Now 0 < θ1 , θ2 , θ3 ≤ 3 and θ1 + θ2 + θ3 = 3 so that using Lagrange multipliers the maximum of p(-1 2 ) is 0 and is obtained for θ1 = θ2 = θ3 = 1. This corresponds to G 0 = M.

Comparison between different models giv-

ing the right Prandtl number.

In this section we want to compare different models of the BGK form giving the right Prandtl number and having an entropy decay. Namely we compare the ES-BGK model, the model of Bouchut and Perthame ( [START_REF] Bouchut | A BGK model for small Prandtl number in the Navier-Stokes approximation[END_REF]) and those of Struchtrup ([12], [START_REF] Mieussens | Numerical comparison of BGK-models with proper Prandtl number[END_REF]). We first recall those models and then derive a general variational principle that can be applied to each of them. Finally we compare them from a theoretical point of view. So let us first recall the model of Bouchut and Perthame ( [START_REF] Bouchut | A BGK model for small Prandtl number in the Navier-Stokes approximation[END_REF]). It writes

R BP (f ) = η λ(V )(M λ -f ) (27)
where η = η(t, x) = η(ρ, T ) and λ(V ) = λ(|V |) is a multiplier acting on

V (t, x, v) = v -u(t, x) [T (t, x)] 1 2
)

Here ρ, u and T are the mass, velocity and temperature of f (3). Finally M λ is a Gaussian distribution

M λ(v) = ρ λ 1 (2πT λ) 3/2 exp - (v -u λ) 2 2T λ (28)
whose parameters (ρ λ, u λ, T λ) are implicitly defined through the relations

R 3 λ(V )(1, v, |v| 2 ) M λdv = R 3 λ(V )(1, v, |v| 2 ) f dv. (29)
λ can be fitted to give the right Prandtl number.

H. Struchtrup ([12]) introduced later on the following relaxation term

R S (f ) = λ1 (V )(Mλ -f ) (30)
where λ1 (V ) is a velocity-dependant collision frequency of the form λ1 = ζ|V | α . Mλ is a Gaussian distribution (see above) whose parameters (ρλ, uλ, Tλ) are implicitly defined through the relations

R 3 λ(V )(1, v, |v| 2 ) Mλdv = R 3 λ(V )(1, v, |v| 2 ) f dv. (31)
Here ζ depends on macroscopic variables and α > 0 can be fitted to give the right Prandtl number. Remark that f does not relax toward Mλ at v = u so that he choose the more appropriate form [START_REF] Mieussens | Numerical comparison of BGK-models with proper Prandtl number[END_REF] λ2

(V ) = ζ(1 + C|V | α ) = ζ λ(V ) (32)
Again constants C, α > 0 can be fitted to give the right Prandtl number. At this stage this last model seems to be just a particular choice of the relaxation model (30) since it suffices to take λ(w) = λ1 (w) or λ(w) = λ2 (w). But it is not since λ(w) satisfies

∃C 1 , C 2 , 0 < C 1 ≤ λ(w) ≤ C 2 < +∞, ∀w ∈ R 3 . ( 33
)
Nevertheless there is some similarity as will be shown in the next section.

Variational principles

As was shown in section 3.1 (Theorem 1), the ES-BGK model enjoys a variational principle. That is G ES is is the solution to the minimization problem

G ES = min g∈C f H(g).
where H(g) is the classical entropy functional (5) and C f is the set of functions satisfying conditions [START_REF] Junk | Maximum entropy for reduced moment problems[END_REF], [START_REF] Schneider | Entropic approximation in kinetic Theory[END_REF] for a given nonnegative function

f ∈ L 1 2 (R 3 ). Now let ω(v) = ω(|v|) be a function such that either ∃ C 1 , C 2 , α > 0 s.t 0 < C 1 ≤ ω(v) ≤ C 2 (1 + |v| α ) ∀v ∈ R 3 . or ω = |v| α
It is clear that for fixed u ∈ R 3 , T > 0 there exists C 2 (u, T ) such that

ω(v)dv = ω((v -u)/ √ T )dv (34) is a measure in L 1 2+α , respectively in V = L 1 2+α ∩ {f ∈ L 1 / |f | |v -u| α < +∞}.
Then one defines the weighted entropy functional

H ω(f ) = R 3 f (ln f -1) ω(v)dv. ( 35
)
One has the following general theorem.

Theorem 3. Let u ∈ R 3 , T > 0 and h = 0 be a nonnegative function in L 1 2+α (respectively in V ) then there exists a unique solution M ω to the minimization problem

M ω = min g∈C h H ω(g). ( 36
)
where C h is the set of function g ≥ 0 s.t. the following relations hold

R 3 ω( v -u T 1 2 )(1, v, |v| 2 ) gdv = R 3 ω( v -u √ T )(1, v, |v| 2 ) hdv. ( 37
)
Moreover this solution is the Gaussian distribution

M ω(v) = ρ ω (2πT ω) 3/2 exp - (v -u ω) 2 2T ω (38)
where (ρ ω, u ω, T ω) are uniquely determined by the relations (37) (replacing g with M ω(v)).

Proof. One can easely extend the proof of Theorem 1 and Remark 1 [START_REF] Schneider | Entropic approximation in kinetic Theory[END_REF] to the case of the weighted entropy H ω.

Now let f ∈ L 1 2+α (respectively in V ) such that ρ = R 3 f dv, u = 1 ρ R 3 vf dv, T = 1 3ρ R 3 |v -u| 2 f dv, (39)
then applying the above theorem with h = f and taking successively ω = λ, ω = λ1 and ω = λ2 proves at the same time the existence of M λ(v) satisfying successively (29), (31) and the variational principle (36). Remark that this simplifies the proof in [START_REF] Bouchut | A BGK model for small Prandtl number in the Navier-Stokes approximation[END_REF] (see Proposition 2 and Theorem 3) for the model R BP (30) since the above variational principle proves at the same time the existence and uniqueness of M λ for f satisfying (39).

Discussion

This subsection is devoted to the comparison between the ES-BGK and the two others models introduced on one hand in [START_REF] Bouchut | A BGK model for small Prandtl number in the Navier-Stokes approximation[END_REF] and on the other hand in [START_REF] Mieussens | Numerical comparison of BGK-models with proper Prandtl number[END_REF], [START_REF] Struchtrup | The BGK model with velocity dependent collision frequency[END_REF].

Let us denote that all the models are solutions to a suitable variational problem. This allows to get the conservation laws and to obtain an Htheorem in the case of [START_REF] Bouchut | A BGK model for small Prandtl number in the Navier-Stokes approximation[END_REF], [START_REF] Struchtrup | The BGK model with velocity dependent collision frequency[END_REF], [START_REF] Mieussens | Numerical comparison of BGK-models with proper Prandtl number[END_REF]. Nevertheless this variational principle -and especially the constraints (37) cannot be used to obtain an H-theorem for the Ellipsoidal Statistical model because of the additional condition [START_REF] Schneider | Entropic approximation in kinetic Theory[END_REF](see Theorem 2). Considering the way to obtain the right Prandtl number one can distinguish two strategies:

• In the ES-BGK model one imposes a relaxation to 0 of the tensor f A(V )dv with a suitable rate λ 1 (23). This leads to the right stress tensor with proper viscosity but not (at first sight) to the right thermal conductivity. Then one compensates this deficiency by modifying the relaxation rate λ to the local equilibrium state (Maxwellian function). Notice at this point two facts. One f tends in reality to an anisotropic Gaussian before this Gaussian becomes the local Maxwellian. Two it is quite remarkable that this anisotropic Gaussian is well defined up to the limit of ν = 1/2 which on one side corresponds to the important physical value of P r = 2/3 and on the other side still gives the Htheorem.

• For the other models the conditions the idea is rather to introduce velocity dependant relaxation coefficients λ, λ1 , λ2 . In [START_REF] Bouchut | A BGK model for small Prandtl number in the Navier-Stokes approximation[END_REF] the relaxation coefficient λ is bounded from above and below (essentially for mathematical reasons) and gives the right Prandtl number. But it is the coefficient η(ρ, T ) that gives the correct viscosity (and henceforth the correct thermal conductivity through P r) -that is the one corresponding to a given interaction law between particles. This is exactly the choice that is made for the Ellipsoidal statistical model. On the other side the functions λ1 , λ2 in [START_REF] Struchtrup | The BGK model with velocity dependent collision frequency[END_REF], [START_REF] Mieussens | Numerical comparison of BGK-models with proper Prandtl number[END_REF] are rather seek on a physical basis, that is to fit with a collision frequency expected to grow up as |v| → +∞. The interest is that from this collision frequency (though not realistic [START_REF] Mieussens | Numerical comparison of BGK-models with proper Prandtl number[END_REF]) one deduces explicit expression for the viscosity and the thermal conductivity.

6 Conclusion and perspectives.

We have constructed a BGK-type operator that replaces the Boltzmann collision operator. This relaxation operator -known as Ellipsoidal statistical model [START_REF] Holway | Kinetic theory of shock structure using an ellipsoidal distribution function[END_REF]-is here presented through a problem of minimization of the entropy under moments constraints. More precisely beside the classical conservation laws one adds a relaxation constraint on the moment of the operator with respect to A(V ) (see [START_REF] Struchtrup | The BGK model with velocity dependent collision frequency[END_REF], [START_REF] Schneider | Entropic approximation in kinetic Theory[END_REF]). This leads to the correct choice of the Prandtl number.

A new proof of entropy dissipation for this model is given (theorem 2) together with a general variational problem for other models [START_REF] Bouchut | A BGK model for small Prandtl number in the Navier-Stokes approximation[END_REF][START_REF] Struchtrup | The BGK model with velocity dependent collision frequency[END_REF][START_REF] Mieussens | Numerical comparison of BGK-models with proper Prandtl number[END_REF] giving the right Prandtl number. We shall address in forthcoming papers the case of polyatomic gases and a possible generalization of this model to more moments constraints so as to fit with higher order(s) in the Chapman-Enskog expansion.

with |v -u| 2 h(α v -u, χ -1 (v -u) )dv

I α = 1 ρ R 3 1 (det(χ))
+ 1 I α R 3 1 (det(χ)) 1 2
|u| 2 h(α v -u, χ -1 (v -u) )dv.

χ being symetric definite and positive such that χ = S 2 . Then,

R 3 |v| 2 g α dv = 1 I α R 3 |S c| 2 h(α|c| 2 )dc + ρ|u| 2 .
But, as the third condition of (18) is satisfied.

1 2 hR 3 3 |v| 2 g α dv = 1 I α R 3 1 (

 123311 (α v -u, χ -1 (v -u) |c 1 | 2 h(|c| 2 )dc R 3 h(|c| 2 )dc .belongs to C f . Indeed the first and the second constraints of (18) on the moments being satisfied let us check that the third relation relation is also satisfied. R

R 3 |S 5 2 R 3 |c 1 |R 3

 352313 c| 2 h(α|c| 2 )dc = tr(χ) α 2 h(|c| 2 )dc, |c 1 | 2 h(|c| 2 )dc R 3 h(|c| 2 )dc.

A proof of lemma 1

For any function h ∈ L 2 1 the function g α defined by


Carlos Carrascal-Manzanares

Alexandre Imperiale

Gilles Rougeron

Vincent Bergeaud

Lionel Lacassagne

A fast implementation of a spectral finite elements method on CPU and GPU applied to ultrasound propagation

Keywords: Spectral Finite Elements, parallelization, vectorization, SIMD, GPU

In this paper we present an optimization of a spectral finite element method implementation. The improvements consisted in the modification of the memory layout of the main algorithmic kernels and in the augmentation of the arithmetic intensity via loop transformations. The code has been deployed on multi-core SIMD machines and GPU. Compared to our starting point, i.e. the original scalar sequential code, we achieved a speed up of ×228 on CPU. We present comparisons with the SPECFEM2D code that prove the good performances of our implementation on similar cases. On GPU, a hybrid solution is investigated.

Introduction

In the context of ultrasonic non-destructive testing (NDT) simulation, a new finite element method based on the spectral finite elements and the domain decomposition method has been developed in the CIVA software [START_REF] Impériale | UT simulation using a fully automated 3D hybrid model: application to planar backwall breaking defects inspection[END_REF]. As the computational cost of such a numerical method is expensive, our goal is to reduce the simulation time by an efficient use of the hardware to reach the speed expected by the NDT engineers who are mainly equipped with standard PCs or workstations.

First, we will perform an overview of the spectral finite elements and the domain decomposition method. Secondly, the equation will be split into workable stages to analyse the performance and which optimizations can be efficient, like memory layout optimizations and loop transformations. Thirdly, an analysis of accuracy, hardware architecture, scalability and vectorization on several multi-cores machines will be presented. We will compare our code performance with SPECFEM2D, a well-known software implementation of SFEM. A feasibility study will be also performed on GPU. Lastly, we will conclude and draw some perspectives.

Overview of the Numerical Method

The ECHO library simulates ultrasounds controls to detect echoes on cracks. It has to model the propagation of waves inside a component several order bigger. For this, it combines the domain decomposition method and the Spectral Finite Element Method (SFEM) to solve the problem.

The SFEM is a powerful tool for solving the wave equation. It relies on high-order finite element spaces while conserving a fully explicit scheme, thanks to the mass lumping technique [START_REF] Cohen | Higher-order finite elements with mass-lumping for the 1d wave equation[END_REF]. For example, applying an explicit time discretization of order two, which is stable under the so-called Courant-Friedrichs-Lewy condition on the time step, the fully discrete equation obtained reads

M → P n+1 -2 → P n + → P n-1 ∆t 2 + K → P n = M → F n (1)
where → P is the unknown, n represents the time step, M is the mass matrix, K is the stiffness matrix and → F is the source. In practice, in order to apply the SFEM, the domain needs to be discretized by a mesh composed of hexahedral elements, which are known to be a major challenge for meshing algorithms [START_REF] Shepherd | Hexahedral mesh generation constraints[END_REF].

This method has been extensively used: to solve transient wave equations [START_REF] Cohen | Higher-Order Numerical Methods for Transient Wave Equations[END_REF], elastodynamic multi-dimensional problems [START_REF] Casadei | A mortar spectral/finite element method for complex 2D and 3D elastodynamic problems[END_REF] or targeting seismic problems [START_REF] Komatitsch | Introduction to the spectral element method for three-dimensional seismic wave propagation[END_REF][7] using also GPU hardware [START_REF] Komatitsch | Porting a high-order finite-element earthquake modeling application to NVIDIA graphics cards using CUDA[END_REF] and distributed computation [START_REF] Komatitsch | Modeling the propagation of elastic waves using spectral elements on a cluster of 192 GPUs[END_REF] [START_REF] Komatitsch | High-order finite-element seismic wave propagation modeling with MPI on a large GPU cluster[END_REF].

The region of interest, where the cracks are situated, is decomposed into deformed cubes, which are called sub-domains. The cracks are situated as boundaries between two or more sub-domains where the wave propagation is altered or interrupted. Each subdomain corresponds to a deformed predefined grid. Every element is associated with a reference spectral finite element. In this reference element, high-order lagrangian functions are defined using Gauss-Lobatto points [START_REF] Durufle | Influence of Gauss and Gauss-Lobatto quadrature rules on the accuracy of a quadrilateral finite element method in the time domain[END_REF]. Using an approximation order q leads to q + 1 symmetrical DoFs per dimension. In order to reduce sub-domains communications we impose one-to-one interfaces and same DoF position within common interfaces. This strategy allows us to solve explicitly interface problems arising from domain decomposition. Hence the computational burden is mostly related to computational steps.

Optimization of the Stiffness Matrix

Local Stiffness Matrix

Since the mass matrix is diagonal, the main computational effort (75% of the total execution time) comes from the product of the stiffness matrix with the previously computed solution vector. Thanks to the SFEM, the structure of the global stiffness matrix is composed of overlapping local stiffness matrices related to every element in the mesh.

In the context of ultrasound propagation modelling, storing the global stiffness is numerically intractable, due to the huge number of elements in the mesh. Hence, we opt for an unassembled strategy where the product between the local stiffness matrices and the finite element vector corresponding to each element is performed on the fly. The only constraint to the unassembled strategy is that the local products within a sub-domain with two elements sharing DoFs (neighbours) cannot be computed at the same time. For this, we use a colouring strategy that gives a color to each element, knowing that two neighbours cannot have the same colour (that is, 4 colors in 2D, 8 colors in 3D). Then, each group of colours can benefit from parallel computing. The colouring strategy is also explicit and does not add heavy complexity. A similar strategy was followed in [START_REF] Komatitsch | Porting a high-order finite-element earthquake modeling application to NVIDIA graphics cards using CUDA[END_REF], with the added complexity of dealing with unstructured meshes. In the following, working on the elements of the same colour is supposed, eliminating all constraints regarding concurrent accesses in the parallelization.

Stage Decomposition

In this paper, we shall address the case of a two-dimensional acoustic wave equation and a propagation inside an isotropic homogeneous material, but other cases, as three dimensions, elastic wave equation and anisotropic materials, have also been treated. The resulting equation is Eq. ( 2), where the input U represents the local finite element vector of size (q + 1) 2 , and the output is the vector with the product. A is the gradient matrix, and A T is its transposition. co(DF) is the cofactor of the linear transformation F and J F the jacobian of F. These terms are related to the transformation and they are applied to each local DoF, a vector of two dimensions. ĉ is the wave speed, which is constant because of the homogeneous isotropic nature of the material. From an algorithmic point of view, this equation can be divided in 5 stages in accordance with their main function.

V Stage 4 = A T Stage 3 ⎛ ⎝ ĉ co(DF) T co(DF) J F ⎞ ⎠ xi Stage 2 A Stage 1 ⋅ U Stage 0 (2)
More specifically, stage 0 finds the degrees of freedom corresponding to a local element on the global vector and then copies them to a local vector. Stage 1 computes the gradient as a product of an invariant matrix by this local vector. Stage 2 consist in two two-dimensional matrix products, adding some scalar products for each DoF. Stage 3 computes the transposed gradient, which are the same calculations as in stage 1. This last one is finally accumulated on the global vector by means of stage 4.

Memory Layout

In order to select which optimization to apply, we first compute the arithmetic intensity (AI) which is the ratio between computations and memory accesses. A small AI usually leads to memory bound algorithm. The original code has an average AI value of 0.5 for approximation order 2 (further details are shown in Table 1). The main target of the following optimizations is to reduce the amount of memory accesses. Let us focus on stage 1 and then apply the optimizations to stage 3, as it is mostly the same calculation. The gradient matrix, A, is a regular sparse matrix. Most of the coefficients are zero, and the non-zero ones follow a regular structure. There is a different matrix A for each dimension. Each (i, j) coefficient shows the relation between the DoF i and the DoF j. Thanks to the finite element of reference, the DoFs elements are symmetrically situated in each dimension with repeated pattern. Considering second order, the reference element size is 9 DoFs and thus there are only 9 different coefficients. As shown in Figure 2, we decided to store once each coefficient, creating a reorganized copy of the local vector to match the new and smaller matrix product. In this example two 9×9 matrix products are reduced to six 3 × 3 matrix products. As Figure 2 shows, the second vector is reorganized. This new memory layout simplifies the code because multiple chained loops with strided accesses become only one loop performing basic small products. We apply the same ideas to stage 3, but using the transposed matrix. Other stages are modified to suit the new layout. Stage 2, being in the middle of this optimizations, benefits largely regarding both memory accesses and computations. New AI values can be found in Table 1. This technique should reduce the number of cache misses.

Loop Transformation

Loop unwinding and scalarization [START_REF] Kennedy | Optimizing Compilers for Modern Architectures: A Dependence-based Approach[END_REF][13] have been applied. It pursues operator fusion and reduces the amount of memory accesses (or increasing the AI). The code is divided in three pieces: memory loads, computations and memory stores. The data is loaded into and stored from scalar variables. The strided accesses are reduced, therefore allowing the compiler to optimize the code freely. Complex codes with multiple reused values can benefit largely, although the resulting code may be longer. After applying this transformation in all stages, the code has an average AI of 13.7 (see Table 1), which means it is computed bound.

Benchmark & Analysis

To benchmark our code, we selected a configuration of a two-dimensional mesh of 64 million global DoFs (that is 8,000 for each dimension), with a polynomial approximation order from 2 to 10. The code solving the stiffness matrix product is intended to be used from order 2 to 5, but the range has been extended to gather more information. As the global points are fixed, increasing the elements size results in decreasing the number of total elements. That allows to evaluate the relation between the order and the optimizations. It should be highlighted that computing each local element separately leads to compute much more DoFs than 64 million because of repeated DoFs in borders. It should also be noted that the bigger the element is, the smaller the border ratio is, and thus less DoFs are computed. 64 million global points on 32-bit floating points leads to a 1.4GB-size problem. The unit of measure selected is cycles per point (cpp), that is total clock cycles (RDTSC instruction) divided per computed DoFs [see Table 2]. All measurements correspond to one time step. The processors used are specified in each case. Turbo-boost and hyperthreading are always disabled. Code is compiled with Intel C++ Compiler 17.0 and highest optimization flag (O3). The source code is parallelized using a # pragma omp parallel OpenMP directive and vectorized using only # pragma ivdep directives.

Accuracy Analysis

Single precision floating point values are used, although the original code used double precision values. This change have been decided based upon the relative error when using single precision compared with a double precision result. To do so, the biggest difference between the points of a defined time-step is compared with the biggest value in that timestep. After 500 time-steps, the relative error is for all approximation orders smaller than 0.01%, and after 5000 time-steps smaller than 0.1%. The error grows linearly with the number of steps and can be limited.

Impact Analysis

Four machines have been used for benchmarking: Westmere (WSM) 2 × 6 cores Intel Xeon Processor X5690 3.46Ghz, Ivy Bridge (IVB) 2 × 4 cores Intel Xeon Processor E3-1245 v2 3.45GHz, Haswell (HSW) 2 × 8 cores Intel Xeon Processor E5-2630 v3 2.4GHz and Broadwell (BDW) 2 × 10 cores Intel Xeon E5-2640 v4 2.4GHz. The following test is performed with the IVB machine. The optimized code is faster in every case than the older code using a scalar mono-thread execution (Figure 3). The cache misses have been reduced in a 50% and the acceleration factor is on average 3. It seems to use less cpp before approximation order 7 and much more after order 8. By using Intel VTune Amplifier XE, an increase of instruction cache misses has been measured from insignificant to around 10%. This is related to the increasing size of the code with the order, from 240 lines for approximation order 2 to 1520 lines for approximation order 10. In order to know if the issue can be overcomed varying the processor architecture, this version has been tested with WSM, IVB and HSW. The L1 cache sizes used are 64KB, 16KB and 32KB respectively. As shown in Figure 4, HSW does not meet the problem until the approximation of order 12, using a smaller L1 cache than WSM. HSW must implement a different algorithm to handle the instructions cache misses, but the information about these algorithms is scarce, so it is just a guess. Regarding why WSM performs better than IVB, the L1 cache sizes difference may be the reason. Then using HSW, the acceleration factor is on average 4.5.

Scaling and Vectorization Efficiency

For clarity, only the first four approximation orders are kept during this section. Then the acceleration factor is on average 5.3. The machine used is the HSW. Our new code benefits of great vectorization and scalability efficiency. The test compares a scalar monothread execution against 2 up to 16 scalar OpenMP threads. Before optimizations, the code had an efficiency of around 22%. Now the efficiency goes from 97% using two threads to 90% using sixteen (Figure 5).

Then we evaluated the impact of different SIMD extensions, comprising SSE4.2 and AVX2 instructions without and with FMA. The theoretical speed up factors are 4, 8 and 16 respectively. Before, the vectorization efficiency with AVX instruction was of 26%. The efficiency now is respectively around 100%, 75% and 35% (Figure 6). Some NUMA issues in bi-processor configuration have been found when combining OpenMP with SSE4.2 and AVX. As HSW has two processors, the first one is used efficiently but the second one suffers a fall of performance. Using SSE4.2 instructions (Figure 7), the first eight threads have an efficiency of around 95%, but then it falls to 50%. On the other hand, AVX performs around 90%, then decreases to 30%, to increase again to 80%. This roller coaster is closely related to our memory layout. The loop is divided between the threads, which process process elements that can be situated in the same line and be 'reasonably' contiguous, or through different lines. Depending on the number of threads and the problem size, this can highly affect the NUMA nodes, or not at all. The code has been tested with different NUMA node sizes, and there are always some strong disturbances using the second one. Finally, the 8-core AVX version using the HSW machine is ×228 faster than the original scalar code.

Comparative Benchmark with SPECFEM2D

SPECFEM2D [START_REF] Tromp | Spectral-element and adjoint methods in seismology[END_REF] is the reference software for seismic wave propagation in twodimensional acoustic and elastic media. Both codes use spectral finite elements for simulation but at different scales: ECHO is coded for NDT and SPECFEM2D is used in geodynamics. Until now, only the product of the stiffness matrix has been discussed. For comparison, the capability to solve the entire wave equation (Eq. ( 1)) has been added, using the same methodology when possible.

The machine used is the BDW. An approximation order of 2 is used, because SPECFEM2D only allows orders 1 and 2. SPECFEM2D uses a square mesh generated with the pre-processing finite element mesh generator Gmsh. The same source and similar variables has been settled between the two codes. SPECFEM2D has been modified to remove all output and related checks it usually does, to compare only computation time. 2,000 time-steps are used. In this case, both codes use a call to CPU time, so the elapsed time is given in seconds, not cycles. The measured scale is time per element. Both of them have been compiled with the Intel compiler, ICC in our code and IFC for SPECFEM2D. OpenMP and AVX2 vectorization flags have been activated. To be fair, we can not compare multi-threading execution. SPECFEM2D is designed to use distributed architectures via MPI, whereas our code is intended to use OpenMP in single workstations. Our code benefits of great scalability but SPECFEM2D suffer with MPI communication overhead, which is pointless using this kind of architecture. We are therefore comparing mono-thread executions. As shown in Table 3, our code has a substantial speed up factor. There are multiple reasons. First, SPECFEM2D aims at treating more general geometries than our code and therefore supports unstructured meshes. This more generic data structure may provoke indirect accesses and more complex memory layouts. Secondly, according to the Intel optimization report, few loops have been vectorized, because vector dependence prevents vectorization.

GPU: a Feasibility Study

Our optimized stiffness algorithm has a high AI and a good parallel efficiency. It should therefore benefit from GPU computational power. Our interest is to target one single GPU, as did Huthwaite [START_REF] Huthwaite | Accelerated finite element elastodynamic simulations using the GPU[END_REF], but it can be done using several GPUs [START_REF] Komatitsch | Modeling the propagation of elastic waves using spectral elements on a cluster of 192 GPUs[END_REF]. The are two main issues regarding the efficient use of GPUs: transfer time and computation time.

The input memory (which is going to be read during stage 0) and the output memory (which is written during stage 4) need to be transferred. The constant coefficients are negligible. Furthermore, these transfers will not fit the GPU global memory in a context of multi-domain problems, where the meshes need dozens of gigabytes [START_REF] Komatitsch | Porting a high-order finite-element earthquake modeling application to NVIDIA graphics cards using CUDA[END_REF]. Our strategy is to load (potentially several) domains in order to fill in the GPU memory and compute one temporal iteration before replacing it with a new domain. These memory transfers can be asynchronous, so they can be hidden. Our tests use a NVIDIA GeForce GTX TITAN, 6GB of global memory, Kepler GK110, 0.8GHz, compute ability 3.5 and 2048 single precision cores. This GPU is connected through a PCI Express x16 2.0. The code is written in CUDA and compiled with NVCC. We consider a two-dimensional mesh of 640 million global DoFs. The problem size is 4.2GB, half input, half output. The bandwidth achieved is 6.29GB/s out of the ideal 6.4GB/s transfer bandwidth given by the tool CUDA Bandwidth Test (that is 749ms). This bandwidth could be higher using a motherboard connection with architecture PCI 3.0, going up to 10.8GB/s.

Regarding the computing time, there are up to 64 warps of 32 threads each and 65536 registers to be distributed equally. Occupancy is fulfilled if all the 64 warps are in use, that is using 2048 threads and having 32 registers for each one. A single thread can not use more than 255 registers. Balancing the warp used against register availability is crucial, because high occupancy is not always better [START_REF] Volkov | Better Performance at Lower Occupancy[END_REF].

The GPU version of the code computes one element using one thread. Shared memory is not used because threads does not access same positions. The cudaOccupancy-MaxPotentialBlocks function helps us to calculate the maximal size of one block considering how many registers the kernel uses (it depends on the approximation order, see Table 4). The code presents two main issues. The first one is that a single thread requires a lot of registers and thus the occupancy decreases. The second issue is about memory layout. As each thread load an element, and the elements are not adjacent, the threads cannot load coalesced memory. The GPU version is slower than our fastest CPU code considering execution time. Our current work is trying to limit the register number to increase the GPU occupancy. We plan to use the maxrregcount flag combined with kernel launch bounds to control the number of registers. By increasing the warps used and the occupancy, the hypothesis is that the performance will be increased.

Even when the GPU is not still fully used, we plan to develop a hybrid computation strategy, following the ideas of Papadrakakis [START_REF] Papadrakakis | A new era in scientific computing: Domain decomposition methods in hybrid CPU-GPU architectures[END_REF]. As all the information regarding the use of different approximation orders is available, we will be able to add an auto-tuning step, as did Dollinger [START_REF] Dollinger | CPU+GPU load balance guided by execution time prediction[END_REF]. We expect to predict the execution time launching a small test-case with different approximation orders and sub-domains sizes, and then decide consequently the placement.

Conclusion and Perspectives

A fast SFEM implementation has been presented. We have performed memory layout optimizations and code transformations to increase the arithmetic intensity and reduce the cache misses. Further loop optimizations lead to a code automatically vectorized and parallelized which is several order faster than the original code. A comparative benchmark of execution time has been performed between SPECFEM2D and our code, show-ing a significant performance improvement. Finally, a feasibility study regarding the use of GPUs has been conducted.

A particular case of the wave equation has been presented. Our next step will be to apply the work done to a wider range of cases: dimension, materials, transformations. Using these code transformations means having one source file for each approximation order (potentially 2 to 10), which size can be considerable. Combining the different binary parameters leads to 144 files for each hardware. We therefore plan to generate automatically these files from another program.

Figure 1 .

 1 Figure 1. From left to right: a 2D surface with a parametric defect; the same surface divided into six sub-domains; a deformed pre-defined grid; a finite reference element and 16 DoFs located by a quadrature rule of order 3.

Figure 2 .

 2 Figure2. Schema of stage 1 calculations. One gradient for each dimension is converted into several small matrix-vector products using the same smaller matrix.

Figure 3 .

 3 Figure 3. Cpp by approximation order+1 comparing a scalar mono-thread execution using a IVB processor.

Figure 4 .

 4 Figure 4. Performance of the code using WSM, IVB and HSW processors. The abscissa has been extended to observe the deferred effect.

Figure 5 .

 5 Figure 5. Scalability of the new code.

Figure 6 .

 6 Figure 6. Vectorization impact of the new code.

Figure 7 .

 7 Figure 7. Scalability of the SSE4.2 code.Figure 8. Scalability of the AVX code.

Figure 8 .

 8 Figure 7. Scalability of the SSE4.2 code.Figure 8. Scalability of the AVX code.

Table 1 .

 1 Computations and memory accesses for the original code, the memory layout optimizations and the loop transformations in each stage for the approximation order 2.

	Computation / IOs	Stage 0	Stage 1	Stage 2	Stage 3	Stage 4	Total AI
	Original	0 18	108 234	198 315	108 225	18 27	432 810 = 0.5
	Memory Layout	0 27	108 216	117 45	108 216	27 36	360 540 = 0.7
	Loop Transformation	0 9	90 0	162 0	90 0	27 18	369 27 = 13.7

Table 2 .

 2 Millions of elements and computed DoFs against approximation order for 64 million points mesh (values are rounded).

	Order	2	3	4	5	6	7	8	9	10
	Elements (M)	16	7	4	2.5	1.8	1.3	1	0.79	0.064
	Computed DoFs (M)	144	114	100	92	87	84	81	79	77

Table 3 .

 3 DoFs, elements and execution time of our code and SPECFEM2D.

		DoFs	Elements	Total seconds	Seconds per element
	SPECFEM2D	130561	32480	1.46 × 10 -2	4.48 × 10 -7
	ECHO	133225	33124	1.76 × 10 -3	5.30 × 10 -8
				Speed up factor	×8.45

Table 4 .

 4 Values related to the GPU algorithm. Each column from left to right: approximation order+1, registers used in each kernel, block size, active warps in GPU at the same time, its proportion against 64 total warps, GPU execution time and CPU execution time (seconds) using the fastest CPU version (HSW mono-socket using AVX instructions).

	Order+1 Registers/thread	Threads/block	Active warps	Occupancy	GPU time	CPU time
	3	71	32x28	28	43.75%	1.57	1.00
	4	101	32x16	16	25.00%	1.82	0.81
	5	135	32x12	12	18.75%	2.20	0.73
	6	181	32x8	8	12.50%	2.98	0.74