domain: spi info

text:
Alexandre Jacquemart
email: alexandre.jacquemart@esi-group.com

Sio-Song Ieng
email: sio-song.ieng@univ-eiffel.fr

Mokrane Hadj-Bachir
email: mokrane.hadj-bachir@esi-group.com

Abdelaziz Bensrhair
email: abdelaziz.bensrhair@insa-rouen.fr

Dominique Gruyer
email: dominique.gruyer@univ-eiffel.fr

3-GWD : A Textural and AI-Based Approach for Real-Time Detection of Weather Disturbances Applied to Autonomous Vehicle

Keywords: weather detection, fog detection, rain detection, co-occurrence matrix, autonomous vehicle, artificial intelligence, ADAS, perception, safety, neural network I

We propose 3-GWD, a new and innovative solution for the particle atmospheric disturbers' detection that can be used to weight the data in a fusion process and adjust the sensor reliability. The proposed approach, based on 3D Gray-Level Co-occurrence Matrices (3D-GLCMs), is generic and could be used for a significant number of sensors (LiDAR, IR camera, Neuromorphic camera...). We demonstrate that this method leads to the identification of very relevant patterns in the symbolic space that shallow Artificial Neural Network (ANN) can classify with a great accuracy. We also discuss the importance of all-environment dataset by creating a realityenriched simulated database with controlled disturbances and numerous environments. Using a single optic camera produces promising results with a high prediction score. We end our study with an attempt to compare our approach with previous works using the DAWN and the MCWCD datasets.

of degraded conditions management and assessment led to research on visibility under degraded weather conditions. Visibility problems induced by these conditions are often caused by the presence of clouds of particles with diverse physical and optical properties. Most of these phenomena deteriorate visibility with the presence of fog. As a result, many research works are focused on fog detection or removal [START_REF] Miclea | Visibility enhancement and fog detection: Solutions presented in recent scientific papers with potential for application to mobile systems[END_REF].

The theoretical works on the characterization of disturbers is exposed in papers such as [START_REF] Koschmieder | Theorie der horizontalen sichtweite[END_REF], [START_REF] Narasimhan | Vision and atmosphere[END_REF] and [START_REF] Garg | When does a camera see rain?[END_REF]. Other works focus on pattern recognition or image processing [START_REF] Vis | Analysis of rain and snow in frequency space[END_REF], [START_REF] Bossu | Rain or Snow Detection in Image Sequences through use of a Histogram of Orientation of Streaks[END_REF]. Finally, a third approach angle consists in removing or mitigating these phenomena on images [START_REF] Garg | When does a camera see rain?[END_REF], [START_REF] He | Single image haze removal using dark channel prior[END_REF], [START_REF] Sun | Icyclegan: Single image dehazing based on iterative dehazing model and cyclegan[END_REF].

Nevertheless, existing works do not yet allow a sufficiently reliable usage in the CAV, and algorithms still need to be developed to assess in real-time the best combination of embedded sensor technologies and the best topology of these sensors. The choice of specific technologies could enable the vehicle to optimize the selection of models and algorithms to use, both increasing reliability and security.

AI approaches on challenging weather have already been proposed, whether it be for image restoration, weather forecast, or labelling for search browser retrieval. All those approaches rely on structures on the image that are relevant of particular weather. [START_REF] Vis | Analysis of rain and snow in frequency space[END_REF] for example, tries to detect rain and snow using a Fourier decomposition, [START_REF] Allamano | Toward the camera rain gauge[END_REF] tries to count droplets and their apparent size to estimate current rainfall amount. Other methods, reviewed in [START_REF] Miclea | Visibility enhancement and fog detection: Solutions presented in recent scientific papers with potential for application to mobile systems[END_REF], use image processing to try and catch relevant clues of weather disturbances' presence: an overview of the methods used reveals a common usage of the Fourier decomposition. Other approaches exist : [START_REF] Asery | Fog detection using glcm based features and svm[END_REF], as an example, uses an approach based on the graylevel co-occurrence matrix (GLCM).

We can draw similarities with texture characterization : Fourier Decomposition is well-known for this purpose, and the GLCM was first described in [START_REF] Haralick | Textural features for image classification[END_REF] as a method of texture classification, alongside derived coefficients.

Our approach is based on the 3D-GLCM [START_REF] Philips | Directional invariance of cooccurrence matrices within the liver[END_REF] [17] algorithm. We found that, if used on an optic RGB camera, meaningful features in symbolic space can be found. We also prove that these features are both sensible and specific to the studied disturbances. Our algorithm is also lightweight, as we believe a security-assessment algorithm should not have a major impact on available resources for critical modules, e.g., detection algorithm or optimal path estimation. The proposed method could also be applied to a whole range of sensors, making further research in this field interesting.

To describe and justify our idea, and then to train our network, we simulated, using the PROSIVIC [START_REF] Gruyer | From virtual to real, how to prototype, test, evaluate and validate adas for the automated and connected vehicle?[END_REF] software, three scenarios : a countryside one, a forest road one, and an urban one. To record the data, we also simulated a vehicle with an embedded camera. Due to the simulation, we can precisely control the parameters of the disturbances. That allowed us to conduct a parametric study on an identical simulation with only the weather disturbance changing. We then added real data from [19] [20] [21] to confirm that the simulation emulates the reality, as well as enriching our neural network with new features.

In this paper, we will start by listing the already developed methods for texture characterization and weather detection. Then, we will present the database we generated and assembled. We will follow by the 3D-GLCM algorithm, and we will show that some patterns on the 3D-GLCM are present or missing depending on the current disturbers. Before assessing our method, we will describe the neuronal classifier used. Finally, we will describe and exploit the results and demonstrate the sensibility and specificity of our method on the studied disturbances.

II. RELATED WORK

A. Texture characterization

There are several ways to describe an image using texture features, which can be used for image segmentation and image classification. [START_REF] Cavalin | A review of texture classification methods and databases[END_REF] divides these approaches into 4 families : Statistical, Structural, Filter-based and Modelbased. For this paper, we will focus on the Statistical family because their methods can easily be abstracted to multidimensional space.

Statistical texture characterization approaches are made up of many algorithms. One of those is the Gray-Level Cooccurrence Matrix (GLCM). A GLCM is defined as a count of couples of pixels separated by a given step. Because of this, there are countless GLCMs for any image. Nevertheless, short steps are usually used, as longer ones do not ensure that we are still in the same image region. The principle was first described in [START_REF] Haralick | Textural features for image classification[END_REF]. This paper also shows an example of its use : the computation of derived coefficients for describing an image. Among the most common are entropy, dissimilarity, and correlation. Another commonly used algorithm is Local Binary Patterns (LBP) [START_REF] Ojala | Performance evaluation of texture measures with classification based on kullback discrimination of distributions[END_REF], based on binarization of the immediate neighborhood and the extraction of an integer coded by the succession of the resulting 1s and 0s.

As for our method, we considered the statistical texture characterization for their explainability and their fast adaptation to different types of data. More precisely, we chose the GLCM, as LBP was already studied, and we found derived methods expanding it to 3 dimensions.

B. Weather detection algorithms

Existing methods for detecting weather conditions can be divided into 2 groups.

The first group is composed of "data driven" algorithms. The principle of these algorithms is to learn the relevant features from the data using a Deep Neural Network (DNN) architecture, usually a CNN. However, this approach suffers from a lack of explainability, that often makes them qualified as "black-box" algorithms. [START_REF] Al-Haija | Detection in adverse weather conditions for autonomous vehicles via deep learning[END_REF] for example, uses already existing architecture, like SqueezeNet, ResNet-50 and EfficientNet CNN layers to extract features, then a single fully connected layer with Softmax activation function to classify the data. [START_REF] Samo | Deep learning with attention mechanisms for road weather detection[END_REF] uses CNNs layers from well-known architectures and transformers to achieve multi-label weather classification. The work presented in [START_REF] Dhananjaya | Weather and light level classification for autonomous driving: Dataset, baseline and active learning[END_REF] is another case of CNN-based weather detection algorithm using the Resnet18 framework. However, this method uses active learning and auto-labelling, which are other ways of training than the classic epoch and batch algorithm.

The second type of algorithm is composed of "model driven" algorithms. These methods rely on deterministic features search. Then, a machine learning classifier is trained and used to predict the label of the current weather situation. In [START_REF] Bossu | Rain or Snow Detection in Image Sequences through use of a Histogram of Orientation of Streaks[END_REF], the method consists in a background subtraction for moving pixels disclosure and a combination of photometric and size rules alongside a Histogram of Oriented Streaks to detect the rain. Another method is proposed in [START_REF] Chu | Camera as weather sensor: Estimating weather information from single images[END_REF] with 543 image characteristic and 2 metadata : month and hour. With a database of the average month by month temperatures on European cities, they achieved better results than with the image characteristics alone. Using a texture-based approach, [START_REF] Khan | Snow detection using in-vehicle video camera with texture-based image features utilizing k-nearest neighbor, support vector machine, and random forest[END_REF] uses contrast, correlation, energy, and homogeneity to extract key features from GLCMs. Features obtained by LBP are also used. The merged features are then fed into a nonneuronal classifier, either a Support Vector Machine (SVM), a k-Nearest Neighbors (k-NN) algorithm or a Random Forest (RF). Lastly, the previously discussed LBP algorithm is used in [START_REF] Khan | Real-time weather detection system with local binary pattern based features using artificial neural network and random forest: An unsupervised learning approach[END_REF] with various machine-learning algorithms to achieve rain detection.

Following our choice of the GLCM for texture characterization, we will follow the model-driven approach family. As such, we will maintain a satisfying level of explainability, however special care must be taken regarding performances.

III. DATASET GENERATION

Regardless of the aimed application: classification, detection, or tracking; all machine learning algorithms require a significant amount of data. It is then essential for obtaining the expected results that this database is of high quality. Therefore, we built a hybrid dataset (simulated and realworld data) with three requirements for the detection of different meteorological disturbances:

1) The characteristics of weather disturbers in the dataset must be configurable according to physical models.

2) The generated data must represent a wide range of values covering the field of application of each type of disturber (different densities of rain drops and fog visibility).

3) The data must present real features to ensure the reliability in reality execution. To meet the first and the second requirements, we have generated a dataset using the PROSIVIC software ( [START_REF] Gruyer | From virtual to real, how to prototype, test, evaluate and validate adas for the automated and connected vehicle?[END_REF]), a road scene and sensor simulator. This simulator enables us to modify both the weather conditions and its parameters, such as particle size or density [START_REF] Gruyer | Modeling and validation of a new generic virtual optical sensor for adas prototyping[END_REF]. This allows us to generate the same scenes with different weather conditions.

The majority of datasets involve urban environment, as the quantity of information to process is a real challenge. However, as stated by the Society of Automotive Engineers, level 4 and 5 AVs will need an all-weather, all environment operation range. So, in accordance with point 2, an emphasis was placed on the diversity of generated environments:

• An urban environment with average distances of view and an occluded horizon line. • A forest road, with close to medium distances of view and an occluded horizon line. Unlike buildings of the city, trees are more scattered with far more contrast. • The last scenes are country roads with long-distance views and a visible horizon line.

Fig. 1 illustrates some data we can generate using PRO-SIVIC. Koschmieder's model [START_REF] Koschmieder | Theorie der horizontalen sichtweite[END_REF] is used for fog rendering, and rain is emulated droplet by droplet. For the latter, the parameters are the number of droplets n in the volume, the dimensions of the box where rain droplets are generated, and the opacity and diameter d of said droplets. These disturbers' model are perfectly fitted to answer point 1. of our requirements. Finally, to fulfil the third requirement, we included in the dataset real data from BeDDE [START_REF] Zhao | Dehazing evaluation: Real-world benchmark datasets, criteria, and baselines[END_REF], O-HAZE [START_REF] Ancuti | Ohaze: A dehazing benchmark with real hazy and haze-free outdoor images[END_REF] and RainCityScape [START_REF] Hu | Depth-attentional features for single-image rain removal[END_REF]. These datasets are well-known and used as reference point for the detection of targets under challenging weather. The weather information included in these datasets are shown in Tab. I.

