Ricky Walsh
email: ricky.walsh@inria.fr

Cédric Meurée
email: cedric.meuree@inria.fr

Anne Kerbrat

Arthur Masson

Burhan Rashid Hussein

Malo Gaubert

Francesca Galassi

Benoit Combès

Expert Variability and Deep Learning Performance in Spinal Cord Lesion Segmentation for Multiple Sclerosis Patients

Keywords: Multiple sclerosis, spinal cord, magnetic resonance imaging, lesion segmentation, inter-rater variability, intrarater variability, deep learning, automated segmentation

Multiple sclerosis (MS) patients often present with lesions in spinal cord magnetic resonance (MR) volumes. However, accurately detecting these lesions is challenging and prone to inter-and intra-rater variability. Deep learning-based methods have the potential to aid clinicians in detecting and segmenting MS lesions, but can also be affected by rater variability.

This study assesses the inter-and intra-rater variability in manual segmentation of spinal cord lesions, and evaluates raters and a state-of-the-art nnU-Net model against a ground truth (GT) segmentation of a senior expert. Four experts segmented twelve spinal cord MR volumes from six patients twice, at a time distance of two weeks.

Considerable inter-and intra-rater variability were observed, with the total number of detected lesions ranging from 28 to 60, depending on the rater. Moreover, the segmented volumes of individual lesions varied substantially between raters. All raters and the model achieved high precision when evaluated against the senior expert GT, but sensitivity was notably lower. These results motivate the need for more sensitive automated methods to aid clinicians in lesion detection, and suggest that consideration should be given to inter-rater variability when training and evaluating automated methods.

I. INTRODUCTION

Multiple Sclerosis (MS) is a chronic inflammatory demyelinating disease of the central nervous system [START_REF] Kutzelnigg | Chapter 2 -Pathology of multiple sclerosis and related inflammatory demyelinating diseases[END_REF]. Magnetic Resonance (MR) imaging of the brain and spinal cord is used to identify MS lesions, providing an essential aid to clinicians in the diagnosis and prognosis of this disease [START_REF] Thompson | Diagnosis of multiple sclerosis: 2017 revisions of the McDonald criteria[END_REF]. In particular, MS lesions in the spinal cord are prognostic of short-term disability progression [START_REF] Leguy | Prognostic value of spinal cord MRI in multiple sclerosis patients[END_REF] and thus can indicate the need for more aggressive treatments.

The identification of MS lesions on a given MR image is a complex and mentally demanding task and can lead to an underestimation of disease activity, even for experienced radiologists and neurologists. The quality of a segmentation is influenced, among other factors, by the characteristics of the lesions, the experience of the clinician examining the image, the amount of information available (e.g., the number of scans or clinical information), and image quality. The quality of MR images of the spinal cord is particularly affected by artifacts from patient movement or respiratory and cardiac cycles [START_REF] Kearney | Spinal cord MRI in multiple sclerosis-diagnostic, prognostic and clinical value[END_REF], which can introduce uncertainty as to whether certain hyper-intense regions in the image are truly lesions or merely artifacts. Even after determining that a lesion exists, delineating the lesion boundaries is also complicated by these artifacts and partial volume effects in the image, and is particularly difficult for diffuse lesions even in the absence of artifacts. Inter-and intra-rater variability in assessing disease involvement in spinal cord MR volumes can arise from these challenges and could have important clinical consequences if it affects the treatment selected for the patient.

Automated methods have been proposed to aid radiologists and neurologists in the challenging task of identifying lesions in the spinal cord and to reduce inter-and intrarater variability [START_REF] Gros | Automatic segmentation of the spinal cord and intramedullary multiple sclerosis lesions with convolutional neural networks[END_REF]. Methods based on deep learning are now the state of the art in this regard. However, these methods rely on learning lesion patterns from a set of images where lesions have been delineated by experts. Inter-rater variability can therefore affect these models in two ways. Firstly, the models may learn to recreate the segmentation patterns of the individual raters rather than the "true" segmentation, leading to sub-optimal results [START_REF] Karimi | Deep learning with noisy labels: Exploring techniques and remedies in medical image analysis[END_REF]. Secondly, to evaluate the model, its predictions are compared to some expert "ground truth" annotation, and variability in these annotations could result in an underestimate or overestimate of model performance. It is therefore important to assess the extent of rater variability, both to explore the potential added value of automated methods and Accepted at CBMS 2023 © 2023 IEEE.

to understand the potential consequences for the training and evaluation of deep learning-based methods.

The main objectives of this study are to:

• assess inter-and intra-rater variability in segmentations of MS lesions in spinal cord MR volumes, and • evaluate the relative performance of experts and assess a state-of-the-art automated method in this context. Several studies have examined inter-rater variability on annotations of MS lesions in brain MR volumes [START_REF] Egger | MRI FLAIR lesion segmentation in multiple sclerosis: Does automated segmentation hold up with manual annotation?[END_REF]- [START_REF]Evaluating White Matter Lesion Segmentations with Refined Sørensen-Dice Analysis[END_REF], and Gros et al. [START_REF] Gros | Automatic segmentation of the spinal cord and intramedullary multiple sclerosis lesions with convolutional neural networks[END_REF] reported the variability in Dice scores of several raters compared to a majority consensus in spinal cord MR volumes. The current study builds upon these previous works to explore and quantify the variability with regard to spinal cord MR imaging. In particular, we include an intrarater analysis, an evaluation of raters with respect to a ground truth adjudicated by a senior expert, and a deeper analysis of variability with respect to lesion-wise sensitivity and precision. Furthermore, we evaluate an automated segmentation method based on the nnU-Net [START_REF] Isensee | nnU-Net: a self-configuring method for deep learning-based biomedical image segmentation[END_REF] framework within the context of this rater variability. Finally, we compare automatically generated consensuses to an adjudicated segmentation of a senior expert. In the absence of an expert consensus or adjudication, studies on inter-rater variability often use automated consensus methods, e.g., Simultaneous Truth and Performance Level Estimation (STAPLE) [START_REF] Warfield | Simultaneous truth and performance level estimation (STAPLE): an algorithm for the validation of image segmentation[END_REF] or taking a majority vote on each voxel. We examine how well these consensus methods approximate the segmentation of a senior expert in the context of spinal cord MS lesions.

II. EXPERIMENTAL SETUP

A. Dataset

The data for this study were gathered from the OFSEP (Observatoire Franc ¸ais de la Sclérose en Plaques -French MS registry) database [START_REF] Vukusic | Observatoire Franc ¸ais de la Sclérose en Plaques (OFSEP): A unique multimodal nationwide MS registry in France[END_REF]. The data consisted of MR volumes of six patients diagnosed with MS. Each patient had one sagittal T2-weighted (T2-w) volume for the upper section of the spinal cord (12-13 vertebrae), and another for the lower spinal cord. Furthermore, three of the six patients also had corresponding Short inversion Time Inversion Recovery (STIR) volumes available from the same acquisition session [START_REF] Krinsky | Nonspecificity of short inversion time inversion recovery (STIR) as a technique of fat suppression: pitfalls in image interpretation[END_REF]. Three subjects were imaged on a Siemens Aera (1.5T) scanner, and the scans of the other subjects came from a Siemens Avanto (1.5T), a Siemens Spectra (3T) and a General Electric Optima MR450w (1.5T).

The data were selected based on an expert visual inspection with the following criteria: the MR volumes should not present strong acquisition artifacts (e.g., noise, motion), and the selected data should be representative of a range of potentially encountered clinical cases, i.e., ranging from a patient with no visible lesions to patients with a high lesion load. Based on the segmentation of the senior expert described in Sec. II-D1, three of the twelve MR volumes contained no lesions (i.e., both volumes of one patient and one volume of another). In the remaining nine volumes, a total of 63 lesions were segmented by the senior expert, with a mean volume of 599.9mm 3 . Each of these images contained between 2 and 15 lesions, with a minimum and maximum volume of 0.4mm 3 and 3056.6mm 3 , respectively. The voxel sizes of the T2-w acquisitions were (3.04 ± 0.56, 0.53 ± 0.14, 0.53 ± 0.14) mm 3 , and (2.90 ± 0.15, 0.61 ± 0.18, 0.61 ± 0.18) mm 3 for the STIR volumes, reported as mean ± standard deviation. The MR volumes were acquired on four different scanners from two vendors.

B. Expert Annotations

Four experts were asked to identify and delineate the lesions on each of the 12 sagittal T2-w volumes, with the aid of STIR volumes when available (3/6 patients). The raters comprised three radiologists and one neurologist, and each had approximately five years of experience. The segmentations were performed using the ITK-SNAP software [START_REF] Yushkevich | User-guided 3D active contour segmentation of anatomical structures: Significantly improved efficiency and reliability[END_REF], all raters having been previously trained to use it for this particular task. The raters were asked to repeat this segmentation task on the same MR volumes in a second session two weeks later, resulting in eight rater annotations per MR volume.

C. Automated Segmentation

A segmentation model was trained using the nnU-Net framework [START_REF] Isensee | nnU-Net: a self-configuring method for deep learning-based biomedical image segmentation[END_REF] on a separate dataset of 62 sagittal T2-w volumes for training and 13 volumes as a validation set (voxel sizes=(2.81 ± 0.20, 0.58 ± 0.12, 0.58 ± 0.12) mm 3 , 44 upper and 18 lower spinal cord volumes, 10 scanner models, 4 brands). The training set contained other volumes from the same four MR scanners on which the six subjects of the current study were imaged. The ground truth (GT) segmentation for each MR volume in the training set was created by one of six experts, including the four experts in the current study. The masks for each MR volume were validated by the senior expert in the current study.

nnU-Net is a self-configuring method based on the U-Net architecture that determines an optimal network topology and hyper-parameters using heuristics and empirical tests during five-fold cross-validation on the training dataset. It has achieved state-of-the-art performance across a variety of medical image segmentation tasks [START_REF] Antonelli | The Medical Segmentation Decathlon[END_REF].

The 3D volumes were pre-processed by zeroing voxels that do not belong to the spinal cord and applying a re-orientation (RPI, i.e. Right-to-left, Posterior-to-anterior, Inferior-to-superior) and a resampling to voxel sizes of (0.5, 0.5, 0.5) mm 3 . Spinal cord masks were obtained using the Spinal Cord Toolbox, as the union of the output masks of the sct deepseg sc function, dilated with a 5×5×5 structuring element, and the sct create mask function, which selects 20mm×20mm axial sections around the spinal cord centreline [START_REF] Leener | SCT: Spinal Cord Toolbox, an open-source software for processing spinal cord MRI data[END_REF]. The union of these two masks was used to avoid occasional errors observed in the individual methods.

