Giuseppe Boccignone
email: giuseppe.boccignone@unimi.it

Donatello Conte
email: donatello.conte@univ-tours.fr

Vittorio Cuculo
email: vittorio.cuculo@unimi.it

Raffaella Lanzarotti
email: lanzarotti@di.unimi.it

AMHUSE: a multimodal dataset for HUmour SEnsing

Keywords: Human-centered computing → Human computer interaction (HCI), Graphical user interfaces, Emotion recognition, Affective computing, Dataset, Multimodal dataset, Amusement evaluation, Annotation tool

come

INTRODUCTION

Emotions lie behind human-human and human-computer interaction [START_REF] Nussbaum | Upheavals of thought: The intelligence of emotions[END_REF][START_REF] Picard | Affective computing[END_REF][START_REF] Reeves | The Media Equation: How People Treat Computers, Television, and New Media Like Real People and Places[END_REF] and the resulting behavioural landscape is complex and hardly decomposable into its affect and cognition dimensions [START_REF] Pessoa | On the relationship between emotion and cognition[END_REF]. However, to bring affective computing [START_REF] Picard | Affective computing[END_REF] into the interaction game, the issue of gathering good affective data is crucial for learning models and subsequent benchmarking [START_REF] Picard | Toward machine emotional intelligence: Analysis of affective physiological state[END_REF].

Early research in affective computing was mostly focused on facial expression recognition and, by and large, datasets were built on a single modality, namely images or videos [START_REF] Calvo | Affect detection: An interdisciplinary review of models, methods, and their applications[END_REF][START_REF] Sariyanidi | Automatic analysis of facial affect: A survey of registration, representation, and recognition[END_REF]. To truly capture the multidimensionality of emotion, in the last fifteen years the number of public repositories has grown larger, where behavioral , , data gathered in realistic and natural setting experiments have been recorded by multiple modalities [START_REF] Calvo | Affect detection: An interdisciplinary review of models, methods, and their applications[END_REF][START_REF] Sidney | A review and meta-analysis of multimodal affect detection systems[END_REF][START_REF] Gunes | Categorical and dimensional affect analysis in continuous input: Current trends and future directions[END_REF][START_REF] Mihalis | Continuous prediction of spontaneous affect from multiple cues and modalities in valence-arousal space[END_REF][START_REF] Picard | Toward machine emotional intelligence: Analysis of affective physiological state[END_REF][START_REF] Poria | A review of affective computing: From unimodal analysis to multimodal fusion[END_REF].

In such perspective, the main contributions of this study can be summarised as follows. First, we present a novel public multimodal dataset focusing on a positive emotion, namely amusement, which has not been covered in other datasets; signals have been gathered from multiple signal sources: RGB image sequences; depth sequences from a Kinect 3D scanner system; physiological signals, namely, temperature, electrodermal activity and blood volume pulse.

Second, we provide DANTE (Dimensional ANnotation Tool for Emotions), a simple and efficient web-based tool for annotating emotions in the arousal (relaxed vs. aroused) and valence (pleasant vs. unpleasant) dimensions [START_REF] Russell | Core affect and the psychological construction of emotion[END_REF].

Data and annotation tool are freely available at http://phuselab. di.unimi.it/resources.php.

In this study, we use the term "amusement" for referring to the affective state evoked by humorous material [START_REF] David R Herring | Coherent with laughter: Subjective experience, behavior, and physiological responses during amusement and joy[END_REF]. The rationale behind the choice of amusement is that a large body of research has investigated negative emotions, but, surprisingly enough, much less is known about positive emotions, markedly their autonomic response. As Shiota et al. [START_REF] Shiota | Feeling good: autonomic nervous system responding in five positive emotions[END_REF] put it , "this lack of attention is consistent with the historical underrepresentation of positive emotions in psychological research, and with the still common perception among theorists that positive emotions have fewer implications for evolutionary fitness, are less differentiated, and have less distinct impact on motivation and behaviour than is true of the negative emotions." However, positive emotions have broad implications for cognition, physiology, behaviour, and more generally for human well-being [START_REF] Barbara L Fredrickson | The value of positive emotions[END_REF][START_REF] David R Herring | Coherent with laughter: Subjective experience, behavior, and physiological responses during amusement and joy[END_REF][START_REF] Shiota | Feeling good: autonomic nervous system responding in five positive emotions[END_REF][START_REF] Sturm | Damage to left frontal regulatory circuits produces greater positive emotional reactivity in frontotemporal dementia[END_REF]. Further, they serve important social interactions, facilitating approach behaviour, motivating social engagement, promoting new social connections, and reversing the physiological activation caused by negative emotions. They represent a hard task, since inducing more subtle autonomic responses, as opposed to negative emotions that are usually characterized by greater physiological activation. Apart from these theoretical points, positive emotions are likely to play a key role for a number of humancomputer applications (e.g, videogames [START_REF] Christian | Gaming well: links between videogames and flourishing mental health[END_REF], HCI [START_REF] Devillers | Multimodal data collection of human-robot humorous interactions in the Joker project[END_REF]), gaining a great interest for a range of fields, notably the gaming industry and advertisement.

To the best of our knowledge, this is the first multimodal dataset focusing on this specific emotion.

Note that the dataset is the result of an experiment concerning amusement elicitation using italian spoken video clips. A previous effort adopting italian language was described in [START_REF] Esposito | The New Italian Audio and Video Emotional Database[END_REF], though in a different setting and without considering physiological signals. Indeed, multimodal datasets, such as the one presented here, allow to explore human emotional experience more in depth, but leave some problems open. Among them, cogently when analysing the continuous dimensions of emotion, there is the necessity of continuous annotations. As summarized in [START_REF] Sariyanidi | Automatic analysis of facial affect: A survey of registration, representation, and recognition[END_REF]: "Labelling data is a challenging and laborious task, particularly for spontaneously displayed expressions and emotions. " Moreover, combining multiple annotations is even harder [START_REF] Metallinou | Annotation and processing of continuous emotional attributes: Challenges and opportunities[END_REF]. In this respect the DANTE tool is friendly and intuitive, allowing the user to seamlessly annotate valence and arousal separately, during the video reproduction of either the stimuli and the subject's recorded behaviour.

RELATED WORKS

In the last decade, several benchmarks for emotions detection have been presented. Yet, new datasets are proposed to cope with novel challenges: multiple data sources, spontaneous facial expressions, continuous emotional space, and so on.

The Kanade et al. dataset was put forward in 2000 [START_REF] Kanade | Comprehensive database for facial expression analysis[END_REF] and subsequently extended in 2010 [START_REF] Lucey | The extended cohn-kanade dataset (ck+): A complete dataset for action unit and emotion-specified expression[END_REF]. It is widely known with the acronym CK (CK+ for the extended version). Its main characteristic is the annotation of images by means of FACS [START_REF] Ekman | What the face reveals: Basic and applied studies of spontaneous expression using the Facial Action Coding System (FACS)[END_REF]. CK+ includes 593 images sequences from 123 subjects. The image sequences vary in duration and incorporate the onset to peak formation of the facial expressions. Full FACS coding of peak frames is provided. These frames are also labeled with subject's impression of each of the 7 basic emotion categories: Anger, Contempt, Disgust, Fear, Happy, Sadness and Surprise. Therefore, this benchmark is used for testing emotion recognition methods measuring their accuracies. The main problem of this dataset is that the images contain unnatural and exaggerated expressions.

A first step towards a more realistic dataset was taken by Pantic et al. [START_REF] Pantic | Web-based database for facial expression analysis[END_REF] in 2005 with the MMI benchmark. It is composed by 2900 videos for 79 subjects. As well as CK+, MMI is annotated in terms of Facial Action Coding System by FACS experts. MMI shows more naturalistic view of faces and subjects are captured in frontal and in profile views. Similarly to CK, images are annotated with emotion labels, and algorithms on this dataset are evaluated with respect to detection accuracy.

In CK+ the actors were expected to mime a requested affective expression, while in MMI the expression was triggered by asking the subjects to watch a video that was supposed to activate the desired emotion. Research in psychology and social sciences has shown that presence or absence of motor mimicry behaviour can serve as an indicator for emotion inference [START_REF] Adolphs | Recognizing emotion from facial expressions: Psychological and neurological mechanisms[END_REF][START_REF] Wood | Fashioning the face: sensorimotor simulation contributes to facial expression recognition[END_REF]. Thus, Bilakhia et al. [START_REF] Bilakhia | The MAHNOB Mimicry Database: A database of naturalistic human interactions[END_REF] proposed a database (MAHNOB mimicry database) suitable for investigation of mimicry and negotiation behaviour. The dataset consists of 54 recordings of face-to-face interactions. It is a multimodal dataset capturing audio and visual data, from different points of view. The data have been fully annotated for 15 out of 54 sessions, in terms of gestures (hand gestures, head movements, etc.).

The evaluation of this database is done in terms of cross-correlation between the annotated and the detected gestures.

All videos of the previous benchmarks were produced in a "labcontrolled" recording environment. In [START_REF] Dhall | Collecting Large, Richly Annotated Facial-Expression Databases from Movies[END_REF] the authors proposed a new facial expression database (AFEW, Acted Facial Expressions in the Wild) consisting in clips of videos obtained from movies (SFEW, Dhall et al. [START_REF] Dhall | Static facial expression analysis in tough conditions: Data, evaluation protocol and benchmark[END_REF], is the static subset of AFEW). By extracting data from public movies, the authors were able to annotate them with dense information about the subjects (theme of the scene, emotion of the actors, information about the context, etc.).

OPEN EmoRec II [START_REF] Rukavina | OPEN EmoRec II -A Multimodal Corpus of Human-Computer Interaction[END_REF] is an open multimodal corpus as the result of an HCI-experiment, concerning the solving of 6 sequences of a mental trainer designed to induce different emotions. The corpus contains sensory signal of video, audio, physiology (SCL, respiration, BVP, EMG Corrugator supercilii and Zygomaticus Major) and facial reaction annotations.

Soleymani et al. [START_REF] Soleymani | A multimodal database for affect recognition and implicit tagging[END_REF] proposed a database, named MAHNOB-HCI, of multimodal recordings of participant responses to affectively stimulating movie excerpts, images and videos. The recordings of this database are precisely synchronized and multimodality allows to study the simultaneous emotional responses using different channels. The videos are annotated by each participant, by means of a self annotation form. In particular annotations concern an emotional label (neutral, anxiety, amusement, sadness, joy, disgust, anger, surprise, and fear), and a measure (on a nine point scale) for arousal, valence, dominance, and predictability ( [START_REF] Fontaine | The world of emotions is not two-dimensional[END_REF]).