TABLE I: Weather disturbances given by real datasets

Name

Clear Weather Fog Rain BeDDE [START_REF] Zhao | Dehazing evaluation: Real-world benchmark datasets, criteria, and baselines[END_REF] Yes Yes No O-HAZE [START_REF] Ancuti | Ohaze: A dehazing benchmark with real hazy and haze-free outdoor images[END_REF] Yes Yes No RainCityScape [START_REF] Hu | Depth-attentional features for single-image rain removal[END_REF] Yes No Yes

The new hybrid database then considers 2 types of data:

• The simulated data (the vast majority of the dataset) with 3 weather conditions (Clear, Fog, Rain). These 3 weather conditions take place in 3 environments: City, countryside, and forest. • The real-world data with 3 weather conditions (Clear, Fog, and rain). The clear weather condition is obtained from the 3 following datasets: BeDDE [START_REF] Zhao | Dehazing evaluation: Real-world benchmark datasets, criteria, and baselines[END_REF], O-HAZE [START_REF] Ancuti | Ohaze: A dehazing benchmark with real hazy and haze-free outdoor images[END_REF], RainCityScape [START_REF] Hu | Depth-attentional features for single-image rain removal[END_REF]. The real fog conditions are taken in 2 datasets: BeDDE and O-HAZE. About the rain condition, the RainCityScape one is used. Each simulated weather class is given a parametric set for quality assessment and feature characterization.

IV. THE TEXTURAL APPROACH

We decide to tackle the problem of weather conditions detection by adapting the 3D approach proposed in [START_REF] Philips | Directional invariance of cooccurrence matrices within the liver[END_REF] to RGB images. The initial method was proposed to process Computed Tomography (CT) scans from a voxel map. The new method starts from the following hypothesis: if the interesting patterns were in a 2D space, it could also be possible to find them in a 3D space that includes this 2D space. Moreover, for a given sensor C, if the data can be mapped into a 3D discrete-finite array, then this approach could be used in the same way. As a result, this approach can be applied not only to RGB images but also to data coming from a large set of sensors with different technologies.

Let M be an RGB image considered as a 3D discrete matrix. Let N be the cardinal of the set of values M , i.e, 256 for optic images. The size of M is (n×m×3). We define our reference frame ρ as (line, column, channel) where :

     0 <= line < n 0 <= column < m 0 <= channel < 3 (1)
(0,0,0) represents the top-left pixel of the shallowest channel.

Let us define the neighborhood Q p of the pixel of interest p (POI). Q p will have to respect 3 rules:

1) ∀q ∈ Q p , L2(p, q) = i (p i -q i ) 2 < t. Let's say that t = 2. This creates our raw neighborhood.

2) The coordinates of any q ∈ Q p must satisfy constraints in Eq. 1. This rule ensure that we don't access out of bound data. 3) Let Q ′ p be the set obtained by mirroring Q p indexes to p, then Q p ∩ Q ′ p = ∅ This equation removes the symmetry from the neighborhood, as described in [START_REF] Haralick | Textural features for image classification[END_REF]. The set we kept, answering those requirements, is shown in the Fig. 2. There are 4 neighboring pixels in the same channel and 9 others in the channels G and R respectively. If p is in the channel R, it only has 4 neighboring pixels in the same channel since others does not satisfy Eq. 1. We did not have represented borders and corners neighborhood, but in these particular cases, rule 2 will reduce Q p 's cardinal.

Let s ∈ Z 3 be the direction vector from p to q ∈ Q p ; the set of all steps S is therefore of size 13. Let M(x) be the value of pixel x. GLCM s are computed following the procedure described in [START_REF] Haralick | Textural features for image classification[END_REF] below: for all pixel P as POI of M: do for all q s ∈ Q p : do GLCM s (M (p), M (q)) ← GLCM s (M (p), M (q))+ 1 GLCM s (M (q), M (p)) ← GLCM s (M (q), M (p))+ 1 end for end for The result of this algorithm is a set of 13 256 × 256-integer GLCMs. Using our parametric dataset, we generated for each image its corresponding 3D-GLCMs, then summed them all together [START_REF] Kovalev | Threedimensional texture analysis of mri brain datasets[END_REF]. The results in Figs. 345shows for each class of weather (Clear, Rain, Fog) specific patterns :

• Clear : The clear weather characterizes itself by the presence of numerous distinct areas of high cardinalities, inside the 3D-GLCM. These little areas, e.g., near (192, 68) in Fig. 3, are due to contrast. • Fog : The foggy weather characterizes itself by the super dominance of the diagonal. This is a due to contrast loss, resulting in the matrix shrinking towards its diagonal. If the depth allows it, we also see dense areas in zones corresponding to the different veils used for our datasets. • Rain : The rainy weather characterizes itself by a dominance of its 3D-GLCM diagonal, with a larger spread than fog. That corresponds to the diffusion caused by the raindrop marks on the image : on these regions, a blur is done, caused by all the optic phenomena a changing of medium implies. We also observe a peak activity in the area where the droplet's light intensity is.

V. OUR NEURONAL CLASSIFIER

An advantage of the GLCM in general is that we can predict its size. For example, for a RBG optic image, the matrix will always be 256×256, since there are 256 distinct values for each pixel. We've also shown that 3D-GLCM behave a certain way depending on the weather studied. Detecting these patterns was done using a single logistical Fig. 3: Original image on clear weather and its associated 3D-GLCM. The threshold allows us to highlight patterns that have medium values, contrasting with some huge diagonal cells. We can notice the presence of little islands of cardinalities, signs of highly contrasted areas in the image. The curves on the bottom right are due to the progressive RGB shift in the sky. Fig. 4: Same frame as in Fig. 3 with rain added and its associated 3D-GLCM. The threshold allows us to notice the loss of the islands that were highlighted in Fig. 3. We can notice the loss of the bottom-right curves, since the sky changed. The droplets streak, bright but translucent, explain the surge of the bottom right activity. Fig. 5: Same frame as in Fig. 3 with heavy fog added and its associated 3D-GLCM. The threshold allows us to notice the loss of the islands that were highlighted in Fig. 3. We also notice a threshold in the 3D-GLCM range due to the dim light of the situation. Due to the contrast loss induced by the weather condition, we notice a compression of the matrix towards its diagonal. perceptron for each class. Let c be the number of classes we want to classify, i.e., the weather situations including the clear weather. Then Fig. 6 is the schematics of our classifier. The fact to use a single stack of perceptron (or 1 layer of fully connected neurons) was already used in [START_REF] Khan | Real-time weather detection system with local binary pattern based features using artificial neural network and random forest: An unsupervised learning approach[END_REF] alongside Softmax activation function. Fig. 6: Schematic of the classifier network. 3D-GLCM is first flattened, then put into a dense layer. Finally, the classifier classifies the presence probability of each weather disturbance using the logistic function. The result is a vector of probability presence.

We do not need CNNs before our classifier, since the feature extraction of those layer is already done in a deterministic and comprehensible way by the 3D-GLCM. A simple perceptron is used, as we need to characterize the relevant areas, and a single weighting should work, given that we demonstrated that there exist areas that looked sensible or specific. The logistical function then follows our perceptron.

By reshaping the weights vector into the original 3D-GLCM matrix shape, we can visualize the relation between these weights and the 3D-GLCM components and thus how the inputs will activate the neurons. We think that this visualization allows a deep understanding of the classifier and this is due to the simple architecture of the classifier. This is a crucial key that deeper network loses in benefit for better space segmentation.

VI. EXPERIMENTS

We trained our model, 3D-GLCMs Weather Detector (3-GWD), on a Dell Precision 7540 equipped with an Intel Core i9-9980HK CPU, 32 GB of RAM and an NVIDIA Quadro T1000.

The training phase was carried out using a L2-regularized Binary-Cross Entropy (BCE). Let y and ŷ be the target (c values in 0; 1) and the prediction (c values in [0; 1]), then BCE is defined as :

BCE(y, ŷ) = - 1 n c i=1 y i •log( ŷi )+(1-y i )•log(1-ŷi ) (2)
The L2 Loss is defined by :

L2(y, ŷ) = c i=1 (y i -ŷi ) 2 (3)
Combining ( 2) and ( 3) we obtain our loss function :

Loss(y, ŷ) = BCE(y, ŷ) + λ • L2(y, ŷ) (4)
Where λ is the regularization hyperparameter. Alongside the L2-regularised BCE Loss, we use a learning rate (lr) scheduler. More specially, we use the "reduce lr on plateau" strategy, which states that lr is multiplied by a learning rate multiplier (0 < lrm < 1) when the training metric has stopped improving.

For our validation, we will use the F 1 -Score, an indicator that measures the quality of the model using its precision and recall.

Precision is the measure of relevance in the positive prediction. Using True Positives (TP), i.e., the correct positive predictions, and False Positives (FP), i.e., incorrect positive predictions, Precision is defined as :

P recision = T P T P + F P (5)
On the other hand, Recall is the measure of retrieval of positive prediction. Using TP and False Negatives (FN), i.e., incorrect negative predictions, Recall is defined as :

Recall = T P T P + F N (6)
Using ( 5) and ( 6), the F 1 -Score is defined as :

F 1 = 2 Recall -1 + P recision -1 (7)

A. Test on our dataset

To test our model, we first trained and validated on our previously described dataset. The architecture was trained using the following hyperparameters : lr = 0.1, lrm = 0.24, batch size = 256, epochs = 6, λ = 0.4. Figs. 789represent the confusion matrices for the different part of the validation set. We can see that real data are more subject to confusion between clear (expected) and rain (true value) while PROSIVIC data are more prone to the opposite. The F 1 score, Precision and Recall measures, are computed for each part of the dataset and summarized in Tab. II. We also pointed out in Sec. V that explainability was a key point of our method. To demonstrate this, we will use the following visualization tool, hereby called "RGB stacking". First, let W be the reshaped weight matrix. Then, we apply a threshold to W between -t and t. t = 3 * std as statistically, value outside this interval can be considered extreme. Then, we will apply Eq. 8 to W , resulting in W G ∈ [0; 1] 256×256 . The G in W G means that the matrix values are in [0; 1], and as such can be considered Gray-scaled.

W G = W + t 2 × t (8)
Let C G , F G and R G be the processed matrix with "Clear weather", "Fog", and "Rain" neuron weights. Let I be a

(3 × 256 × 256)matrix such that I[0] = C G , I[1] = F G and I[2] = R G .
I is then a matrix which values are in [0; 1] and, as such, is alike an RGB image, with R channel related to "Clear Weather" neuron weights, G channel to "Fog" and B channel to "Rain". As such, if we plot I as an image, we can set visual explanation of different areas in the 3D-GLCM.

• Red : Red areas represent areas where only the "Clear Fig. 9: Confusion matrix on real data from validation set Weather" neuron have high positive weights. These are then statically sensible to clear weather.

• Green: Same explanation as Red but for sensibility to Fog • Blue : Same explanation as Red but for sensibility to Rain • Yellow (Red + Green) : Yellow areas represent area where "Clear weather" and "Fog" neurons gives credit.

On the opposite, that means the Blue channel is close to 0. As such, "Rain" neuron is the only one that refuse a positive classification. As such, these are statically specific to "Rain". • Magenta (Red + Blue) : Same explanation as Yellow but for specificity to Fog • Cyan (Green+Blue) : Same explanation as Yellow but specificity to Clear weather Fig. 10 is the RGB staking of our obtained model. We immediately notice the presence of Red, Green, and Blue areas : our architecture is then statistically sensible to each condition alone. We also can see small areas of Yellow near (60,60), Magenta areas in bottom right and dim cyan area away from the diagonal. Therefore, we can then conclude that our architecture is statistically specific, but we have to nuance, since those areas are dim or small.