The size of input patches to the model was chosen by the nnU-Net framework to be 40×128×496 voxels. Training batches comprised two patches from a single MR volume, where at least one of the patches contained lesion voxels, except for MR volumes with no lesions. The U-Net model chosen by nnU-Net consisted of 5 convolution levels, with

D. Ground Truth

1) Senior Expert Adjudication: Following the rater annotations and model predictions, a total of nine annotations were available for each volume (two annotations per rater and one model prediction). These nine annotations were merged into a single mask of "plausible lesion" voxels for each volume, i.e., any voxel labelled as a lesion voxel in any of the nine annotations was kept as a lesion voxel. A senior neurologist with 10 years of experience was asked to adjudicate on this mask, keeping valid lesion voxels and removing non-lesion voxels. This adjudicated segmentation was used as GT for quantitative evaluations.

Fig. 1 shows an example of a segmented lesion in the dataset, which was detected in seven of the eight rater annotations. The segmentation mask created by the senior expert during the adjudication phase is slightly smaller than accepting all voxels labelled by any of the raters.

2) Automated Consensus: We assessed how well several consensus generation techniques approximated the senior expert adjudication. These included accepting all lesion voxels marked in any of the eight annotations or by the model (Any Rater), accepting a lesion voxel if it was segmented by at least five of the nine segmentations per volume (Majority), and applying the STAPLE algorithm [START_REF] Warfield | Simultaneous truth and performance level estimation (STAPLE): an algorithm for the validation of image segmentation[END_REF] to the nine segmentations for each volume. STAPLE iteratively estimates the GT by weighting each annotation by its estimated sensitivity and precision, and it is one of the most popular automated consensus generation methods [START_REF] Lampert | An Empirical Study Into Annotator Agreement, Ground Truth Estimation, and Algorithm Evaluation[END_REF]. One of the assumptions of STAPLE is that the annotations are independent of each other, which is not the case in our study as we have two annotations per rater. However, as the number of annotations per rater are equal across raters, it is less serious than if one senior rater were outnumbered by junior raters with correlated annotations, the example given in the original paper on STAPLE [START_REF] Warfield | Simultaneous truth and performance level estimation (STAPLE): an algorithm for the validation of image segmentation[END_REF].

We evaluated these consensus segmentation masks against the senior expert segmentation in the same way as when evaluating the raters and model, which will be described in Sec. II-F. We also assessed the impact on reported Dice similarity scores of the raters and model when using these automated consensuses as GT.

E. Variability Evaluation

The inter-and intra-rater variability was assessed by summarising the number of lesions and lesion volumes segmented by each rater in each session. The Cohen kappa value was computed to evaluate the intra-rater agreement in lesion count between the two sessions for each rater, while the Fleiss kappa calculation was used to assess the inter-rater agreement in each session.

F. Performance Evaluation

Evaluation metrics were based on those presented in Gros et al. [START_REF] Gros | Automatic segmentation of the spinal cord and intramedullary multiple sclerosis lesions with convolutional neural networks[END_REF] and the MICCAI16 brain MS segmentation challenge [START_REF] Commowick | Objective Evaluation of Multiple Sclerosis Lesion Segmentation using a Data Management and Processing Infrastructure[END_REF]. Voxel-wise sensitivity, precision and Dice similarity coefficient (F1-score) were used to measure the overlap between the senior expert GT and the manual or automatic segmentations. The detection and count of lesions is vital in monitoring MS patients undergoing disease-modifying treatments. In this context, lesion-wise metrics have more value for clinicians than voxel-wise metrics. Therefore, sensitivity, precision and F1-score were also defined at the lesion level.

A GT lesion was treated as correctly detected for the purposes of sensitivity if more than 10% of its voxels were labelled as lesion voxels by a manual or automatic segmentation, i.e., if the voxel-wise sensitivity for this GT lesion exceeded 10%. The 10% threshold follows that used in [START_REF] Commowick | Objective Evaluation of Multiple Sclerosis Lesion Segmentation using a Data Management and Processing Infrastructure[END_REF]. Lesionwise sensitivity is then the ratio of the number of correctly detected GT lesions to the total number of GT lesions.

A predicted lesion was considered to be precise if over 30% of its voxels are truly lesion voxels in the GT, i.e., if the voxelwise precision of this predicted lesion exceeded 30%. Again, this threshold value follows [START_REF] Commowick | Objective Evaluation of Multiple Sclerosis Lesion Segmentation using a Data Management and Processing Infrastructure[END_REF]. To avoid a segmentation being considered more precise if it splits one GT lesion into many predicted lesions, we counted precise predicted lesions only once per GT lesion when calculating precision. Lesionwise precision is then the ratio of precise predicted lesions to the total number of predicted lesions.

Defining lesion-wise sensitivity and precision separately in terms of only voxel-wise sensitivity and precision, respectively, allows a more straightforward interpretation of the results, i.e., how sensitive and how precise the segmentations were at the lesion level. We further used these lesion-wise sensitivity and precision values in the calculation of lesionwise F1-score. 1) Lesion Counts: Inter-and intra-rater variability can first be observed in the number of lesions segmented on each MR volume. Differences between raters, as well as between sessions for the same rater, can be observed in the lesion counts in Fig. 2. Concerning intra-rater variability, the number of lesions varied between sessions for all raters. For example, rater 1 segmented six lesions in the upper spinal cord of subject 5 in session 1 but segmented ten lesions in the same MR volume in session 2 (brown blocks in Fig. 2). Comparing the number of lesions in each MR volume across the two sessions for each rater, Cohen kappa values of 0.61, 0.24, 0.50 and 0.50 were obtained for raters 1 to 4, respectively, indicating a high variability from poor to substantial agreement between both sessions, depending on the rater. Inter-rater variability is also apparent, as the total number of lesions ranges from 28 to 60 between raters 3 and 2, respectively. Fleiss kappa values of 0.26 and 0.15 were obtained for sessions 1 and 2, respectively, indicative of the low agreement between raters.

III. RESULTS

A. Inter-and Intra-Rater Variability

2) Lesion Volumes: Inter-and intra-rater variability is also noticeable in the volume of segmented lesions. Fig. 3 highlights this aspect by presenting the distribution of volumes of individual lesions. Table I complements this figure with summary statistics of the segmented lesion volumes by rater. Raters 2 and 4 tended to delineate smaller regions, with median lesion volumes of 52.6mm 3 and 64.5mm 3 , respectively. Raters 1 and 3, on the other hand, had a median lesion volume over double that of rater 2, with values of 113.3mm 3 and 110.3mm 3 , respectively.

While Fig. 3 and Table I provide a quantitative evaluation of segmented lesion volumes from the whole dataset (12 MR volumes), a finer analysis at the image level provides nuances to these observations. Indeed, comparisons between raters do not always hold when comparing across MR volumes (Fig. 4). For example, rater 4 segments the highest total lesion volume for subject 3 but the lowest volume for subject 5. Fig. 4 also shows that a higher number of segmented lesions does not necessarily correspond to a higher segmented volume. For instance, rater 2 delineates ten lesions in the upper spinal cord of subject 5 (image 5-U), compared to the six lesions of rater 1. However, the total volume segmented by rater 1 is over three times that of rater 2, indicating a higher disease involvement. This demonstrates that although the number of lesions is a commonly used clinical indicator, it is important to also take the lesion volume into account.

3) Lesion-level Agreement: Fig. 5 shows the level of agreement between the raters' segmentations on individual lesions that appear in the GT adjudicated by the senior expert. 13 GT lesions were detected by different numbers of annotations (four raters with two sessions each). A GT lesion is counted as detected if the annotation has segmented over 10% of its voxels. For each level of rater agreement, the mean volume of lesions is shown in orange. For example, 13 GT lesions were present in all eight annotations, and these lesions had a mean volume of 956mm 3 . lesions were detected by all four raters and in both sessions, suggesting obvious lesions on which all raters agree, although there may still be disagreement on the exact delineation. On the other side of the figure, 16 GT lesions were detected in only one annotation, meaning that even the same rater did not detect them across sessions. Moreover, two lesions in the senior expert GT were not present in any of the raters' annotations; these two lesions were detected by the automated method. Importantly, 34 (54%) of the 63 GT lesions are present in only one to four annotations. Therefore, if a majority vote was used to create a consensus GT from these eight annotations, more than half of the valid lesions would be removed. Lesions with a higher level of rater agreement tended to have higher volumes, as the mean volume of lesions detected in all eight annotations was 956.0mm 3 compared to 53.1mm 3 for those present in only one annotation.

B. Performance against Ground Truth

Most of the voxels predicted to be lesion voxels by the raters and model were validated by the senior expert, as seen in the voxel-wise precision distribution in Fig. 6. In particular, for 7 of the 9 MR volumes with lesions, the automated segmentation method has a voxel-wise precision of over 98%. On the other hand, the voxel-wise sensitivity scores are significantly lower. Rater 1 has the highest sensitivity, followed by the other three raters who all have median voxel-wise sensitivity between 20-30%. Finally, the median voxel-wise sensitivity for the model is 14%. The patterns of Dice score across raters are similar to voxel-wise sensitivity, albeit slightly higher because of the influence of voxel-wise precision on the Dice similarity.

Lesion-wise sensitivity is derived from voxel-wise sensitivity (see Sec. II-F), so we might expect the distributions of the two metrics to be similar. However, different patterns arise the distribution of lesion-wise sensitivity. For example, rater 1 has a higher voxel-wise sensitivity and Dice coefficient in general, but this is not reflected in a generally higher lesionwise sensitivity. This suggests that rater 1 is not detecting more lesions than the other raters, but that rater 1 segments comparatively larger areas for the lesions they detect, which also aligns with the segmented lesion volume statistics presented in Sec. III-A2.

A further difference when comparing voxel-wise and lesionwise sensitivity is that the automated method has a more similar distribution to the four raters in lesion-wise sensitivity, with its median sensitivity even surpassing that of rater 3. The different patterns in voxel-wise and lesion-wise sensitivity demonstrate the added value of lesion-wise metrics in addition to voxel-wise metrics.

The results in Fig. 6 present an important overall picture. Most of the lesions that were segmented by the raters were accepted by the senior expert in the adjudication phase, which can be seen by the high precision. However, many lesions were missed by each rater, leading to an underestimate of disease involvement. In particular, none of the raters achieved a median lesion-wise sensitivity of 60%.