Nowadays, we are moving more and more towards a continuous emotional space [START_REF] Fontaine | The world of emotions is not two-dimensional[END_REF] and the measure of emotion is often shaped in terms of arousal and valence. The Audio Visual Emotion Recognition Challenge (AVEC [START_REF] Ringeval | Av+ ec 2015: The first affect recognition challenge bridging across audio, video, and physiological data[END_REF]) has been proposed, focusing on affect analysis as a regression problem. The RECOLA database deployed for this challenge [START_REF] Ringeval | Introducing the RECOLA multimodal corpus of remote collaborative and affective interactions[END_REF] addresses multimodal data recordings in the context of remote collaborative work. Emotion has to be detected in terms of continuous time and continuous valued dimensional affect in the two dimensions of arousal and valence; the Concordance Correlation Coefficient (CCC) [START_REF] Lawrence | A concordance correlation coefficient to evaluate reproducibility[END_REF] has been chosen as evaluation measure.

In the same vein, our multimodal AMHUSE dataset is conceived following such research trend, but focusing on the positive emotion of amusement. We also cope with the labelling problem, by exploiting DANTE to produce video annotations continuously and not on a frame-by-frame basis. The goal is to provide a suitable tool for continuous affect state modelling and benchmarking.

PROPOSED DATASET

Participants. AMHUSE collects the data of 36 different subjects, who agreed with the scientific use of the recorded material. The participants were 9 females and 27 males, with an age varying from 18 to 54 years old (µ = 26.7 and σ = 8.8). All participants were Italian speaking and did not receive any monetary contribution for the experiment.

Procedure. The experiment was conducted on different days, with the same acquisition protocol and in the same room, namely the PHuSe Laboratory at the department of Computer Science of the Università degli Studi di Milano, Italy. Each participant was welcomed and did receive instructions on the experiment without

Dataset

Subjects Audio RGB Depth EEG ECG BVP EDA EMG R T G FP Pose AU MAHNOB-HCI [START_REF] Soleymani | A multimodal database for affect recognition and implicit tagging[END_REF] 30

✓ ✓ - ✓ ✓ - ✓ - ✓ ✓ ✓ - - - RECOLA [44] 46(18) ✓ ✓ - - ✓ - ✓ - - - - - ✓ ✓ OPEN_EmoRec_II [45] 30 ✓ ✓ - - - ✓ ✓ ✓ ✓ - - - - - AMHUSE 36 - ✓ ✓ - - ✓ ✓ - -✓ - ✓ - ✓ Table 1:

Dataset

Annotator Type Emotion space MAHNOB-HCI [START_REF] Soleymani | A multimodal database for affect recognition and implicit tagging[END_REF] S D PAD + 9 tags RECOLA [START_REF] Ringeval | Introducing the RECOLA multimodal corpus of remote collaborative and affective interactions[END_REF] E(6) + S C VA + 5 tags OPEN_EmoRec_II [START_REF] Rukavina | OPEN EmoRec II -A Multimodal Corpus of Human-Computer Interaction[END_REF] E(4) + S D VA + 6 tags AMHUSE E(4) + S C VA Table 2: Emotional annotations provided in each corpus. S = Self report, E = External. In brackets the number of external annotators. C = Continuous, D = Discrete. PAD = Pleasure, Arousal, Dominance, VA = Valence, Arousal specific details on the purpose of the study, in order to avoid biased reactions. In particular, they were told that they would have to watch four videos, for the total duration of about 10 minutes. At the end of each video, they were asked to annotate the emotional state they felt when watching the movie. Each participant was asked to sit in front of a screen equipped with stereo speakers and with their backs toward a neutral background, as shown on the left of Fig. 1. Physiological data were acquired via e-Health Sensor Platform together with an Arduino UNO, connected to three wired sensors worn by participants. Sensors were placed all on the left hand, leaving the option to annotate the videos with the right hand. Specifically the sensors are: 1) a body temperature sensor, placed on the little finger; 2) a Galvanic Skin Response sensor (GSR), recording the electrical resistance between the medial phalanges of the middle and index fingers; 3) a pulse oximeter sensor, placed on the ring finger. A detail of such setting is shown in Fig. 1, on the right. We encountered artefacts, especially for the GSR signal, in some subjects affected by "cold hands": hence, we waited for the subject's adaptation. A few subjects were excluded since post-analysis revealed artefacts due to excessive movement or adjustment of the device.

In addition to the physiological sensors, a Microsoft Kinect v2 camera, was placed on top of the screen to record the participant reactions in terms of color and depth streams, together with facial landmarks both in two and three dimensional spaces. Once inserted anonymous information such as sex, age and nationality, the stimulus video started playing together with the data recording.

Stimuli. The stimuli (video clips) have been chosen being aware that people have different sense of humour. Thus, aiming at stimulating amusement in each subject at least once, each clip features a different kind of comicality. All videos were in Italian language, all participants were Italian speaking, and no one had never seen the stimuli before.

The four videos chosen as stimuli and shown consecutively are:

(1) A fragment of a documentary1 made by an Italian TV program on the topic of Hammond organ instruments. This video was chosen to induce in the subjects a neutral state and to create a baseline for the following measurements (duration: 56 seconds). (2) A fake movie trailer 2 by the Italian comedian Marcello Macchia. This video is a parody of "The Sixth Sense" drama movie. Marcello Macchia's parodies are always very biting and are mostly appreciated among a younger audience. For this reason, the target of this clip are people between 18 and 30 years old (duration: 87 seconds). (3) A satirical gag of Maurizio Crozza 3 , one of the most appreciated comedians in Italy, specialized in political satire. The target of this video are people who enjoy political satire (duration: 25 seconds). (4) A snippet of a sketch of Aldo, Giovanni & Giacomo 4 . This trio of comedians exploits a quite classical type of comicality, and they are appreciated by people of all ages. For this reason, the target of this video is the widest of all three (duration: 117 seconds).

Gathered data. As mentioned earlier, several spontaneous responses of the participants were gathered during the experiment. One is their facial expressions. For this purpose we adopted a low cost Microsoft Kinect Sensor for Xbox One, equipped with the Kinect Adapter for Windows. This sensor presents a color camera with a full HD resolution of 1920 × 1080 pixels, a frequency of 30 fps and a FOV of 84.1 × 53.8 degrees resulting in an average of about 22 × 20 pixels per degree (Fig 2a). The depth camera, instead, is able to record images with a resolution of 512 × 424 pixels, a frequency of 30 fps and a field of view (FOV) of 70.6 × 60 degrees, resulting in about 7 × 7 pixels per degree. Depth frames are stored in binary format, each value is a 16-bit unsigned integer (distance in mm). Nevertheless, also a grayscale version of each frame is provided (cfr. Fig. 2b).

Even if the camera captures only the frontal view of the participant, the depth stream permits the extraction of 3D head and body movements. Moreover, the wide angle lens of the color camera captures part of the body, arms and hands, which also carry important cues about the affective state of the participant. Furthermore, the Microsoft Kinect SDK extracts in real-time and for each frame a pool of 1347 face points. Such points rely on the camera's 3D space, as shown in Fig. 2c, resulting in a very dense mesh of participant's face. The landmarks are translated in color and depth camera space, making it easy to work with such spaces. As to data compression, we aggregated all the color frames in a video for each session, varying the frame rate and the video duration to pair the stimulus within an AVI container using MPEG-2 codec. The videos were cropped based on the white background, in order to give more focus on subject's reactions, eventually resulting in a resolution of 1024 × 768 pixels. We also extracted Action Units (AU) activations, relying on a freely available AU detector [START_REF] Baltrušaitis | Cross-dataset learning and person-specific normalisation for automatic action unit detection[END_REF]. The AU detector provides, at each frame, the activation level of the following AUs AU k , k = 1, 2, 4, 5, 6, 7, 9, 10, 12, 14, 15, 17, 20, 23, 25, 26, 45, plus the binary presence of AU 28.

Along with the visible cues, a set of physiological signals was acquired during the experiment. Electrodermal activity (EDA) was measured via a Galvanic Skin Response (GSR) sensor, gauging conductance between two points. The level of conductivity is provided in terms of microSiemens (µS) and Volts (V), as well as its reciprocal resistance, expressed in Ohms (Ω). Put simple, a GSR sensor measures the electrical skin resistance in presence of sweat produced by the body: when a high condition of sweating occurs, the electrical skin resistance drops down. Emotions with a prominent presence of positive or negative arousal, such as excitement, stress or fear can induce fluctuations of skin conductivity [START_REF] Peter | Looking at pictures: Affective, facial, visceral, and behavioral reactions[END_REF][START_REF] Nakasone | Emotion Recognition from Electromyography and Skin Conductance[END_REF]. Skin temperature (in Celsius degrees) was also recorded since it changes in different emotional states [START_REF] Plutchik | The psychophysiology of skin temperature: A critical review[END_REF]. Eventually, a blood volume pulse (BVP) sensor was used as a non invasive mean to obtain an indirect measure of the heart rate, via the arterial oxygen saturation of hemoglobin. This signal is strictly correlated with the heart rate measured via ECG and increases in presence of pleasant stimuli [START_REF] Selvaraj | Assessment of heart rate variability derived from finger-tip photoplethysmography as compared to electrocardiography[END_REF]. Heart rate values are provided as the number of contractions of the heart per minute (bpm). For completeness, values of blood oxygenation are also included in the dataset, though there is not a direct relationship between this kind of signal and emotional states.

All the physiological signals (Fig. 3) were recorded at 40 Hz during each session, and provided in CSV format. Acquisition software was developed adopting MATLAB Parallel Computing Toolbox. Each worker is responsible for a task: video playback, Kinect recording and physio recording. All workers are kept in sync and share messages by relying on the MATLAB built-in synchronization mechanism.

Annotations. At the end of each stimulus, all subjects were asked to evaluate their levels of pleasure, arousal and dominance (controlling vs controlled feeling, e.g. anger vs. fear). The annotations were made by using the AffectButton [START_REF] Broekens | AffectButton: A method for reliable and valid affective self-report[END_REF]. The latter was chosen for its intuitive usage: an emoticon representation of emotion, as shown in Fig. 4, where no prior training is required. The Affect Button measures the emotional states in three dimensions: pleasure, arousal and dominance, following the PAD emotional state model [START_REF] Mehrabian | An approach to environmental psychology[END_REF]. To employ AffectButton, the user moves the mouse cursor over the button and the emoticon-like face changes expression accordingly. Every position over the button maps to a corresponding point in the PAD space, with values ranging from -1 to 1, for each dimension. The user clicks when the face expression corresponds to his/her emotional state and the chosen value is saved in the dataset.