B. Tests on DAWN+MCWCD datasets

To validate the study about the proposed classifier, 3-GWD, and in order to quantify its performances relatively to existing classifiers, a comparison process has been carried out by using two annotated datasets DAWN [START_REF] Kenk | Dawn: Vehicle detection in adverse weather nature dataset[END_REF] and MCWCD [START_REF] Oluwafemi | Multi-class weather classification from still image using said ensemble method[END_REF]. In these datasets, the authors of [START_REF] Al-Haija | Detection in adverse weather conditions for autonomous vehicles via deep learning[END_REF] applied data augmentation process to extend the base by using cropping, rotation, translation, resizing among others.

Unfortunately, after obtaining the results of the comparison shown in Table III, it appeared that these datasets were TABLE III: Accuracy comparison on the dataset as given by [START_REF] Al-Haija | Detection in adverse weather conditions for autonomous vehicles via deep learning[END_REF] Name Accuracy Technology SqueezeNet [START_REF] Al-Haija | Detection in adverse weather conditions for autonomous vehicles via deep learning[END_REF] 96.05% CNN (4,8M+ param.) ResNet-50 [START_REF] Al-Haija | Detection in adverse weather conditions for autonomous vehicles via deep learning[END_REF] 98.50% CNN+DNN (25M+ param.) EfficientNet-b0 [START_REF] Al-Haija | Detection in adverse weather conditions for autonomous vehicles via deep learning[END_REF] 97.78% CNN+DNN (5,6M+ param.) [START_REF] Oluwafemi | Multi-class weather classification from still image using said ensemble method[END_REF] 86.50% LBP, HSV, Grad. + Rand. Forest 3-GWD 83.32% GLCM + Percept. (196K param.) Fig. 11: Sample of images from MCWCD that have no relevance for our use unsuitable and inappropriate to our domain of application (Road scenes seen from an embedded camera). Indeed, a detailed analysis of the datasets show important defects in the type of images and annotations provided in these datasets.

• In the MCWCD [START_REF] Oluwafemi | Multi-class weather classification from still image using said ensemble method[END_REF] dataset, the data includes very poor quality images and some include watermarks and no road scenes. Many scenes represent the beach or the sea. The data augmentation procedure proposed by [START_REF] Al-Haija | Detection in adverse weather conditions for autonomous vehicles via deep learning[END_REF] brings eventually more confusion especially when the generated data have the identical label to the original image. For instance, some cropped and resize images from an image of beach with sun represent only water, but they are annotated as "clear weather". These images should not be generated for the weather classification. These examples of training set used in [START_REF] Al-Haija | Detection in adverse weather conditions for autonomous vehicles via deep learning[END_REF] were not in the training set of the 3-GWD because they are not road scenes or urban scenes. • In the DAWN [START_REF] Kenk | Dawn: Vehicle detection in adverse weather nature dataset[END_REF] dataset, the annotations provided are clearly ambiguous and sometimes wrong relatively to the classical definitions of rain, fog, snow, dust and sand. Indeed, many examples of the dataset represent highway scenes with wet ground, traffic vehicles raising water drops. For this type of data, we labelled these images with "rain". On the other hand, in the used datasets, the labels of these are mostly "fog". We think this limitation is due to the specificity of this kind of weather that induces behaviors that can be both classified as "rain" and "fog", whereas the dataset is a single-label one.

VII. CONCLUSION

With the acceleration on the development of AI-based target detection architecture, the need for external and internal validation algorithms is vital for real-time and after-crash risk assessment. Fig. 12: Sample of images from DAWN that are symptomatic of the label problem. In the dataset, they are labelled "Fog" while we would classify them as "Rain"

To address this issue, research on visibility under degraded weather conditions has led to the development of several intelligent weather detection systems, including the proposed 3-GWD. This system is a GLCM-based (Grey Level Cooccurrence Matrix) weather detection system designed to support AVs. It has been trained and evaluated on a custom hybrid dataset involving both simulated and real data, and tested on three different road situations under three controlled weather disturbances.

The GLCM space has proven to be useful for feature extraction for weather detection, and a set of rules for generating the hybrid dataset has been proposed. The proposed method has been evaluated using various metrics, including precision, recall and F1-score, all of which have shown promising results.

To ensure explainability, a visualization tool called RGBstacking has been developed, which allows users to understand how the classifier segments the symbolic space of the GLCMs. This visualization tool helps to ensure that the results of the 3-GWD system are easily interpretable.

Future studies will aim to extend the 3-GWD system to new sensor technologies, such as LiDAR, IR, and eventbased cameras. These studies will confirm the presence of significant patterns inside the 3D-GLCM and help to identify the best sensors for detecting weather disturbances in different situations.

Another important study is the comparison of the 3-GWD system with existing works, such as [START_REF] Al-Haija | Detection in adverse weather conditions for autonomous vehicles via deep learning[END_REF], which provides a comprehensive survey and evaluation of many methods over the DAWN [START_REF] Kenk | Dawn: Vehicle detection in adverse weather nature dataset[END_REF] and MCWCD [START_REF] Oluwafemi | Multi-class weather classification from still image using said ensemble method[END_REF] datasets. However, our work is not conclusive due to the problems mentioned in the section VI-B but underlined the importance of the quality of the data and their annotation in Machine Learning approach. As such, a new comparison on the relabelled dataset have to be conducted.

In conclusion, the development of intelligent weather detection systems such as the 3-GWD system is a crucial step in the development of AVs that can operate safely and efficiently in all weather conditions. The 3-GWD system, with its innovative approach to feature extraction, hybrid dataset generation, and visualization tools, shows great potential in this regard and will undoubtedly inspire further research in this exciting field.

Fig. 1 :

 1 Fig. 1: Sample of the simulated dataset on the 3 different scenes : the city, the countryside, and forest.

Fig. 2 :

 2 Fig. 2: Valid neighbors for each channel of an RGB image, depending on the channel. Gray pixel is the Pixel of Interest (POI) called p, lighter pixels around represents valid neighborhood Q and darker pixels are incorrect neighbors

Fig. 7 :Fig. 8 :

 78 Fig. 7: Confusion matrix on all validation set

Fig

  Fig. RGB stacking of the trained classifier

TABLE II :

 II Precision, Recall and F1-score on the different dataset parts

	Part	Precision	Recall	F 1 -Score
	Whole dataset	0.9800	0.9801	0.9800
	Simulated images	0.9930	0.9931	0.9930
	Real images	0.9670	0.9670	0.9670

This work has been funded by ANRT organism and ESI group.


Matthieu Dolbeault
email: dolbeault@igpm.rwth-aachen.de

Moulay Abdellah Chkifa
email: abdellah.chkifa@um6p.ma

Randomized least-squares with minimal oversampling and interpolation in general spaces

Keywords: Least squares, Interpolation, Christoffel function MSC 2020. 65D15, 41A65, 65D05, 41A81, 41A05

In approximation of functions based on point values, least-squares methods provide more stability than interpolation, at the expense of increasing the sampling budget. We show that near-optimal approximation error can nevertheless be achieved, in an expected L 2 sense, as soon as the sample size m is larger than the dimension n of the approximation space by a constant ratio. On the other hand, for m = n, we obtain an interpolation strategy with a stability factor of order n. The proposed sampling algorithms are greedy procedures based on [BSS09] and [LS18], with polynomial computational complexity.

Introduction and main results

Let (X, A, µ) be a probability space. We consider the problem of estimating an unknown function f : X → C from evaluations of f at chosen points x 1 , . . . , x m ∈ X. We assess the error between f and its estimator f either in the L 2 (X, µ) norm

g L 2 := X |g(x)| 2 dµ(x) 1/2
, or in the uniform norm g L ∞ (X,µ) . Given a subspace V n ⊂ L 2 (X, µ) of dimension n, we would like the estimator f to belong to V n and to perform almost as well as the best approximation of f in V n , that is, its orthogonal projection

P n f = arg min v∈Vn f -v L 2 .
As we only have access to point-wise observations, we cannot explicitly compute P n f in general. In this context, a classical approach consists in considering a solution to the weighted least-squares problem

P m n f ∈ arg min v∈Vn m i=1 s i |f (x i ) -v(x i )| 2 ,
where we may use some weights s 1 , . . . , s m > 0. We would like this problem to admit a unique solution, and therefore require m n. Similar to P n for the L 2 norm, the operator P m n is the orthogonal projector onto V n with respect to the empirical norm

g m := m i=1 s i |g(x i )| 2 1/2 .
The approximation accuracy is inherently related to the points x i and the weights s i . For the sake of illustration, consider the setting where V n = P n-1 is the space of algebraic polynomials of degree less than n, restricted to the interval X = [-1, 1], and choose m = n, so that P n n is the Lagrange interpolation operator associated with {x 1 , . . . , x n }. For equally spaced points x i , this corresponds to interpolation on a uniform grid, which is known to be highly unstable, failing to converge towards f , even when f is infinitely smooth. This is the so-called Runge phenomenon, see e.g. [START_REF] Mastroianni | Interpolation Processes: Basic Theory and Applications[END_REF]. The phenomenon persists with equally spaced points even for m = rn with r > 1 constant, as observed in [START_REF] Boyd | Divergence (Runge phenomenon) for least-squares polynomial approximation on an equispaced grid and Mock-Chebyshev subset interpolation[END_REF] and theoretically explained in [START_REF] Platte | Impossibility of fast stable approximation of analytic functions from equispaced samples[END_REF]. In fact, the Runge phenomenon occurs for any points that do not cluster quadratically like Chebyshev points, see [START_REF] Adcock | Optimal sampling rates for approximating analytic functions from pointwise samples[END_REF]. It is however defeated by interpolation on Chebyshev type points.

As far as interpolation in spaces V n is concerned, there exists no systematic choice of points that prevents all instabilities. A generic choice is that of Fekete points, which ensure that

P n n g L ∞ n g L ∞ in the strong uniform norm g L ∞ = sup x∈X |g(x)| 1 , resulting in the stability inequality f -P n n f L ∞ (n + 1) inf v∈Vn f -v L ∞ ,
see for instance Proposition 1.2.5 in [START_REF] Novak | Deterministic and stochastic error bounds in numerical analysis[END_REF]. This guarantees a good approximation P n n f of f , provided the latter is sufficiently smooth and V n is well chosen. However, for general domains X and spaces V n , the computation of Fekete points can be intractable.

For polynomial interpolation over compact domains X ⊂ R or C, Leja sequences are greedy alternatives to Fekete points. For instance, for interpolation by polynomials in P n-1 , restricted to union of closed intervals of R, they provably yield P n n g L ∞ n 13/4 g L ∞ , see the recent paper [START_REF] Andrievskii | A simple upper bound for Lebesgue constants associated with Leja points on the real line[END_REF]. We note that numerical evidence shows that n 13/4 can be replaced by n for intervals. Similar polynomial growths, with a factor n 3 in the estimate, hold for a certain type of polynomial interpolation over tensor product domains, see [START_REF] Chkifa | High-dimensional adaptive sparse polynomial interpolation and applications to parametric PDEs[END_REF][START_REF] Chkifa | On the stability of polynomial interpolation using hierarchical sampling[END_REF].

On the other hand, near optimal approximation error in expected L 2 sense

E( f -P m n f 2 L 2 ) C inf v∈Vn f -v 2 L 2 ,
can be attained by taking m larger than n. In the case of uniformly distributed points and equal weights s i = 1, this usually requires m to scale polynomially in n, see [CDL13, CCM + 15], but a logarithmic oversampling is achievable if one considers a different sample distribution [START_REF] Cohen | Optimal weighted least-squares methods[END_REF].