A comparison of the lesion-wise sensitivities between the two manual segmentation sessions (Fig. 7) indicates that although some intra-expert differences exist (e.g., for rater 2), the differences between experts are more substantial. This aligns with the observations of variability in lesion counts and volumes in Sec. III-A.

Finally, specificity was calculated with respect to the MR volumes. According to the senior expert GT, there were no lesions in three of the 12 MR volumes in this study. Subject 2 had no GT lesions in either the upper (2-U) nor lower (2-L) volumes, whereas subject 1 had no lesions in the upper (1-U) spinal cord, but had two lesions in the lower (1-L) spinal cord. Indeed, raters 3 and 4 correctly segmented no lesions in the three volumes with no GT lesions. Rater 1, however, segmented one lesion in volume 2-U in both sessions. Rater 2 segmented one lesion in 2-L in both sessions, one lesion in 2-U in the first session, and one lesion in 1-U in the second session. The model segmented five lesions in 2-U, but returned no lesions in 2-L nor 1-U. In other words, raters 1 to 4 had volume-level specificities of 66.7%, 33.3%, 100%, and 100%, respectively, while the model had a specificity of 66.7%.

C. Potential Added Value of the Automated Method

To simulate the potential added value of the model to clinical practice, the number of lesions detected by the automated method but not by individual raters was calculated and is presented in Table II. A mean number of 4.9 (min=3, max=8) lesions could have been potentially retained with the help of the model, which could lead to a mean lesion-wise sensitivity increase of 7.8% (min=4.8%, max=12.7%) for each rater.

D. Evaluation of Automated Consensuses

Using a majority vote on the nine segmentation masks per MR volume (one automated and eight manual segmentations) returned a total of 39 lesions with a mean volume of 112.1mm lesions, leaving 59 lesions with a mean volume of 216.0mm 3 . A total of 75 lesions with a mean volume of 272.1mm 3 were present after taking a union of all nine masks, i.e., accepting any voxel marked in any of the nine masks as a lesion voxel (Any Rater).

We calculated the same evaluation metrics for the consensus segmentations as before for the individual rater segmentations, comparing each consensus to the GT adjudicated by a senior expert. These results are presented in Table III. As might be expected after considering the high precision of individual raters observed in Sec. III-B, the 'Any Rater' consensus best approximated the final senior expert GT, with a median Dice similarity of 83.7%. Although the majority consensus was highly precise, with a median voxel-wise precision of 97.9%, it discarded many valid lesion voxels, resulting in a low voxel-wise sensitivity of 28.2% and Dice of 44.0%. Finally, STAPLE kept more valid lesion voxels, with a median voxelwise sensitivity of 63.2% and Dice of 74.5%. As a result, STAPLE outperformed the median rater in voxel-wise Dice and lesion-wise F1-score, whereas majority voting did not. The previous analysis focused on how well automated consensuses approximated the senior expert GT, but we also examine how reported performance metrics would change by using an automated consensus as GT. For this, we calculated the median Dice scores by taking the median first over the nine MR volumes and then over the raters and sessions. The median Dice score for the raters against the senior expert GT was 48.9%, but when evaluated against the STAPLE and majority consensuses, this rose to 51.8% and 67.3%, respectively, whereas reported performance would drop to 40.2% if the 'Any Rater' mask was used as GT. The median Dice of the model changes in a similar way, with scores of 16.4%, 24.8%, 26.4% and 42.3% using the 'Any Rater', adjudicated, STAPLE or majority consensuses, respectively, as GT.

IV. DISCUSSION

A. Inter-rater Variability

Considerable inter-rater variability was observed in this study. The total number of lesions segmented by each rater in all 12 MR volumes ranged from 28 to 60, and only 13 of the 63 lesions validated by the senior expert were segmented by all four raters in both sessions. Even when raters agreed on lesions, the segmented lesion volumes varied significantly. For example, the median lesion volumes segmented by raters 1 and 3 were over double that of rater 2, and the maximum lesion volume segmented by rater 1 (2549mm 3 ) was over five times that of rater 4 (449mm 3 ). This inter-rater variability leads to very different depictions of the total lesion load and the disease involvement for individual patients.

Clinical patient stratification based on the number of spinal cord lesions can inform patient care and treatment decisions. The disagreements observed between raters, and between the two segmentations of the same rater, could therefore lead to different patient care pathways. For instance, for subject 1, rater 3 segmented no lesions in one session, whereas all other manual annotations showed a lesion of approximately 100mm 3 for this subject. This single lesion could change the prognosis for the patient and may thus affect therapeutic decisions.

The raters were instructed to segment lesions only if they had a high confidence that they were indeed lesions. This may have led to increased variability, as it is not just a question of detection of hyper-intense areas, but also a subjective assessment of the probability those areas represented lesions. Furthermore, the raters were instructed that segmentation boundaries should include only hyper-intense voxels. However, areas of hyper-intensity can also arise from artifacts and partial volume effects with cerebrospinal fluid, for example.

Moreover, what appears as hyper-intense is influenced by the manual adjustment of image contrast by each rater using the ITK-SNAP viewer, which introduces further subjectivity. These factors partly explain the observed disagreements on lesions and lesion volumes.

B. Rater Performance

Evaluated against a ground truth (GT) of a senior expert adjudicating on all possible lesion voxels, both the raters and the automated segmentation method tended to be "precise". Specifically, when the raters or model found and delineated a lesion, the senior expert tended to agree that a lesion existed in that location. On the other hand, the median lesion-wise sensitivity of each rater ranged from 28% to 60%, which demonstrates that detecting lesions in spinal cord MR volumes remains a challenging task and can lead to underestimates of disease activity for patients.

In this study, we considered the segmentations of the senior expert to be the "ground truth" and a lower lesion sensitivity was interpreted to mean that some lesions were missed. However, it is possible that the raters detected and examined the same hyper-intense regions but disagreed as to whether they represented lesions. Indeed, for certain lesions it is challenging to conclusively determine if it is truly a lesion or not based on a single spinal cord MR volume.

Comparing to previous studies, Combès et al. [START_REF] Combès | A Clinically-Compatible Workflow for Computer-Aided Assessment of Brain Disease Activity in Multiple Sclerosis Patients[END_REF] reported a Fleiss kappa coefficient of 0.47 in an inter-rater study of new lesion segmentation in brain MR volumes, compared to 0.26 and 0.15 for the two annotation sessions here, indicating a lower agreement between raters in the current study. The same study reported a range of lesion-wise sensitivity of the three raters between 60-66%. The median Dice score was 48.9% in the current study, which is lower than other studies on MS lesion segmentation in brain and spinal cord MR volumes, albeit with differing methodologies. Gros et al. [START_REF] Gros | Automatic segmentation of the spinal cord and intramedullary multiple sclerosis lesions with convolutional neural networks[END_REF] reported a median Dice of 60.7% comparing seven raters to a majority consensus in spinal cord MR volumes, Carass et al. [START_REF]Evaluating White Matter Lesion Segmentations with Refined Sørensen-Dice Analysis[END_REF] found a Dice of 66-67% in brain MR volumes, comparing two raters to a STAPLE consensus of the raters and 14 algorithms, while Egger et al. [START_REF] Egger | MRI FLAIR lesion segmentation in multiple sclerosis: Does automated segmentation hold up with manual annotation?[END_REF] observed a mean Dice of 66% over pairwise comparisons between three raters on brain MR volumes. However, caution must be used when comparing these figures, as if a majority vote had been used in the current study in place of an adjudication by a senior expert, the reported median Dice would rise from 48.9% to 67.3%. Moreover, with the small sample size used in this study and by Gros et al. [START_REF] Gros | Automatic segmentation of the spinal cord and intramedullary multiple sclerosis lesions with convolutional neural networks[END_REF] (ten MR volumes), these figures may not be representative of the performance over a wider population.

C. Model Performance

The variability among raters and the potential of missing lesions and underestimating disease activity is a motivating factor for the creation and adoption of automated methods as aids to radiologists and neurologists. Indeed, the automated method tested in this study could have increased sensitivity by between 4.8 and 12.7 percentage points. However, taking the union of the segmentations of a clinician and a model is not how an automated method would be used in practice, so future work will conduct a more in-depth study of the usefulness of an automated method in clinical practice, similar to [START_REF] Combès | A Clinically-Compatible Workflow for Computer-Aided Assessment of Brain Disease Activity in Multiple Sclerosis Patients[END_REF].

The high precision and low sensitivity of the algorithm indicate that a different threshold on the model outputs may yield a better balance between sensitivity and precision. Model output probability maps were thresholded at the default value of 0.5. However, we observed a low model output score (< 0.01) for the majority of voxels labelled as lesions in the GT, so adjusting the binarisation threshold could increase sensitivity, as found in previous studies [START_REF] Combès | A Clinically-Compatible Workflow for Computer-Aided Assessment of Brain Disease Activity in Multiple Sclerosis Patients[END_REF], [START_REF] Hussein | A study on loss functions and decision thresholds for the segmentation of multiple sclerosis lesions on spinal cord MRI[END_REF]. The default value of 0.5 was kept for this study as the segmentation masks provided to the senior expert to create an adjudicated GT had been binarised using this threshold (see Sec. II-D1).

The high precision and low sensitivity of the model could also be influenced by differences between the GT masks used for training the model and those used for evaluation in the current study. In both cases, the masks were validated by a senior expert; however, the initial mask for evaluation was based on the union of eight manual segmentations and the model predictions, whereas for training each mask was created by just one expert. As a result, it is possible that not all lesions were detected in the segmentation masks for training, leading the model to reproduce a lower lesion sensitivity. The high precision and low sensitivity of individual experts in the current study can inform future choices when training deep learning models to segment spinal cord lesions. Moreover, the high variability observed demonstrates that the choice of GT used to evaluate a model could have a substantial effect on reported results.

D. Automated Consensuses

A consensus created with a voxel-wise majority vote poorly approximated the adjudication of the senior expert, with a Dice similarity of 44.0% when taking the median over the nine MR volumes containing GT lesions. STAPLE improved on this approximation with a median Dice of 74.5%, while the best approximation of the GT in this study was given by a union of all the manual segmentations and model output, giving a median Dice of 83.7%. This result is influenced by the high precision of both individual raters and the model; majority voting may still be useful in scenarios where raters are expected a priori to have a high sensitivity and low precision.

E. Limitations and Future Work