Along with self reports, a team of four annotators (3 males and 1 female, with previous experience) was engaged to annotate in the continuous valence/arousal space the reactions of all subjects to each stimulus. The annotators first performed the annotation of a couple of sequences to become familiar with the annotation interface, then started with the real data, observing the participant expressions only. Given the amount of data to be annotated, the annotation process was performed remotely via web browser, in order to split the workload in time. The subject videos were annotated along the arousal and valence dimensions, separately and timecontinuously by using the annotation tool 'DANTE', developed for this purpose and presented in the next section. The frequency of annotations is 25 Hz, provided in CSV format, with values ranging from -1 to 1 and a step of 0.001. DANTE -Dimensional ANnotation Tool for Emotions. Several emotional annotation tools were taken into account and tested. Some of these, such as 'Feeltrace' [START_REF] Douglas-Cowie | A new emotion database: considerations, sources and scope[END_REF] and its successor 'Gtrace' were excluded since the beginning because they rely on software intended to be installed locally and do not allow remote annotations. Aiming at web-based annotation, we considered the 'Valence/Arousal Online Annotation Tool' released together with the AFEW-VA dataset [START_REF] Kossaifi | AFEW for Valence and Arousal estimation In-The-Wild[END_REF] and the 'ANNEMO' [START_REF] Ringeval | Introducing the RECOLA multimodal corpus of remote collaborative and affective interactions[END_REF] annotation tool.

The first requires people to annotate video clips frame by frame, providing continuous annotations for valence and arousal within the range [-10, 10]. We find that this tool could be adopted exclusively for very short videos, as occurs in the dataset AFEW-VA (average length of 50 frames). Moreover, we noticed that the process of per-frame annotation, even if it overcomes the problems related to the delays between the annotation and the video, introduces a bias resulting in sharp annotation signals that hardly follow the dynamics of an expression. In AFEW-VA it has been pointed out another drawback concerning continuous annotation tools: annotators could have a lapse of concentration and inaccuracy due to the sensitivity of the slider. However, we observe that this concern decreases as the annotator training increases.

The 'ANNEMO' tool, instead, is the one which best matches the desired features: separated annotations while playing video, and remote web-based framework. Nevertheless, we experienced some limitations, which may interfere with the annotation process and the usability of the platform. Some of the drawbacks are: i) the absence of indication for the videos already annotated and those to be done; ii) lack of recording with fixed rate; iii) possibility to save the annotations only on text files; iv) missing administration interface to manage annotators and videos; v) unable to differentiate annotators in groups, providing different videos. For these and other minor UI-related reasons, we decided to develop DANTE (Dimensional ANnotation Tool for Emotions), by taking advantage of widely adopted programming languages such as PHP, JavaScript and HTML, backed with a mySQL database.

DANTE allows to create a new annotator account via a dedicated administration page and to assign him/her a random unique identifier corresponding to a private URL for accessing the personal annotation page. The web interface (Fig. 5) presents two main parts: a sidebar on the left, which lists all the videos assigned to the specific annotator, marked with an icon to distinguish between the already annotated videos and those to be done. The center of the page is dedicated to annotation itself. Besides showing the video to be annotated and a sliding bar (with values ranging from -1 to +1 and a step of 0.001), we included a SAM (Self Assessment Manikin) visualization specific for the selected affective dimension (arousal or valence), to help annotators. The actual recording of annotations occurs when the video reaches the end, with a fixed rate specified in the configuration file (default is 25) and saved in CSV format files or directly in the database.

BASELINE RESULTS

Self reports As to the self-report annotations of the AMHUSE dataset, we collected discrete PAD annotations from all the 36 subjects over the four stimuli, as previously described. Table 3 summarises the mean and standard deviation of the annotations, for each stimulus and each emotional dimension. The ICC(3,k) figure of merit (intraclass correlation coefficient, two-way randomized) was used to calculate the agreement between the subjects, where k = 36 is the number of the only raters of interest. The coefficient leads to ICC = 0.844, obtained as the mean of the ICC for pleasure (= 0.664), arousal (= 0.921) and dominance (= 0.948). This shows that on average, all the subjects had remarkable agreement on their ratings, validating stimulus' effectiveness. The visualization in Fig. 6 shows, as expected, a neutral state with very low arousal for the first stimulus, and a general high amusement state for the other videos. Interestingly enough, as shown in Fig. 7 and differently from what expected, the satirical video clip (Stimulus 3) received most annotations with a value of arousal and pleasure higher than 0.5. We surmise that selected subjects were politically interested and, on average, on the same wing. External annotations. In addition to the time-continuous annotations, we derived an aggregated annotation by means of an estimator that weights the annotations of each rater by his/her respective agreement with the others, defining an individual evaluator confidence score. Such technique is called Evaluator Weighted Estimator (EWE) [START_REF] Grimm | Evaluation of natural emotions using self assessment manikins[END_REF]. The individual evaluator score for the annotator k and emotional dimension i ∈ {V , A} over N annotations is

r i k = N n=1 x i n,k -µ i k x M LE,i n -µ M LE,i N n=1 x i n,k -µ i k 2 N n=1 x M LE,i n -µ M LE,i 2 , (1)
where µ i k is the mean annotation of the evaluator k, µ M LE,i is the mean value considering all the evaluators. This score considers noise in the individual annotation as well as his/her grade of experience, with r i k = 0 interpreted as a completely unreliable evaluator. Using these measures as weights, we obtain for each annotation x i n , the corresponding "gold standard":

xi n = 1 K k =1 r i k K k=i r i k x i n,k . (2)
Figure 8 shows an example of annotation together with the weights for each annotator and the resulting gold standard. In order to show the temporal dynamics of annotations, in Fig. 9 we visualize the gold standard of a subject in terms of valence, arousal and time. Consistently with theoretical predictions, it shows a behaviour typical of a mean-reverting random walk [START_REF] Kuppens | Feelings change: accounting for individual differences in the temporal dynamics of affect[END_REF][START_REF] Oravecz | A hierarchical latent stochastic differential equation model for affective dynamics[END_REF].

Other analyses were performed on the annotations to assess the agreement between annotators. Namely, the Cronbach's α [START_REF] Lee | Coefficient alpha and the internal structure of tests[END_REF], the mean Pearson's correlation coefficient and the Concordance Correlation Coefficient (CCC) [START_REF] Lawrence | A concordance correlation coefficient to evaluate reproducibility[END_REF]. Results indicate a good interrater reliability for the valence, whilst, as expected, a poorer value for arousal is observed. Indeed, it is well known that the level of arousal is more difficult to distinguish than valence, resulting in a lower agreement between the annotators.

In particular, the Cronbach's α is an estimate of the consistency between annotations. The valence's value is 0.842 and falls in the range 0.8 ≤ α < 0.9, which is to be considered a good internal consistency, while is equal to 0.325 for the arousal. The mean correlation coefficient, instead, is equal to 0.310 for arousal and 0.742 for valence. The CCC reflects such trend, giving 0.093 for the arousal and 0.509 for the valence.

CONCLUSIONS AND PERSPECTIVES

Two are the main contributions of this study: i) AMHUSE, a multimodal dataset of induced amusement emotional states, and ii) DANTE, a novel web-based annotation tool for valence and arousal continuous values. Both are made available to the research community. Recordings include multimodal data (video, depth sequence, facial landmarks and AUs, EDA, BVP and temperature) gathered from 36 participants who watched and rated their emotional response to 4 stimuli (videos). The subjects' reactions were also annotated by 4 annotators on all recorded sequences in terms of valence and arousal. Significant correlation were found between the participant self-reports. Note that, since there is a gender imbalance (3:1 males to females) both in terms of participants and annotators this could be a possible issue to take into account when using the dataset. fMRI studies suggest that women show greater correlates to humour appreciation than men; however, Chan's study [START_REF] Chan | Neural Correlates of Sex/Gender Differences in Humor Processing for Different Joke Types[END_REF] more subtly identified gender differences in amusement that were specific to particular types of jokes. In this terms, gender differences should be averaged, because of the different kinds of humour stimuli. More controversial are physiological measures: men and women respond electrodermally similarly to pleasant stimuli, except for erotic, where men show larger SCRs than women [START_REF] Bradley | Emotion and motivation II: sex differences in picture processing[END_REF].

Amusement is important from a theoretical standpoint. Its distinction with respect to joy [START_REF] David R Herring | Coherent with laughter: Subjective experience, behavior, and physiological responses during amusement and joy[END_REF], the subtle autonomic behaviour with respect to that of negative emotions, the involvement of key expressive components such as smiling [START_REF] Paula M Niedenthal | The Simulation of Smiles (SIMS) model: Embodied simulation and the meaning of facial expression[END_REF], all make its analysis and modelling a challenging task. For instance, smiles themselves can be either simple or complicated things [START_REF] Paula M Niedenthal | The Simulation of Smiles (SIMS) model: Embodied simulation and the meaning of facial expression[END_REF]. They can be triggered by positive emotion, positive social motives, but also exploited to communicate and maintain social status. Further, coping with such challenge, might require to go beyond the current paradigm of the pattern recognition "pipeline" [START_REF] Sariyanidi | Automatic analysis of facial affect: A survey of registration, representation, and recognition[END_REF] and engage with embodied and simulation-based accounts [START_REF] Goldman | Simulationist models of face-based emotion recognition[END_REF][START_REF] Paula M Niedenthal | The Simulation of Smiles (SIMS) model: Embodied simulation and the meaning of facial expression[END_REF][START_REF] Vitale | Affective facial expression processing via simulation: A probabilistic model[END_REF]. Clearly, in such case, the availability of multimodal data involving physiological signals is mandatory. Eventually, positive emotions serve important social functions [START_REF] Barbara L Fredrickson | The value of positive emotions[END_REF] and there is evidence that human-computer interaction is natural and social too [START_REF] Reeves | The Media Equation: How People Treat Computers, Television, and New Media Like Real People and Places[END_REF]. We surmise that such emotions and amusement to a great extent, are likely to play a key role for a number of human-computer applications.

  Recorded signals in each dataset. In brackets the number of complete data. EEG = electroencephalography, ECG = electrocardiogram, BVP = blood volume pulse, EDA = electrodermal activity, R = respiration, T = temperature, G = eye gaze tracking. Extracted visual features in each dataset. FP = Fiducial Points, AU = Action Units