Numerical methods [START_REF] Haberstich | Boosted optimal weighted least-squares[END_REF] and a theoretical solution [START_REF] Cohen | Optimal pointwise sampling for L 2 approximation[END_REF] have been proposed to reduce the sample size m to a constant multiple of n. A discrete version of the above bound can also be found in [START_REF] Chen | Active regression via linear-sample sparsification[END_REF], in the context of statistical machine learning.

In the last five papers, the sample points are drawn at random according to a prescribed measure, and the error bounds are presented in expectation. This setting does not require any additional assumption: indeed, point evaluations of a function f ∈ L 2 are defined in an almost sure sense, and so is the weighted least-squares projection P m n f . This will also be the main framework for the present article.

Our main theorem, stated below, provides new bounds on the L 2 approximation error, depending on the ratio between m and n.

Theorem 1.1. Let m n.

• The conditional weighted least-squares approximation f ∈ V n defined in (15), using m evaluations of f at points selected by Algorithm 3.1, with inputs ε = r -1/4 and γ = r 1/2 -r 1/4 , satisfies

E( f -f 2 L 2 ) 1 + 1 r(1 -r -1/4 ) 7 min v∈Vn f -v 2 L 2 , (1)
where r = m+1 n > 1.

• In turn, the weighted least-squares estimator f = P m n f , based on points selected by Algorithm 4.1, with inputs δ = r -1/2 and κ ∈ [0, 1], simultaneously satisfies

E( f -f 2 L 2 ) 1 + 1 1 -κ 1 (1 -1/ √ r) 2 min v∈Vn f -v 2 L 2
(2)

if κ < 1 and f -f L 2 1 + 1 √ κ 1 1 -1/ √ r min v∈Vn f -v L ∞ a.s. ( 3
) if κ > 0.
Here r = m n-1 > 1, and we assume that V n ⊂ L ∞ in the last bound.

The first two bounds are of the same kind, the main difference being in the constant factor on the right-hand side, which is of optimal order 1 + 1/r + o(1/r) in (1) as r tends to infinity, but behaves better in (2) when r gets close to 1. The parameter κ is a free input parameter selected by the user. The choice κ = 0 gives the best constant in (2), while κ = 1 is best suited for (3). One can also construct an estimator f = P m n f satisfying both estimates at the same time, by running Algorithm 4.1 with an intermediate value of κ.

The uniform bound (3) follows the approach developed in [LT22, Tem21, PU22, DT24, BSU23], slightly improving the constants when compared to Theorem 1.1 in [START_REF] Temlyakov | On optimal recovery in L 2[END_REF] and Theorem 6.3 in [START_REF] Bartel | Constructive subsampling of finite frames with applications in optimal function recovery[END_REF], and linking it to the approach in expectation (2) by the use of a common algorithm. In a third approach, similar deterministic bounds can be proved by assuming more regularity on f through a nested sequence of approximation spaces (V n ) n∈N , in a Hilbert space setting [KU21a, KUV21, MU21, NSU22, BSU23, GW24] or in more general Banach spaces [START_REF] Krieg | Function values are enough for L 2 -approximation: Part II[END_REF][START_REF] Dolbeault | A sharp upper bound for sampling numbers in L 2[END_REF][START_REF] Krieg | Sampling recovery in the uniform norm[END_REF].

Taking m = n in the last two estimates, and observing that 1/ √ r = 1 -1/n 1 -1/2n, we immediately obtain the following result on interpolation in L 2 . Corollary 1.2. For m = n, the interpolation f = P n n f ∈ V n of f at random points x 1 , . . . , x n selected by Algorithm 4.1 with input δ = 1 -1/n achieves the accuracy bounds

E( f -f 2 L 2 ) (4n 2 + 1) min v∈Vn f -v 2 L 2 if κ = 0,
and, assuming

V n ⊂ L ∞ , f -f L 2 (2n + 1) min v∈Vn f -v L ∞ a.s. if κ = 1.
Remark 1.3. In the case of uniform approximation, if each function in V n is uniformly bounded on X and not just essentially bounded, one can replace L ∞ by the strong supremum norm L ∞ , and remove the almost sure restriction, by considering deterministic samples satisfying the constraints of our algorithms. A similar observation can be found in Remark 3 of [START_REF] Krieg | Sampling recovery in the uniform norm[END_REF].

Finally, one can combine the second estimate in Corollary 1.2 with an inverse inequality between L 2 and L ∞ in V n .

Corollary 1.4. If the functions in V n are bounded and X is compact, Algorithm 4.1, applied with the optimal measure µ * from [KW60], provides n points for which the interpolation

f * = P n n f ∈ V n of f satisfies f -f * L ∞ (1 + 2n √ n) min v∈Vn f -v L ∞ .
Note that in the case m = n, the values of the weights s i have no importance, since the minimum in the definition of P m n is zero. In fact, if measure µ * is known, we exhibit a constructive set of points such that the Lebesgue stability constant

L n = max f ∈C 0 (X)\{0} P n n f L ∞ f L ∞ is at most 2n √ n.
Although Fekete points achieve L n n, their computational complexity is exponential in n. In our approach, the main computational challenge is to find an optimal measure µ * as defined in [START_REF] Kiefer | The equivalence of two extremum problems[END_REF], see also [START_REF] Bos | Some remarks on the Fejér problem for Lagrange interpolation in several variables[END_REF], which is in general as difficult as looking for the Fekete points.

However, using a different measure µ, we still obtain a result similar to Corollary 1.4 however with the factor √ n replaced by a larger power of n. For example, if X = [-1, 1] d and V n = span{x ν : ν ∈ Λ} is a space of multivariate polynomials indexed by a lower set Λ ⊂ N d 0 , the uniform measure yields a factor n instead of √ n and an estimate L n 2n 2 , while the tensor product arcsine measure achieves a factor min(2 d/2 √ n, n log 3/ log 4 ) instead of √ n, resulting in the estimate L n 2n min(2 d/2 √ n, n log 3/ log 4 ). Note that the Lebesgue constant does not only determine the convergence of the approximation, but also its robustness to numerical errors and noise in the measurements, see Section 4 in [START_REF] Cohen | Optimal weighted least-squares methods[END_REF], as well as [START_REF] Adcock | Optimal sampling rates for approximating analytic functions from pointwise samples[END_REF][START_REF] Platte | Impossibility of fast stable approximation of analytic functions from equispaced samples[END_REF].

The above discussion demonstrates that good interpolation points can be constructed for multivariate polynomial approximation. The estimate 2n 1+log 3/ log 4 on the Lebesgue constant is moderate. It outperforms the estimate n 3 which was established in [START_REF] Chkifa | On the stability of polynomial interpolation using hierarchical sampling[END_REF] using constructions based on -Leja sequences. We note however that the new procedure is not hierarchical.

Remark 1.5. One can also use an inverse inequality between L 2 and L ∞ in the least-squares regime m > n. More precisely, combining (3) with [START_REF] Kiefer | The equivalence of two extremum problems[END_REF] as in Corollary 1.4, we obtain

f -P m n f L ∞ 1 + √ n 1 -1/ √ r min v∈Vn f -v L ∞ ,
with r = m/(n -1). Taking the supremum for f in a compact class of functions, and optimizing on both sides over V n , this implies that the uniform sampling numbers are bounded by O( √ n) times the uniform Kolmogorov n-widths, improving on the O(n) factor of Fekete points, if we allow for a constant oversampling r > 1. This can already be seen by applying Corollary 5.3 in [START_REF] Pozharska | A note on sampling recovery of multivariate functions in the uniform norm[END_REF] together with the optimal density from [START_REF] Kiefer | The equivalence of two extremum problems[END_REF]. The same result is obtained in Theorem 3 of [START_REF] Krieg | Sampling recovery in the uniform norm[END_REF], where error bounds are also derived in any L p norm. We additionally refer to [START_REF] Geng | On the power of standard information for tractability for L ∞ approximation of periodic functions in the worst case setting[END_REF], where implications in the field of Information Based Complexity are drawn, in a specific Hilbert space setting. The two very recent papers [START_REF] Krieg | Sampling recovery in the uniform norm[END_REF][START_REF] Geng | On the power of standard information for tractability for L ∞ approximation of periodic functions in the worst case setting[END_REF] rely on the pioneering works [START_REF] Krieg | Function values are enough for L 2 -approximation[END_REF][START_REF] Krieg | Function values are enough for L 2 -approximation: Part II[END_REF], and on the infinite-dimensional adaptation [START_REF] Dolbeault | A sharp upper bound for sampling numbers in L 2[END_REF] of the result from [START_REF] Marcus | Interlacing families II: Mixed characteristic polynomials and the Kadison-Singer problem[END_REF], see also [START_REF] Freeman | The discretization problem for continuous frames[END_REF].

The rest of the paper is organized as follows. In Section 2, we inspect the weighted leastsquares projection, and point out why different strategies should be used when r 1 and r ≈ 1. Sections 3 and 4 introduce and analyze the sampling Algorithms 3.1 and 4.1. They are independent from one another, up to the shared use of a few formulas, and the various estimates of Theorem 1.1 are proved separately. Finally, we discuss some numerical aspects of the presented algorithms in Section 5, and provide numerical illustrations in Section 6.

Least-squares

Let ϕ = (ϕ 1 , . . . , ϕ n ) be an orthonormal basis of V n in L 2 (X, µ). For any x ∈ X, we consider ϕ(x) as a vector in C n , denote ϕ(x) * its conjugate transpose, and |ϕ(x)| 2 = ϕ(x) * ϕ(x) its squared euclidian norm. Observe that, given any n × n matrix M ,

X ϕ(x) * M ϕ(x)dµ(x) = Tr M X ϕ(x)ϕ(x) * dµ(x) = Tr(M ). ( 4
)
Adopting the formalism from [KUV21, MU21, NSU22]), we let f = ( √ s i f (x i )) 1 i m ∈ C m be the measurement vector, and

L = ( √ s i ϕ j (x i )) 1 i m,1 j n ∈ C m×n
the collocation matrix. Then the weighted least-square estimator writes

P m n f = n j=1 a j ϕ j , a = L + f ,
where L + = (L * L) -1 L * stands for the Moore-Penrose pseudo-inverse of L. In the sequel, we will make sure that L has full column rank n, so that the Gram matrix

A m := L * L = m i=1 s i ϕ(x i )ϕ(x i ) *
is indeed invertible. Denoting g = f -P n f the optimal residual error and g = ( √ s i g(x i )) 1 i m the associated vector, the least-squares error decomposes as

f -P m n f 2 L 2 = f -P n f 2 L 2 + P n f -P m n f 2 L 2 = g 2 L 2 + P m n g 2 L 2 . ( 5
)
There are two possible strategies for bounding P m n g 2 L 2 : either we use

P m n g 2 L 2 = |A -1 m L * g| 2 A -1 m 2 2→2 |L * g| 2 = λ min (A m ) -2 | ϕ, g m | 2 , (6)
where ϕ, g m := ( ϕ j , g m ) 1 j n ∈ C n , or we bound it by

P m n g 2 L 2 = |L + g| 2 L + 2 2→2 |g| 2 = A -1 m 2→2 g 2 m = λ min (A m ) -1 g 2 m . (7)
The first approach is expected to give better estimates when m is much larger than n, since in that case the discrete inner product •, • m weakly converges to the continuous inner product

•, • L 2 , and as g is orthogonal to V n , the vector ϕ, g m should be small. On the other hand, when m is close to n, L may be ill-conditioned, leading to small values of λ min (A m ), thus favouring the second approach.

In both situations, one should choose the points x i and weights s i in order to control the smallest eigenvalue of A m from below. This can be ensured through a greedy selection of points, based on the effective resistance [START_REF] Spielman | Graph sparsification by effective resistances[END_REF] ϕ(x) * (A -I) -1 ϕ(x) of a point x ∈ X with respect to a hermitian matrix A, with a lower bound on the eigenvalues of A, and I the identity matrix. Intuitively, if the points x 1 , . . . , x i-1 are already fixed, the effective resistance of x i with respect to the partial Gram matrix