Further analysis did not indicate significant differences in rater performance where both STIR and T2-w acquisitions were available compared to T2-w alone. However, the added value of STIR to lesion sensitivity has been discussed in the literature [START_REF] Nayak | A comparison of sagittal short T1 inversion recovery and T2-weighted FSE sequences for detection of multiple sclerosis spinal cord lesions[END_REF], and it is one of the recommended sequences for spinal cord imaging [START_REF] Wattjes | MAGNIMS-CMSC-NAIMS consensus recommendations on the use of MRI in patients with multiple sclerosis[END_REF]. We note that no STIR volumes were used by the model, which could be another factor partly explaining the lower voxel-wise sensitivity relative to the raters. Future work will use both modalities as the inputs of an automated segmentation method, and assess the potential of such a method to further aid spinal cord MS lesion detection in clinical practice.

Future work will also include exploring a better balance between sensitivity and precision in the nnU-Net model, and assessing the impact of the definition of lesion-wise metrics. For example, Gros et al. [START_REF] Gros | Automatic segmentation of the spinal cord and intramedullary multiple sclerosis lesions with convolutional neural networks[END_REF] determined a GT spinal cord lesion to be correctly detected if 25% of its voxels had been segmented, whereas we opted for the lower threshold of 10% used by Commowick et al. [START_REF] Commowick | Objective Evaluation of Multiple Sclerosis Lesion Segmentation using a Data Management and Processing Infrastructure[END_REF] for brain lesions. This decision was motivated by the large variability observed between raters in volumes of segmented lesions. However, this threshold and the metrics in general need to be assessed as to whether they are clinically meaningful.

A limitation of this study is the small sample size of six patients. We make no claim that the reported statistics represent the true sensitivity and precision of experts in clinical practice, nor even for the specific experts in this study. Moreover, while we have demonstrated that significant variability can exist between raters, more large-scale studies are required to determine the factors that influence this variability.

Fig. 1 .

 1 Fig. 1. Example of a lesion in the dataset (at vertebra T6). The bright vertical bands are cerebrospinal fluid surrounding the spinal cord. Left: subsection of sagittal T2-w volume. Centre Left: contrast adjusted to improve lesion visibility. Centre Right: voxel-wise agreement between the eight annotations (four raters with two sessions each). One of the raters did not detect this lesion in one of the sessions. Right: predicted segmentation of nnU-Net model (blue) and adjudicated segmentation by senior expert (green).

Fig. 2 .

 2 Fig. 2. Number of lesions segmented by each rater in each session. The blocks in each column correspond to individual MR volumes.

Fig. 3 . 3 )Fig. 4 .

 334 Fig.3. volumes segmented by each rater in each session. Each point represents one lesion. This plot is cropped to 600 mm 3 , preventing some outliers to appear.

3 )Fig. 5 .

 35 Fig.5. Rater agreement and lesion volume. The blue bars show how many GT lesions were detected by different numbers of annotations (four raters with two sessions each). A GT lesion is counted as detected if the annotation has segmented over 10% of its voxels. For each level of rater agreement, the mean volume of lesions is shown in orange. For example, 13 GT lesions were present in all eight annotations, and these lesions had a mean volume of 956mm 3 .

Fig. 6 .Fig. 7 .

 67 Fig.6. Lesion-wise and voxel-wise metrics for each rater and the algorithm predictions in the first annotation session evaluated against the senior expert GT. Metrics are first calculated for each of the nine MR volumes with GT lesions (individual points), and boxplot distributions over these nine values are shown.

TABLE I SUMMARY

 I STATISTICS OF VOLUMES (MM 3 ) OF LESIONS DELINEATED BY EACH RATER ACROSS BOTH ANNOTATION SESSIONS.

	Rater	Min	Max	Mean Std. Dev. Median
	Rater 1	1.1	2548.5	246.6	397.3	113.3
	Rater 2	2.2	1633.9	131.2	248.1	52.6
	Rater 3	0.4	1298.2	203.3	260.8	110.3
	Rater 4	5.3	448.5	86.8	72.4	64.5

  3 . A STAPLE consensus removed fewer of the

		1.0
		0.8
	Metric Value	0.4 0.6					Rater 1 Rater 2 Rater 3 Rater 4 nnU-Net
		0.2
		0.0
		Sensitivity	Precision	Dice	Sensitivity	Precision	F1-score
			Voxel-wise			Lesion-wise

TABLE II NUMBER

 II OF GT LESIONS DETECTED BY THE NNU-NET MODEL ANDMISSED BY A RATER, FOR BOTH SEGMENTATION SESSIONS.

	Rater	Session 1 Session 2
	Rater 1	3	3
	Rater 2	7	5
	Rater 3	8	5
	Rater 4	3	5

TABLE III EVALUATION

 III OF CONSENSUS METHODS AGAINST SENIOR EXPERT ADJUDICATED GT. SHOWN ARE THE MEDIAN SCORES ACROSS THE NINE MR VOLUMES WITH GT LESIONS. THE MEDIAN RATER PERFORMANCE IS GIVEN FOR COMPARISON. (Sens=Sensitivity, Prec=Precision)

			Voxel-wise			Lesion-wise
	Consensus Sens Prec Dice	Sens	Prec	F1
	Majority	28.2	97.9	44.0	40.0	100.0	57.1
	STAPLE	63.2	84.4	74.5	60.0	100.0	75.0
	Any Rater	99.3	72.3	83.7	100.0	100.0 100.0
	Raters	33.1	88.9	48.9	50.0	100.0	66.7

This study was partly funded by the French doctoral program in artificial intelligence (reference: ANR-20-THIA-0018). The study was supported by the French National Research Agency (ANR) within the France 2030 program, 3rd PIA (reference ANR-21-RHUS-0014). Data collection was supported by a grant provided by the French State and handled by ANR within the France 2030 program (reference ANR-10-COHO-002 OFSEP).


q q q q q q q q q q q q q q q q q q q q q q q q q q q q All data fit Working set (MB) GFlop/s q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q HFP HFP only HFP with load balance hMETIS hMETIS + HFP q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q B = M A + B = M GFlop/s max 0 2500 5000 7500 10000 12500 500 1000

Working set (MB) GFlop/s q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q HFP HFP only ready DMDAR MST HFP only ready and flip HFP only flip RCM Eager q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q B = M A + B = M All data fit GFlop/s max 0 3000 6000 9000 2500 5000 7500

Working set (MB)

GFlop/s (a) Performance.

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q B = M (b) Amount of data transfers.

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q HFP HFP only Belady DMDAR MST

HFP only ready HFP only flip RCM Eager

Maxime GONTHIER -Samuel THIBAULT -Loris MARCHAL maxime.gonthier@ens-lyon.fr samuel.thibault@inria.fr -loris.marchal@ens-lyon. Initial intuition: Group tasks with common data into packages of maximum size M = GPU's memory. The goal is to minimize the number of packages. This will maximize data reuse Hierarchical Fair Packing step by step: 0 Initially each task is an elementary package 1 Compute the number of common data on all pairs of packages 2 Consider packages P i with the fewest tasks 3 Merge couples (P i , P j ) which reach max number of common data such that Weight(P i P ∪ P j ) ≤ M 4 Repeat from step 1 until there are no more possible merge 5 Consider that M = ∞ and repeat from step 1 as long as there is more than 1 package q q q q q q q q q q q q q q q q q q q q q q q q GFlop/s max 0 200 400 600 500 1000

Working set (MB) GFlop/s q q q q q q q q q q q q q q q DMDAR Eager Dynamic outer We use a set of independant tasks sharing input data. We want to produce a task order σ(t) and an eviction policy ν(t) The processing order is the order in which vertices are added to the spanning tree. Schedules tasks where their completion times is expected to be minimal. Uses a ready strategy to favor tasks whose data has already been loaded into memory.

▸ Eager Eager

Process tasks in the order in which it receives them. Serves as a baseline. Experiments on a tesla V100 GPU Experiments on a tesla V100 GPU as well as simulations on Simgrid using the same performance model.

as well as simulations on Simgrid using the same performance model. ▸ 6000 MB/s PCI bandwidth to represent the bandwidth share typically available for a given GPU in a multi-GPU 6000 MB/s PCI bandwidth to represent the bandwidth share typically available for a given GPU in a multi-GPU platform platform ▸ GPU memory limited to 500 MB in order to better distinguish the performance of different strategies even on small GPU memory limited to 500 MB in order to better distinguish the performance of different strategies even on small datasets datasets The dotted horizontal black line represents the maximum Gflop/s that the GPU can achieve when processing elementary matrix product and is our asymptotic goal.

Applications

The red dotted vertical line denotes the situation when the GPU memory can fit exactly only one of the two input matrices, and the orange line denotes the situation when it can accommodate both input matrices.

The black diagonal dottedline is the maximum number of transfers that can be done during the minimum time for computation, thus the hard limitation induced by the PCI bus bandwidth.

▸ Our contibutions: Our contibutions: ▸ We provided a formal model of the optimization problem of ordering independent tasks sharing input data in order to minimize the ▸ amount of data transfers ▸ HFP allows on average a 8.3% (resp. 11% and 11.6 %) improvement for a 2D matrix multiplication (resp. 3D matrix multiplication and Cholesky).

▸ HFP is good in distributing data transfer over time to increase transfer/computation overlap ▸ HFP is well suited to multi GPU Dynamic outer is a dynamic scheduler Dynamic outer is a dynamic scheduler based on processing task whose data based on processing task whose data are already loaded. We adapted it to are already loaded. We adapted it to choose to load (resp. evict) data allowing choose to load (resp. evict) data allowing to process the most task (resp. the least to process the most task (resp. the least task) task)

Results:

Results: ▸ Dynamic outer achieve performance Dynamic outer achieve performance similar to HFP with a smallest complexity similar to HFP with a smallest complexity HFP's packing allow us to easily adapt it to multi GPU HFP's packing allow us to easily adapt it to multi GPU cases. Load balancing was added in order to equilibrate cases. Load balancing was added in order to equilibrate the predicted computation time on each GPU. the predicted computation time on each GPU. We also added a partitionner, hMETIS, and used it in We also added a partitionner, hMETIS, and used it in combination with HFP. combination with HFP.

Results:

Results: ▸ HFP with load balance reaches good performances HFP with load balance reaches good performances ▸ hMETIS coupled with HFP also achieve satisfying hMETIS coupled with HFP also achieve satisfying performance performance

GFlop