Figure 1 :

 1 Figure 1: Experimental setup (left) and detail of the physiological sensors used during the acquisition process (right).

Figure 2 :

 2 Figure 2: Examples of video-based acquired data: color camera frame, depth camera frame and 3D facial landmarks.

Figure 3 :Figure 4 :

 34 Figure 3: Visualization of raw data of the considered physiological signals (one session)

Figure 5 :

 5 Figure 5: An annotation session using DANTE tool.

Figure 6 :

 6 Figure 6: Distribution of the self-annotations in terms of pleasure, arousal and dominance for each of the four stimulus. The red filled points represent the mean value for each of the three dimensions.

Figure 7 :

 7 Figure 7: Histogram of frequencies for the self-annotations, in terms of pleasure and arousal for each of the four stimuli.

Figure 8 :

 8 Figure 8: Annotations of valence from 4 annotators and the corresponding gold standard in dashed green line. In brackets it is shown the evaluator score, namely the inter-rater agreement.

Figure 9 :

 9 Figure 9: Visualization of gold standard annotation in terms of valence, arousal and time, binned with 2 seconds overlapped moving windows.

Table 3 :

 3 Mean and standard deviation of the self report annotations for each stimulus and each emotional dimension.

	P		A		D
	µ	σ	µ	σ	µ	σ
	Stim 1 0.47 0.40 -0.31 0.80 -0.02 0.44
	Stim 2 0.63 0.50 0.49 0.62 0.54 0.53
	Stim 3 0.76 0.36 0.50 0.64 0.69 0.36
	Stim 4 0.67 0.39 0.47 0.59 0.64 0.43

https://youtu.be/-CZWXgPFj6A (from 1:17 to

2:13).[START_REF] Baltrušaitis | Cross-dataset learning and person-specific normalisation for automatic action unit detection[END_REF] https://youtu.be/Nwc1kRJdtYw (full length).[START_REF] Bilakhia | The MAHNOB Mimicry Database: A database of naturalistic human interactions[END_REF] https://youtu.be/C7xZYfxC8k4 (from 0:29 to end).[START_REF] Bradley | Emotion and motivation II: sex differences in picture processing[END_REF] https://youtu.be/tb_2TJjsq4g (from 4:53 to 6:50).

ACKNOWLEDGMENTS

This research was carried out as part of the project "Interpreting emotions: a computational tool integrating facial expressions and biosignals based shape analysis and bayesian networks", supported by the Italian Government, managed by MIUR, financed by the Future in Research Fund.


Nicolas Tollenaere

Auguste Olivry

Guillaume Iooss

Hugo Brunie

Albert Cohen

P Sadayappan

Fabrice Rastello

Efficient convolution optimisation by composing micro-kernels

come

Introduction

Tensor computations are at the core of many applications in scientific computing, data analytics and machine learning. Their optimized implementation is therefore of considerable interest. The current options for optimizing the implementation of a tensor operator, such as a 2D convolution, are:

• Polyhedral compilers like Diesel [START_REF] Elango | Diesel: Dsl for linear algebra and neural net computations on gpus[END_REF], Polly [START_REF] Grosser | Polly -performing polyhedral optimizations on a low-level intermediate representation[END_REF], Pluto [START_REF] Bondhugula | A practical automatic polyhedral program optimization system[END_REF], PPCG [START_REF] Verdoolaege | Polyhedral parallel code generation for cuda[END_REF], Tensor Comprehensions, [START_REF] Vasilache | Tensor comprehensions: Framework-agnostic high-performance machine learning abstractions[END_REF], Tiramisu [START_REF] Baghdadi | Tiramisu: A polyhedral compiler for expressing fast and portable code[END_REF] automatically generate multi-level tiled code for any affine loop computation such as 2D convolutions. However, a significant limitation is that none of them can directly optimize across tile sizes, which is critical for efficient CNN implementations. • Vendor libraries like oneDNN [START_REF]Intel. oneAPI deep neural network library (oneDNN)[END_REF] and cuDNN [START_REF]Cudnn: Gpu accelerated deep learning[END_REF] have been manually optimized by expert HPC and software engineers. While these implementations use JIT optimization, they cannot fully adapt to every given tensor extents of a CNN layer in a DNN pipeline.

• Autotuning can be performed by systems like AutoTVM [START_REF] Chen | Learning to optimize tensor programs[END_REF] or AutoScheduler [START_REF] Zheng | Ansor: Generating high-performance tensor programs for deep learning[END_REF] both part of the DNN specific compiler TVM [START_REF] Chen | TVM: An automated end-to-end optimizing compiler for deep learning[END_REF]. An expert user provides a multi-level tiled loop structure along with a specification of a search space including permutation among subsets of loops and parametric tile sizes. A search process guided by a dynamically constructed machine learning model [START_REF] Chen | Learning to optimize tensor programs[END_REF] iterates through tiled loop configurations, where code is generated, compiled and executed on the target platform. AutoTVM has been demonstrated to outperform polyhedral compilers [START_REF] Chen | Learning to optimize tensor programs[END_REF], and AutoScheduler to outperform AutoTVM [START_REF] Zheng | Ansor: Generating high-performance tensor programs for deep learning[END_REF], but both of them still do not match reference libraries like oneDNN. • Analytical modeling and optimization. Recent research has shown that a comprehensive characterization and optimization across all possible tiled loop configurations for CNNs is feasible [START_REF] Li | Analytical characterization and design space exploration for optimization of cnns[END_REF].

The approach is semi-automatic: manual reasoning to build analytical cost models for data movement and pruning, in conjunction with the automated resolution of nonlinear optimization problems to optimize tile sizes. The innermost loops use a manually created microkernel, similar to the BLIS one [START_REF] Field | Blis: A framework for rapidly instantiating blas functionality[END_REF] and vendor libraries. Over a set of 20+ CNN stages from 2 DNN pipelines, performance was shown to be consistently higher than state-of-the-art autotuning [START_REF] Chen | TVM: An automated end-to-end optimizing compiler for deep learning[END_REF] and comparable or better than the state-of-the-art oneDNN library [START_REF]Intel. oneAPI deep neural network library (oneDNN)[END_REF]. The best performing implementations in the state-of-the-art are based on a single microkernel. In this paper, we challenge that assumption, by claiming that limiting ourselves to a single microkernel for all problem sizes is too restrictive. Indeed, we show that there is a collection of well performing microkernels, and some of them might better fit the considered problem sizes.

We also claim that we should use combinations of microkernels of different sizes to cover exactly the whole space, as an alternative of partial tiles and the padding technique. This is especially critical for computation whose sizes are small: indeed the effect of padding or of a partial tile is significant on the performance. Therefore, combining microkernels allows us to have more consistant and performant implementation. In particular, such situation happens for the convolutional layers of neural networks.

Our main contribution is an end-to-end compiler flow for the optimization of tensor operations, with these properties: No partial tiles: It is possible to compose microkernels of different shapes to cover the iteration space, obviating the need for mixing full tiles (optimized) with partial tiles (suboptimal or unoptimized). Hybrid compilation strategy: We compose techniques with complementary strengths and weaknesses, such as domain-specific compilation, autotuning and analytical modeling, in order to use them effectively while compensating their limitations.

• To provide an example of the limitations of analytical modeling, some recent efforts [START_REF] Li | Analytical characterization and design space exploration for optimization of cnns[END_REF][START_REF] Olivry | IOOpt: Automatic derivation of I/O complexity bounds for affine programs[END_REF] developed analytical models for optimizing CNNs, establishing tight data-movement bounds in an idealized model of computation. However, their experimental results on real processors show that it is not as effective. A long-time barrier to accurate performance modeling is that architectural features like out-of-order execution or hardware prefetchers are virtually impossible to incorporate into a model for code optimization. • Another example of limitation, concerning autotuning techniques, is the explosion of the size of the search space: the full space of possible tiled loop configurations, choice of tile sizes and microkernels is enormous and only a tiny fraction is feasible to search. Autotuners resort to manually pruning the search-e.g. scheduling templates in AutoTVM [START_REF] Chen | TVM: An automated end-to-end optimizing compiler for deep learning[END_REF]-and/or restricting its expressiveness at the expense of missing high-performing variants-e.g. narrowing to a single microkernel. Our hybrid strategy overcomes these limitations by combining them: we use experimental evaluation of microkernels to prune the search space early, and we also use analytical modeling to guide autotuning. Decoupling of the search space: Both the modeling work of Olivry et al. [START_REF] Olivry | IOOpt: Automatic derivation of I/O complexity bounds for affine programs[END_REF] and the AutoTVM autotuning framework [START_REF] Chen | TVM: An automated end-to-end optimizing compiler for deep learning[END_REF] take a uniform or "flat" perspective of the tiled iteration space of a computation. In contrast, we use a "split" approach where a small set of innermost loops is decoupled from the remaining outer tile loops for the purpose of optimization. The inner set of loops constitute a microkernel that we separately optimize using extensive autotuning to generate a collection of high-performance variants. We then use a combination of analytical modeling and much smaller scale autotuning to optimize the outer tile loops.

While our approach applies to a wider class of tensor computations, we focus on CNNs with 2D convolutions in the evaluation. The main contributions of the paper are as follows:

• Instead of a single fixed (hand-coded) microkernel used by prior approaches [START_REF] Field | Blis: A framework for rapidly instantiating blas functionality[END_REF][START_REF] Li | Analytical characterization and design space exploration for optimization of cnns[END_REF], we conduct a comprehensive search among the space of possible microkernel configurations, identifying a set of microkernel instances as base building blocks for use in generating optimized code for every specific convolution instance. • For a specific convolution and microkernel, a judicious mix of automated analytical modeling along with limited autotuning allows to optimize the tile loops enclosing the microkernel. • For a specific convolution, a combination of several microkernels can be considered using imperfectly nested loops-or "beyond perfect" tile loops.