A i-1 = i-1 ι=1 s ι ϕ(x ι )ϕ(x ι ) *
quantifies how close ϕ(x i ) is to the eigenvectors of A i-1 with small eigenvalues, helping the lower potential Tr((A -I) -1 ) to decrease when going from A i-1 to A i . The key idea in [START_REF] Batson | Twice-ramanujan sparsifiers[END_REF] is to increase the lower barrier at each step, without increasing the potential more than it had decreased when adding sample x i . Keeping the potential bounded ensures that all updates of are of the same size δ, and its trace form is chosen to allow rank-one matrices s i ϕ(x i )ϕ(x i ) * to decrease the potential, independently of the distribution of eigenvalues of A i-1 . In this way, we obtain the desired bound λ min (A m ) > m at the end of the algorithm.

As in [START_REF] Batson | Twice-ramanujan sparsifiers[END_REF], it is also possible to ensure an upper bound λ max (A m ) < u m on the eigenvalues of A m , by considering the upper potential Tr((uI -A) -1 ). Although λ max (A m ) does not explicitly appear in the least-squares formulation, it is worth noticing that the ratio u m / m controls the condition number λ max (A m )/λ min (A m ) of A m , which determines the computational cost of solving the least-squares problem via an iterative solver, as well as its robustness to numerical error. This is discussed in, e.g., Section 5.3.4 of [START_REF] Adcock | Sparse polynomial approximation of high-dimensional functions[END_REF].

Nevertheless, we chose not to include this upper potential, because it would have prevented us from achieving the critical sampling budget m ≈ n. Heuristically, when we use both upper and lower potentials, the vector ϕ(x i ) should increase the smallest eigenvalues of A i-1 , while staying far from the eigenvectors with largest eigenvalues. Due to the second condition, the lower barrier can only increase half as much as it could without it, resulting in the constraint m 2n. In addition, we still achieve an a posteriori upper bound, see (21), which is not particularly sharp, but sufficient for our purposes. Lastly, a full inversion of matrix A m , with complexity O(n 3 ), remains reasonable in view of the moderate values of n imposed by the sampling costs, see Section 5.

Random sampling by effective resistance

We first consider an approach combining (5) and ( 6) to obtain the first estimate (1) in Theorem 1.1. To observe the appropriate decay as r → ∞, we impose that the vector

ϕ, g m = m i=1 s i ϕ(x i )g(x i )
is an unbiased estimator of ϕ, g L 2 = 0, by taking weights s i inversely proportional to the sampling density of x i . If the points x i are independent, the euclidian norm | ϕ, g m | 2 is therefore bounded in expectation by the sum of the variances of each term, of the form

E(s 2 i |ϕ(x i )| 2 |g(x i )| 2 ) ∝ X s i |ϕ(x i )| 2 |g(x i )| 2 dµ(x i ) s i |ϕ(x i )| 2 L ∞ g 2 L 2 .
As a result, to control m from below and these variances from above, the appropriate sampling density for point x i is a combination of the effective resistance mentioned earlier and the socalled Christoffel function |ϕ(x)| 2 = n j=1 |ϕ j (x)| 2 . Algorithm 3.1 is inspired by the deterministic procedure from [START_REF] Batson | Twice-ramanujan sparsifiers[END_REF] and its randomization presented in [START_REF] Tat | Constructing linear-sized spectral sparsification in almostlinear time[END_REF], which has already been applied to least-squares recovery in [START_REF] Chen | Active regression via linear-sample sparsification[END_REF]. It takes as inputs the parameters

ε ∈ (0, 1), η = ε 1 -ε ∈ (0, ∞), and γ 0,
which respectively influence the size of the barrier increments δ i = ii-1 , the weights s i , and the balance between effective resistance and Christoffel function in the sampling density.

Algorithm 3.1: Random sampling by effective resistance Input: probability measure µ, parameter ε ∈ (0, 1), and γ 0.

A 0 = 0 ∈ R n×n , 0 = -n for i = 1, . . . , m do Let Y i = (A i-1 -i-1 I) -1 Update i = i-1 + δ i , where δ i := ε/Ψ i with Ψ i = Tr(Y i ) + γ Let Z i = (A i-1 -i I) -1 Let ρ i : x → ϕ(x) * (Z i + γI/n)ϕ(x) Draw x i from probability measure ρ i (x) Ξ i dµ(x), where Ξ i = Tr(Z i ) + γ Let s i = η/ρ i (x i ) Update A i = A i-1 + s i ϕ(x i )ϕ(x i ) *
end for Output: sample {x 1 , . . . , x m }, weights {s 1 , . . . , s m }, and matrix A m .

For the sake of the analysis, we also define

Y m+1 = (A m -m I) -1 , Ψ m+1 = Tr(Y m+1 ) + γ, δ m+1 = ε Ψ m+1 and m+1 = m + δ m+1 ,
as would have been done at iteration i = m + 1.

Lemma 3.1. The algorithm is almost surely well defined and A i i+1 I for all 0 i m, where stands for the Loewner order between positive semi-definite matrices.

Proof. We proceed by induction on i. At initialization A 0 -nI = 0 I. For i ∈ {1, . . . , m} fixed, assume that A i-1 i-1 I. This implies that Y i is well-defined, positive definite, and

δ i I = εI Tr(Y i ) + γ εI λ max (Y i ) = λ min (A i-1 -i-1 I) εI ≺ A i-1 -i-1 I. ( 8
)
Hence A i-1 ( i-1 + δ i )I = i I, so Z i is well-defined, positive definite, and by identity (4),

X ρ i (x) dµ(x) = Tr(Z i + γI/n) = Ξ i ,
proving that ρ i dµ/Ξ i is indeed a probability density. Finally s i > 0 almost surely, so A i A i-1 i I. This completes the induction, and a last iteration shows that A m m+1 I, concluding the proof.

Contrarily to most variations on the algorithm from [START_REF] Batson | Twice-ramanujan sparsifiers[END_REF], no upper potential is used to bound the eigenvalues of A i from above. Instead, parameter γ provides a lower bound on ρ i in terms of the Christoffel function ϕ(x) * ϕ(x), which turns into an upper bound on the norm of rank-one terms s i ϕ(x i )ϕ(x i ) * , and therefore on the norm of their sum A m .

To bound the eigenvalues of A i from below, we use the lower barrier i+1 , which should increase at each step, at a speed δ i controlled by the lower potential Tr(Y i ). As the densities ρ i are positive, the sample (x 1 , . . . , x m ) could be in any part of X m , so there cannot be any positive deterministic lower bound on λ min (A m ). However, the following lemma proves a monotonicity property in expectation, similar to Lemmas 4.4 in [START_REF] Tat | Constructing linear-sized spectral sparsification in almostlinear time[END_REF] and [START_REF] Tat | An SDP-based algorithm for linear-sized spectral sparsification[END_REF], themselves inspired by Lemmas 3.3 and 3.4 in [START_REF] Batson | Twice-ramanujan sparsifiers[END_REF].