▸

  The communication speed is limited The communication speed is limited between the GPU and the main memory between the GPU and the main memory which induce communication cost which induce communication cost ▸ Reverse-Cuthill-McKee (RCM): Reverse-Cuthill-McKee (RCM): Permutes a sparse matrix into a minimal band matrix. Vertices sharing an edge are at most k edges away. Thus all tasks sharing a data are processed in the interval [t -k ; t + k].

Figure 2 :

 2 Figure 2: Example with a memory holding at most M=3 data. The graph of input data dependencies is shown on the left. The schedule on the right corresponds to processing the tasks in the natural order with the following eviction policy: ν(1) = ν(2) = ø, ν(3) = {1}, ν(4) = {2}, ν(5) = {3, 4}. This results in 7 loads.

Figure 1 :

 1 Figure 1: Model representing the link between main memory and local memory of the processor unit

Figure 4 :Figure 3 :

 43 Figure 4: Exemple of a maximum spanning tree on a graph

Figure 5 :

 5 Figure 5: Visualisation of DMDAR's ordering on a 2D matrix multiplication

Figure 6 :Figure 7 :

 67 Figure 6: Package flipping. The pair of sub-packages (P i end , P j end ) is the one with the most shared data. Thus P j is reversed before merging packages.

▸Figure 8 :Figure 9 :Figure 10 :

 8910 Figure 8: GFlop/s on 2D matrix multiplication using a Tesla V100 GPU

▸

  We adapted two heuristic from the literature and designed a new heuristic based on gathering tasks with similar data into packages ▸ We implemented all three heuristics into the StarPU runtime ▸ We implemented a new scheduler : Dynamic Outer Results: Results:

Figure 11 :

 11 Figure 11: Gflop/s on a 2D matrix multiplication in simulation with 3 GPUs

Figure 12 :

 12 Figure 12: GFlop/s on a 2D matrix multiplication in simulation


Maxime Gonthier -Samuel

Thibault -Loris Marchal

Memory-Aware Scheduling Memory-Aware Scheduling Of Tasks Sharing Data On Multiple GPUs Of Tasks Sharing Data On Multiple GPUs

Figure: Example with a memory of size 2 data. The graph of

Schedules tasks where their completion times are expected to be minimal. Uses a ready strategy to favor tasks whose data has already been loaded into memory.

▸ EAGER EAGER

Process tasks in the order in which it receives them our baseline ⇒ our ba ▸ Tesla V100 GPUs Tesla V100 GPUs ▸ 12000 MB/s PCI bandwidth 12000 MB/s PCI bandwidth ▸ GPU memory limited to 500 MB to better distinguish the performance of different strategies GPU memory limited to 500 MB to better distinguish the performance of different strategies even on small datasets even on small datasets ▸ Pathological matrix size due to LRU's eviction policy when the memory is a constraint (after the red Pathological matrix size due to LRU's eviction policy when the memory is a constraint (after the red dotted line) for EAGER, hMETIS and DARTS dotted line) for EAGER, hMETIS and DARTS

▸ DMDAR suffers from a conflict between prefetches and evictions DMDAR suffers from a conflict between prefetches and evictions The dotted horizontal black line represents the maximum Gflop/s the GPUs can achieve when processing elementary matrix product and is our asymptotic goal. The red dotted vertical line denotes the situation when the cumulated GPUs memory can fit exactly only one of the two input matrices, and the orange line denotes the situation when it can accommodate both input matrices.

The black diagonal dotted line on the figure on the right is the maximum number of transfers that can be done during the minimum time for computation, thus the hard limitation induced by the PCI bus bandwidth.

Motivation: Extract peak performance from GPUs Motivation: Extract peak performance from GPUs

To get peak performance we need to minimize the number of access to the main memory ▸ Motivations : Motivations :

• High-performance computing applications, like weather and climate forecasting, have an increasing demand in computer power • GPUs have a large computation power but have a limited memory. This, coupled with the share of a single PCI express bus by multilple GPUs induce a large amount of data transfers which reduces performance Limiting data movement is crucial to extract the most out of GPUs ⇒ our ba ⇒ our ba We provide a new algorithm (DARTS) and a new eviction policy (LUF) focused on data locality ▸ Results : Results : DARTS+LUF achieves good performance because it : :

• Limits data transfers thanks to the finding of an optimal data and an adapted eviction policy

• Overlaps communication and computations by distributing transfers over time

• Can be used with a reduced complexity ▸ Areas for improvement: Areas for improvement:

• Consider tasks with dependencies

Flow of tasks Flow of tasks

• (1). DARTS fills plannedTasks (1). DARTS fills plannedTasks • (2). Some tasks are popped by (2). Some tasks are popped by taskBuffer and start to be taskBuffer and start to be prefetched prefetched • (3). If the GPU is idling, a task

(3). If the GPU is idling, a task is popped from taskBuffer is popped from taskBuffer

▸ Strategy Strategy

• (1).Find the optimal data D such that the number of tasks (1).Find the optimal data D such that the number of tasks depending on D and on other data in memory is maximum depending on D and on other data in memory is maximum

▸

  Figure: Gflop/s (left) and data transfers (right) on a tiled (1 x n) Matrix Multiplication in outer product order with 2 Tesla V100 GPUs

•Figure :

 : Figure:A hyper-graph with 2 partitions (grey and white)

queues of tasks 2 queues of tasks

	▸ Scientific application Scientific application uses a huge quantity of uses a huge quantity of data data
	GPUs can calculate ▸ GPUs can calculate quickly but their memory quickly but their memory are limited are limited
	The PCI bus becomes a ▸ The PCI bus becomes a bottleneck when large bottleneck when large amounts of data need to amounts of data need to be loaded be loaded
	▸ Independant tasks Independant tasks sharing data sharing data
	Tasks in plannedTasks can be • Tasks in plannedTasks can be removed. removed. Tasks in taskBuffer are being • Tasks in taskBuffer are being prefetched and can't be removed prefetched and can't be removed • Tasks in GPU are being computed Tasks in GPU are being computed
	▸

• Consider inter-GPU communications (NVLink) • Manage multiple MPI nodes ▸ Same size data Same size data ▸ Multiple GPUs with Multiple GPUs with limited memory limited memory Figure: Task flow when using DARTS ▸ 2

input data dependencies is shown on the left. The figure on the right corresponds to the partition and schedule produced


Sarah L Thomson
email: s.l.thomson@stir.ac.uk

Nadarajen Veerapen
email: nadarajen.veerapen@univ-lille.fr

Gabriela Ochoa
email: gabriela.ochoa@stir.ac.uk

Daan Van Den Berg
email: d.van.den.berg@vu.nl

Randomness in Local Optima Network Sampling

Keywords: • Mathematics of computing → Graph algorithms, Combinatorial algorithms, • Theory of computation → Evolutionary algorithms Fitness Landscapes, Quadratic Assignment Problem (QAP), Local Optima Networks (LONs), Iterated Local Search

published or not. The documents may come

L'archive ouverte

INTRODUCTION

Local optima networks (LONs) [START_REF] Ochoa | A study of NK landscapes' basins and local optima networks[END_REF] model local optima connectivity in fitness landscapes and can serve as a vehicle for understanding or predicting algorithm-problem interactions [START_REF] Bożejko | Local optima networks in solving algorithm selection problem for tsp[END_REF][START_REF] Herrmann | Communities of local optima as funnels in fitness landscapes[END_REF]. Although several studies use small search spaces so that exhaustive enumeration of a LON can be conducted [START_REF] Ochoa | A study of NK landscapes' basins and local optima networks[END_REF][START_REF] Ochoa | Understanding phase transitions with local optima networks: number partitioning as a case study[END_REF][START_REF] Ochoa | Local optima networks: A new model of combinatorial fitness landscapes[END_REF][START_REF] Potgieter | A Local Optima Network Analysis of the Feedforward Neural Architecture Space[END_REF][START_REF] Sarah L Thomson | Comparing communities of optima with funnels in combinatorial fitness landscapes[END_REF][START_REF] Verel | Local optima networks of NK landscapes with neutrality[END_REF], there has been an emergent appetite for larger spaces and therefore LON sampling in recent years [START_REF] Melike D Karatas | Towards population-based fitness landscape analysis using local optima networks[END_REF][START_REF] Ochoa | Perturbation strength and the global structure of QAP fitness landscapes[END_REF][START_REF] Marcondes Pavelski | Local Optima Network Sampling for Permutation Flowshop[END_REF][START_REF] Tari | Local search pivoting rules and the landscape global structure[END_REF][START_REF] Sarah L Thomson | Inferring future landscapes: sampling the local optima level[END_REF][START_REF] Verel | Sampling local optima networks of large combinatorial search spaces: The qap case[END_REF][START_REF] Whitley | Local optima organize into lattices under recombination: an example using the traveling salesman problem[END_REF]. The sampling algorithms contain elements of randomness such as starting location; however, it is not known how much of an impact this randomness has. Previous literature shows that evolutionary algorithm (EA) performance can be affected (sometimes to a substantial degree) by randomness [START_REF] Cárdenas-Montes | Sensitiveness of evolutionary algorithms to the random number generator[END_REF][START_REF] Krömer | Randomness and chaos in genetic algorithms and differential evolution[END_REF][START_REF] Maucher | Search heuristics and the influence of non-perfect randomness: examining Genetic Algorithms and Simulated Annealing[END_REF][START_REF] Mark | Randomness and GA performance, revisited[END_REF]; it therefore follows that LON construction -which is done by evolutionary and heuristic search -may partially depend on random factors. We argue that this is an important consideration, especially as the popularity of LONs for understanding and predicting search is growing. The purpose of this paper is to investigate the variation seen in LONs when providing different random seeds to an iterated local search algorithm used for constructing them. Instances from the quadratic assignment problem (QAP) library are used. Although the QAP is a case study here, the consideration of randomness in LONs would be relevant for other problem domains. We consider the nature of the LONs themselves, visualising and analysing several extracted features. Additionally, the predictive potency of LON metrics used in algorithm performance models is studied, with iterated local search (ILS) runtime as the response variable. Throughout the study, the effect of randomisation is given particular attention. The contributions of the work are as follows:

(1) An preliminary, exploratory investigation into the effect of randomness in LON construction (2) Insights about LON metrics which show varying degrees of stability across different random seeds [START_REF] Cárdenas-Montes | Sensitiveness of evolutionary algorithms to the random number generator[END_REF] The finding that different algorithm configurations can be recommended, depending on the seed provided to LON construction

BACKGROUND 2.1 Fitness Landscapes