• Experiments with all 20+ CNN layers from 2 ML inference models (ResNet-18 and Yolo-9000) matches or outperforms the state-of-the-art oneDNN library [START_REF]Intel. oneAPI deep neural network library (oneDNN)[END_REF] and TVM [START_REF] Chen | TVM: An automated end-to-end optimizing compiler for deep learning[END_REF] autotuning frameworks. We observe that we are able to automatically match or outperform the performance of vendor libraries without copying and packing. This goes against the established methodologies [START_REF] Field | Blis: A framework for rapidly instantiating blas functionality[END_REF][START_REF]Intel. oneAPI deep neural network library (oneDNN)[END_REF] and raises open questions about the impact of this optimisation on performance. The rest of the paper is organized as follows: Sec. 2 provides a high-level overview of our approach. Sec. 3 shows the importance of the divisibility constraint, which motivates the choices made in the design of the search space in Sec. 4. Sec. 5 details the search strategies and Sec. 6 describes our code generator. Sec. 7 reports experimental comparisons against state of the art frameworks and libraries. Sec. 8 discusses related work before the conclusion in Sec. 9.

Overview of the Approach

Let us now present a high-level overview of our optimization approach. Its rationale derives from the following observations based on the advances made as well as the challenges faced by prior optimization efforts:

• Accurate fine-grained performance modeling of a multi-level tiled loop computation down to the innermost levels is extremely challenging and has not been achieved for complex computations like CNNs by any prior effort, but coarse-grained analytical modeling of the memory hierarchy has been shown to be effective, e.g., the work of Li et al [START_REF] Li | Analytical characterization and design space exploration for optimization of cnns[END_REF]. Their use of a microkernel made the analytical modeling coarse-grained and sufficiently accurate, but by using a fixed microkernel they limited the design space they explored. • Autotuning is an effective approach to overcome the challenges of accurate performance modeling.

The main blocker to its application to multi-level tiled loops is the enormous explosion of the optimization space and its effective traversal. Adding degrees of freedom with the ability to select among multiple microkernels, and to combine them, makes the problem even more challenging. AutoTVM does find good solutions but only by using expert-engineered scripts that limit the search space. • A key insight behind the optimization framework presented in this paper is that we can explore the for ( n = 0 ; n < N ; n+= 1 )

for ( k = 0 ; k < K ; k+= 1 ) for ( c = 0 ; c < C ; c+= 1 ) for ( h = 0 ; h < H ; h+= 1 ) for ( w = 0 ; w < W ; w+= 1 ) for ( r = 0 ; r < R ; r+= 1 ) for ( s = 0 ; s < S ; s+= 1 ) full design space effectively using a split two-level strategy: 1) Develop a collection of microkernels for a given target hardware platform via extensive one-time autotuning. The set of microkernels is agnostic to the actual sizes of CNN layers to be optimized. At this level, accurate performance modeling is infeasible, but the optimization space is quite manageable through autotuning via the execution of all potential variants of interest (detailed in Sec. 5.1). We then combine them in order to exactly fit a problem size without any partial tiles (detailed in Sec. 5.3).

O[k, h, w] = K[k, c, r, s] * I[c, h + r, w + s]
2) Use a combination of autotuning and analytical modeling to dramatically prune the space of possible tile loop permutations and degenerate loops (explained shortly). This is detailed in Sec. 5.2 to Sec. 5.4. We illustrate our integrated approach to optimized code generation for convolutions on Fig. 1 and Fig. 2. A 2D convolution is a 7-dimensional nested loop. Its optimized implementation requires multi-level tiling. Given a d-dimensional nested loop (d = 7 here), and a 5-level memory hierarchy (main-memory, L3, L2, L1 caches and registers), the total number of nested loops for tiling at all levels is 5d (35 here). This is illustrated in Fig. 2(a) as a set of outermost d tile-loops that step through L3-level tiles. Each L3-level tile has d tile-loops to step through a set of L2-level tiles, and so on, with the register-level tiles being marked as a microkernel. In practice, efficient tiled implementations will only have a small subset of active tile loops at a level, while the remaining ones are degenerate with a range of a single iteration and hence removed from the code. However, we cannot know a priori which tile-loops in a band are active versus degenerate. We first identify the microkernels that are performing well in isolation. Then we pick, among those, one or two combined microkernels that fit the considered problem size. Then, we use analytical modeling to identify the active loops in each band and the permutation within active loops in a band. This modeling is similar to other recent efforts on analytical model-driven optimization of CNN (e.g. the work of Olivry et al. [START_REF] Olivry | IOOpt: Automatic derivation of I/O complexity bounds for affine programs[END_REF] or that of Li et al. [START_REF] Li | Analytical characterization and design space exploration for optimization of cnns[END_REF]), so they will just be summarized in this paper. Finally, we use autotuning to search across the dramatically pruned space of outer-level tile configurations.

Fig. 2(b) shows the code generated by our optimization framework on one sample convolution for a target platform with a vector size of 16 elements. It uses two microkernels, one corresponding to a slice of the convolution iteration space with tile extents [H : 8, W : 1, C : 1, K : 2 × 16], and another with tile extents [H : 13, W : 1, C : 1, K : 2 × 16]. The L1-level tile (color coded blue) spans the full range of 128 iterations along C, which covers the full problem extent along C. An L2-level tile (color-coded red) spans a range of 8 + 2 × 13 = 34 along H and a range of 2 × 32 = 64 along K (which is the full problem extent). An L3-level tile (color-coded green) spans 4 × 34 = 136 along H, 136 along W. At this point the full problem extents have been covered and therefore the outer-most level of tiling loops (color-coded black in Fig. 2(a)) are degenerate. This example illustrates how combining two well-performing microkernels can be used to divide a problem size. It also shows that only a subset of the tile loops at any level are non-degenerate, thus demonstrates the importance of having an analytical model to prune the search space, before exploring it using autotuning.

Divisibility constraint and microkernels

In this section, we demonstrate the importance of combining microkernels instead of relying on (suboptimal) partial tiles. We consider the multiplication of very small matrices, such that the data footprint fits inside the L1 cache, and we measure performance for a continuous range of problem sizes.

If the microkernel sizes divide exactly the problem sizes, then it fits perfectly, and we observe a peak in performance. If the microkernel sizes does not divide exactly, the classical options are (i) to have a partial tile, smaller than the microkernel, that finishes the coverage of the iteration space; or (ii) to pad the space in order to continue using the microkernel one last time, at the cost of additional computation. In this paper, we take a third route: (iii) to combine two of the best performing microkernels to cover the space without partial tiles. The method to determine the best performing microkernel will be described in Section 5.1, and the selection algorithm is explained in Section 5.3.

Figure 3 compares the sequential performance of small matrix multiplication implementations, for problem sizes J = K = 128 and 8 ≤ I ≤ 49, on a Intel Xeon Gold 6230R CPU (Cascade Lake-SP, with AVX512). The performances are shown as percentage of the absolute peak performance, corresponding to the maximal utilisation of the two vectorized FMA units of the architecture.

MKL, Blis and libxsmm report the performance of these libraries. Notice the peak every 8 elements of I for MKL and a peak every 12 elements for BLIS. This gives us an indication about the size of their microkernel along the i dimension. Libxsmm also considers combination of microkernels, but restricted to predefined sizes such as multiples of 2 along the i dimension. Our experiment shows that this is not enough to obtain consistent performance for all problem sizes.

"Single microkernel, partial tile" is the performance of code generated by our framework, but only using the BLIS microkernel, with an unrolled partial tile. We observe a fluctuation of periodicity 12 in its performance. Notice that for values of I with a low modulo 12, the performances are worse than for the high modulo 12, because of the low performance of the partial tile. "Single microkernel, padded " is also the performance of the code generated by our framework, but using a padding strategy instead of a partial tile. We assumed that the padding overhead is free. As expected, the performance for low modulo are quite low, due to the significant additional amount of computation performed. However, this penalty decreased with the size of I.

Finally, "Combination of microkernels" corresponds to our microkernel combination strategy. The performances are more stable for any value of I.

This shows the importance of using all the microkernels available and to combine them, to avoid loss of performance due to padding or partial tiles. This is particularly important for some convolution benchmarks, such as Yolo9000, which have small problem sizes along most dimensions, which amplifies the penalty due to a partial tile, and which can have uncooperative divisors, (such as 34 = 2 × 17 for Yolo9000-12). Therefore, we build our optimisation space around this constraint, as described in Sec.4.

Defining the Optimization Space

In the following, a microkernel refers to an efficient region of code composed of a (large) basic block resulting from the full unrolling of innermost parallel loops, enclosed into zero or more perfectly nested reduction loops. It is generally written in assembly language or using vector intrinsics, aiming for the following objectives: (i) effective utilization of vector ALUs; (ii) effective reuse of data in (vector) registers across iterations through unrolling and register promotion; (iii) adequate Instruction-Level Parallelism (ILP) to hide the latency of pipelined functional units (multiply-and-add).

The iteration space is the set of integer vectors taken by the loop indices enclosing a given computational statement. Tiling [START_REF] Rivera | A comparison of compiler tiling algorithms[END_REF][START_REF] Coleman | Tile size selection using cache organization and data layout[END_REF] is a loop transformation that partitions the iteration space into sets, for ( it = 0 ; it < I ; it += 6 )

for ( jt = 0 ; jt < J ; jt += 32 ) for ( k = 0 ; k < K ; k += 1 ) µkernel gemm 6,32 (C, A, B, it, jt, k) called tiles and executed atomically. We only consider programs with rectangular iteration spaces, and rectangular tiling. Tiled code has additional loops compared to the original code: loops over tiles, and loops inside a tile. This partitioning allows us to control the amount of data accessed per tile, a.k.a. footprint, to make sure it does not exceed a given cache capacity. Fig. 4 shows a tiled matrix multiplication kernel as an illustrative example. It relies on an (inline) fully-unrolled and vectorized microkernel of size 6 × 32.

We use lowercase for problem dimensions (i, j, k), i.e. loop iterations, and uppercase to name the (possibly symbolic) upper bound on each dimension (resp. I, J, K), a.k.a. problem size. We also assume that any dimension is either parallel-i and j-or a reduction-k and all dimensions are permutable (loop interchange). While associativity can be used to parallelize a reduction, we do not exploit it.

In the class of computations we consider, a tensor may be accessed multiple times but always with the same subscript expressions, which are affine functions of surrounding loop iterators. For example, tensor A of shape {i, k | 0 ≤ i < I, 0 ≤ k < K} may be subscripted by [i, k], corresponding to the access function (i, j, k → i, k). We also assume that a loop index cannot appear twice inside an access function: for example E[i, i] is forbidden. These conditions are satisfied by all tensor contractions and convolutions, including strided variants.