Lemma 3.2. The sequence (E(Tr(Y i )) 1 i m+1 is non-increasing.

Proof. Let 1 i m, and denote

u i = √ s i ϕ(x i ). Applying the Sherman-Morrison formula yields Y i+1 = (Z -1 i + u i u * i ) -1 = Z i - Z i u i u * i Z i 1 + u * i Z i u i . ( 9
)
By design of the sampling density and weights,

u * i Z i u i = s i ϕ(x i ) * Z i ϕ(x i ) s i ρ i (x i ) = η, (10)
so by definition of η = ε/(1 -ε), which implies 1 + η = η/ε, we obtain

Y i+1 -Z i - Z i u i u * i Z i 1 + η = - ε η s i Z i ϕ(x i )ϕ(x i ) * Z i . (11)
On the other hand, Y i and Z i commute because their inverses do, and we can write

Z i -Y i = (Y -1 i -Z -1 i )Y i Z i = δ i Y i Z i .
Combining this identity to the previous estimate and taking traces gives

Tr(Y i+1 ) -Tr(Y i ) = δ i Tr(Y i Z i ) - ε η s i ϕ(x i ) * Z 2 i ϕ(x i ).
Observe that Y i and Z i only depend on the samples x 1 , . . . , x i-1 . As a result, if these samples are fixed, taking the expectation

E x i = E( • |x 1 , . . . , x i-1
) with respect to point x i , we arrive at

E x i Tr(Y i+1 ) -Tr(Y i ) = δ i Tr(Y i Z i ) - ε η X η ρ i (x i ) ϕ(x i ) * Z 2 i ϕ(x i ) ρ i (x i ) Ξ i dµ(x i ) = ε Tr(Y i Z i ) Tr(Y i ) + γ - Tr(Z 2 i ) Tr(Z i ) + γ , (12)
where we have used identity (4). Since

Z i = Y i + δ i Y i Z i and Z 2 i = Y i Z i + δ i Y i Z 2 i , (13)
we have by Cauchy-Schwarz inequality

Tr(Y i Z i ) Tr(Z i ) -Tr(Y i ) Tr(Z 2 i ) = δ i Tr(Y i Z i ) 2 -Tr(Y i ) Tr(Y i Z 2 i ) 0. (14) Moreover, Tr(Y i Z i ) -Tr(Z 2 i ) = -δ i Tr(Y i Z 2 i ) = -δ i Tr(Z i Y i Z i ) < 0.
Hence the right-hand side in (12) is negative. Taking an expectation over the previous samples x 1 , . . . , x i-1 in turn implies that E Tr(Y i+1 ) is smaller than E Tr(Y i ) , which concludes the proof.

One main difference with [START_REF] Tat | Constructing linear-sized spectral sparsification in almostlinear time[END_REF] and [START_REF] Tat | An SDP-based algorithm for linear-sized spectral sparsification[END_REF] is the choice of different matrices Y i and Z i for updating the lower barrier i and selecting a new point x i . This comes at the expense of the refined analysis (13), ( 14) to show that (12) is non-positive. Moreover, in view of (8), this changes the sampling density at most by a factor 1/(1 -ε).

However, if we had used Y i instead of Z i in the sampling density ρ i , in other words if we had defined ρ i (x) = ϕ(x) * (Y i + γI/n)ϕ(x) and Ξ i = Tr(Y i ) + γ, we would have obtained (10), where the first inequality comes from the fact that Y -1 i -Z -1 i εY -1 i in view of (8). In order to recover the same right-hand side as before in inequality (11), one would then need to assume

u * i Z i u i u * i Y i 1 -ε u i s i ρ i (x i ) 1 -ε = η 1 -ε instead of
1 1 + η 1-ε ε η ⇐⇒ 1 + η 1 -ε η ε ⇐⇒ η 1 1 ε -1 1-ε = ε(1 -ε) 1 -2ε .
So this simpler approach requires ε < 1 2 , which is not possible unless r > 2, and even in that case our approach gives slightly better estimates.

The control on the lower potentials in expectation from Lemma 3.2 yields a lower bound on the eigenvalues of A m in probability.

Proposition 3.3. Let m n, p ∈ (0, 1), and define r = (m + 1)/n. The random matrix A m generated by Algorithm 3.1 satisfies

A m αI, α = n εr 1/p + γ -1 ,
with probability at least 1 -p.

Proof. In view of the monotonicity property established in Lemma 3.2, and thanks to the initialization 0 = -n, we have

E Tr(Y i ) Tr(Y 1 ) = 1, 1 i m + 1.
By Markov's inequality

P m+1 i=1 Tr(Y i ) > m + 1 p p E( m+1 i=1 Tr(Y i )) m + 1 p.
As Ψ i = Tr(Y i ) + γ, we deduce that with probability at least 1 -p, it holds

m+1 i=1 Ψ i (m + 1) 1 p + γ .
Together with the inequality

m+1 i=1 Ψ i m+1 i=1 Ψ -1 i (m + 1) 2 , this implies m+1 i=1 Ψ i -1 m + 1 1/p + γ = rn 1/p + γ ,
and therefore

m+1 = 0 + m+1 i=1 ε Ψ i -n + rnε 1/p + γ = α,
with probability at least 1 -p. Since A m m+1 I, see Lemma 3.1, the proof is complete.

Remark 3.4. The lower bound A m αI can be rewritten as a frame inequality

c * A m c = m i=1 s i | c, ϕ(x i ) | 2 α|c| 2 , c ∈ C n ,
or as a Marcinkiewicz-Zygmund inequality: for all v ∈ V n , v 2 m α v 2 L 2 . All three versions have been extensively used in the literature for emphasizing the relations with subsampling of frames and discretization of continuous norms, see [FS19, NOU13, NOU16, LT22]. Note that our approach is not performing a subsampling in general, although it can be interpreted as such in the context of Remark 5.1.

We are now ready for the proof of the first statement of Theorem 1.1. Let E be the event of probability at least 1 -p where A m αI, for p ∈ (0, 1) and α > 0 given by Proposition 3.3, and define

f = P m n f |E (15)
as the weighted least-squares projection, conditioned to event E. In practice, this amounts to relaunching Algorithm 3.1 until E occurs, and computing the weighted least-squares projection of f with the last sample.

Proof of Theorem 1.1, equation (1). As V n ⊂ L 2 , we can assume that f is also in L 2 , otherwise the right-hand side would be infinite. Similar to (5), by Pythagoras theorem,

E( f -f 2 L 2 ) = f -P n f 2 L 2 + E( f -P n f 2 L 2 ) = g 2 L 2 + E( P m n f -P n f 2 L 2 |E) = g 2 L 2 + E( P m n g 2 L 2 |E).
where g = f -P n f . Combining this with equation ( 6), together with the definition of E, yields

E( f -f 2 L 2 ) g 2 L 2 + α -2 E(| ϕ, g m | 2 |E).
We now use Proposition 3.3:

E(| ϕ, g m | 2 |E) = E(| ϕ, g m | 2 χ E ) P(E) 1 1 -p E(| ϕ, g m | 2 ).
Develop the discrete inner products •, • m as

E(| ϕ, g m | 2 ) = E n k=1 | ϕ k , g m | 2 = n k=1 E   m i=1 s i ϕ k (x i )g(x i ) m j=1 s j ϕ k (x j )g(x j )   .
When i > j, we have that

E x i s i ϕ k (x i )g(x i ) = η Ξ i X ϕ k (x)g(x)dµ(x) = 0,
since g is orthogonal to V n , so the corresponding term is zero. The same holds true for i < j, by exchanging i and j. We are only left with the diagonal terms

E(| ϕ, g m | 2 ) = E m i=1 s 2 i |ϕ(x i )| 2 |g(x i )| 2 nη γ E m i=1 s i |g(x i )| 2 = nη γ mη Ξ i g 2 L 2 ,
where we have first used the uniform bound s i |ϕ(x i )| 2 nη/γ, then integrated over x i using s i ρ(x i ) = η, both identities following from the definition of s i in Algorithm 3.1. Since Ξ i γ, m rn and η = ε/(1 -ε), we deduce that

E( f -f 2 L 2 ) 1 + 1 α 2 (1 -p) rn 2 ε 2 γ 2 (1 -ε) 2 g 2 L 2 = 1 + r 1 -p 1 β(ε, γ) 2 g 2 L 2 , (16)
where

β(ε, γ) = αγ(1 -ε) nε = γ rε 1/p + γ -1 1 ε -1
given the definition of α in Proposition 3.3. For p fixed, maximizing β(ε, γ) over ε ∈ (0, 1) and γ > 0 yields conditions

1 p + γ = rε 2 and 1 p + γ 2 = rε p ,
which are equivalent to rε 3 = 1/p and γ = rε 2 (1 -ε), and the maximum value of β(ε, γ) is r(1 -ε) 3 . Letting p = r -1/4 ∈ (0, 1), and thus ε = r -1/4 ∈ (0, 1) and γ = r 1/2 -r 1/4 > 0, we obtain the upper bound in (1). Note that the first optimality condition above automatically implies α = n(ε -1 -1) > 0, making the proof consistent.

Notice that, with this last choice p = ε = r -1/4 , γ = (1 -ε)/ε 2 and α = n(1 -ε)/ε. In particular, as r tends to infinity, the probability of failure 1 -p in Proposition 3.3 goes to zero, and γ goes to infinity, meaning that the sampling density gets closer to the Christoffel function, as in [START_REF] Cohen | Optimal weighted least-squares methods[END_REF].

Remark 3.5. In the example where V n consists of piecewise constant functions on a fixed partition of X into pieces of equal measure, consider f such that f -P n f is a realisation of a white noise of variance 1. It is easily seen that the optimal approximation error is the Monte-Carlo rate 1+1/r+o(1/r), with f made of averages over r random samples in each piece. In that sense, Algorithm 3.1 achieves the optimal decay rate of the error, similar to [CM17, HNP22], but with a remainder independent of n.

Refined randomized sampling algorithm

We now seek to optimize the approximation strategy in the regime r ≈ 1. Although Algorithm 3.1 works as soon as r 1, several improvements can be performed for small values of r. First, using (7) instead of (6) reduces the factor η 2 /α 2 to η/α in bound (16). Secondly, bounding the weights s i by a multiple of 1/γ becomes a crude estimate when r → 1, since γ → 0 in that case. Lastly, the use of Markov's inequality in Proposition 3.3 comes at the expense of a factor 1/(1 -p), which grows as r gets close to 1; to avoid the last issue, we would like the lower barrier to grow in a steady, deterministic fashion. This prevents the sample points x i from being drawn anywhere in the domain X.

For all these reasons, we consider Algorithm 4.1, with inputs δ ∈ (0, 1), the size of lower barrier increments, and κ ∈ [0, 1], a parameter balancing the L 2 and L ∞ settings. Note that we will always fix m beforehand, and set δ = 1/ √ r, where r = m/(n -1). Again, we define Y m+1 = (A mm I) -1 , as we would have done at iteration i = m + 1.

Algorithm 4.1: Random sampling with fixed increments of the barrier

Input: probability measure µ, parameter δ ∈ (0, 1), and κ ∈ [0, 1].

A 0 = 0 ∈ R n×n , 0 = -n for i = 1, . . . , m do Let Y i = (A i-1 -i-1 I) -1 Update i = i-1 + δ Let Z i = (A i-1 -i I) -1 Let W i = Tr(Z i ) -Tr(Y i ) -1 Z 2 i -Z i and w i : x → ϕ(x) * W i ϕ(x) Define R i (x) = w i (x)1 w i (x) κ 1-δ δ Draw x i from probability density R i (x) Γ i dµ(x), where Γ i = X R i (x)dµ(x) Let s i = 1/R i (x i ) = 1/w i (x i ) Update A i = A i-1 + s i ϕ(x i )ϕ(x i ) *
end for Output: Sample {x 1 , . . . , x m }, weights {s 1 , . . . , s m }, and matrix A m .

Lemma 4.1. The algorithm is well-defined and Tr(Y i ) = 1 for 1 i m + 1.

Proof. We use an induction on i. First, note that Y 1 = I/n hence Tr(Y 1 ) = 1. Let i ∈ {1, . . . , m}, and assume that Y i is positive definite with Tr(Y i ) = 1. We have

λ min (Y -1 i ) = λ max (Y i ) -1 (Tr(Y i )) -1 = 1 > δ, Therefore Z i = (Y -1 i -δI) -1
is well defined and positive definite, and W i is also well defined since Z i Y i . We then need to show that R i (x) is not zero a.e., in other words that w i (x) κ(1 -δ)/δ on a set of positive mesure. To this end, we will rely on an averaging argument by showing that X w i dµ, which is equal to Tr(W i ) according to (4), is larger than (1 -δ)/δ. Recalling (13) and ( 14), we observe as in Claim 3.6 of [START_REF] Batson | Twice-ramanujan sparsifiers[END_REF] Tr

(W i ) = Tr(Z 2 i ) Tr(Z i ) -Tr(Y i ) -Tr(Z i ) = 1 δ Tr(Z 2 i ) Tr(Y i Z i ) - Tr(Z i ) Tr(Y i ) + 1 -δ δ Tr(Z i ) > 1 -δ δ , (17)
where we used the induction hypothesis Tr(Y i ) = 1, and its immediate consequence Tr(Z i ) > 1.

Finally, using the Sherman-Morrison formula as in (9), with

u i = √ s i ϕ(x i ), Tr(Y i+1 ) = Tr(Z i ) - u * i Z 2 i u i 1 + u * i Z i u i = Tr(Z i ) - ϕ(x i ) * Z 2 i ϕ(x i ) w i (x i ) + ϕ(x i ) * Z i ϕ(x i ) = Tr(Y i ) = 1,
which concludes the induction.

Remark 4.2. The symmetric matrix W i is not necessarily positive semi-definite, but we have a framing on its maximal eigenvalue

0 < Tr(W i ) n λ max (W i ) λ max (Z 2 i ) Tr(Z i ) -Tr(Y i ) Tr(Z 2 i ) δ Tr(Y i Z i ) 1 δ(1 -δ) , since Z i Y i (1 -δ)Z i
. This property will be useful in the discussion of Section 5.

As i = iδ -n, we easily get a lower bound on the eigenvalues of A m . A m αI a.s., with α = (n -1)(δr -1).

Proof. After the last iteration of the algorithm, we have

λ min (A m ) = m + (λ max (Y m+1 )) -1 mδ -n + (Tr(Y m+1 )) -1 = (δr -1)(n -1).
We now have all the tools for the proof of the main theorem.

Proof of Theorem 1.1, equation (2). Introducing the set

E i = x ∈ X, w i (x) κ 1 -δ δ , observe that Γ i = E i R i dµ = X w i (x)dµ(x) - E c i w i (x)dµ(x) Tr(W i ) - X κ 1 -δ δ dµ (1 -κ) 1 -δ δ . (18)
Taking the expectation E x i over point x i , with x 1 , . . . , x i-1 fixed, we compute

E x i (s i |g(x i )| 2 ) = E i 1 R i (x) |g(x)| 2 R i (x) Γ i dµ(x) g 2 L 2 Γ i 1 1 -κ δ 1 -δ g 2 L 2 ,
where g = f -P n f ∈ L 2 (X, µ). The same holds for all 1 i m, and consequently

E g 2 m m 1 -κ δ 1 -δ g 2 L 2 ,
which together with (7) and Proposition 4.3 gives

E( P m n g 2 ) α -1 E P m n g 2 m α -1 E g 2 m r m(δr -1) m 1 -κ δ 1 -δ g 2 L 2 .
Taking δ = 1/ √ r and adding (5), we arrive at (2).

Remark 4.4. For κ = 0, the density of point x i is proportional to w i (x) + = max(w i (x), 0). A sampling density proportional to the positive part of a function also occurs in [START_REF] Tat | An SDP-based algorithm for linear-sized spectral sparsification[END_REF], equation (6). There, it is used to select points which increase the lower potential without increasing too much the upper potential.

Proof of Theorem 1.1, equation (3). As V n ⊂ L ∞ , we can assume that f is also in L ∞ , otherwise the right-hand side would be infinite. We introduce the notation

P ∞ n f ∈ arg min v∈Vn f -v L ∞
for a Chebyshev projection of f onto V n with respect to the Banach norm • L ∞ , and take

g = f -P ∞ n f the associated residual. As 1/s i = w i (x i ) κ(1 -δ)/δ, it holds g 2 m = m i=1 s i |g(x i )| 2 m i=1 s i g 2 L ∞ mδ κ(1 -δ) g 2 L ∞ a.s.
By a triangular inequality

f -P m n f L 2 f -P ∞ n f L 2 + P m n f -P ∞ n f L 2 = g L 2 + P m n g L 2 .
We conclude by bounding the operator norm P m n L ∞ →L 2 through

P m n g 2 L 2 α -1 P m n g 2 m α -1 g 2 m r δr -1 δ κ(1 -δ) g 2 L ∞ . ( 19
)
As a result, Algorithm 4.1 with δ = 1/ √ r implies inequality (3).

Remark 4.5. In the randomized setting, µ could be any positive measure on X. On the contrary, in the uniform setting, we need µ to be a probability measure to bound L 2 norms by L ∞ norms. Similar results would hold if µ was a positive measure of finite mass, with all L ∞ norms multiplied by a factor µ(X), see for instance [START_REF] Limonova | On sampling discretization in L 2[END_REF][START_REF] Temlyakov | On optimal recovery in L 2[END_REF][START_REF] Krieg | Sampling recovery in the uniform norm[END_REF].

We turn to the proof of the last result of the introduction, which is stated in strong uniform norm to avoid technicalities.

Proof of Corollary 1.4. Denote L ∞ the set of bounded functions, and f L ∞ := sup x∈X |f (x)| the associated norm. As V n ⊂ L ∞ , we can assume that f is also in L ∞ , otherwise the right-hand side would be infinite. Taking g = f -P ∞ n f the L ∞ -residual, we write as above

f -P m n f L ∞ f -P ∞ n f L ∞ + P m n f -P ∞ n f L ∞ = g L ∞ + P m n g L ∞ .
For any probability measure µ on X, consider (ϕ 1 , . . . , ϕ n ) an orthonormal basis of V n in L 2 = L 2 (X, µ). Then the supremum of the Christoffel function is, by a Cauchy-Schwarz inequality,

K n := sup x∈X |ϕ(x)| 2 = sup x∈X max c∈R n \{0} | n j=1 c j ϕ j (x)| 2 |c| 2 = max v∈Vn\{0} v 2 L ∞ v 2 L 2
.

By simply considering (19) in the case m = n, κ = 1, and δ = 1/ √ r 1 -1/2n, we derive the inequality P n n g L 2 (1 -δ) -1 g L ∞ 2n g L ∞ . As P n n g ∈ V n , we obtain

P n n g L ∞ K n P n n g L 2 2n K n g L ∞ .
We now rely on the main corollary of the paper [START_REF] Kiefer | The equivalence of two extremum problems[END_REF]. Since X is compact, there exists a measure µ * for which K n = n. Using such a measure in our framing, we obtain by application of the above

P n n g L ∞ 2n
√ n g L ∞ , which concludes the proof.

Let us briefly explain the result in [START_REF] Kiefer | The equivalence of two extremum problems[END_REF] to our setting. Given φ = (φ 1 , . . . , φ n ) a basis of V n and µ a probability measure over X, we consider the n × n symmetric positive semi-definite matrix M (µ) := X φφ * dµ. Assuming this matrix is non-singular, ϕ := M (µ) -1/2 φ is also basis of V n and is orthonormal with respect to µ. In particular ϕ(x) * ϕ(x) = φ * M (µ) -1 φ is the associated Christoffel function. The optimal probability measures from [START_REF] Kiefer | The equivalence of two extremum problems[END_REF] are characterized by µ * ∈ arg min µ∈Prob(X)

sup x∈X φ(x) * M (µ) -1 φ(x) = arg max µ∈Prob(X)
det (M (µ)) .

These two equivalent extremum problems are difficult to solve for general compacts X and spaces V n , thus making the approach from Corollary 1.4 unpractical. Nevertheless, it gives a theoretical bound on the approximation error achievable by this method. We refer to [START_REF] Bos | Some remarks on the Fejér problem for Lagrange interpolation in several variables[END_REF] for an alternative proof of the result in [START_REF] Kiefer | The equivalence of two extremum problems[END_REF], with applications to interpolation, in particular to the so-called Fejer problem. We also note that [START_REF] Kiefer | The equivalence of two extremum problems[END_REF] answers Remark 5.5 in [START_REF] Pozharska | A note on sampling recovery of multivariate functions in the uniform norm[END_REF].

By inspection of the proof of Corollary 1.4, we observe that any measure µ used in Algorithm 4.1 yields a bound of the form

f -P n n f L ∞ (1 + 2n K n ) f -P ∞ n f L ∞ . ( 20
)
Corollary 1.4 can also be stated for m n, as done in Remark 1.5. For instance, Algorithm 4.1 for m = 2n, hence with parameters δ = 1/ √ r 1/ √ 2, κ = 1, and any measure µ, yields

f -P m n f L ∞ (1 + 4 K n ) f -P ∞ n f L ∞ .
The above remarks has to be combined with upper bounds on K n , for which we give a few examples in polynomial approximation.

We first consider the space V n = P n-1 of algebraic polynomials of degree less than n, restricted to X = [-1, 1]. By considering the uniform measure dµ(x) = dx/2, we have

K n = n-1 j=0 (2j + 1) = n 2 , while with the arcsine measure dµ(x) = dx π √ 1-x 2 , we obtain K n = 1 + n-1 j=1 2 = 2n -1,
where we used uniform bounds on the Legendre and Chebyshev polynomials, respectively.

Similar results are also available for spaces of multivariate polynomials indexed in lower set, see [CCM + 15, Lemma 3.1]. We recall that a lower set, or downward closed set, is a set of multi-indices Λ ⊂ N d 0 such that ν ∈ Λ implies ν ∈ Λ for all ν ν, the last inequality being understood component-wise. For X = [-1, 1] d , we consider polynomial approximation spaces of the form

V n = Span{x ν 1 1 . . . x ν d d , ν
∈ Λ}, with Λ a lower set of cardinality n. Among the associated spaces V n , one can find polynomials of bounded degree in each variable, of bounded total degree, as well as hyperbolic cross spaces, which are of important use in high-dimensional approximation [START_REF] Dung | Hyperbolic cross approximation[END_REF]. The following estimates hold:

• for µ the uniform measure, K n = ν∈Λ d j=1 (1 + 2ν j ) is bounded by n 2 ,
• for µ the tensorized arcsine measure, K n = ν∈Λ 2 |ν| 0 is bounded by n 2 d , respectively. Such estimates might be much smaller, for instance in the second setting when d is small. Combining these bounds with (20) implies the results given at the end of the introduction.

In the particular case of polynomials of fixed total degree, invariance by affine transformations allows to treat more general cases.

Theorem 4.6 ([CD21], Theorems 5.4 and 5.6). For X ⊂ R d a compact domain, dµ(x) = dx/|X| the uniform measure, and V n a set of polynomials of fixed total degree, it holds K n C X n 2 if X has a Lipschitz boundary, and

K n C X n 1+1/d if X has a smooth boundary.
It is worth noting that the bound K n C X n 2 also holds for arbitrary lower sets when the domain is a union of rectangles of fixed volume, see Theorem 6.5 of [START_REF] Adcock | Approximating smooth, multivariate functions on irregular domains[END_REF]. All these results show that, in many instances, we have K n = O(n 2 ) with the uniform measure. Therefore Algorithm 4.1 provides a feasible way to sample points for which the Lebesgue constant is

L n = O(n 2 ).

Numerical aspects

Our two algorithms have a polynomial complexity in n, since one has to solve linear systems of size n at each iteration. In fact, Tr(Y i+1 ) can be computed without the knowledge of Y i+1 , by taking a trace in the Shermann-Morrison formula (9), i.e.

Tr(Y i+1 ) = Tr(Z i ) - ϕ(x i ) * Z 2 i ϕ(x i ) s i + ϕ(x i ) * Z i ϕ(x i ) ,
so we only need a matrix inversion for Z i . Moreover, this can be done efficiently by updating the eigenvalues and eigenvectors of A i from one iteration to the next, see [START_REF] Tat | Constructing linear-sized spectral sparsification in almostlinear time[END_REF] for details.

In [START_REF] Tat | Constructing linear-sized spectral sparsification in almostlinear time[END_REF] and [START_REF] Tat | An SDP-based algorithm for linear-sized spectral sparsification[END_REF], polynomial and exponential modifications of the lower potential are also proposed, allowing to draw the samples by batches, and therefore reducing the complexity to a linear expression in n up to logarithmic factors, in the particular case of subsampling a graph laplacian. The proofs are more involved, produce larger constants and require a large value of r, so we did not adapt them to our setting.

The essential remaining difficulty is to draw each point x i according to its prescribed density ρ i (or R i ). We start with the following observation.

Remark 5.1. An important application is the case where X = {x 1 , . . . , x M } is a finite set, and µ is the uniform measure on X. The orthonormality of the basis ϕ = (ϕ 1 , . . . , ϕ n ) is encoded by the following decomposition of the identity into a sum of rank-one matrices:

M i=1 v i v * i = I, v i := 1 √ M (ϕ 1 (x i ), . . . , ϕ n (x i )) .
In this situation, it suffices to evaluate ρ i at all points of X, and to select x i = x with probability ρ i (x)/ x ∈X ρ i (x ). Our algorithms 3.1 and 4.1 can also work as a sample reduction technique, see for example [START_REF] Cohen | Optimal pointwise sampling for L 2 approximation[END_REF], where the evaluation points are extracted from a large initial sample x 1 , . . . , x M . Indeed, it suffices to apply our algorithms to the set {x 1 , . . . , x M }, equipped with the empirical measure µ = 1 M M i=1 δ x i .

When drawing points directly from an infinite set X, one may use an acceptance-rejection strategy. Notice that ρ i in Algorithm 3.1 and w i in Algorithm 4.1 are of the form ϕ(x) * M ϕ(x), with either M = Z i + γI/n or M = W i . Moreover, some routines have been developed [CM17, AC20, ACD23, Mig21, Mig19, CD21] for drawing points from the Christoffel measure 1 n |ϕ(x)| 2 dµ(x) in relevant multivariate settings.

In the case of Algorithm 3.1, as M = Z i + γI/n is positive semi-definite, it suffices to draw candidate points x from the Christoffel measure, and to accept them with probability

ϕ(x) * M ϕ(x) λ max (M )|ϕ(x)| 2 1.
The probability of accepting a point is

x∈X ϕ(x) * M ϕ(x) λ max (M )|ϕ(x)| 2 |ϕ(x)| 2 n dµ(x) = Tr(M ) nλ max (M ) 1 n ,
therefore we need in average nλ max (M )/ Tr(M ) n draws from the Christoffel measure to find x i .

In the case of Algorithm 4.1, M = W i may have negative eigenvalues, so the bound λ max (M ) Tr(M ) may not hold. In order to draw a sample from the probability density R i (x)/Γ i , we can still rely on acceptance-rejection from the Christoffel measure, and accept points with probability

R i (x) λ max (W i )|ϕ(x)| 2 ∈ [0, 1].
The acceptance probability is then

x∈X R i (x) λ max (W i )|ϕ(x)| 2 |ϕ(x)| 2 n dµ(x) = Γ i nλ max (W i ) (1 -κ)(1 -δ) 2 n
in view of Remark 4.2, and we need to at most n(1 -κ) -1 (1 -δ) -2 candidate points in average. When κ = 1, inspection of ( 18) and ( 17) respectively shows that

Γ i Tr(W i ) - 1 -δ δ and Tr(W i ) 1 -δ δ Tr(Z i ), so that δ 1 -δ Γ i Tr(Z i ) -1 = δ Tr(Y i Z i ) δ Tr(Y 2 i ) δ n Tr(Y i ) 2 = δ n ,
hence we can replace the factor 1 -κ in the acceptance probability by a factor δ/n.

Remark 5.2. In the first iterations 1 i < n, the maximal eigenvalue of M has multiplicity at least n -i, so the bound on the acceptance probability is n -i times higher. This implies for instance that there are in average less than 2 (respectively, less than 2(1 -κ) -1 (1 -δ) -2 ) rejections for i = 1, . . . , n/2 in Algorithm 3.1 (respectively, Algorithm 4.1), which we observe in practice.

Even with such methods, the computational complexity remains dominated by the sampling step: Algorithm 3.1 runs in O(mn 3 t) time, where the factor m comes from the number of iterations, a factor n 2 from the matrix-vector multiplication in the evaluation of ρ i , the last factor n from the average number of rejections, and t is the time needed to generate a point with the Christoffel density. The linear algebra takes O(n 3 ) time per iteration, hence O(mn 3 ) time in total. Similarly, Algorithm 4.1 runs in O mn 3 t(1 -κ) -1 (1 -δ) -2 time, the only difference coming from the bound on the acceptance ratio. Therefore, applying our sampling strategies may prove quite challenging for large values of n.

Remark 5.3. In the deterministic setting, it is not necessary to draw each point x i with a density proportional to ρ i (x): the only requirement is that w i (x i ) κ 1-δ δ for all 1 i m. The idea of using the algorithm from [START_REF] Batson | Twice-ramanujan sparsifiers[END_REF], still in a deterministic context, but with an infinite set X, seems to originate in [DPS + 21], which investigates Marcinkiewicz-type discretization theorems. However, it seems that searching for x i by rejection sampling works better in practice than sorting X and looking for the first point that achieves the above condition [START_REF] Bartel | Constructive subsampling of finite frames with applications in optimal function recovery[END_REF].

We next comment on the relatively limited importance of the weights. In Algorithm 3.1, the sampling density decomposes as

ρ i (x) Ξ i = (1 -p) ϕ(x) * Z i ϕ(x) Tr(Z i ) + p |ϕ(x)| 2 n ,
with p = γ/Ξ i . In others words, we are sampling from a mixture of effective resistance and Christoffel density. Here we only know that s i |ϕ(x i )| 2 nη/γ, which by the way yields

λ max (A m ) Tr(A m ) = m i=1 s i |ϕ(x i )| 2 mnη γ . ( 21
)
But there is no a priori uniform bound on the weights. If we desire such a bound, it suffices to add a third, constant term in the sampling density:

ρ i (x) Ξ i = Tr(Z i ) Ξ i ϕ(x) * Z i ϕ(x) Tr(Z i ) + γ Ξ i |ϕ(x)| 2 n + γ ∞ Ξ i
for a new parameter γ ∞ > 0, where we redefine Ξ i as Tr(Z i ) + γ + γ ∞ . To choose a point according to this density, one can simply draw it from dµ(x) with probability γ ∞ /Ξ i , from the Christoffel measure with probability γ/Ξ i , and from the effective resistance measure otherwise. This choice immediately yields s i = η/ρ i (x i ) η/γ ∞ , however taking γ ∞ > 0 deteriorates a bit the other estimates, which is why we did not include it earlier.

Having a uniform bound on the weights s i is particularly important in case of noisy samples, where one can only observe f (x i ) + ε i for some random ε i . Then the estimates of Theorem 1.1 hold with an additional noise term of the form

L + ( √ s i ε i ) 1 i m 2 2 λ min (A m ) -1 η γ ∞ (ε i ) 1 i m 2 2 ,
due the definition of the weighted least-squares estimator in Section 2. We refer to Theorems 5.3 and 5.19 of [START_REF] Adcock | Sparse polynomial approximation of high-dimensional functions[END_REF] for details.

Concerning Algorithm 4.1, in the deterministic setting, one can replace each weight s i by its upper bound δ κ(1-δ) , resulting in an unweighted discrete norm

g 2 m = 1 m m i=1 |g(x i )| 2 ,
without changing our estimates. We refer to [START_REF] Bartel | Constructive subsampling of finite frames with applications in optimal function recovery[END_REF] for earlier results on subsampling of frames with unweighted discrete norms. In any case, in numerical experiments, the weights seem to have lesser importance than the position of the points.

Numerical experiments

We conclude with a few illustrations in the space of multivariate polynomials, focusing on a function whose best approximation is explicitly known. Our code is available at the following address: github.com/Belloliva/minimal-oversampling-for-multivariate-polynomial-least-squares.

In the univariate setting, the Legendre polynomials (L k ) k 0 are orthonormal in L 2 ([-1, 1], dx/2), and satisfy

L k L ∞ = L k (1) = √ 2k + 1.
We recall that their generating function is given by

g y (x) = ∞ k=0 L k (x) L k (1) y k = 1 1 -2xy + y 2 , x ∈ [-1, 1], y ∈ (-1, 1). Now, on the domain X = [-1, 1] d with d ∈ N, one can consider tensorized Legendre polynomials L k (x) := d j=1 L k j (x j ), x = (x 1 , . . . , x d ) ∈ X, k = (k 1 , . . . , k d ) ∈ N d 0 ,
which are orthonormal for the uniform measure dµ(x) = dx/2 d over X. Fixing y = (y 1 , . . . , y d ) ∈ (0, 1) d , we would like to approximate the multivariate function

x → g y (x) = d j=1 g y j (x j ) = k∈N d 0 c k L k (x), c k = d j=1 y k j j 2k j + 1 .
This is a typical representer of holomorphic functions on X with anisotropic dependance on x parameterized by y. Its L 2 (X, µ) norm is analytically given by

g y 2 L 2 = d j=1 1 -1 |g y j (x j )| 2 dx j 2 = d j=1 1 2y j ln 1 + y j 1 -y j .
Moreover, the best approximation space of dimension n for g y w.r.t.

• L 2 is V n = Span{L k , k ∈ Λ n },
where Λ n is the set of n multi-indices k with largest values of c k . Note that Λ n is lower and is easily computed, since the coefficients c k are decreasing for the component-wise order, and that the approximation error is

min v∈Vn g y -v 2 L 2 = g y - k∈Λn c k L k 2 L 2 = g y 2 L 2 - k∈Λn c 2 k . (22)
Remark 6.1. When looking at approximations in the uniform norm, although V n and k∈Λn c k L k may not be optimal, one can still compute an upper bound

min v∈Vn g y -v L ∞ g y - k∈Λn c k L k L ∞ = d j=1 1 1 -y j - k∈Λn d j=1 y k j j ,
since the maximum is attained at (1, . . . , 1) ∈ X.

In Figure 1, we take parameters y = (0.9, 0.8, 0.7, 0.6) in spatial dimension d = 4, use Legendre expansions of size n = 128, and draw m = 2n = 256 points x 1 , . . . , x m ∈ X according to one of the following strategies: We immediately observe that the last two methods achieve smaller condition numbers, often between 30 and 40 for Algorithm 3.1, and between 20 and 30 for Algorithm 4.1, compared to i.i.d points, which generally give values larger than 40, and sometimes much larger, especially in the case of uniform points. Heuristically, we expect the condition number of Algorithm 4.1 to be in average a bit larger than the eigenvalues of

E(A m ) λ min (A m ) (1 -κ) -1 mδ(1 -δ) -1
(n -1)(δr -1) I = 2 (1 -1/ √ r) 2 I 22.9 I, where we used arguments from the proof of Theorem 1.1, equation (2), together with the fact that r = 256 127 . Therefore, there is a very good agreement between the theoretical and numerical spectral properties of A m .

It should be mentioned that in Algorithm 3.1, the random event from Proposition 3.3 was always realized, resulting in usual least-squares instead of the conditioned version (15). Nevertheless, redrawing the whole sample when the final condition number is large, as proposed in [START_REF] Haberstich | Boosted optimal weighted least-squares[END_REF] for Christoffel points, remains a good option in practice for all sampling schemes, since it amounts to truncating the tails of the above histograms.

In Figure 2, we plot the number of rejections as a function of the iteration index i = 1, . . . , m, averaged over the 400 runs, for Algorithms 3.1 and 4.1. In accordance with Remark 5.2, few rejections are observed in the first iterations i < n. However, the number of rejections remains moderate in the second half n i m, staying below 15 in the case of Algorithm 3.1, which is much better than the pessimistic upper bound n = 128 from Section 5. Algorithm 4.1 incurs about 7 times more rejections, which is again smaller than the factor (1 -κ) -1 (1 -δ) -2 23 encountered in the analysis.

Finally, in Figure 3, we keep the same parameters d = 4, y = (0.9, 0.8, 0.7, 0.6) and n = 128, and plot the normalized error E emp ( g y -P m n g y L 2 )/E * as a function of the number of samples m, for m = n + 2p and p = 0, . . . , 40. Here E emp stands for an empirical average over 100 runs, and E * 0.402882 is the error of best approximation, computed by (22). For the sake of clarity, we only display the results for the Christoffel measure Figure 2: Average over 400 runs of the number of rejections in each iteration and our algorithms. We again observe that the latter perform better, especially in the regime m ≈ n, and that they are within a factor 2 of the optimal error as soon as m n + 30. Remark 6.2. In view of the discussion from Section 5, the last three schemes necessitate sampling from Christoffel measure, which can be efficiently implemented by first drawing k uniformly in Λ n , and then drawing x from the measure

|L k (x)| 2 dµ(y) = d j=1 |L k j (x j )| 2 dx j 2 .
For the second step, exploiting the product structure of L k , it suffices to draw each component x j independently from |L k j (x)| 2 dx 2 . When k j = 0, this is just the uniform measure over [-1, 1]. Otherwise, we rely once more on acceptance/rejection from the arcsine measure, that is from the cosine of uniform points. Thanks to the so-called Berstein inequality for Legendre polynomials

|L k (x)| 2 2 2k + 1 k 1 π √ 1 -x 2 ,
the acceptance probability is at least k j /(2k j + 1) 1/3 for k j 1, so there are at most 2 supp(k) 2 min(n, d) rejections in average. As a conclusion, the time needed to generate a point from the Christoffel probability measure is t = O(d min(n, d)). We refer to Section 5 of [START_REF] Cohen | Optimal weighted least-squares methods[END_REF], as well as [START_REF] Cohen | Optimal sampling and christoffel functions on general domains[END_REF] and [START_REF] Migliorati | Multivariate approximation of functions on irregular domains by weighted least-squares methods[END_REF], for an overview of efficient sampling strategies on more general domains.

Aknowledgement: The authors would like to thank Albert Cohen for insightful feedback and discussions all along the elaboration of the paper, David Krieg, Mario Ullrich and Tino Ullrich for their enriching questions and comments, and the reviewers for their careful reading and valuable suggestions.

Proposition 4. 3 .

 3 Let m n and define r = m/(n -1). The random matrix A m generated by Algorithm 4.1 verifies

  established in [CCM + 15]. We note also the plain linear estimates C Λ n where C Λ = max ν∈Λ d j=1 (1 + 2ν j ) and C Λ = max ν∈Λ 2 |ν| 0

  (a) The points are i.i.d according to the uniform measure dµ(x) (b) The points are i.i.d according to the tensor product arcsine measure d j=1 dx j π 1 -x 2 j (c) The points are i.i.d according to the Christoffel measure 1 n k∈Λn |L k (y)| 2 dµ(x) (d) The points and weights are generated by Algorithm 3.1 with input parameters ε = r -1/4 and γ = r 1/2 -r 1/4 , where r = (m + 1)/n (e) The points and weights are generated by Algorithm 4.1 with input parameters δ = r -1/2 and κ = 1/2, where r = m/(n -1).

Figure 1 :

 1 Figure 1: Histograms of condition number of A m for various random sampling strategies

Figure 3 :

 3 Figure 3: Ratio between the least-squares error and the best approximation error for n = 128 and m ranging between 128 and 208, when the points and weights are i.i.d according to the Christoffel measure, or generated by Algorithms 3.1 and 4.1.

We assume here that Vn is included in L ∞ , the space of bounded functions.

The author acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Project number 442047500 through the Collaborative Research Center Sparsity and Singular Structures (SFB 1481).