A fitness landscape [START_REF] Peter | Fitness landscapes. Biological Evolution and Statistical Physics[END_REF] is composed of three parts: (𝑆, 𝑁 , 𝑓 ) : S is the full set of possible solutions; 𝑁 : 𝑆 -→ 2 𝑆 is the neighbourhood function, which assigns a set of adjacent solutions 𝑁 (𝑠) to every 𝑠 ∈ 𝑆; and f is a fitness function 𝑓 : 𝑆 -→ R that provides a mapping from solution to associated fitness. That fitness can be conceptualised as the solution height within the landscape metaphor.

The local optima network (LON) model [START_REF] Ochoa | A study of NK landscapes' basins and local optima networks[END_REF] was introduced as a tool for studying the connectivity of local optima in a fitness landscape, and has subsequently shown proficiency in helping with explaining metaheuristic search behaviour [START_REF] Bożejko | Local optima networks in solving algorithm selection problem for tsp[END_REF][START_REF] Herrmann | Communities of local optima as funnels in fitness landscapes[END_REF]. In this work, we consider LONs which are sampled using iterated local search. We now formally define this type of LON.

Local optima. We assume a search space 𝑆 with a fitness function 𝑓 (𝑆) and a neighbourhood function 𝑁 (𝑠). A local optimum, which in the QAP is a minimum, is a solution 𝑙 such that ∀𝑠 ∈ 𝑁 (𝑙), 𝑓 (𝑙) ≤ 𝑓 (𝑠).

Monotonic perturbation edges.

There is an edge from local optimum 𝑙 1 to local optimum 𝑙 2 , if 𝑙 2 can be obtained after applying a random perturbation (in the present study, this is 𝑘-exchange) to 𝑙 1 followed by local search, and 𝑓 (𝑙 2 ) ≤ 𝑓 (𝑙 1 ). In LON terminology, these are called escape edges [START_REF] Verel | Local optima networks with escape edges[END_REF]. The edges are called monotonic because they record only non-deteriorating, directed transitions between local optima. Edges are weighted with the frequency of transition. The set of edges is denoted by 𝐸.

Compound local optima.

A compound local optimum is a set of connected LON nodes with the same fitness value. Two nodes are connected if there is a monotonic perturbation edge between them. The set of connected optima with the same fitness, denoted by 𝐶𝐿, corresponds to the set of nodes in the compound monotonic LON model.

Compound monotonic LON. Is the directed network CMLON = (CL, CE), where nodes are compounded local optima 𝐶𝐿, and the edges 𝐶𝐸 are aggregated from the monotonic edge set 𝐸 by summing up the edge weights. Monotonic sequence. A monotonic sequence is a path of connected nodes 𝑀𝑆 = {𝑐𝑙 1 , 𝑐𝑙 2 , . . . , 𝑐𝑙 𝑠 } where 𝑐𝑙 𝑖 ∈ 𝐶𝐿. By definition of the edges, 𝑓 (𝑐𝑙 𝑖 ) ≤ 𝑓 (𝑐𝑙 𝑖 -1 ). There is a natural end to every monotonic sequence, 𝑐𝑙 𝑠 , when no improving transitions can be found. This node, 𝑐𝑙 𝑠 , is called a sink as it does not have improving outgoing edges 1 The global optimum will always be a sink, and there can be sub-optimal sinks as well. Note that 𝑐𝑙 𝑠 can also be a compounded node; that is, the combination of more than one connected nodes of equivalent fitness.

Funnel. A funnel in the CMLON comprises the aggregation of all monotonic sequences ending at the same sink (which can be a compound node) [START_REF] Ochoa | Understanding phase transitions with local optima networks: number partitioning as a case study[END_REF] and is a subset of the CMLON. Conceptually, funnels can be viewed as being basins of attraction at the level of local optima. The funnel whose sink is the global optimum can be called the "optimal" funnel (there can be more than one optimal funnel); all other funnels are sub-optimal.

LON Construction

For a few years after local optima networks were proposed, exhaustive enumeration of the complete network was standard [START_REF] Ochoa | A study of NK landscapes' basins and local optima networks[END_REF][START_REF] Ochoa | Local optima networks: A new model of combinatorial fitness landscapes[END_REF][START_REF] Verel | Local optima networks of NK landscapes with neutrality[END_REF]; indeed, some later works also took this approach [START_REF] Ochoa | Understanding phase transitions with local optima networks: number partitioning as a case study[END_REF][START_REF] Potgieter | A Local Optima Network Analysis of the Feedforward Neural Architecture Space[END_REF][START_REF] Sarah L Thomson | Comparing communities of optima with funnels in combinatorial fitness landscapes[END_REF]. In these works, the node set is constructed by mapping every solution to its local optimum through local search. Edges are then defined according one of two approaches: basin transition edges [START_REF] Ochoa | A study of NK landscapes' basins and local optima networks[END_REF], which capture the probability of moving between two basins of attraction under the chosen neighbourhood operator, and escape edges [START_REF] Verel | Local optima networks with escape edges[END_REF]: these connect nodes if the destination local optimum can be reached from the source by carrying out a defined operation (usually perturbation) followed by local search. The nature of exhaustive enumeration necessarily constrains the studied instances to small sizes. In response to this, sampling methodologies have been proposed to construct partial local optima networks; these implement the escape edge model. Algorithms include iterated local search (ILS) [START_REF] Ochoa | Perturbation strength and the global structure of QAP fitness landscapes[END_REF]. ILS is a metaheuristic which repeatedly combines large and random mutations (perturbations, for the purpose of diversification) with local search (intensification); it navigates through the local optima space, and is therefore suited for extracting a LON. Other LON sampling approaches include snowball sampling [START_REF] Verel | Sampling local optima networks of large combinatorial search spaces: The qap case[END_REF], recombination approaches [START_REF] Melike D Karatas | Towards population-based fitness landscape analysis using local optima networks[END_REF][START_REF] Whitley | Local optima organize into lattices under recombination: an example using the traveling salesman problem[END_REF], and local search based methods [START_REF] Bożejko | Local optima networks in solving algorithm selection problem for tsp[END_REF][START_REF] Iclanzan | Data-driven local optima network characterization of QAPLIB instances[END_REF]. Literature has suggested that snowball-sampled LON metrics correlate well with those of a complete LON [START_REF] Verel | Sampling local optima networks of large combinatorial search spaces: The qap case[END_REF]. For LON sampling by ILS, the number of nodes -and to a lesser degree, the number of edges -were found to be related to the exhaustively-enumerated quantity [START_REF] Sarah L Thomson | Inferring future landscapes: sampling the local optima level[END_REF]. The issue with this kind of comparative analysis is that it is necessarily constrained to problem sizes small enough that LONs can be completely enumerated; it is not clear how the perceived correlations scale to larger problem dimension.

In terms of connection between LONs and their predictive potency, authors have shown that LONs sampled using ILS seem to be able to explain search performance better than enumerated LONs [START_REF] Sarah L Thomson | Inferring future landscapes: sampling the local optima level[END_REF]. This might be in line with the argument from network sampling literature that sampling algorithms can provide a kind of "positive" bias [START_REF] Arun | Benefits of bias: Towards better characterization of network sampling[END_REF], in that they may encourage the samples to contain representative information. In general, questions surrounding sampling bias and the effect of randomness in LON construction algorithms are still open. Iterated local search based sampling seems to be the most prevalent in the literature [START_REF] Ochoa | Neural Architecture Search: A Visual Analysis[END_REF][START_REF] Marcondes Pavelski | Local Optima Network Sampling for Permutation Flowshop[END_REF][START_REF] Tari | Local search pivoting rules and the landscape global structure[END_REF][START_REF] Sarah L Thomson | The fractal geometry of fitness landscapes at the local optima level[END_REF][START_REF] Veerapen | Modelling genetic improvement landscapes with local optima networks[END_REF]; we therefore focus our attention on it for the purposes of this study.

Randomness and Sampling

Statistical randomness, otherwise known as pseudorandomness, is obtained from a pseudorandom number generator [START_REF] Maurice | Randomness and random sampling numbers[END_REF]. This behaviour does not imply or guarantee true randomness. We consider the effect of statistical randomness in this study, and refer to it with the term "randomness" for the purpose of concision. In statistical language, a partial local optima network is essentially a sample from a larger population: the exhaustively-enumerated sets of nodes and edges. For random samples, it has been argued that 30-50 are sufficient [START_REF] Christopher Z Mooney | Bootstrapping: A nonparametric approach to statistical inference[END_REF]; this aligns with the Central Limit Theorem, which stipulates that from around 30 observations, sample means begin to resemble a normal distribution [START_REF] Gyu | Central limit theorem: the cornerstone of modern statistics[END_REF].

Literature has indicated that evolutionary algorithm performance can be affected by the choice of random number generator [START_REF] Cárdenas-Montes | Sensitiveness of evolutionary algorithms to the random number generator[END_REF][START_REF] Krömer | Randomness and chaos in genetic algorithms and differential evolution[END_REF][START_REF] Maucher | Search heuristics and the influence of non-perfect randomness: examining Genetic Algorithms and Simulated Annealing[END_REF][START_REF] Mark | Randomness and GA performance, revisited[END_REF]. It has been argued that it is possible to remove randomness from EAs while maintaining algorithm proficiency [START_REF] Zelinka | Do evolutionary algorithms indeed require random numbers? Extended study[END_REF]. Some authors have proposed a "quasi-random" approach as an alternative to random population instantiation [START_REF] Heikki Maaranen | Quasi-random initial population for genetic algorithms[END_REF]; a quasi-random generator places points in a space such that the points are maximally distant from each another. The construction of sampled LONs is influenced by randomness in the heuristic algorithms used to identify the nodes and edges. Authors use varying numbers of independent algorithm runs for iterated local search LON construction to try and mitigate the effects of randomness. There seems to be no agreed value: as few as 30-50 runs have been employed [START_REF] Ochoa | Neural Architecture Search: A Visual Analysis[END_REF][START_REF] Marcondes Pavelski | Local Optima Network Sampling for Permutation Flowshop[END_REF][START_REF] Tari | Local search pivoting rules and the landscape global structure[END_REF]; other works use in the region of 100-200 [START_REF] Ochoa | Local optima network analysis for MAX-SAT[END_REF][START_REF] Sarah L Thomson | The fractal geometry of fitness landscapes at the local optima level[END_REF][START_REF] Sarah L Thomson | Inferring future landscapes: sampling the local optima level[END_REF][START_REF] Treimun-Costa | Modelling parameter configuration spaces with local optima networks[END_REF]; as many as 1000 have also been considered [START_REF] Veerapen | Modelling genetic improvement landscapes with local optima networks[END_REF]. This lack of agreement is one reason why an assessment for the influence of randomness on LON construction is important.