As mentioned previously in Sec. 3, high-performance libraries, such as BLIS, TCCG, oneDNN, rely on the use of a single microkernel with some fixed tile sizes within the microkernel, e.g., 6 and 32 in the example of Fig. 4. When tile sizes do not divide tensor shapes, the traditional approach involves conditional execution or padding to manage partial tiles. We consider a broader optimization space, using a collection of microkernels so that their combination eliminates the need for partial tiles. We relax the divisibility constraint that must be satisfied in order to avoid partial tiles, enabling the ability to compose multiple, fully-optimized microkernels.

Our code generator is driven by a so-called optimization scheme. Conceptually, it can be seen as as a specialized abstraction, higher level than TVM schedules. An optimization scheme is a list of specifiers that describe the layered structure of the generated code, from the outermost loop inwards:

• [START_REF] Baghdadi | Tiramisu: A polyhedral compiler for expressing fast and portable code[END_REF][START_REF] Elango | Diesel: Dsl for linear algebra and neural net computations on gpus[END_REF], [START_REF] Abadi | Tensorflow: A system for large-scale machine learning[END_REF][START_REF] Feautrier | Some efficient solutions to the affine scheduling problem: I. one-dimensional time[END_REF]] fulfills the divisibility constraint (no partial tiles) while involving high-performance microkernels of size 11 and 12 along y.