The random mechanism which is used during LON construction in this work is the pseudorandom number generator, Mersenne Twister (MT) [START_REF] Matsumoto | Mersenne twister: a 623dimensionally equidistributed uniform pseudo-random number generator[END_REF]. The quality of a pseudorandom number generator can be defined by its period, 𝑝 -how many numbers before the generator repeats itself. Mersenne Twister has 𝑝=2 19337 -1. Randomness is found in a few places in the LON construction algorithm. The most consequential of these may be the starting solution, but there is also the order of neighbour exploration in the first-improvement local search and the random perturbations to consider.

METHODOLOGY 3.1 Quadratic Assignment Problem

Definition. A solution to the QAP is generally written as a permutation 𝑠 of the set {1, 2, ..., 𝑛}, where 𝑠 𝑖 gives the location of item 𝑖. Therefore, the search space is of size 𝑛!. The cost associated with a permutation 𝑠 is a quadratic function of the distances between the locations, and the flow between the facilities, 𝑓 (𝑠) = 𝑛 𝑖=1 𝑛 𝑗=1 𝑎 𝑖 𝑗 𝑏 𝑠 𝑖 𝑠 𝑗 , where 𝑛 denotes the number of facilities/locations and 𝐴 = {𝑎 𝑖 𝑗 } and 𝐵 = {𝑏 𝑖 𝑗 } are the distance and flow matrices, respectively.

Instances. We consider instances from the QAPLIB2 [START_REF] Burkard | QAPLIB -A Quadratic Assignment Problem Library[END_REF] with between 25 and 50 facilities; these are of moderate size, and yet are not always trivial to solve. There are 47 instances in this size range, but we remove the esc group (which are real-world QAP occurrences for the testing of sequential circuits in computer hardware): their local optima networks are uninteresting to study because there is a very high degree of LON neutrality [START_REF] Sarah | On funnel depths and acceptance criteria in stochastic local search[END_REF]. The reduction results in 40 QAPLIB instances for the study, ranging in size between 25 and 50 facilities and locations. Eleven of the instances in this group have not been solved to optimality; for those, we use their best-known fitness as the surrogate global optimum. In the rest of this paper, for simplicity we refer to these as the global optimum.

LON Metrics

Funnel-like organisation of local optima has been shown to be related to search difficulty for the travelling salesman problem [START_REF] D R Hains | Revisiting the big valley search space structure in the TSP[END_REF], number partitioning problem [START_REF] Ochoa | Understanding phase transitions with local optima networks: number partitioning as a case study[END_REF] and also for the QAP [START_REF] Sarah | On funnel depths and acceptance criteria in stochastic local search[END_REF]. We therefore concentrate our focus towards measurements relating to funnels. Unless otherwise specified, all metrics are computed on the CMLON. To calculate funnel depth measurements, we first compute all finite shortest paths in the CMLON from origin nodes (local optima at the beginning of ILS runs) towards sink nodes (termination points of ILS runs). A sink is the bottom of a funnel. This process is completed separately for optimal and suboptimal sinks. From the resultant sets of shortest paths (which are monotonic in nature) we can extract metrics. For both optimal and suboptimal funnels, we consider the mean depth and the maximum depth of them; that is, the mean and maximum of the shortest paths which reach the end of the funnel from an origin point.

Other LON metrics are included too: the incoming strength to optimal sinks, and the relative size, in nodes, of the optimal funnel(s). The former is calculated as the weighted in-degree to optimal sinks in CMLON -as a proportion of the total weighted incoming degree to all sinks. The size of optimal funnel(s) is defined to be the number of nodes which can reach the optimal sinks using a monotonic path; the size is normalised by the total number of nodes. Also considered are the number of global optima and local optima, respectively, in the un-compounded LON; the number of nodes in the CMLON; the mean size of compounded nodes; and the cardinality for optimal and sub-optimal funnels. Abbreviated names for LON features are used in the results section; we therefore provide a list with shortnames and descriptions now, in the order they first appear in a Figure .   (1) cnodes_size_avg: the mean size of compounded LON nodes (nodes which are a combination of more than one connected nodes of equal fitness) (2) depth_gfunnel_avg: the mean depth (in edges) of optimal funnels (3) depth_gfunnel_max: the maximum depth (in edges) of optimal funnels (4) depth_sofunnel_avg: the mean depth (in edges) of sub-optimal funnels (5) depth_sofunnel_max: the maximum depth (in edges) of suboptimal funnels (6) num_cnodes: number of compounded nodes [START_REF] Iclanzan | Data-driven local optima network characterization of QAPLIB instances[END_REF] num_gfunnel: number of optimal funnels (8) num_global: number of global optima (9) num_nodes: number of local optima [START_REF] Knowles | Instance generators and test suites for the multiobjective quadratic assignment problem[END_REF] num_sofunnel: number of sub-optimal funnels (11) size_gfunnel: proportional size, in nodes, of optimal funnel(s) (12) strength_gfunnel: proportional strength of optimal funnel(s), measured as weight of incoming edges to the funnel sink(s)

EXPERIMENTAL SETUP

Throughout the experiments, a random swap of facilities in the QAP permutation is considered as the neighbourhood structure. All experiments are run on an 8-core MacBook Pro equipped with the Apple M1 chip and using OS X Ventura. We use Stützle's iterated local search (ILS) for both gathering performance data and as the foundation of LON construction [START_REF] Stützle | Iterated local search for the quadratic assignment problem[END_REF]. The local search stage uses a first improvement hill-climbing, and the perturbation strength is 𝑘 = 𝑁 2 . This setting was selected because large perturbation strengths of around 50% of 𝑁 are generally advantageous in this problem domain [START_REF] Ochoa | Perturbation strength and the global structure of QAP fitness landscapes[END_REF][START_REF] Stützle | Iterated local search for the quadratic assignment problem[END_REF].

LON Construction

For each instance, 30 LONs are constructed -with a different random seed provided to the Mersenne Twister (MT) mechanism each time. The C implementation of MT is used. For each instance and seed pair, 100 runs of ILS are amalgamated to form the LON. Only local optima with equal or improving fitness are accepted. Worsening local optima are never accepted. Runs terminate after 10 000 ILS iterations with no improvement. This termination condition has been used in the literature to estimate sinks and funnel structure [START_REF] Ochoa | The multi-funnel structure of TSP fitness landscapes: a visual exploration[END_REF].

Algorithm Performance

We consider runtime over 100 independent executions (with different random seeds from those used for the LON construction runs) to summarise ILS performance on the instances. Each of the three ILS variants have the same seeds provided to them. Runs terminate when either the known best fitness is found or after 10 000 ILS iterations with no improvement. To give an indication of how many iterations this may result in: the instance tai50a was associated with 12 178 iterations in an example run. Performance metrics were computed using three distinct local optima acceptance strategies for the ILS, detailed now.

Better or equal. Only local optima which have improved or equal fitness to the current local optimum are accepted. Worsening local optima are never accepted.

Annealing. Improving and equal local optima are always accepted. Worsening local optima are accepted according to probability. This is governed by an annealing-like cooling schedule termed Large Step Markov Chain (LSMC) in the work proposing ILS for QAP. We set the LSMC parameters to mirror those used in Stützle's work [START_REF] Stützle | Iterated local search for the quadratic assignment problem[END_REF]: that is, an initial temperature of 2.5% of the current solution fitness -a solution which is 2.5% lower fitness is allowed with the probability 1 𝑒 ; the number of steps at a temperature is 10. Temperature is reduced according to a scaling factor, 𝛼 = 0.9.

Restart. Improving and equal local optima are always accepted. If iterations since an improvement in local optima quality have exceeded 3𝑁 (where 𝑁 is the problem dimension), a restart happens.

Regression Model Setup

LON metrics serve as the independent variables in the algorithm performance models. The dependent variable is runtime of iterated local search; the modelling is therefore a regression task. As the learning algorithm, Random Forest in R is used. For each individual model, the observations include metrics from LONs constructed using one of the 30 random seeds, combined with ILS runtime for the problem instance as the dependent variable. To be clear, the same ILS performance data is used regardless of the LON random seed: for example, the runtime value for bur26a is mapped to all 30 of its LONs. Because we use one row per QAP instance to build the models, there are limited number of observations in the datasets: 40. The one-in-ten rule [START_REF] Harrell | Regression modelling strategies for improved prognostic prediction[END_REF] stipulates that roughly ten observations are required per independent variable. Owing to data splitting during cross-validation, our training sets can be as small as 32. We compute 12 LON metrics and then conduct recursive feature elimination to obtain the top three features. The feature selection bases its decisions on the minimisation of RMSE and executes 100 bootstrap replications. One hundred bootstraps is a quantity backed up by literature [START_REF] Nicholas D Pattengale | How many bootstrap replicates are necessary[END_REF]. Feature selection is conducted separately for each of the 30 set of LONs (one set per seed), and the elite features stored. For the final modelling stage, features which appeared the most at the top three positions were included. Each feature selection process includes 100 bootstrapping iterations as cross-validation.

For both feature selection and final modelling, default hyperparameters for the Random Forests are employed; these are: 500 trees, 𝑛𝑝 3 variables sampled [where 𝑛𝑝 is the number of predictors], subsampling of observations with replacement, and a sample size of 𝑛 [where 𝑛 is the number of the whole sample]. We conduct Monte Carlo cross-validation (random repeated sub-sampling) to estimate model quality metrics. This is done with an 80-20 trainingvalidation split and for 100 iterations; we consider the 𝑅 2 and root mean squared error (RMSE) as model metrics. These are computed on the validation sets and averaged over 100 cross-validation iterations.

RESULTS

Network Comparison

Figure 1 shows, for each QAP instance, the distribution of LON features across different random seeds provided to the LON construction process. There are 30 seeds for each instance. Each horizontal line on a box is for a particular feature, as indicated on the left of the overall Figure . Recall the descriptions for features from Section 3.2. Colours represent the normalised value of the feature: dark is closer to zero, with yellow closer to one. Instance names are the labels for individual boxes.

Notice from the Figure that particular features, such as the metrics size_gfunnel and strength_gfunnel, are often consistent across seeds. This can be observed with the solid blocks of either yellow or dark blue at the bottom of several of the instance boxes. Other features which display little variation for many of the instances include size_cnodes_avg, num_global, and num_sofunnel. There are, however, features which show less stability across seeds. We notice that metrics relating to funnel depth, particularly depth_gfunnel_avg and depth_gfunnel_max, often depend on the seed. This can be seen by checking for the rows labelled with those metric names and observing the variation in colour there.

With all that said, there seems to be no universal rule about which metrics are dependable across seeds. Even size_gfunnel and strength_gfunnel display lack of uniformity in some instances (for example, kra30a and kra30b); we therefore stipulate that stability of metrics depends on the problem instance.

Figure 2 presents distributions of LON metrics across seeds as parallel coordinate plots. Notice that the 𝑥-axis now measures the normalised value for features. Each line traces feature values for a particular seed. Colour is used to identify an instance as belonging to an instance group, as labelled in the legend. Surveying the Figure ,  we can see that the eight instances comprising the bur category are each associated with LONs which are remarkably constant across different seeds. This can be seen by considering the tightness of the line groups in these instances when compared to instances from other instance types. In the tainna category, there is a lot of variation in metrics relating to depth of the optimal funnel(s): depth_gfunnel_avg and depth_gfunnel_max. The number of local optima, as well as features which capture neutrality at the local optima level -cnodes_size_avg and num_cnodes -show stability across seeds almost universally for problem instances (note the lowest part of the individual boxes). The tainnb category each display the same "skeleton" of feature patterns. The first four instances have extremely similar plots, while tai50b shows more diversity of funnel depth metrics across different seeds.

To take a closer look into possible differences between QAP instance groups, we show the distributions for selected features as boxplots in Figures 3 and4. The diagrams represent the features num_nodes and depth_gfunnel_avg, respectively. Those two were selected because one of them displays relative consistency across groups (and across random seeds), while the other shows different behaviour: the amount of variation across seeds seems to be somewhat related to instance group. Due to space constraints, we cannot show these plots for every feature. tured by LON sampling is usually quite consistent even with differing random seeds. We can see this from the narrow distributions throughout most of the plot. These are present in every instance group, although some instances within real have more local optima and slightly more varied distributions. Turning our attention to Figure 4, the pattern is unlike that of Figure 3. Here, there are several distributions which exhibit high variation. In some cases, this seems like it might be related to instance group: compare grid and tainna to the other categories. For those two, the sampled depth of optimal funnel is often highly diverse depending on seed. For the real instances, however, the depth is quite stable across seeds.

Feature Selection

For feature selection, the set of observations used to build a model comprise metrics computed on a set of 40 LONs -one per QAP instance -constructed using a particular random seed. There are therefore sets 30 models, differentiated by the seed provided to the LONs. Each of the three ILS variants serve as a separate response variable, so each has an associated group of 30 model configurations. We investigate in this Section how consistent the top-selected feature was for these models across the different random seeds. Figure 5, presents the number of times particular features ranked first place for the three response variables: ILS with default, annealinglike, and restart acceptance. Notice that some features from the full set (described in Section 3.2) do not appear in Figure 5; this means that they did not appear as the top feature in any of the models.

Comparing the distributions in the three sub-plots, we can see that the top feature for the default and annealing ILS variants was rather consistent across different LON seeds. In the case of default, the strength_gfunnel feature ranked as first for 17 of the 30 possible seeds (observe the lowest bar of the left-most plot). The most-common top feature for annealing ILS, depth_gfunnel_avg, is even more dominant across seeds: around 22 out of 30. With restart ILS as response variable, it is a different story: a few different features have several occurrences, but none prevail as being distinctly the most salient. These results suggest that when it pertains to explaining and predicting default and annealing ILS, the LON features which are most potent are consistent across random seeds provided to the LON construction algorithm. This does not appear to be the case for restart ILS, however.

Modelling

The final models contain the top three selected features for each ILS variant. Different LON seeds sometimes resulted in different features ranking as the best. We took a pooling approach to feature choice for the final models: selecting the features which most often ranked at places one, two, and three across the 30 seeds.

For default acceptance, these were strength_gfunnel, size_gfunnel, and depth_sofunnel_avg; in the case of annealing ILS, they were depth_gfunnel_avg, num_cnodes, and depth_sofunnel_avg; and for restart acceptance: strength_gfunnel, depth_gfunnel_avg, and size_gfunnel. Figures 6 and7 present the distributions across 30 LON seeds of model 𝑅 2 and RMSE, respectively, on validation sets -averaged over 100 Monte Carlo cross-validation iterations. One outlier has been removed for visual clarity: a model in the restart set which had an 𝑅 2 of -13.44. For random forest regression, a negative 𝑅 2 like this can occur when the model prediction is less accurate than simply predicting the mean of the response variable. Looking at Figure 6 and comparing the three boxes, it can be observed that models with default ILS as response variable have the narrowest range of 𝑅 2 across seeds. Even so, there are values as far apart as ∼0.5 and ∼0.97. This is a substantial difference. The ranges for annealing and restart are even more vast; this is particularly the case for the latter. Depending on the random seed provided to the LON construction algorithm, the subsequent model 𝑅 2 can vary from below zero -indicating a model which is worse than predicting the mean response variable value -to ∼0.96 (almost all variance explained). Recalling the variation of top feature across seeds during the feature selection process (Figure 5, it is possible that the low 𝑅 2 values arise because the most powerful metrics of a LON constructed with a specific seed were not used (we used the same three features to build models regardless of seed). 4: For each instance, the distribution of the average depth of optimal funnel(s) in the LON across different random seeds provided to LON construction narrowest distribution, arises from models with default ILS as the response variable. This communicates that these models are the strongest, and also that the dispersion of model error is relatively low across LON seeds. The annealing and restart models have larger errors and wider distributions. The consistency of model error depending on LON seeds seems quite low.

Algorithm Configuration.

We now report an algorithm configuration decision experiment. The intention was to investigate whether different recommendations might be made depending on the seed provided to LON construction.

For all three ILS variants, we set their predictor features to be the best three for the default ILS. This decision was made for simplicity and consistency. Thereafter, we obtain the predicted runtimes for each QAP instance and for the three algorithm configurations. These are produced by executing 100 iterations of Monte Carlo cross-validation, with the relevant instance always excluded from the training set and included in the validation. The models predict the runtime. Then each instance is mapped to the ILS variant with associated with the lowest predicted runtime.

In Table 1 we report statistics describing the consistency of algorithm configuration recommended across different seeds: the column labelled #config. contains values between 0 and 1 -the number of recommended algorithm configurations as a proportion of the available number (three). In column freq. prevailing config. the proportions are out of 30, and represent for how many seeds the most prevalent algorithm configuration was recommended. Let us consider the #config. column initially. The minimum proportion is 0.667, which means that at least two algorithm configurations were recommended for each instance. Both the 3rd quartile and the maximum have proportional values of 1, which tells us that some instances had all three ILS variants recommended. In the freq. prevailing config. column, notice that the minimum proportional representation of the dominant recommended configuration is 0.5. The maximum is 0.967, which represents 29 out of the 30 seeds agreeing on the best algorithm variant. A mean and median of 0.700 and 0.633 respectively indicate that, for most instances, the dominant variant is recommended for a reasonable majority of LON seeds.

CONCLUSIONS

We have conducted a preliminary study on statistical randomness in a local optima network (LON) construction algorithm: iterated local search (ILS). LONs were extracted for a set of quadratic assignment problem (QAP) instances and with a range of different random seeds provided to the construction algorithm. Metrics were computed on the obtained LONs. These were analysed with the aid of visualisation. Algorithm performance prediction models were built using the metrics as features and performance of ILS as the response variable. The results showed that some LON metrics seem consistent across seeds, while others vary substantially. Additionally, the quality of algorithm performance models using LON metrics as predictors can differ depending on randomness. Finally, LON metrics associated with separate seeds can lead to different algorithm configuration recommendations for the same instance. For future work, we will consider a larger QAP instance set. In particular, an instance generator (such as [START_REF] Knowles | Instance generators and test suites for the multiobjective quadratic assignment problem[END_REF]) could be used to obtain an expansive set of similar problems. This should remove the variation seen between instances from different classes. We would also like to further remove noise from the performance modelling pipeline; consider the study of other problems, including pseudoboolean spaces; and include other heuristics or alternative LON sampling methods such as snowballing [START_REF] Verel | Sampling local optima networks of large combinatorial search spaces: The qap case[END_REF] and those relating to recombination [START_REF] Melike D Karatas | Towards population-based fitness landscape analysis using local optima networks[END_REF][START_REF] Whitley | Local optima organize into lattices under recombination: an example using the traveling salesman problem[END_REF]. The data from this work are available at https://doi.org/10.5281/zenodo.7900303.

  tai35a

Figure 1 :

 1 Figure 1: For each instance, the distribution of LON metrics across different random seeds provided to LON construction

Figure 2 :

 2 Figure 2: For each instance, the parallel coordinate representation of the distribution of LON metrics across different random seeds provided to LON construction

Figure 7 contains

 7 RMSE values. RMSE is in the same units as the response variable, which is ILS runtime, measured in iterations. The Figure indicates that the lowest RMSE values, and indeed the

Figure 3 :

 3 Figure 3: For each instance, the distribution of the number of local optima in the LON across different random seeds provided to LON construction

Figure

  Figure 4: For each instance, the distribution of the average depth of optimal funnel(s) in the LON across different random seeds provided to LON construction

Figure 5 :

 5 Figure 5: Distribution across 30 seeds of the top feature selected with each ILS variant as response variable

Figure 6 :Figure 7 :

 67 Figure 6: Distribution across 30 LON seeds of model 𝑅 2 on validation data (mean over 100 cross-validation iterations) with three ILS variants as the response variable. Predictors are the top three LON metrics obtained from feature selection

Table 1 :

 1 Algorithm configuration distributions: #config. captures how many different ILS variants were recommended for the instance, depending on LON seed (as a proportion of three total variants); freq. prevailing config. is how often the most common recommendation for an instance appears (presented as a proportion of 30 models)

	statistic	#config. freq. prevailing config.
	minimum	0.667	0.500
	1st quartile	0.667	0.567
	median	0.667	0.633
	mean	0.825	0.700
	3rd quartile	1.000	0.833
	maximum	1.000	0.967

In directed networks, a node without outgoing edges is called a 𝑠𝑖𝑛𝑘.

https://coral.ise.lehigh.edu/data-sets/qaplib/qaplib-problem-instances-andsolutions/