Example The naive implementation of a matrix multiplication would be represented as [R i , R j , R k ]. An implementation for higher performance, based on the BLIS [START_REF] Field | Blis: A framework for rapidly instantiating blas functionality[END_REF] microkernel for floats (f32

) on AVX2 is: [R j , R k , R i , T nc 16 ,j , T mc 6 ,i , T n k ,k , U 6,i , U 2,j , V j ]
The generated code contains a microkernel of size (i = 6, j = 16, k = n k ) known to be quite efficient as it requires only 15 vector registers and exposes enough ILP (12 independent multiply-add instructions issued between two accumulation steps) [START_REF] Field | Blis: A framework for rapidly instantiating blas functionality[END_REF]. Above it, loops i and j induce a 2D tile of size (m c , n c ). One may immediately notice that this approach assumes that I is a multiple of m c , itself being a multiple of 6 (similar constraints apply for j and k). State of the art libraries rely on fixed-size microkernels and tuned tiles sizes, and thus introduce partial non-optimized tiles to cope with arbitrary problem sizes that do not fulfill such a divisibility constraint. Assume for example a matrix-multiplication of size I × J × K = 128 × 128 × 64. 128 is not divisible by 6, but 128 = 12 × 6 + 8 × 7, and efficient code can be obtained using the following scheme: [R j , λ seq i α. [ [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem: I. one-dimensional time[END_REF][START_REF] Bondhugula | A practical automatic polyhedral program optimization system[END_REF], (8, 7)] , T n k ,k , U α,i , U 2,j , V j ] which leads to the loop structure shown in Fig. 5.

Optimization space exploration

We propose a novel optimization algorithm capable of exploring the expressive optimization space introduced in the previous section. This algorithm uses a combination of analytical modeling and autotuning; it is split into the offline optimization of the microkernel (register level), before specializing on a particular convolution and deciding on a tiling structure (cache level). Among all the possible combinations for a tiled convolution, we apply the following pruning strategy:

• we only use the best performing microkernels: the performance of this portion of the code limits the performance of the whole application, so we need to ensure that it runs as efficiently as possible;

• we forbid partial tiles, which is a source of slowdown, especially on the innermost levels of the generated code; because the divisibility constraint would be too strict for some problem sizes, we consider combinations of microkernels. To speed up the search, we use a metric that focuses on the variations which maximize the sizes of the reduction loop above the microkernel, and the operational intensity of the computation.

Finally, we consider multiple parallelisation strategies, allowing to collapse and parallelize across multiple non-reduction loops.

The full algorithm is summarized by Fig. 6:

Step 1. We measure the performance of many microkernels in isolation. The set of microkernel candidates is formed of the best-performing ones (Sec. 5.1). This phase is problem size agnostic but specific to each target architecture.

Step 2. For each microkernel candidate, accounting for the problem size and cache sizes, we determine the best loop permutation (list of loop dimensions) enclosing the microkernel, that is, the main structure of the loop nest(which dimensions are tiled at each level). Tile sizes are not determined at this point. This permutation is obtained through operational research, using an analytical model of the footprint, data movement and reuse across tiles (Sec. 5.2).

Step 3. From the selected microkernels and loop permutations, we generate the space of optimization schemes (Sec. 5.3). In general, we require any tile size picked at a given dimension to be a multiple of the size of its sub-tiles and to divide the (full) problem size along that dimension. Except for microkernels where such a restriction would be impractical: e.g. Yolo9000-8 has a problem size H = 17, which is a prime number and all microkernels of size 17 have terrible performance. We leverage the λ seq specifier to compose microkernels of size 5 and 6 to solve the problem: microkernels of size {µS = 3, µH = 5, µK = 4} and {µS = 3, µH = 6, µK = 4} can be combined by sequencing two iterations of the latter before one iteration of the former, fulfilling the divisibility constraint without compromising performance (Sec. 5.3).

Step 4. At this point, the search space is already tractable for autotuning (a few hours per convolution), but we can do better. A simple metric (detailed in Sec. 5.4) sorts all the resulting schemes, from which we can retrieve any number of candidate implementations. We pick the top 200.

Step 5. Optionally, we may integrate the top-performing scheme as a sub-tree in a TVM program. We use this ability to facilitate performance comparisons, and to automatically upgrade the optimal sequential scheme into a parallel implementation (see Sec. 5.5).

Microkernel definition and evaluation

We identified 4 different unrolling schemes to form efficient microkernels for a 2D convolution. Each one is suitable for a different class of tensor shape and convolution operator. In all 4 schemes, dimension k is selected for vectorization because it contains the simplest access pattern among w, h and k. Unrolling takes place along k and h, and optionally along the convolution kernel (stencil) dimensions r and s to form the 4 schemes:

• [T 512,c , U β,h , U α,k , V k ] • [T 512,c , U 3,s , U β,h , U α,k , V k ] • [T 512,c , U 3,r , U β,h , U α,k , V k ] • [T 512,c , U 3,r , U 3,s , U β,h , U α,k , V k ]
where α and β are the sizes of the microkernel along the k and the h dimension, respectively.

To evaluate performance, we repeat the resulting unrolled basic block many times along the c dimension (T 512,c ) and run the microkernel on a matching problem size. The results for the first scheme on We observe that the graph is roughly convex with some local fluctuations. Many microkernels are near-optimal for a given unrolling scheme. We select the ones above 85% optimal and partition them into classes of fixed sizes for C, S, R, and where H belongs to an interval defining the class. For example, {[U β,h , U 2,k , V k ], 8 ≤ β < 15} is one of the class of microkernels that is selected for AVX-512, as shown with the leftmost red vertical rectangular contour on Fig. 7.

This step is problem size agnostic and needs to be done only once per target architecture.

Loop permutation above the microkernel

The next step is to find a suitable permutation of tiling loops, for fixed values of parameters and for each candidate microkernel. For example, for Yolo9000 layer 0, the permutation that is chosen for microkernel [U 8,h , U 2,k , V k ] is (from outer to inner): K, H, W, H, H, S, R, W, C . This permutation is found using an analytical model of the data movement across all cache levels, similar to the one taken in recent work by Li et al. [START_REF] Li | Analytical characterization and design space exploration for optimization of cnns[END_REF] and Olivry et al. [START_REF] Olivry | IOOpt: Automatic derivation of I/O complexity bounds for affine programs[END_REF]. For a given permutation, we derive an analytical expression of the data movement volume at each cache level, as a function of tiling loop extents and cache sizes. This is done by computing the footprint of each array at each level of the loop nest, as well as the level at which the total memory footprint exceeds the cache size. First, a pre-processing step prunes the space of all possible loop permutations. Indeed, we have (7!) loop permutation for each of the 4 remaining level of memory, and many permutations are equivalent in term of data movement cost, or can be shown to have worse data reuse than others. Following the methodology proposed by Olivry et al. [START_REF] Olivry | IOOpt: Automatic derivation of I/O complexity bounds for affine programs[END_REF], we manage to reduce the number of considered permutations from 7! to only 6 (for each level of the memory hierarchy).

Then, for given values of array sizes, the analytical model is fed to a non-linear problem optimizer which select the best permutation overall. The result is actually richer than that, as the solver also produces (generally non-integral) tile sizes that minimize the overall data movement [START_REF] Olivry | IOOpt: Automatic derivation of I/O complexity bounds for affine programs[END_REF]. Yet we observed that these tile sizes are not necessarily useful, as finer-grained performance considerations come into play when getting closer to the peak performance of the machine. We thus only retain the loop permutation and delay tile size selection to a later autotuning step. Interestingly, we observe that the permutation remains stable when varying the unroll factor of microkernels along the h dimension. This is easily understood, considering that the loop permutation mostly relates to cache-level optimization while microkernels operate at the register level. As a result, we use the same permutation and only perform the analysis once per unrolling factor along k, r and s. As a fortunate side-effect, this means that when combining two microkernels, since they only differ in their unrolling factor along h, we can be sure that the loop order above both of them is the same.

Space of the valid optimization schemes

Microkernels and combination We consider each class of microkernels and select those whose sizes divide the problem sizes. Then, we look for the combination of two microkernels differing only along the h dimension that allow to cover the size of the h dimension. For all pairs of microkernel in the same class, of sizes h 1 and h 2 , and given a problem size H, we look for a number of repetitions a and b of these microkernels such that: (a × h 1 + b × h 2 ) divides H. For example, if H = 34 and considering two microkernels of the same class of size h 1 = 11 and h 2 = 12, one may combine two microkernels of size 11 followed by a microkernel of size 12, for a total of 34.

If no single microkernel or combination of microkernels would be found with this process, the fallback would be to use a suboptimal microkernel, which is what we wanted to avoid by ruling out partial tiles. Fortunately, this situation never happens on the AVX512 microarchitectures we considered. Indeed, the classes of microkernels identified in Section 5.1 are large enough to accommodate for any possible size through the combination of two microkernels (as long as the size of the convolution dimension of interest is greater than or equal to the smallest microkernel in the selected class). For example, for the microkernel class {[U β,h , U 2,k , V k ], 8 ≤ β < 16}, assuming that the problem size along the K dimension is even, all problem sizes H above 8 can be obtained by a linear combination of two integral elements h 1 and h 2 from the interval [START_REF] Chen | TVM: An automated end-to-end optimizing compiler for deep learning[END_REF][START_REF] Goto | High-performance implementation of the level-3 blas[END_REF] (e.g. 17 = 8 + 9).

Completing the scheme Given (i) a single microkernel that divides the problem sizes, or a microkernel combination that divides the problem sizes, and (ii) a loop permutation, completion into an optimization scheme proceeds as follows:

• Set the optimization scheme (the list of specifiers) to the one corresponding to the chosen microkernel or combination of microkernels. • For each dimension, consider the divisors of the problem size divided by the microkernel size. These divisors need to be allocated to the different strip-mined occurrences of the dimension across the whole permutation. There are multiple solutions and the algorithm considers all of them. Notice that there are (by construction) at most 4 occurrences of a dimension in a loop permutation, which limits the amount of possibilities. • For each element d in the loop permutation, consider the product π of the divisors allocated to this occurrence and insert the corresponding specifier T π,d to the left of the current scheme. • When considering a sequence of two microkernels at dimension d with the combination a×h 1 +b×h 2 , insert λ seq d α.

[l] at any occurrence of dimension d in the scheme, and consider all possibilities of placement of this insertion. The value of the list is [(a, h 1 ), (b, h 2 )].

Example Consider the range of microkernels:

{[U β,h , U 2,k , V k ], 8 ≤ β ≤ 15}
for an AVX512 architecture, and the Yolo9000-13 problem sizes (K, C, H/W, R/S) = (512, 256, 34, 3). Assume that the corresponding loop permutation found was:

[[K, H], [W, H], [H], [S, R, W, C]]
The problem size H on dimension h is 34 = 2×17, hence there is no single microkernel from the considered class that matches one of its divisors. Next, we consider combinations of 2 microkernels from that class; 2 × 11 + 12 is one such combination. Now, we need to distribute the multiples of the other dimensions across the different levels of tiling described by the loop permutation:

• The k dimension is trivial: there is only one loop above the microkernel and the microkernel size along this dimension has a footprint of 32. Thus, we need to tile by a factor of 16 on the outer loop to reach 512. • The c dimension is also trivial: the only loop need to be tiled by a factor of 256. Likewise for r and s dimensions, the tiling factor should be 3 for both. • The h dimension is already managed by the combination of microkernels. We have 3 locations where we can place the λ seq h β. [ [START_REF] Baghdadi | Tiramisu: A polyhedral compiler for expressing fast and portable code[END_REF][START_REF] Elango | Diesel: Dsl for linear algebra and neural net computations on gpus[END_REF], [START_REF] Abadi | Tensorflow: A system for large-scale machine learning[END_REF][START_REF] Feautrier | Some efficient solutions to the affine scheduling problem: I. one-dimensional time[END_REF]] compositon of two microkernels. Let us consider the outer one. • The w dimension has 34 to be distributed across 2 level of tiling. There are 4 combinations: 34 × 1, 17 × 2, 2 × 17 and 1 × 34. Let us consider the second one. The resulting optimization scheme (among many) is:

[T 16,k , T 1,h , λ seq h β. [(2, 11), (1, 12)] , T 17,w , T 1,h , T 1,h , T 3,s , T 3,r , T 2,w , T 256,c , U β,h , U 2,k , V k ]

Further pruning of the optimization space

At this stage, we have a space of sequential scheme candidates, whose size range from a few hundreds to tens of thousands of potential schemes, depending on both the size of the problem and the architecture. Our objective is to prune this space by combining two criteria. The first criterion is to maximize the Operational Intensity, that is, the number of operations divided by the volume of data movement. The second criterion is to reduce the overhead of control flow, control misprediction, and data misprefetching. In particular, the larger the size of the reduction loop above the microkernel (along dimension c) is, the better the scheme is in relation to this criterion.

So, we propose the following metric, to focus on a region where there is at least one of the best performing schemes:

• First, select 40% of the schemes that have the largest reduction sizes (on c) above the microkernel.

• Then, sort these schemes according to their volume of data movement. The smallest volume is the best candidate, and we select up to 200 candidates.

Integration and parallelization with TVM

So far, we focused on optimising sequential performance. We have interfaced our tool with TVM to facilitate comparisons with the state of the art and to generate parallel code. Because TVM does not have currently a notion of microkernel, we use its tensorize construct to import a generated C code (see Sec. 6) corresponding to the innermost loops inside TVM's Python intermediate representation of a convolution. These innermost loops are the ones up to dimensions h and w, which are usually part of the L1 tiling level. By construction, these loops will include the microkernel and some of the L1-resident loops, including the innermost reduction loop (on c) enclosing the microkernel. This allows us to express completely a sequential scheme in the TVM framework. Then, we consider the loops above the tensorized section, and we consider several strategies to introduce parallelism:

• We consider the biggest band of non-reduction loops (to have enough iterations), collapse it and parallelize the resulting loop. • We regroup the reduction loops at the innermost level and collapse/parallelize the rest of the loops.

• We regroup the reduction loops at the outermost level and collapse/parallelize the rest of the loops. We add these strategies as a part of the optimization space, and consider all variations when evaluating the performance of the candidates.

Code generation

We now describe how to generate C code from a computation specification, problem size and the associated optimization scheme. Generating a loop requires to know the size of the sub-tiles, so our code generator proceeds from innermost outwards. Calling a sub-scheme the suffix of an optimization scheme, at a given step the already generated code (that corresponds to inner levels) is fully specified by the corresponding sub-scheme. In the following, the size of a sub-scheme refers to the size of the corresponding (parameterized) sub-iteration space. Taking the example from Sec. 4, the sub-scheme of the BLIS microkernel (including the reduction loop on k) is:

S µkernel = [T n k ,k , U 6,i , U 2,j , V j ].
Its size along i, j, and k is respectively 6, 16 and n k .

Overview of the code generation algorithm

Our code generator iterates right to left on the optimization scheme in a single pass. At every level, we keep track of the following information: (i) the size of the loops that are already generated; (ii) for each dimension, the name of the last index used by a for loop (to handle tiling).

Before applying our code generation algorithm, we apply a preprocessing step to get rid of the λ seq specifier and its parameter α. We introduce a new specifier Seq that corresponds to the sequential composition of a list of strategies. In our case, the list of the λ seq specifier is always of size 2. The corresponding rewriting rule is:

[ . . . , λ seq d α. [(i 1 , v 1 ), (i 2 , v 2 )] , S ] ⇒ [ . . . , Seq([T i1,d , S[α/v 1 ]], [T i2,d , S[α/v 2 ]]) ]
where S is the sub-scheme following the λ seq specifier and S[α/v] is this sub-scheme where α was substituted by the value v. We now have a tree of specifiers instead of a list of specifiers, on which we can still iterate from the leaves (innermost loops) to the root of the tree (outermost loops).

Code generation rules

Let us now survey the different specifiers and how code generation operates for each one:

• Sequence Seq: Combine sequentially the generated code corresponding to the sub-schemes inside the Seq. • Vectorization V d : Considering the definition of the computation described in Sec. 4, one may determine which operations should be vectorized by traversing the graph starting from the loads:

read (T, f ) is vectorized if d appears in the access function f .

-Op (x, y) is vectorized if one of its operands (x or y) is vectorized. If one of them is a scalar, it is broadcasted. -write (v, T, f ) is vectorized if v is a vector and d appears in the access function f . These conditions must be both true or false, else this is an error. The C code uses Intel intrinsics to manipulate vectors.

• Unroll U k,d : Unroll the computation over the d dimension k times by duplicating the generated code of its sub-scheme, while updating the value of the loop index on the d dimension in each duplication. • Tiling T k,d or R d : Add a loop over the generated code of its sub-scheme that iterates k times, and whose value is increased by the value of the sub-scheme. In the case of R d , one may deduce the correct number of iterations by comparing the size of the sub-scheme with the problem sizes. This changes the current loop index in use over the d dimension.

Performance results

In this section, we evaluate the performance of the code generated by TTile with several state-of-theart approaches: oneDNN [START_REF]Intel. oneAPI deep neural network library (oneDNN)[END_REF] (Intel library, V2.3), AutoTVM [START_REF] Chen | TVM: An automated end-to-end optimizing compiler for deep learning[END_REF] (autotuning, v0.8.dev0 of June 2021), Autoscheduler [START_REF] Zheng | Ansor: Generating high-performance tensor programs for deep learning[END_REF] (Ansor, autotuning, in TVM) and Mopt [START_REF] Li | Analytical characterization and design space exploration for optimization of cnns[END_REF] (analytical modeling).

Setup

The experiments were carried out on two architectures: (a) a 18-core Intel Xeon Gold 5220 Cascade Lake processor (frequency set to 2.2 GHz) with with 1 socket and 1 AVX512 fused multiply-add unit per core; and (b) a 32-core Intel Xeon Gold 6130 Skylake processor (frequency set to 2.3 GHz) 2 The OS is Debian GNU/Linux with kernel version 4.19, monitoring hardware counters using PAPI version 5.7.0. All codes were compiled using gcc with the flags -O3 -march=native -fno-align-loops. For AutoTVM and AutoScheduler, we used the recommended template conv2d NCHWc from the TVM library. AutoTVM's machine learning model ran over 1000 trials in order to find its best configuration.

We evaluate performance over the convolutions of two networks: Yolo-9000 [START_REF] Redmon | YOLO9000: better, faster, stronger[END_REF] and ResNet-18 [START_REF] He | Deep residual learning for image recognition[END_REF].The sizes of their convolution layers can be found in Fig. 8. To measure performance with AutoTVM and TVM+TTile, we used the TVM function evaluator, with parameters repeat and number set to 1. Runs of oneDNN use the input tensor layout nchw. Every convolution is run 20 times on 20 different copies of the tensors (to avoid reusing data in the last-level cache for repeated experiments). We eliminated the 5 best and the 5 worst results and take the median of the 10 remaining results. We also consider a cold cache and we flush it between every run.

Performance measurement Fig. 9 presents the performance results, reported as a fraction of the peak performance of the multicore CPU. A few MOpt results are missing due to reproducibility issues with the author's artifact. We reported such situations with a performance at 0% of the machine peak. We also report the geometric mean for each tool and in the sequential/parallel cases, across all the available benchmarks. We notice that TVM+Ttile is well above the other state-of-the-art tools in average, which demonstrates the viability of our approach. The performance results for oneDNN seem surprisingly low, despite our efforts to explore the relevant configuration settings for oneDNN. However, we have confirmed with the developers that these performance were coherent with their expectation.

In term of total compilation time to find the best implementation found, over all benchmarks, AutoTVM and AutoScheduler both take around 5h. MOpt takes around 1h30 and TTile takes around 2h. OneDNN is performing JIT compilation thus does not have an offline exploration phase.

Microkernel variations

We report in Fig. 10 the microkernel used by our best performing schemes. We notice a variety of 16 different microkernels used for 23 different convolutions. For example, the size of the microkernel picked for ResNet18-1 is "[14H2K]", i.e., 14 unrolled iterations along the h dimension and 2 iterations along the k dimension. In comparison, Yolo9000-4 uses a combination of two microkernels "[3R10H] + 2[3R12H]". Both microkernels are of size 3 along the r dimension, but differ in the size along the h dimension. The second microkernel is also used twice as often than the first microkernel. As another example illustrating this result, the best scheme found for Yolo9000-5 was introduced as the example of Sec. 2.

Interpretation We can draw several observations from these experiments. First, the divisibility constraint is the key restriction of our optimization space, that we have slightly relaxed with the possibility to combine microkernels. Despite this restricted, we still manage to find an implementation with excellent performance. This confirms that the divisibility constraint is effective at aggressively pruning implementations (including good ones) while preserving enough of the best-performing ones.

Also, the diversity of the microkernels available and their combination is the key point that differentiates our optimization space with those of the state-of-the-art tools. The fact (i) that we managed to find better implementations, sometimes by a factor of 2, and (ii) the variability of microkernels used by the best implementation we found, support the claims made in Sec. 3. In order to reach the best performance, we should not restrict ourselves to a single microkernel. Also, combining microkernels happens to outperform competing alternatives such as partial tiles (limited or no unrolling when hitting tensor boundaries) or padding.

As a final observation, copying and packing the footprint of register tiles in microkernels is generally assumed to be essential to performance (on both tensor contractions and convolutions) [START_REF] Field | Blis: A framework for rapidly instantiating blas functionality[END_REF][START_REF]Intel. oneAPI deep neural network library (oneDNN)[END_REF]. Our results challenge this assumption by obtaining an efficient implementation without considering this transformation in our search space. This shows that this transformation should not be applied systematically for all convolution kernels. We plan to study when the packing transformation should be applied in our future work.

Related work

Optimization of affine programs To optimize affine programs, some methods are based on analytical models and operation research. This is the main approach used by polyhedral based compilers [START_REF] Grosser | Polly -performing polyhedral optimizations on a low-level intermediate representation[END_REF][START_REF] Bondhugula | A practical automatic polyhedral program optimization system[END_REF][START_REF] Verdoolaege | Polyhedral parallel code generation for cuda[END_REF][START_REF] Vasilache | Tensor comprehensions: Framework-agnostic high-performance machine learning abstractions[END_REF][START_REF] Elango | Diesel: Dsl for linear algebra and neural net computations on gpus[END_REF][START_REF] Baghdadi | Tiramisu: A polyhedral compiler for expressing fast and portable code[END_REF] that leverage parametric integer linear programming. Although such approaches are well suited to expose parallelism [START_REF] Feautrier | Some efficient solutions to the affine scheduling problem: I. one-dimensional time[END_REF][START_REF] Feautrier | Some efficient solutions to the affine scheduling problem. part ii. multidimensional time[END_REF] and coarse grain locality [START_REF] Bondhugula | A practical automatic polyhedral program optimization system[END_REF], we believe it may not be the right formalism for tile size selection or register level optimizations.

On the other hand, the ability to count points in a polyhedra [START_REF] Barvinok | A polynomial time algorithm for counting integral points in polyhedra when the dimension is fixed[END_REF] allows to automatically generate (non-linear) cost models. We use the Barvinok [3] library for this purpose, generalizing the approach of Li et al. [START_REF] Li | Analytical cache modeling and tilesize optimization for tensor contractions[END_REF] for the selection of a permutation scheme.

Cloog [START_REF] Bastoul | Code generation in the polyhedral model is easier than you think[END_REF] is a powerful algorithm to automatically generate imperative code for scanning a union of polyhedra. Polyhedral compilers leverage such code generation capabilities but face the challenge of dealing with a very general class of imperfect nests and transformations. It is difficult in such a broad context to compete with domain-specific optimizations. Our code generator involves simple polyhedron scanning algorithms, and the divisibility constraint allows to generate high-quality compiler friendly code without heroic efforts [START_REF] Grosser | Polyhedral ast generation is more than scanning polyhedra[END_REF].

Optimization of machine learning programs There exist many compilers specialized for machine learning: PlaidML [START_REF] Chen | Plaidml-he: Acceleration of deep learning kernels to compute on encrypted data[END_REF] using polyhedral techniques, XLA [START_REF] Google | Xla : optimiser le compilateur pour le machine learning[END_REF] for TensorFlow [START_REF] Abadi | Tensorflow: A system for large-scale machine learning[END_REF], Halide [START_REF] Ragan-Kelley | Halide: a language and compiler for optimizing parallelism, locality, and recomputation in image processing pipelines[END_REF], or TVM [START_REF] Chen | TVM: An automated end-to-end optimizing compiler for deep learning[END_REF]. TVM, as opposed to most approaches does not rely on numerical libraries. Its strategy is to select the best schedule using autotuning with an ML-based performance model. Contrary to our approach that decouples the search into microkernel optimization and loop tiling/permutation search, the TVM search space is flat. In TVM, optimizations related to strength reduction and register tiling are left to the compiler. TVM has been extended with FlexTensor [START_REF] Zheng | Flextensor: An automatic schedule exploration and optimization framework for tensor computation on heterogeneous system[END_REF] and Ansor (a.k.a AutoScheduler) [START_REF] Zheng | Ansor: Generating high-performance tensor programs for deep learning[END_REF]. We also compare to AutoTVM, the auto optimiser of TVM. Telamon [START_REF] Beaugnon | Optimization space pruning without regrets[END_REF] tackles this problem by building a very large, flat search space where optimization choices are tied together by dependency constraints. Then the exploration combines an elaborate performance model to prune the search space with feedback from actual executions.

Linear algebra and CNN libraries Frameworks such as TBLIS [START_REF] Matthews | High-performance tensor contraction without transposition[END_REF] or TCCG [START_REF] Springer | Design of a high-performance GEMM-like Tensor-Tensor Multiplication[END_REF] aim at creating portable optimized code for BLAS or tensor contraction kernels. These frameworks implement an efficient predefined scheduling scheme which is very effective, in particular for matrix multiplication [START_REF] Goto | High-performance implementation of the level-3 blas[END_REF]. These frameworks take advantage of advanced optimizations: tensor transposition, tensor blocking or subviewing, data prefetching, vectorization, block scheduling, unrolling and register promotion. The register tile shape is predefined using expert knowledge on instruction level and register pressure. Thanks to aggressive autotuning and JIT/AoT code versioning, MKL [START_REF] Wang | Intel Math Kernel Library[END_REF] and oneDNN [START_REF]Intel. oneAPI deep neural network library (oneDNN)[END_REF] are the best available Intel libraries which implement all these techniques.

Conclusion

We presented TTile, a tool which generates high-performance code for tensor computation, and we evaluated it on static-shape 2D convolutions. The generated code exploits a large set of performant microkernels and combines them in order to avoid partial tiles. The methodology itself is a combination of different kind of state-of-the-art methods, including analytical modeling, autotuning and experimentation. We present performance results that demonstrate higher performance than oneDNN in the vast majority of cases, and consistently higher performance than the state-of-the-art TVM tensor compiler and autotuner and MOpt. We plan to release publicly a version of our tool in the next few months.

Future work includes extension of the search space to include additional transformations such as packing or prefetching, in order to carefully study their impact on the performance and determine when they are needed. We also plan to generalize the method to a wider set of tensor operators and layout configurations, using MLIR [START_REF] Lattner | Mlir: Scaling compiler infrastructure for domain specific computation[END_REF]. This will also facilitate the exploration of interactions with graph-level optimizations, leveraging optimization opportunities currently accessible only to JIT compilers such as XLA.

Figure 1 :

 1 Figure 1: 2D Convolution (unit stride).

••

  Once per HW platform • Agnostic to layer specification { Per-convolution autotuning (a) Multi-level tiling and split two-level optimization.

Figure 2 :

 2 Figure 2: Code generation sketch using microkernel composition. Convolution sizes are K = 64, C = 128, H = W = 136 and R = S = 1. Note that 136 = (8 + 13 × 2) × 4. Loop colors in (b) match the cache level they fit into.

Figure 3 :

 3 Figure 3: Performance of small matrix multiplication kernels, for J = K = 128 and 8 ≤ I ≤ 50.

Figure 4 :

 4 Figure 4: Tiled sgemm with microkernel.

Figure 5 :

 5 Figure 5: Microkernel composition example.

Figure 6 :

 6 Figure 6: Flow of the optimization algorithm.

Figure 7 :

 7 Figure 7: Performance of microkernels in isolation for AVX512 in percentage of the machine peak, for R = S = 1. Microkernel sizes-α along the k dimension (horizontal axis) and β along the h dimension (vertical axis)-vary between 1 and 15. Only the upper-left triangle was evaluated. Red-bordered microkernels are the ones selected (offline) for our algorithm.

Figure 8 :

 8 Figure 8: Convolution benchmarks and sizes. The kernels marked with a * are stride 2, else stride 1. Dimension k of Yolo9000-23 was padded to 28272 (a multiple of 16) to vectorize it on AVX512.

Figure 9 :

 9 Figure 9: Comparison with AutoTVM, AutoScheduler, oneDNN, Mopt for AVX512 (Intel Xeon Gold 5220 and 6130), shown as percentage of machine peak. The left side of each line is the sequential performance and the right side is the parallel performance for the same convolution sizes. The last entry is the geometric mean for each tool, for the sequential and parallel case.

Figure 10 :

 10 Figure10: Microkernels of the best parallel scheme found by TTile on an Intel Xeon Gold 6130. The corresponding unrolling factors and split factors of (combined) microkernels are reported with a compact notation where "2[3R11H]+[3R12H]" stands for the combination of 2 microkernels: the first is used twice as much as the second; "[3R11H]" stands for a microkernel where dimension r is unrolled 3 times and h is unrolled 11 times.