Reda Khoufache
email: reda.khoufache@uvsq.fr

Mustapha Lebbah
email: mustapha.lebbah@uvsq.fr

Hanene Azzag
email: azzag@univ-paris13.fr

Etienne Goffinet
email: etienne.goffinet@tii.ae

Djamel Bouchaffra
email: djamel.bouchaffra@gmail.com

Distributed Collapsed Gibbs Sampler for Dirichlet Process Mixture Models in Federated Learning

Keywords: Federated learning, Distributed computing, Dirichlet process mixture models, Markov Chain Monte Carlo, Bayesian non-parametric

Dirichlet Process Mixture Models (DPMMs) are widely used to address clustering problems. Their main advantage lies in their ability to automatically estimate the number of clusters during the inference process through the Bayesian non-parametric framework. However, the inference becomes considerably slow as the dataset size increases. This paper proposes a new distributed Markov Chain Monte Carlo (MCMC) inference method for DPMMs (DisCGS) using sufficient statistics. Our approach uses the collapsed Gibbs sampler and is specifically designed to work on distributed data across independent and heterogeneous machines, which habilitates its use in horizontal federated learning. Our method achieves highly promising results and notable scalability. For instance, with a dataset of 100K data points, the centralized algorithm requires approximately 12 hours to complete 100 iterations while our approach achieves the same number of iterations in just 3 minutes, reducing the execution time by a factor of 200 without compromising clustering performance. The code source is publicly available at https://github. com/redakhoufache/DisCGS.

Introduction

Clustering is an unsupervised learning method that aims to partition data into clusters such that elements of the same cluster are similar while those of different clusters are dissimilar. It plays a crucial role in various domains, such as data analysis, pattern recognition, and data mining. To address clustering problems, mixture models have emerged as popular generative probabilistic models. They assume the existence of a mixture of distributions over the observation space, where each cluster is associated with a latent component distribution. The flexibility of mixture models allows researchers and practitioners to analyze and interpret data effectively.

Dirichlet Process Mixture Models (DPMMs) represent an extension of mixture models into Bayesian non-parametric models. DPMMs assume a mixture model with an infinite number of latent components and make a prior distribution over the model parameters. In literature, two common processes are used to represent the Dirichlet Process (DP). The first is the Chinese Restaurant Process (CRP) [START_REF] Antoniak | Mixtures of Dirichlet Processes with Applications to Bayesian Nonparametric Problems[END_REF], which induces a distribution over the space of partitions. The second representation is known as the Stick-Breaking (SB) process [START_REF] Sethuraman | A constructive definition of dirichlet priors[END_REF], which provides a construction of the Dirichlet Process distribution. In [START_REF] Miller | An elementary derivation of the chinese restaurant process from sethuraman's stickbreaking process[END_REF], the authors proved that each of these representations could be derived from the other.

The Gibbs Sampling algorithm [START_REF] Neal | Markov chain sampling methods for dirichlet process mixture models[END_REF] is a famous MCMC method that performs the inference of DPMM's parameters. This algorithm iteratively updates the membership vector and the parameters associated with each cluster. When the prior conjugacy assumption holds, it becomes possible to integrate out the parameters. Consequently, the inference process only requires sampling the memberships, leading to a variant known as collapsed Gibbs sampling.

While the DPMM has the advantage of automatically estimating the number of clusters and discovering new structures, its inference process becomes prohibitively slow and non-scalable when the number of observations is large. This poses a significant inconvenience for handling massive datasets and limits the applicability of the DPMM in real-world scenarios.

Distributed computing consists in distributing data across multiple workers; the workers could be components in the edge network, cores in a cluster, computers in a network, processors in a computer, etc. This distribution enables them to execute parallel computations independently. It presents an opportunity to accelerate computations and circumvent memory limitations, making it particularly suitable for handling large datasets.

Federated Learning (FL), a pioneering concept introduced by [START_REF] Konečný | Federated optimization: Distributed machine learning for on-device intelligence[END_REF], revolutionized large-scale distributed learning by harnessing the machine intelligence residing in local devices. Federated learning is used to train machine learning models directly on edge devices or computing nodes. In the following, we will refer to a component of the edge network or a network node as a "worker". One of the defining characteristics of FL is that the data distributions may possess very different properties across the workers. Hence, any potential FL clustering method is explicitly required to be able to work under the heterogeneous data setting. Therefore, the availability of DPMM models tailored to these distributions becomes an indispensable requirement in our case. However, due to its highly decentralized nature, FL poses numerous statistical and computational challenges. In addition, the complexity of the clustering problem increases significantly when applied in a federated context. Furthermore, in many real-world applications, data is actually stored at the edge, spanning nodes, phones, and computers.

In this paper, we address a specific DPMM problem that arises in federated learning, contributing to the advancement of this exciting field. We propose a new distributed MCMC-based approach for DPMM inference (DisCGS) using sufficient statistics. It is important to emphasize that our approach focuses on a specific inference algorithm for the DPMM, the Collapsed Gibbs Sampler proposed in [START_REF] Neal | Markov chain sampling methods for dirichlet process mixture models[END_REF]. The scalability of this algorithm has not been addressed yet.

We highlight our key contributions as follows: (1) In the federated context, we have developed a new formulation of distributed MCMC inference of the DPMM using the Master/Worker architecture. The data is evenly distributed among the workers, ensuring a balanced workload. The workers operate independently and do not share information with each other.

Instead, they communicate solely with the master, exchanging minimal information. (2) Each worker executes a local collapsed Gibbs sampler, which enables the discovery of local clusters and the inference of a local DPMM. Then, the sufficient statistics associated with each local cluster are sent to the master. (3) At the master level, a global collapsed Gibbs sampler is performed using only sufficient statistics. This allows the synchronization and estimation of the global DPMM and global clustering structure without accessing the content of each cluster. Theoretical background and computational details are provided. The general workflow of our model is illustrated in Figure 1. (4) We conduct several experiments to validate our method's merits and scalability with the number of nodes on large datasets. For a dataset of 100K data points, the centralized collapsed Gibbs sampler takes 12 hours to perform 100 iterations, whereas our approach does in just 3 minutes. (5) Distributing the collapsed Gibbs sampler for DPMMs is a crucial intermediate step towards the distribution of Bayesian non-parametric coclustering [START_REF] Meeds | Nonparametric bayesian biclustering[END_REF] and Multiple coclustering [START_REF] Goffinet | Multivariate time series multicoclustering. application to advanced driving assistance system validation[END_REF].

Related Work

Clustering problems in federated learning have drawn attention in recent studies. For instance, [START_REF] Kumar | Federated k-means clustering: A novel edge ai based approach for privacy preservation[END_REF] introduced an approach for federated K-Means clustering, [START_REF] Stallmann | Towards federated clustering: A federated fuzzy c-means algorithm (FFCM)[END_REF] proposed a Federated Fuzzy c-Means algorithm. In [START_REF] Jichan | Federated unsupervised clustering with generative models[END_REF], the authors presented a federated clustering method using probabilistic models.

Several parallelized and distributed approaches have been proposed to address the slow inference of DPMMs. [START_REF] Lovell | Clustercluster: Parallel markov chain monte carlo for dirichlet process mixtures[END_REF] introduced a reparameterization of the Dirichlet Process that facilitates the learning of clusters representing the data and super-clusters defining the granularity of parallelization. [START_REF] Williamson | Parallel Markov chain Monte Carlo for nonparametric mixture models[END_REF] incorporated auxiliary variables into the Dirichlet Process and Hierarchical Dirichlet Process to establish a conditional independence structure, enabling a parallel Gibbs sampler without the need for approximations. However, [START_REF] Gal | Pitfalls in the use of parallel inference for the dirichlet process[END_REF] demonstrated that these approaches are impractical due to extremely unbalanced data distributions. Another parallelized MCMC sampler was proposed by [START_REF] Chang | Parallel sampling of dp mixture models using sub-cluster splits[END_REF], which combines a non-ergodic, restricted Gibbs sampler in the workers with split/merge proposals [START_REF] Jain | A split-merge markov chain monte carlo procedure for the dirichlet process mixture model[END_REF] at the master level to ensure ergodicity. However, this method is implemented only for a single machine. This approach is extended in [START_REF] Dinari | Distributed mcmc inference in dirichlet process mixture models using julia[END_REF] to a distributed version across multiple machines using a distributed-memory model.

In [START_REF] Ge | Distributed inference for dirichlet process mixture models[END_REF], the authors proposed a distributed inference algorithm for DPMMs based on a slice sampler developed by [START_REF] Walker | Sampling the dirichlet mixture model with slices[END_REF]. They addressed the scalability by developing a distributed architecture that estimates the Dirichlet Process and Hierarchical Dirichlet Process mixture models. [START_REF] Wang | Scalable estimation of dirichlet process mixture models on distributed data[END_REF] presented a scalable estimation method for DPMMs in a distributed architecture. Their approach allows the creation of new components locally in individual computing nodes and proposes a probabilistic consolidation scheme for merging the created components. In the [START_REF] Meguelati | Dirichlet process mixture models made scalable and effective by means of massive distribution[END_REF], the authors introduced another distributed Markov Chain Monte Carlo (MCMC) inference method based on the Gibbs sampler. Their approach performs a local Gibbs sampler to infer new local clusters and a new Gibbs sampler on the means of each cluster to synchronize the clustering. However, this method assumes a Dirichlet Process Gaussian mixture model with known variances. This assumption emphasizes known variances within and between clusters due to the prior, which limits the method's applicability.

The concept of using sufficient statistics is wellknown and has been widely used to distribute the DPMM inference in many works, including [START_REF] Dinari | Distributed mcmc inference in dirichlet process mixture models using julia[END_REF][START_REF] Ge | Distributed inference for dirichlet process mixture models[END_REF][START_REF] Meguelati | Dirichlet Process Mixture Models made Scalable and Effective by means of Massive Distribution[END_REF][START_REF] Wang | Scalable estimation of dirichlet process mixture models on distributed data[END_REF]. However, it has not been applied to distribute the Collapsed Gibbs Sampler, which constitutes the main focus of this paper. Moreover, none of the existing methods is exploitable to enhance the scalability and the distribution of Bayesian non-parametric coclustering [START_REF] Meeds | Nonparametric bayesian biclustering[END_REF] and Multiple coclustering [START_REF] Goffinet | Multivariate time series multicoclustering. application to advanced driving assistance system validation[END_REF]. Moreover, to the best of our knowledge, no prior work has been proposed on federated clustering using Dirichlet Process models.

Background

Model definition

Let n and d be two positive integers, x = (x 1 , • • • , x n ) T ∈ R n×d the observed data, where (•) T denotes the transpose operator. Let z = (z 1 , • • • , z n ) be the membership vector, where z i is a latent variable such that z i = k means that the observation x i belongs to the cluster k. The DPMM assumes that the observed data are generated according to the following model:

x i | {z i = k, θ k } i.i.d. ∼ f (x i , θ k ) , ∀i ∈ {1, • • • , n} θ k i.i.d. ∼ G 0 , ∀k ∈ {1, 2, • • • } z i i.i.d. ∼ Mult (π) , ∀i ∈ {1, • • • , n} π ∼ SB(α).
Under this model assumption, an observation x i is generated by first drawing z i from the Multinomial distribution parameterized by the vector of weights π = (π k ) ∞ k=1 (also called the mixture proportions), and then x i is sampled from f (x i , θ zi ), where f (•, θ k ) is the parameterized component distribution associated to the cluster k. The components parameters θ k follow a prior distribution G 0 (also called the base distribution); in the multivariate Gaussian case, we have θ k = (µ k , Σ k ). We let Θ = {θ k , k > 1} be the set of components parameters. The mixture proportions follow the Stick-Breaking process [START_REF] Sethuraman | A constructive definition of dirichlet priors[END_REF] parameterized by a concentration parameter α > 0. We recall that the Stick-Breaking process is defined as follows:

π k = v k k-1 k =1 (1 -v k ) , v k i.i.d. ∼ Beta (1, α) .
We have π k > 0 for every k > 1 and

∞ k=1 π k = 1.
It is proved by [START_REF] Sethuraman | A constructive definition of dirichlet priors[END_REF] that the distribution G = ∞ k=1 π k δ θ k (where θ k are sampled i.i.d according to G 0 and δ is the indicator function) follows the Dirichlet Process with a concentration parameter α and a base distribution G 0 , noted DP (α, G 0 ). In the following, we denote Ω = (α, G 0 ) the hyper-parameter set.

Inference.

The Gibbs Sampling algorithm [START_REF] Neal | Markov chain sampling methods for dirichlet process mixture models[END_REF] is a popular algorithm used for the DPMM inference based on Monte-Carlo sampling. It alternates between updating the membership vector and updating the parameters associated with each cluster. In the membership update step, the Gibbs sampling simulates the posterior distribution p(z | x, Θ, Ω) by sampling each z i from the conditional distribution p(z i | x, z -i , Θ, Ω), where z -i = {z l , l = i} is the set of the remaining memberships. In the parameter update step, the component parameter of each cluster k is updated by sampling according to the posterior distribution p(θ k |x k , G 0 ), with x k = {x i , z i = k} the set of observations that belong to cluster k. The computing of such distributions is analytically tractable when G 0 is the conjugate prior to the density f . The collapsed Gibbs sampler, which corresponds to the third algorithm proposed in [START_REF] Neal | Markov chain sampling methods for dirichlet process mixture models[END_REF] skips the parameter sampling step because under the conjugacy assumption, it is possible to integrate out the parameters θ k which allows to directly compute the predictive (prior and posterior) without the parameter's values. Each membership z i is sampled according to p(z i | z -i , x, Ω) ∝:

n k p(x i | z i = k, x k , G 0 ), existing cluster k,(1) αp(x i | Ω), new cluster, (2)
where n k is the cardinal of cluster k. The posterior predictive p(x i | z i = k, x k , G 0 ) and prior predictive p(x i | Ω) can be obtained by integrating over θ.

In the multivariate Gaussian case with unknown mean and variance, we choose the Normal Inverse Wishart [START_REF] Gelman | Bayesian data analysis[END_REF] (NIW) to be the conjugate prior G 0 , with hyper-parameters λ 0 = (µ 0 , κ 0 , Ψ 0 , ν 0 ), where µ 0 ∈ R d is the prior mean, κ 0 ∈ R is the number of prior measurements, ν 0 is the degrees of freedom and Ψ 0 ∈ R d×d is the precision matrix. In this case, the probabilities in equations 1 and 2 can be computed analytically [START_REF] Murphy | Conjugate bayesian analysis of the gaussian distribution[END_REF].

To summarize, the collapsed Gibbs sampler is an efficient MCMC method because it avoids component parameter sampling. This technique, known as Rao-Blackwellization, is due to the Rao-Blackwell theorem [START_REF] Blackwell | Conditional Expectation and Unbiased Sequential Estimation[END_REF]. This theorem ensures that the estimator's variance obtained by integrating out θ is always lower than or equal to the one of a direct estimator. This theorem remains true in collapsed Gibbs sampling [START_REF] Liu | Covariance structure of the gibbs sampler with applications to the comparisons of estimators and augmentation schemes[END_REF].

The proposed method

The collapsed version of the Gibbs sampling inference inspires our approach. The main objective of DisCGS is to make the inference scalable while keeping the MCMC's precision and the DPMM's flexibility. Our approach consists in distributing the data evenly over the workers and alternating at each iteration between two steps: Worker level and Master level. Below, we provide a theoretical detailed of the collapsed Gibbs samplers executed at each level. The workflow of our approach is described in figure 1 4.1 DisCGS at worker level We denote by x j = {x j 1 , • • • , x j n j } the set of observations assigned to the j-th worker, n j the cardinal of x j , and z j = {z j 1 , • • • z j n j } the local membership vector such that z j i = k means that the observation x j i (assigned to the worker j) belongs to the local cluster k.

In this level, the local memberships are updated one by one using the collapsed version of Gibbs sampling detailed in section 3. Each z j i is updated by sampling from p(z j i | z j -i , x j , Ω) ∝

n j k p(x j i | z j i = k, x j k , G 0 ), existing cluster k, (3) αp(x j i | Ω), new cluster, (4)
where n j k is the size of cluster k in worker j and x j k is the contents of cluster k in worker j. The posterior and prior predictive distributions are computed analytically using the same global prior G 0 , as detailed in [START_REF] Murphy | Conjugate bayesian analysis of the gaussian distribution[END_REF]. After updating the local membership vector, we compute the sufficient statistics [START_REF] Silvey | Statistical Inference[END_REF] associated with each cluster. In the multivariate Gaussian case, the sufficient statistics (T j k , S j k ) for a cluster x j k are given by [START_REF] Gelman | Bayesian data analysis[END_REF]:

T j k = 1 n j k x∈x j k x ∈ R d , (5)
S j k = x∈x j k (x -T j k )(x -T j k ) T ∈ R d×d . (6)
Finally, sufficient statistics and the sizes of each cluster are sent to the master. The DisCGS inference process at the worker level is described in Algorithm 1.

Algorithm 1 DisCGS inference at worker level 1: Input: Dataset x j , concentration parameter α and prior G 0 . 2: For i ← 1 to n j do:

3:

Remove x j i from its local cluster.

4:

Compute p(z j i | z j -i , x j , Ω) as defined by eq. 3 and 4.

5:

Sample z j i .

6:

Add x j i to its new cluster. 7: For k ← 1 to K do:

8:
Compute the sufficient statistics (T j k , S j k ) using eq. 5 and 6. 9: Output:

Sufficient statistics {(T j 1 , S j 1 ) • • • , (T j K , S j K )} and cluster sizes (n j 1 , • • • , n j K ).

DisCGS at master level

The master receives from each worker the sample size of each cluster and its associated sufficient statistics. The objective is to estimate the membership vector z = (z 1 , • • • , z n ) and update the prior hyperparameters of each cluster.

In this level, the observations are assigned by batch; a batch corresponds to a set of observations that belong to the same cluster. In fact, instead of assigning the observations one by one to their clusters, we assign a group of observations that already share the same local cluster (i.e., at the worker level) to a global cluster at the master level. Hence, the observations assigned to the same global cluster will share the same label. We sample the global membership z j h of the cluster x j h (the local cluster h of worker j) according to p(z

j h | z -j -h , x, Ω) ∝ n k p(x j h | z j h = k, x k , G 0 ), existing cluster k, (7) αp(x j h | G 0 ), new cluster, (8)
In practice, the joint posterior predictive and the joint prior predictive distributions (equations 7 en 8 respectively) are computed analytically by only using sufficient statistics, i.e. without having access to the content of cluster x j h . In fact, we have:

p x j h | Ω = π -n j h d 2 • κ d/2 0 κ j h d/2 • Γ d ν j h /2 Γ d (ν 0 /2) • |Ψ 0 | ν0/2 Ψ j h ν j h /2
where | • | is the determinant, and the hyper-parameter values (µ j h , κ j h , Ψ j h , ν j h ) are obtained as follows:

µ j h = κ 0 µ 0 + n j h T j h κ j h , κ j h = κ 0 + n j h , ν j h = ν 0 + n j h , Ψ j h = Ψ 0 + S j h + κ 0 n j h κ j h µ 0 -T j h µ 0 -T j h T ,
where T j h and S j h are the sufficient statistics obtained from the workers. Moreover, we have p(

x j h | z j h = k, x k , G 0 ) = π -dn j h 2 • κ d/2 k κ j h d/2 • Γ d ν j h /2 Γ d (ν k /2) • |Ψ k | ν k /2 Ψ j h ν j h /2
and the posterior distribution parameters (µ k , κ k , Ψ k , ν k ) associated to the global cluster k, updated from the prior:

µ k = κ 0 µ 0 + n k T k κ k , κ k = κ 0 + n k , ν k = ν 0 + n k , Ψ k = Ψ 0 + S k + κ 0 n k κ k (µ 0 -T k ) (µ 0 -T k ) T .
With T k and S k , the aggregated sufficient statistics obtained when local clusters are assigned to the same global cluster k are computed as follows:

T k = 1 n k j,h| z j h =k n j h • T j h , S k = j,h| z j h =k S j h + j,h| z j h =k n h j • T j h • T j h T -n k • T k • T T k .
The inference process at the master level is described in Algorithm 2.

Algorithm 2 DisCGS inference at master level 1: Input: Sufficient statistics, cluster sizes, α and prior G 0 . 2: For each (j, h) do:

3:
Remove x j h from its global cluster.

4:

Compute p(z j h | z -h -i , x j h , Ω) ∝ as defined by eq. 7 and 8.

5:

Sample z j h .

6:

Add x j h to its new global cluster. 7: Update the membership vector z. 8: Output: Membership vector z.

Collapsed Gibbs sampler in federated learning

In federated learning, the observations are distributed on different workers "components". This corresponds to the horizontal decomposition of the dataset. In this context, each worker is initialized with the same global model, which is identical to the model present on the server. Then, each worker updates its model using its private data through the collapsed Gibbs sampler detailed in section 4.1; this step allows the estimation of the local clusters and the local model. Then, the sufficient statistics and cluster sizes associated with each local cluster are computed and transmitted to the master "server". The master proceeds to update the global model and to estimate the global clustering structure without having access to the data. This process is achieved using the collapsed Gibbs sampler detailed in section 4.2. These updates are then shared with each worker, which allows them to update the local model with the global model. This iterative process continues alternating between the worker and master steps until the global model is fully estimated.

Experiments

To evaluate the effectiveness of our approach, we conduct three types of experiments on synthetic and realworld datasets. Firstly, we compare our distributed algorithm with other state-of-the-art clustering algorithms in terms of clustering performance and convergence rate. Secondly, we compare the execution time and clustering performance of our distributed algorithm DisCGS and the centralized CGS (Collapsed Gibbs Sampler from [START_REF] Neal | Markov chain sampling methods for dirichlet process mixture models[END_REF]) on synthetic datasets of different sizes. Lastly, we investigate the scalability of DisCGS by increasing the number of nodes while keeping the number of observations fixed. For this purpose, we execute our distributed algorithm on 10 6 data points multiple times, varying the number of cores.

Implementation settings and distributed environment

In the following experiments, we use an uninformative prior NIW for both CGS and DisCGS algorithms. Therefore, both methods are implemented by setting the NIW hyper-parameters as follows: µ 0 and the matrix precision Ψ 0 are respectively set to be empirical mean vector and covariance matrix of all data, as we want them as uninformative as possible. κ 0 and ν 0 represents our confidence in µ 0 and Ψ 0 , are set to their lowest values, which are 1 and d+1, respectively, where d is the dimension of the observation space. The initial state is a one-cluster partition. Finally, we have executed the distributed algorithms using the Neowise machine (1 CPU AMD EPYC 7642, 48 cores/CPU) hosted by the cluster grid50001 . The centralized algorithm is executed on the same machine using one core.

Clustering performance

To evaluate the clustering performance of our algorithm, we compare our approach with two distributed algorithms for the DPMM inference: M-R2 [8] and SubC3 [START_REF] Dinari | Distributed mcmc inference in dirichlet process mixture models using julia[END_REF], and two parallelized clustering algorithms: Kmeans and GMM (for both methods, we have used the Spark Mllib implementation). These two parametric methods require the number of clusters. To ensure a fair comparison of the clustering quality with our approach, we set the number of clusters in these methods equal to the one inferred by DisCGS. It is important to mention that SubC and M-R are the only approaches of distributed inference for DPMM for which opensource working code is available.

All the executions are performed on the Neowise machine by distributing the data evenly on 32 cores. In this experiment, we use 8 different datasets described in table 1. Due to the smaller size of the EngyTime dataset, we distributed the data only on two workers. The image datasets are encoded using a variational auto-encoder, and each image is encoded into an 8dimensional vector.

To evaluate the clustering performance, we compute the three clustering metrics: Adjusted Rand Index (ARI) [START_REF] Hubert | Comparing partitions[END_REF], Normalized Mutual Information (NMI) [START_REF] Strehl | Cluster ensembles -a knowledge reuse framework for combining multiple partitions[END_REF], and clustering accuracy (ACC) [START_REF] Yang | Image clustering using local discriminant models and global integration[END_REF].

Table 2 reports the mean and the standard deviation of the clustering metrics, ARI, NMI, and ACC, achieved by each method on each dataset over 10 trials. The results show that our proposed method (DisCGS) outperforms other methods or achieved the secondbest score in almost all the datasets for the three clustering metrics. It is important to note that in this experiment, we only focus on comparing the clustering performance of the different methods; we do not compare the execution times of the three methods because the other approaches proposed an inference algorithm that differs from the collapsed Gibbs sampler.

Figure 2 illustrates the estimated clusters obtained by our distributed approach on the synthetic 10K dataset. This figure represents the best partition achieved over the ten trials. The purpose of using this dataset is to evaluate DisCGS's performance when confronted with complex datasets that exhibit overlapping clusters. The results demonstrate that our method performs remarkably well even when dealing with such datasets. DisCGS achieved an ARI, NMI, and ACC score of 0.85, 0.89, and 0.91, respectively. It estimated seven clusters, while the number of clusters is six.

Convergence

In this experiment, we examine the convergence of both likelihood and ARI score of three methods: Dis-CGS, CGS (Collapsed Gibbs Sampler), and SubC. These evaluations are performed on three datasets: Synthetic 100K, Fashion-mnist, and Balanced. It's important to note that only CGS and SubC share the same model assumption as DisCGS, resulting in a comparable likelihood. Handwritten digits [START_REF] Cohen | Emnist: Extending mnist to handwritten letters[END_REF]. UrbanGB 360177 3 469 Coordinates (longitude and latitude) of road accidents [START_REF] Dua | UCI machine learning repository[END_REF].

Table 1: The description of the datasets used to evaluate the clustering performance of our distributed inference approach. n denotes the size, d the dimension, K the true number of clusters. Figure 3 illustrates the evolution of the loglikelihood and ARI score at each iteration. We observe that our algorithm converges almost at a similar rate as the centralized CGS algorithm and is much faster than SubC, which takes more iterations to converge. This is because our algorithm is able to discover new local clusters inside each worker. Whereas, in SubC, the number of clusters is fixed when performing the restricted Gibbs sampler in each worker, and new components are only discovered at the master level during the split step. Thus, more iterations are required to generate enough components to model the data. This phenomenon is also observed in [START_REF] Wang | Scalable estimation of dirichlet process mixture models on distributed data[END_REF]. Overall, our algorithm converges really fast and maintains a stable ARI score over the iterations. Whereas the CGS and SubC may downgrade their ARI score after some iterations, as can be observed for the Fashion-Mnist dataset. In this experience, we compare the execution time and clustering performance of the distributed collapsed Gibbs sampler (DisCGS) and the centralized collapsed Gibbs sampler (CGS) from [START_REF] Neal | Markov chain sampling methods for dirichlet process mixture models[END_REF], we execute both algorithms on synthetic datasets of different sizes (from n = 20K to n = 100K) generated from K = 10 Gaussian components of dimension 2. The centralized version is too slow; running over 100K observations would take too much time. We use 32 cores for this experiment. Figure 4 illustrates, using a logarithmic scale, the execution time for M = 100 iterations of both the distributed and centralized inference methods. The results show that our distributed algorithm significantly outperforms the centralized approach. For instance, when considering 100K data points, the centralized algorithm takes approximately 12 hours to complete 100 iterations, whereas our approach achieves the same number of iterations in just 3 minutes, reducing the execution time by a factor of 200. Table 3 presents the clustering metrics (ARI, NMI, and ACC) obtained by each algorithm on each dataset. The results indicate that our approach consistently achieves high scores and outperforms the centralized algorithm in almost all cases. We observe that using the dataset of size 40K, the centralized algorithm obtained slightly higher ARI and NMI scores than the distributed algorithm. This is not surprising since both approaches sample the memberships from an approximation of the posterior distribution, resulting in noisy inferred partitions. These samples can be aggregated after a given number of burn-in iterations with a consensus partition estimation. Overall, the findings confirm that the distributed inference does not compromise the clustering performance while considerably reducing the execution time. ) and centralized (Cen.) inference on synthetic datasets of different sizes.

Distributed algorithm scale-up

In this experiment, we use n = 10 6 data points generated from K = 10 components of two-dimensional Gaussian. We run our distributed inference several times by increasing the number of cores from 8 to 48. Figure 5 represents the running time as a function of the number of cores. We observe that the running time decreases when the number of cores increases, showing that our algorithm scales efficiently with the number of workers. Figure 6 shows the clusters (labels) in- ferred by our approach on 10 6 data points distributed across 32 cores. As depicted, our approach successfully identifies meaningful and coherent clusters, resulting in high ARI, NMI, and ACC scores of approximately 0.98. Similar results were observed when distributing the data on different numbers of cores. Additionally, it has inferred 14 clusters while the ground truth is 10. However, only 10 of them are significant clusters, and the 4 others are only "outliers". Moreover, the clustering performance and the number of clusters depend on the auto-encoder. Also, the number of clusters might be influenced by the concentration parameter α [START_REF] Murphy | Machine learning : a probabilistic perspective[END_REF].

Conclusion and perspectives

This article presents a novel distributed MCMC inference method, called DisCGS, for DPMMs. DPMMs are highly useful for clustering problems, especially when the number of clusters is unknown. However, the inference process of DPMMs tends to become significantly slow as the dataset size increases. To overcome this limitation, our proposed DisCGS is specifically designed to handle distributed data across independent and heterogeneous machines, making it suitable for horizontal federated learning scenarios, i.e., when the workers are the different components of the edge network. The experimental results showed highly promising outcomes. The proposed method significantly reduces the inference time while maintaining accurate results. Our ongoing research exploits the DisCGS approach to massively distribute the Non-Parametric Latent Block Model (NPLBM) [START_REF] Meeds | Nonparametric bayesian biclustering[END_REF] and the multiple Coclustering model.

Figure 1 :

 1 Figure 1: DisCGS workflow. Workers have only access to their own local data. Each worker infers a local DPMM and discovers local tables (clusters). Sufficient statistics are sent from each worker to the master, which infers the global DPMM and estimates the global partition. The global model is shared to each worker.

Figure 2 :

 2 Figure 2: Unlabeled (left) and labeled data (right), after 100 iterations of DisCGS on a synthetic dataset with overlapped clusters, ARI = 0.85, NMI = 0.89, and ACC = 0.91. The number of estimated clusters is 7.

  ∓ 0.06 0.31 ∓ 0.04 0.63 ∓ 0.00 0.54 ∓ 0.00 ACC 0.46 ∓ 0.01 0.07 ∓ 0.02 0.10 ∓ 0.02 0.44 ∓ 0.02 0.30 ∓ 0.01 Digits ARI 0.55 ∓ 0.01 0.28 ∓ 0.14 0.74 ∓ 0.05 0.46 ∓ 0.02 0.36 ∓ 0.01 NMI 0.71 ∓ 0.01 0.51 ∓ 0.14 0.79 ∓ 0.02 0.71 ∓ 0.01 0.63 ∓ 0.00 ACC 0.58 ∓ 0.00 0.38 ∓ 0.09 0.81 ∓ 0.05 0.44 ∓ 0.03 0.34 ∓ 0.01 UrbanGB ARI 0.63 ∓ 0.01 0.12 ∓ 0.05 0.09 ∓ 0.00 0.67 ∓ 0.05 0.49 ∓ 0.06 NMI 0.68 ∓ 0.01 0.21 ∓ 0.07 0.24 ∓ 0.00 0.77 ∓ 0.01 0.81 ∓ 0.01 ACC 0.45 ∓ 0.01 0.29 ∓ 0.02 0.29 ∓ 0.00 0.53 ∓ 0.02 0.54 ∓ 0.02 Table 2: The mean and the standard deviation of the three clustering metrics, ARI, NMI, and ACC, over 10 runs on different datasets. The best result within each row is marked as bold, and the runner-up is underlined.

Figure 3 :

 3 Figure 3: Log-likelihood and ARI score every iteration 5.4 Comparison of the distributed and centralized collapsed Gibbs sampler

Figure 4 :

 4 Figure 4: Running time (in logarithmic scale) for 100 iterations of DisCGS and CGS inference for DPMM on synthetic datasets of different sizes.

Figure 5 :

 5 Figure 5: Running time (hours) as a function of the number of cores of DisCGS on 10 6 data points.

Figure 6 :

 6 Figure 6: Unlabeled data (left) and labeled data (right), after M = 100 iterations of DisCGS on 10 6 data points, using 32 cores, ARI = 0.98, NMI = 0.98, and ACC = 0.98. The number of inferred clusters is 14.

  ∓ 0.07 0.47 ∓ 0.04 0.87 ∓ 0.00 0.73 ∓ 0.13 0.54 ∓ 0.20 NMI 0.92 ∓ 0.07 0.42 ∓ 0.02 0.79 ∓ 0.00 0.67 ∓ 0.09 0.56 ∓ 0.11 ACC 0.96 ∓ 0.04 0.75 ∓ 0.04 0.97 ∓ 0.00 0.85 ∓ 0.11 0.68 ∓ 0.20 Synthetic 10K ARI 0.80 ∓ 0.06 0.10 ∓ 0.01 0.38 ∓ 0.07 0.80 ∓ 0.07 0.75 ∓ 0.03 NMI 0.86 ∓ 0.04 0.12 ∓ 0.01 0.61 ∓ 0.08 0.88 ∓ 0.04 0.85 ∓ 0.01 ACC 0.88 ∓ 0.06 0.29 ∓ 0.01 0.38 ∓ 0.08 0.85 ∓ 0.08 0.76 ∓ 0.05 Mnist ARI 0.72 ∓ 0.01 0.20 ∓ 0.07 0.66 ∓ 0.04 0.39 ∓ 0.02 0.26 ∓ 0.01 NMI 0.74 ∓ 0.00 0.38 ∓ 0.08 0.79 ∓ 0.01 0.69 ∓ 0.01 0.62 ∓ 0.00 ACC 0.79 ∓ 0.01 0.30 ∓ 0.06 0.71 ∓ 0.03 0.38 ∓ 0.02 0.23 ∓ 0.01

	Dataset	DisCGS	M-R	SubC	GMM	Kmeans
	EngyTime 0.94 Fashion-Mnist ARI ARI 0.45 ∓ 0.02 0.35 ∓ 0.02 0.40 ∓ 0.02 0.41 ∓ 0.02 0.37 ∓ 0.02 NMI 0.60 ∓ 0.01 0.54 ∓ 0.01 0.60 ∓ 0.01 0.59 ∓ 0.02 0.57 ∓ 0.01
		ACC 0.55 ∓ 0.02 0.45 ∓ 0.03 0.48 ∓ 0.01 0.52 ∓ 0.04 0.48 ∓ 0.01
	Letter	ARI NMI 0.56 ∓ 0.01 0.23 ∓ 0.05 0.47 ∓ 0.06 0.61 ∓ 0.00 0.52 ∓ 0.00 0.30 ∓ 0.01 0.07 ∓ 0.03 0.23 ∓ 0.05 0.30 ∓ 0.01 0.20 ∓ 0.00
		ACC 0.41 ∓ 0.01 0.14 ∓ 0.04 0.31 ∓ 0.06 0.33 ∓ 0.01 0.24 ∓ 0.01
	Balanced	ARI NMI 0.59 ∓ 0.00 0.17 0.35 ∓ 0.00 0.02 ∓ 0.01 0.05 ∓ 0.02 0.35 ∓ 0.01 0.22 ∓ 0.00

Table 3 :

 3 Clustering metrics ARI, NMI, and ACC obtained by the distributed (Dis.

	Dataset size	ARI Dis. Cen. Dis. Cen. Dis. Cen. NMI ACC

20K 0.99 0.89 0.99 0.96 0.99 0.89 40K 0.96 0.99 0.97 0.99 0.99 0.97 60K 0.91 0.89 0.92 0.96 0.92 0.89 80K 0.94 0.89 0.96 0.96 0.91 0.89 100K 0.91 0.89 0.94 0.89 0.91 0.89

https://www.grid5000.fr/w/Grid5000:Home

The code source is taken from: https://github.com/ wangruohui/distributed-dpmm

The code source is taken from: https://github.com/ BGU-CS-VIL/DPMMSubClusters.jl/tree/master

Acknowledgments

This work has been supported by the Paris Île-de-France Région in the framework of DIM AI4IDF.


Mary E Savino

Céline Lévy-Leduc

A NOVEL VARIABLE SELECTION METHOD IN NONLINEAR MULTIVARIATE MODELS USING B-SPLINES WITH AN APPLICATION TO GEOSCIENCE

Keywords: Variable selection, nonparametric regression, B-splines, Group Lasso. Andra

In this article, we introduce a novel data-driven variable selection approach in a multivariate nonparametric regression model designed to capture only the variables on which the regression function depends. The core concept of our method consists in approximating the underlying function by a linear combination of B-splines of order M and their pairwise interactions. The coefficients of this linear combination are estimated by minimizing the standard least-squares criterion penalized by the sum of the 2-norms of the partial derivatives with respect to the different variables on which the function depends. We demonstrate that our proposed method can be formulated as a Group Lasso problem, aiming to discard irrelevant variables for which the corresponding coefficients are close to zero. We validate our approach through numerical experiments varying the number of observations, the noise level and the total number of variables and compared it to two other state-ofthe-art methods. An application to a geochemical system based on calcite precipitation is also explored. In these different contexts, our approach exhibits better performance than the others. Our completely data-driven method is implemented in the absorber R package which is available on the Comprehensive R Archive Network (CRAN).

Introduction

The simulation of geochemical models that incorporate precipitation and dissolution reactions of minerals coupled to other physical processes represents a challenging task. Reactive transport modeling (RTM) serves as an illustration, striving to simultaneously consider geochemical reactions, fluid flow, heat transfer and solute transport, see [START_REF] Steefel | Reactive transport at the crossroads[END_REF] for various applications. To achieve feasible computational times for these simulations, oversimplifications of the model are usually required. Despite significant improvements in the past few decades, solving three-dimensional-large-scale modelling of complex reactive transport over many time steps remains nearly impossible using standard computers. This challenge has led to the development of Machine Learning (ML)-based approaches aimed at estimating real solutions for full simulation models through the use of surrogate models. The main idea here consists in solving the transport equations explicitly and approximating solutions for geochemical reactions at equilibrium using surrogate models at each time step. A wealth of reviews and surveys on surrogate models for RTM is available in the works of [START_REF] Razavi | Review of surrogate modeling in water resources[END_REF][START_REF] Asher | A review of surrogate models and their application to groundwater modeling[END_REF][START_REF] Jatnieks | Data-driven surrogate model approach for improving the performance of reactive transport simulations[END_REF][START_REF] Lary | Machine learning in geosciences and remote sensing[END_REF]. Among these models, Artificial Neural Networks (ANN) have gained prominence, see for instance [START_REF] Guérillot | Geochemical equilibrium determination using an artificial neural network in compositional reservoir flow simulation[END_REF]Prasianakis et al., 2020;[START_REF] Laloy | Speeding up reactive transport simulations in cement systems by surrogate geochemical modeling: deep neural networks and k-nearest neighbors[END_REF][START_REF] Demirer | Improving the performance of reactive transport simulations using artificial neural networks[END_REF]. These surrogate models can provide highly accurate approximations and with optimized hyperparameters, they can outperform other surrogate methods [START_REF] Laloy | Emulation of cpu-demanding reactive transport models: a comparison of gaussian processes, polynomial chaos expansion, and deep neural networks[END_REF]. However, the computational efficiency of ANN in reducing simulation time comes at the cost of requiring a large dataset and often demands extensive CPU times for training and tuning the hyperparameters [START_REF] Karpatne | Machine learning for the geosciences: Challenges and opportunities[END_REF]. Furthermore, an approach based on an on-demand training algorithm allows to train the model at runtime to iteratively build the training dataset [START_REF] Leal | Ultra-fast reactive transport simulations when chemical reactions meet machine learning: chemical equilibrium[END_REF]. Analogously, an active learning approach has also been introduced to RTM [START_REF] Savino | An active learning approach for improving the performance of equilibrium based chemical simulations[END_REF] to drastically diminish the dataset size while insuring good approximation accuracy. Finally, a novel approach based on B-splines and on an adaptive knot selection was proposed in [START_REF] Savino | A novel approach for estimating functions in the multivariate setting based on an adaptive knot selection for B-splines with an application to a chemical system used in geoscience[END_REF] to improve the approximation accuracy while having only a few parameters to tune.

Another approach to improve the surrogate model accuracy while reducing the CPU times is to reduce the number of input variables to consider in the model. This can be reformulated as a variable selection problem in the following framework.

Let us consider that we have n observations satisfying the following nonparametric regression model:

Y i = f (x i ) + ε i , x i = (x (1)
i , . . . , x

(p) i ) ∈ R p , 1 ≤ i ≤ n
(1) where f is an unknown real-valued function and where the ε i 's are i.i.d centered random variables of variance σ 2 . We will also assume that f actually depends only on d variables instead of p, with d < p, which means that there exists a real-valued function f such that f (x) = f ( x), where x ∈ R p and x ∈ R d . Variable selection consists in identifying the components of x.

Efficient methodologies have been devised over the last few decades particularly when the variables x i s and Y i s in (1) are linearly related. Notable examples include the Lasso regression formulated by [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] and one of its variant the Elastic Net defined by [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF]. However, dealing with the nonlinearity of the relationship between the x i s and the Y i s in (1) poses a greater challenge.

In their paper, [START_REF] Yamada | High-dimensional feature selection by feature-wise kernelized lasso[END_REF] introduced a feature-wise kernelized Lasso method tailored for variable selection in nonlinear models: HSIC-Lasso and its variation NOCCO-Lasso. This approach employs the kernel trick within the regular Lasso problem in combination with kernel-based independence measures to discern and selectively choose relevant variables. While these methods offer scalability to high-dimensional variable selection problems, they exhibit sensitivity to the number of observations. Tree-based methods such as Random Forests introduced by [START_REF] Breiman | Random forests[END_REF] are widely used for variable selection in regression models since they are well-suited to describe nonlinearity in (1). Numerous applications and challenges associated with its use for variable selection are discussed in [START_REF] Genuer | Variable selection using random forests[END_REF].

More recently, ANN have gained interest for variable selection and regression. We will present just a few examples of them. For instance, [START_REF] Liang | Bayesian neural networks for selection of drug sensitive genes[END_REF] developed a method based on a Bayesian neural network architecture to select variables for which the marginal inclusion probability exceeds a predefined threshold. Furthermore, regularized approaches with different ANN architectures were proposed as seen in the work of [START_REF] Li | Deep feature selection: theory and application to identify enhancers and promoters[END_REF] where a ridge regularization approach is considered for the weights of the first and hidden layers. Similarly, [START_REF] Feng | Sparse-input neural networks for high-dimensional nonparametric regression and classification[END_REF] introduced the SPINN method which is a single-layer neural network with a sparse group lasso regularization to shrink the weights corresponding to the units of the irrelevant variables. They later proposed another approach with deeper neural networks in [START_REF] Feng | Ensembled sparse-input hierarchical networks for highdimensional datasets[END_REF]. [START_REF] Ye | Variable selection via penalized neural network: a drop-out-one loss approach[END_REF] adapted the SPINN method by adding a greedy elimination algorithm to iteratively drop one variable at a time and determine if the empirical loss decreases. Finally, [START_REF] Lemhadri | Lassonet: A neural network with feature sparsity[END_REF] presented LassoNet, a residual feed-forward neural network architecture as introduced in [START_REF] He | Deep residual learning for image recognition[END_REF] which incorporates a regularization approach based on Lasso regression to selectively use a subset of the features in the network. Similar work on variable selection with ANN can be found in [START_REF] Chen | Nonlinear variable selection via deep neural networks[END_REF][START_REF] Lu | Deeppink: reproducible feature selection in deep neural networks[END_REF][START_REF] Zhu | Deep-gknock: nonlinear group-feature selection with deep neural networks[END_REF]. All these methods display interesting results especially for high-dimensional problems but often suffer from high training CPU times and a large number of hyperparameters to tune.

Another research direction focuses on developing flexible and interpretable methods using splines for piecewise polynomial fitting such as Multivariate Adaptive Regression Splines introduced by [START_REF] Friedman | Multivariate adaptive regression splines[END_REF]. This method enables the description of interactions and nonlinear relations by automatically pruning the most irrelevant terms. In the same vein, [START_REF] Lin | Component selection and smoothing in multivariate nonparametric regression[END_REF] developed COSSO, a regularization method for component selection and smoothing splines where the penalty term is the sum of the component norms. This approach can be considered as a more generalized form of the Lasso approach with a Reproducing Kernel Hilbert Space (RKHS) constraint. Compared to MARS, it shows better results except for higher dimensional cases with small datasets. Similarly, [START_REF] Ravikumar | Sparse additive models[END_REF] proposed sparse additive models (SpAM) based on a generalized additive model with a 2 -norm regularization. A few articles have proposed considering a sparse additive model using a linear combination of B-splines of order M (M ≥ 1), introduced by [START_REF] De Boor | A practical guide to splines[END_REF] in Chapter 9. Their ability to approximate complex functions without being significantly altered by the presence of noise has made them very attractive in the past few decades. As an illustration, [START_REF] Huang | Variable selection in nonparametric additive models[END_REF] approximated the underlying function f in (1) using an additive B-spline estimator and subsequently, employed an adaptive group lasso approach for variable selection. In their study, they presented both numerical applications and theoretical results regarding the selection consistency of their proposed method. [START_REF] Antoniadis | Variable selection in additive models using P-splines[END_REF] leveraged the benefits of Bsplines by incorporating a penalized version known as P-splines, introduced by [START_REF] Eilers | Flexible smoothing with B-splines and penalties[END_REF], and compared their results with various adaptations of COSSO. Additional references on variable selection with P-splines can be found in the review by [START_REF] Gijbels | Variable selection using P-splines[END_REF]. While these approaches have proven to be efficient for high-dimensional nonparametric additive models, they fall short in describing interactions that may exist in real datasets. Therefore, [START_REF] Radchenko | Variable selection using adaptive nonlinear interaction structures in high dimensions[END_REF] extended the SpAM approach to consider both single and interaction terms, aiming to construct a more interpretable approximation of f . This presented approach, known as VANISH, strongly penalizes interaction terms to simplify the model as much as possible and has demonstrated efficiency for small datasets. In parallel, [START_REF] Rosasco | A regularization approach to nonlinear variable selection[END_REF] proposed a novel method for variable selection based on a regularized least-square estimator penalizing large values of the partial derivatives to select the most relevant variables in a multivariate nonlinear regression model with a RKHS constraint.

In this article, we propose a novel method for variable selection motivated by (Radchenko and James, 2010) using a multivariate nonparametric regression model to retrieve the d relevant variables on which f in (1) truly depends. Our approach involves approximating f using a linear combination of B-splines and their pairwise interactions. Additionally, drawing inspiration from the methodology of [START_REF] Rosasco | A regularization approach to nonlinear variable selection[END_REF], the coefficients of the linear combination are estimated by minimizing the usual least-squares criterion penalized by the sum of the 2 -norms of the partial derivatives with respect to the different variables on which f depends. We will demonstrate that our proposed method can be formulated as a Group Lasso problem defined by [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF] and thus, can be easily implemented. Two different approaches to choose the penalization parameter will be presented to the reader.

This paper is organized as follows. Section 2 presents the methodology that we propose for variable selection in nonlinear models. Section 3 investigates the performance of our approach through numerical experiments. Finally, in Section 4, we apply our method to a real geochemical application that motivated this study.

Methodology

2.1. Approximation of f using B-splines. Let us first recall how the B-spline basis associated to a given dimension among the p, the th for instance, is defined.

Let t = (t ,1 , . . . , t ,K ) be a set of K points called knots and let S be a compact subset of R. Following (De Boor, 1978, p. 89-90) and (Hastie et al., 2009, p. 160), the augmented knot sequence τ is defined as follows: k,m the kth B-spline basis function of order m with m ≤ M for the knot sequence τ and for the dimension , B-splines are defined by the following recursion:

τ ,1 = . . . = τ ,M = x ( ) min , τ ,j+M = t ,j , j = 1, . . . , K, τ ,K+M +1 = . . . = τ ,K+2M = x ( ) max , τ = (τ ,1 , . . . , τ ,K+2M ) = x ( ) min , . . . ,
B ( ) k,1 (x ( ) ) = 1 if τ ,k ≤ x ( ) < τ ,k+1 0 otherwise for k = 1, . . . , K + 2M -1, (2)
and for 2 ≤ m ≤ M ,

B ( ) k,m (x ( ) ) = x ( ) -τ ,k τ ,k+m-1 -τ ,k B ( ) k,m-1 (x ( ) ) + τ ,k+m -x ( ) τ ,k+m -τ ,k+1 B ( ) k+1,m-1 (x ( ) ), (3)
for k = 1, . . . , (K + 2M -m).

Inspired by [START_REF] Radchenko | Variable selection using adaptive nonlinear interaction structures in high dimensions[END_REF], we propose approximating the function f (x (1) , . . . , x (p) ) appearing in (1) by a linear combination of B-splines of each variable x (1) , . . . , x (p) and of pairwise interaction of them as follows:

F (x (1) , . . . , x (p) ) = p =1 K+M k=1 β ( ) k B ( ) k (x ( ) ) + p-1 =1 p j= +1   K+M k=1 K+M q=1 β ( ,j) k,q B ( ) k (x ( ) )B (j) q (x (j) )   , (4) where B
( ) 2) and ( 3) and where β

k = B ( ) k,M is defined in (
( ) k and β ( ,j)
k,q are unknown coefficients. Observe that the column vector (F (x 4) can be rewritten as follows:

(1) i , . . . , x (p) i )) 1≤i≤n (
p =1 Ψ β + p-1 =1 p j= +1 Φ j β ,j . ( 5
)
where Ψ is a n×(K +M ) matrix such that its ith row is equal to (B ( )

1 (x ( ) i ), . . . , B ( ) K+M (x ( ) i )) and β = β ( ) 1 . . . β ( ) K+M T for 1 ≤ ≤ p, A T denoting the transpose of the matrix A. Moreover, Φ j is an n × (K + M ) 2 matrix such that its ith row satisfies (Φ j ) i,• = ((Ψ ) i,• ⊗ (Ψ j ) i,• ), ⊗ denoting the Kronecker product, (Ψ ) i,• denoting the ith row of Ψ and β ,j = β ( ,j) 1,1 β ( ,j) 1,2 . . . β ( ,j) K+M,K+M T for 1 ≤ < j ≤ p.
2.2. Description of our variable selection method. Inspired by the methodology of [START_REF] Rosasco | A regularization approach to nonlinear variable selection[END_REF], we propose selecting the variables on which f depends by estimating the coefficients β and β ,j appearing in (5) by minimizing the following regularized criterion:

β 1 (λ ), . . . , β p (λ ), β 1,2 (λ ), . . . , β (p-1),p (λ ) = argmin (β 1 ,...,β p ) (β 1,2 ,...,β (p-1),p )   Y - p =1 Ψ β - p-1 =1 p j= +1 Φ j β ,j 2 2 + λ p =1 n i=1 ∂ F (x i ) 2   , where Y = (Y 1 , . . . , Y n ), the Y i 's being defined in (1), ∂ F (x i ) denotes the th partial derivative of F defined in (4) at some observation point x i = (x (1) i , . . . , x (p) i ) and ||y|| 2 2 = n i=1 y 2 i .
Note that the idea underlying this criterion is that when a function does not depend on a variable its partial derivative with respect to this variable is equal to zero.

Using the definition of F given in ( 5) the criterion can be rewritten as follows:

β 1 (λ ), . . . , β p (λ ), β 1,2 (λ ), . . . , β (p-1),p (λ ) = argmin (β 1 ,...,β p ) (β 12 ,...,β (p-1)p )   Y - p =1 Ψ β - p-1 =1 p j= +1 Φ j β ,j 2 2 +λ p =1 Ψ β + p j= +1 (∂ Φ j ) β ,j + 1≤j< (∂ Φ j ) β j, 2   , ( 6
)
where

Ψ is the n × (K + M ) matrix such that (Ψ ) i,k = B ( ) k (x ( ) i ), B ( ) k denoting the first derivative of B ( ) k . The ith row of (∂ Φ j ) (resp. (∂ Φ j )) is defined by (∂ Φ j ) i,• = ((Ψ ) i,• ⊗ (Ψ j ) i,• ) (resp. (∂ Φ j ) i,• = ((Ψ j ) i,• ⊗(Ψ ) i,• )). By denoting (∂ Φ • ) = (∂ Φ ( +1) ) . . . (∂ Φ p ) , (∂ Φ • ) = (∂ Φ 1 ) . . . (∂ Φ ( -1) ) , β • = β ,( +1) . . . β ,p and β • = β 1, . . . β ( -1
), , the penalty term can be written as:

λ p =1 Ψ β +(∂ Φ • ) β • +(∂ Φ • ) β • 2 =: λ p =1 (∂ Θ ) γ 2 , (7)
where

γ = β T β T , +1 ... β T ,p β T 1, ... β T -1, T . Using that p-1 =1 p j= +1 Φ j β ,j = p j=2 j-1 =1 Φ j β ,j = p =2 -1 j=1 Φ j β j, ,
the least-squares term can be rewritten as follows:

Y - p =1 Ψ β - p-1 =1 p j= +1 Φ j β ,j 2 2 = Y - p =1 Ψ β - 1 2   p-1 =1 p j= +1 Φ j β ,j + p =2 -1 j=1 Φ j β j,   2 2 =: Y - p =1 Θ γ 2 2 . (8)
Equation ( 8) comes by setting

Θ 1 = Ψ 1 1 2 Φ 1• and Θ p = Ψ p 1 2 Φ •p , where Φ • = (Φ ( +1) . . . Φ p ) and Φ • = (Φ 1 . . . Φ ( -1)
). Combining ( 7) and ( 8), ( 6) can be rewritten as:

γ 1 (λ ), . . . , γ p (λ ) = argmin (γ 1 ,...,γ p ) Y - p =1 Θ γ 2 2 + λ p =1 (∂ Θ ) γ 2 . ( 9
)
By defining α = (∂ Θ ) γ and X = Θ (∂ Θ ) + , A + being the Moore-Penrose inverse of matrix A, (9) can be rewritten as:

( α 1 (λ ), . . . , α p (λ )) = argmin (α 1 ,...,αp) Y - p =1 X α 2 2 + λ p =1 α 2 . ( 10
)
The last formulation of our variable selection criterion (10) can be seen as a group lasso problem introduced by [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF], where the size p of each group belonging to {1, . . . , p} is equal to n. This approach is implemented in numerous R packages such as the most recent one sparsegl developed by [START_REF] Liang | sparsegl: An r package for estimating sparse group lasso[END_REF] that we used in the numerical experiments section. For a fixed number of parameters λ , this package provides a set of penalization parameters Λ and the coefficients α (λ ) for λ belonging to Λ. The coefficients γ (λ ) are thus obtained as follows

γ (λ ) = (∂ Θ ) + α (λ ). ( 11
)
We then define the active variables for each λ in Λ as follows:

V λ =    , k≥1 | γ ,k (λ )| = 0    , (12)
where γ ,k (λ ) is the kth coefficient of γ (λ ). We also introduce the set V f of the indices of the d relevant variables on which f in (1) actually depends that we seek to select among the p variables and the set V f of the indices of the irrelevant variables on which f does not depend.

2.3. Choice of K. Our method relies on the definition of the B-spline basis for each in {1, . . . , p} and thus on the choice of the set of knots t used for defining them. For simplifying this choice, we considered evenly spaced knots in the interval [0, 1]. For regularity purposes, we use quadratic B-splines with M = 3. Thus, we are only interested in optimizing the number of knots K. To find the best value of K, we use two sensitivity measures. Firstly, for each λ belonging to Λ, we computed the True Positive Rate (TPR) and the False Positive Rate (FPR), defined as:

TPR(λ ) = TP(λ ) d = |V λ ∩ V f | d and FPR(λ ) = FP(λ ) p -d = V λ ∩ V f p -d ,
where d < p, |A| is the cardinality of the set A, TP(λ ) and FP(λ ) are the number of true selected variables and the number of false selected variables for λ , respectively. V f , V f and V λ are introduced in the previous section. We can then draw the ROC curve where each point has as coordinates (FPR(λ ), TPR(λ )) for λ belonging to Λ.

In order to have an idea of the quality of our variable selection procedure, we calculate the Area Under Curve (AUC) of the ROC curves as well as a complementary indicator that we want to maximize:

max(TPR -FPR) = max λ ∈Λ (TPR(λ ) -FPR(λ )).
To assess the quality of our variable selection procedure according to K, we define two functions on which our method is applied:

f 1 (x (1) , x (2) , x (3) , x (4) , x (5) ) = 2B (3) 2 (x (3) )B (5) 4 (x (5) ) + 2B (5) 4 (x (5) ) -5B (3) 2 (x (3) ), (x (1) , . . . , x (5) ) ∈ [0, 1] 5 , (13)
f 2 (x (1) , . . . , x (10) ) = 1.8 cos (x (1) ) sin (x (7) + 1) -5 ln (x (3) + 1) -0.9

x (10) 2 + 1 , (x (1) , . . . , x (10) ) ∈ [0, 1] 10 . (14)
In ( 13), the B-spline bases are defined using t = (0.2, 0.5, 0.6, 0.75, 0.8) for each belonging to {1, . . . , 5}. Here, V f 1 = {3, 5} and V f 2 = {1, 3, 7, 10}. Results for the two metrics defined above, AUC and max(TPR-FPR), are shown for f = f 1 and f = f 2 for 10 random samplings of the observation set and for σ = 0 and σ = 0.25 in Figure 1 and in Figure 3 for n = 700 and n = 2000, respectively. Firstly, we can clearly see for n = 700 in Figure 1 that for f 1 , all the values of K are satisfying as they allow us to get max(TPR -FPR) = 1 and AUC = 1. However, for f = f 2 , we do not have necessarily AUC = 1 when max(TPR -FPR) = 1 for instance for K = 1, which does not imply that this method does not select properly the relevant variables. To illustrate this idea, one can relate to Figure 2 in which the ROC curves for f = f 2 are drawn for each value of K belonging to {1, . . . , 10}. Here, we can observe that a good variable selection method will not necessarily lead to AUC = 1 since FPR(λ ) < 1 for every λ belonging to Λ, which indicates that the even smallest value of λ will not select all the irrelevant variables. These phenomena are even more visible for n = 2000 in Figure 3 for f 2 . The results for n = 700 in Figure 1 allow us to discriminate a value of K which gives good selection for both σ = 0 and σ = 0.25 and for both functions f 1 and f 2 . Firstly, for f 2 without any noise in the observation set (σ = 0) we can see that only the true variables are selected for K = 1 and K = 2 since AUC = 0 and max(TPR -FPR) = 1. Moreover, these two values of K are the only cases where max(TPR -FPR) = 1 for all the different samplings of the observation set when σ = 0.25. Higher values of K drastically deteriorate the AUC and the max(TPR-FPR), especially for noisy observation sets. Choosing a small value for K offers the advantage of simplifying our model and speeding up computational executions. For all these reasons, we decide to only focus on using our method with K = 1 knot in the B-spline basis. The first one is based on the percentage of variable selection defined for each variable belonging to {1, . . . , p} by:

σ = 0 σ = 0.25 f 1 f 2 K=1 K=2 K=3 K=4 K=5 K=6 K=7 K=8 K=9 K=10 K=1 K=2 K=3 K=4 K=5 K=6 K=7 K=8
P = 100 |Λ| λ ∈Λ 1{ ∈ V λ } , ( 15
)
where |Λ| is the total number of parameters in Λ, 1{A} = 1 if the event A holds and 0 if not and V λ is defined in (12).

Results for the percentage of selection of variables are displayed in Figure 4 (resp. Figure 5) for f 1 (resp. for f 2 ). Firstly, we obtain a high percentage of selection for the relevant variables V f 1 = {3, 5} since they are selected for more than 62.5% of the λ s belonging to Λ. Moreover, we can observe a noticeable gap between the frequency for the relevant (V f 1 = {3, 5}) and the irrelevant (V f 1 = {1, 2, 4}) variables. This gap is amplified as we increase the number of observations from n = 700 to n = 2000. The noise of the observations does not seem here to deteriorate the results. For f 2 , the percentage of the relevant variables (V f 2 = {1, 3, 7, 10}) are lower than the previous function (25% for variable 7 and σ = 0 and 35% for σ = 0.25). However, we can see smaller percentages for the irrelevant variables (V f 2 = {2, 4, 5, 6, 8, 9}) and a clear gap for the unnoisy observation sets (σ = 0) starting with only n = 700 as V f 2 = {2, 4, 5, 6, 8, 9} are never selected. The noise has here an influence on the quality of selection but by increasing the number of observations we can circumvent this issue as the gap is visible for σ = 0.25 with n = 2000. We encourage the user to add known fake variables in order to know which threshold of percentage of selection to use. All the variables having a percentage of selection close to the one of the added fake variables can then be considered as irrelevant as we can see in Figures 4 and5 which are visualizations of the output of our method for one sample of the observation set for f 1 and f 2 , respectively. We also propose another method to automatically choose λ . Among several existing criteria for model selection, we initially assessed a cross-validation criterion using the mean-square error as a loss function. However this approach was not satisfactory as it tended to overestimate the number of relevant variables. Hence, the suggested method leverages the popular Akaike Information Criterion (AIC) introduced in [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF] and defined by:

AIC(λ ) = n ln RSS(λ ) n + 2T λ , (16)
where T λ is the number of terms appearing in (5) by keeping only the variables selected with λ and RSS(λ ) is the residual sum of squares defined as follows:

RSS(λ ) = Y -Y(λ ) 2 2 , with Y(λ ) = p =1 Θ γ (λ ), (17)
where γ (λ ) is defined in (11). Then, the chosen λ = λ AIC is such that:

λ AIC = argmin λ ∈Λ (AIC(λ )). ( 18
)
The TPR(λ ) and FPR(λ ) obtained with λ = λ AIC for both functions f 1 and f 2 are displayed in Figure 8 for n = 700 and n = 2000 and for σ = 0 or σ = 0.25. We can observe for f 1 that even with a small number of observations (n = 700) this criterion allows us to get TPR(λ AIC ) = 1 while having FPR(λ AIC ) = 0 which means that the relevant variables are selected and not the irrelevant ones. The noise in the observation set has a stronger influence on the detection of the relevant variables of f 2 than for f 1 . With σ = 0.25, we indeed have TPR(λ AIC ) < 1 and FPR(λ AIC ) = 0. By increasing the value of n from 700 to 2000, the value of TPR(λ AIC ) is increased and thus the number of relevant selected variables. Moreover, for unnoisy set of observations the relevant variables are recovered from n = 700. Since we are interested in geochemical applications where the noise in the observation sets is very small, we will not be concerned by this issue.

Numerical experiments

In this section, we will assess the robustness of our method called ABSORBER implemented in the absorber R package when the variance of the noise σ 2 increases as well as the number of observations n. We will also study how this novel method behaves when the number of variables p grows. To demonstrate its efficiency, we will compare it to two state-of-the-art methods for feature selection: LassoNet introduced in [START_REF] Lemhadri | Lassonet: A neural network with feature sparsity[END_REF] and the widely used Random Forests (RF) introduced by [START_REF] Breiman | Random forests[END_REF] and used in [START_REF] Genuer | Variable selection using random forests[END_REF] for variable selection.

LassoNet is an open-source package available on GitHub implemented in Python under the name lassonet. The algorithm generates a grid of penalization parameters and the corresponding selected variables. We can thus calculate the percentage of selection for each variable as defined in (15). In the following, the residual neural network is built by taking a hidden layer with 10 neurons as proposed in the notebook example given by this package which focuses on a 26-dimensional regression problem. Despite an increase of the number of neurons in the hidden layer and an extension of the maximal number of epochs, no significant improvement in the selection of relevant variables was observed. This is the reason why we did not explore more complex neural network architectures.

In order to apply Random Forests to our data, we used the R package randomForest, implemented by [START_REF] Liaw | Classification and regression by randomforest[END_REF], with 500 trees. This package provides the percentage of increased mean square error for each variable as the model excludes them one-by-one. This metric is then converted into a percentage of selection for each variable ensuring comparable results for the three methods.

3.1. Influence of n and σ on the quality of variable selection. In the following, we explore the impact of both the number of observations and the noise level in the observation sets on the efficiency of the three previously introduced variable selection methods.

As a result, the percentage of selection for each variable is calculated as described in (15) for 10 samplings of the observation set and the results are shown in Figure 9 (resp. Figure 10) for f = f 1 (resp. f = f 2 ) with n belonging to {700, 1000, 2000} (resp. {700, 1000}). The relevant variables are displayed in green and the irrelevant variables are displayed in red. Additional results for f 2 are presented in the Appendix in Figure 16. As we can see in Figure 9, for a given n, the noise does not seem to have a significant influence on the percentage of variables of f 1 selected by our method ABSORBER. However, as we increase the value of n of the corrupted observation set with σ = 0.25, the percentage of irrelevant variables selected with our method drops from nearly 37.5% for n = 700 to 25% with n = 2000. As observed in these figures, LassoNet tends to select irrelevant variables since the variable selection percentage of the variables belonging to V f 1 is close to that of variable 5 belonging to V f 1 , with a high selection rate of 62.5%. It has to be noticed that Random Forests selects variable 5 with only 37.5% of selection whereas our method (resp. LassoNet) selects it with a percentage of 65% (resp. 73%).

Let us now study the application of these methods to f 2 . The corresponding results are displayed in Figure 10 and in Figure 16 of the Appendix. The noise has an effect on our method ABSORBER in this case since there is no selection of irrelevant variables regardless of the value of n with unnoisy observation sets against 18.5% of selection when σ = 0.25 for the smallest values of n. However, increasing the value of n in this case allows us to reduce the percentage of selection for irrelevant variables to 10% while maintaining the minimum percentage of relevant variables to 32%. Here, using a threshold at 25% allows us to discriminate the relevant variables from the irrelevant ones.

The two other methods appear to be unaffected by changes in both σ and n. Nevertheless, as observed previously with f = f 1 , 50% of the penalization parameters of LassoNet select irrelevant variables. As a consequence, there is no distinct gap between these and relevant variable 7 since its percentage is very close to 50% as well. This statement suggests that using a high threshold on the percentage of selection obtained with LassoNet can result in omitting a relevant variable even for large n. Conversely, a low threshold includes irrelevant variables. In opposition to these two methods, the Random Forests approach tends to fail in detecting the relevant variables as variables 7 and 10 are selected nearly 0% and 5% of the time, respectively, regardless of σ and n. The same conclusion as with LassoNet can be drawn here, emphasizing that our method ABSORBER outperforms those two methods for variable selection while requiring only a few parameters to choose. 3.2. Influence of p on the quality of variable selection. In this section, we seek at studying the behavior of our method when the total number of variables p increases. To do so, we define two additional functions f 3 and f 4 such that: 1) , . . . , x (5) = 1.8 cos (x (1) ) sin (x (3) + 1) -5 ln x (3) + 1 -0.9 (x (4) ) 2 + 1

f 3 x (
x (1) , . . . , x (5) ∈ [0, 1] 5 , (19)
f 4 x (1) , . . . , x (25) = 1.8 cos (x (1) ) sin (x (7) + 1) -5 ln (x (3) + 1) -0.9 1) , . . . , x (25) ∈ [0, 1] 25 .

x (10) 2 + 1 , x (
(20)

We apply all three variable selection methods to f 2 , f 3 and f 4 with observation sets of varying sizes n, all corrupted with the same noise levels as assessed in the previous section. Next, we compute the AUC and max(TPR-FPR) as defined in Section 2.3. The results for these comparisons are displayed in Figure 11 for n = 700 and n = 2000.

We can see from this figure that our method is affected by p but only when the observation set is corrupted with significant noise (σ = 0.25) and has a reduced size (n = 700). Specifically, the value of max(TPR -FPR) equals 1 for p ≤ 10 and is less than 1 for p = 25. In contrast, the efficiency of the two other methods is impacted by the value of p even for σ = 0 as the metrics presented in Figure 11A are deteriorated for p ≥ 10. For instance, with LassoNet, max(TPR -FPR) = 1 and AUC = 1 with p = 5, regardless of σ. However, max(TPR -FPR) becomes less that 1 for p = 10 and it keeps decreasing as the values of p and σ increase. Even in the easiest case where p = 5 and σ = 0, Random Forests is less competitive than the other approaches and its performance remain unchanged regardless the values of n, σ and p. For both of these two methods and for σ = 0, the AUC being close to 1 indicates a high FPR, suggesting that these methods select multiple irrelevant variables while omitting one or more relevant ones. Increasing the number of observations up to n = 2000 improves the results in Figure 11B for our method ABSORBER and Lassonet. Nevertheless, Random Forests continue to display less satisfactory results. Our method is the only one showing satisfactory results as max(TPR -FPR) = 1, regardless of p and σ. Globally, our method outperforms LassoNet and RF in this case and appears to be more robust when facing with higher dimensions p with or without noisy observation sets.

In Figure 12, the impact of the value p on our variable selection procedure using AIC is assessed. We can observe that p has little effect on the efficiency of our method, in contrast to σ which has a more pronounced negative impact on the metrics. Specifically, the TPR values are affected by σ as it necessitates an increasing value of n as σ grows to achieve TPR(λ AIC ) = 1. However, for smaller noise levels (σ < 0.25) and with n ≥ 700 our method enables TPR(λ AIC ) = 1 while maintaining FPR(λ AIC ) = 0. This means that no irrelevant variables are chosen. Since the number of relevant variables is d = 4 even for p = 5, the selection of just one irrelevant variable leads to FPR(λ AIC ) = 0.25, thereby explaining the observed results for n = 2000. However, such trends are not visible for high value of p, demonstrating the efficiency of our variable selection procedure.

3.3. Numerical performance. The goal of this section is to investigate the computational times of our variable selection approach implemented in the absorber R package. Our variable selection method is applied to f 2 , f 3 and f 4 (defined in ( 14), ( 19) and (20), respectively) for p = 5, p = 10 and p = 25, respectively. The timings were obtained on a workstation with 31.2GB of RAM and Intel Core i7 (3.8GHz) CPU. The log 10 -transformed average computational times and their standard deviation obtained from 20 independent executions are displayed in Figure 13. We can see from this figure that it only takes 57 seconds to perform variable selection on a function with p = 25 variables from n = 2000 observations. It has to be noticed that the execution times reported are mainly due to the use of the sparsegl package.

Application to a multidimensional geochemical system

In this section, we apply our variable selection method to real geochemical systems. We start by defining the geochemical system derived from the calcite dissolution and precipitation study in [START_REF] Kolditz | Thermo-hydro-mechanical-chemical processes in porous media: benchmarks and examples[END_REF]. In the following, the observation sets are generated using the geochemical solver PHREEQC as in [START_REF] Parkhurst | Description of input and examples for PHREEQC version 3: a computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations[END_REF]. The corresponding thermodynamic data for aqueous species and minerals are available in the Phreeqc.dat distributed the purposes of this paper, we specifically focus on the calcite precipitation/dissolution: Calcite Ca 2+ + CO 2- 3 , logK = -8.48. The previous equation represents the dissolution reaction along with its corresponding log 10transformed equilibrium constant value logK. The amount of calcite can be here computed with PHREEQC as a function of the total elemental concentrations (C, Ca), the pH (as -log(H + )) and the amount of calcite initially present. The pH is here fixed at 9.8 and we introduce two additional concentrations (K, Cl) which do not participate in the calcite precipitation. This allows us to assess our variable selection on real datasets. We can thus formulate the problem as:

Calcite = f 5 (C * , Ca * , K * , Cl * , Calcite * ) = f 5 (C * , Ca * , Calcite * ),
where C * , Ca * , K * , Cl * , Calcite * are the normalized concentrations and quantities initially present of C, Ca, K, Cl and Calcite, respectively. The normalization of each variable is done by taking into account the minimal and the maximal bound of the values so that each variable belongs to [0, 1]. We also define another function f 6 which takes into account known fake variables to study the behavior of all three methods: Calcite = f 6 (C * , Ca * , K * , Cl * , Calcite * , x (6) , x (7) , x (8) , x (9) , x (10) )

= f 5 (C * , Ca * , Calcite * ), where x (6) , x (7) , x (8) , x (9) , x (10) are synthetic irrelevant variables obtained through uniform sampling between 0 and 1. Hereafter, V f 5 = V f 6 = {C * , Ca * , Calcite * }, V f 5 = {K * , Cl * } and V f 6 = {K * , Cl * , x (6) , x (7) , x (8) , x (9) , x (10) }.

The results for the application of our method, LassoNet and RF to f 5 and f 6 are displayed in Figure 14. Here, our method consistently selects the relevant variables belonging to V f 5 65% of the time and the irrelevant ones belonging to V f 5 are almost never selected. By increasing n, the percentage of selection for these irrelevant variables reaches 0%. In contrast, LassoNet selects all the variables and fails to discriminate the relevant from the irrelevant ones. Random Forests, on the other hand, tends to detect only one variable of V f 5 which is the calcite quantity. Furthermore, it tends to select K and Cl (5%) more than the relevant variables C and Ca (0%), even with n = 2000.

Adding fake variables does not deteriorate our method which continues to select only the relevant variables in V f 6 . However, it does not improve the results for LassoNet which continues to select all ten variables from V f 6 and V f 6 . Interestingly, it helps the Random Forests approach in detecting the relevant variables C and Ca resulting in an increase in the percentage of selection up to 5%. Nevertheless, this emphasizes the efficiency of our method which outperforms the other two in this geochemical case.

Furthermore, we used the AIC to select the parameter λ and to automatically choose the relevant variables. The corresponding results are shown in Figure 15. This statistical criterion proves to be highly efficient as evidenced by TPR(λ AIC ) = 1 and FPR(λ AIC ) = 0, regardless of n and p.

Figure 1 .Figure 2 .Figure 3 .

 123 Figure 1. AUC and max(TPR -FPR) calculated for an increasing value of K for f 1 (top) and f 2 (bottom) with noise (right) or without noise (left) in the observation set Y = (Y 1 , . . . , Y 700 ). 10 random samplings of Y were used to obtain these results. The empty bullets inside the boxplots represent the mean value and the plain bullets outside the boxplots are the extreme values.

Figure 4 .Figure 5 .Figure 6 .Figure 7 .

 4567 Figure 4. Percentage of selection of each variable for f 1 with n = 700 (left) and n = 2000 (right) and for σ = 0 or σ = 0.25. The green (resp. red) variables indicate the true relevant (resp. irrelevant) variables for f 1 . 10 random samplings of Y were used to obtain these results. The empty bullets inside the boxplots correspond to the mean value and the plain bullets outside the boxplots are the extreme values.

Figure 8 .

 8 Figure 8. TPR(λ ) and FPR(λ ) values by choosing λ = λ AIC for f 1 (top) and f 2 (bottom) with an unnoisy (σ = 0) or noisy (σ = 0.25) set of observations. 10 random samplings of Y were used to obtain these results. The empty diamonds inside the boxplots correspond to the mean value and the plain bullets outside the boxplots are the extreme values.

Figure 9 .Figure 10 .

 910 Figure 9. Percentage of selection of each variable of f 1 with three different methods: ABSORBER, LassoNet and RandomForests (RF) with an increasing number of observations n (left to right) and of the value of σ (top to bottom). 10 random samplings of Y were used to obtain these results. The empty diamonds inside the boxplots correspond to the mean value and the plain bullets outside the boxplots are the extreme values.

Figure 11 .Figure 12 .Figure 13 .

 111213 Figure 11. AUC and max(TPR-FPR) calculated for three different variable selection methods: ABSORBER, LassoNet and RandomForests (RF) applied to three functions (top to bottom) f 3 , f 2 and f 4 with n = 700 (A) and n = 2000 (B) and an increasing value of σ (left to right). 10 random samplings of Y were used to obtain these results. The empty diamonds inside the boxplots correspond to the mean value and the plain bullets outside the boxplots are the extreme values.

Figure 14 .

 14 Figure 14. Percentage of selection of each variable of f 5 (top) and f 6 (bottom) with three different methods: ABSORBER, LassoNet and RandomForests (RF) with an increasing number of observations n (left to right). 10 random samplings of Y were used to obtain these results. The empty diamonds inside the boxplots correspond to the mean value and the plain bullets outside the boxplots are the extreme values.

Figure 15 .Figure 16 .

 1516 Figure 15. TPR(λ ) and FPR(λ ) values by choosing λ = λ AIC for f 5 (left) and f 6 (right). 10 random samplings of Y were used to obtain these results. The empty diamonds inside the boxplots correspond to the mean value and the plain bullets outside the boxplots are the extreme values.

  x

	( ) min	, t ,1 , . . . , t ,K	, x ( ) max , . . . , x ( ) max	,
	M times	t	M times
	where x Denoting by B ( ) min and x ( ) ( ) max are the lower and upper bounds of S , respectively.


Clément Delmas
email: delmas.clement@chu-toulouse.fr

Laurent Bonello

François Roubille

For the best management, please ask for assistance!

may come

For the best management, please ask for assistance!

Cardiogenic shock (CS) is defined as a profound impairment of cardiac function, leading to severe and life-threatening systemic organ malperfusion. Long considered to be largely secondary to acute myocardial infarction (AMI), it is now recognized as being of nonischaemic (NI) origin in more than one out of twocases.1Over the past years, despite extensive research efforts, the CS-associated mortality risk has remained unacceptably high and even probably under-evaluated.2 A recent multidisciplinary international impetus has enlightened numerous observational and then interventional studies in CS, however the vast majority of these studies have focused on AMI-CS and neglected shocks occurring in other contexts. Future and ongoing studies also focus on AMI-CS (Table1) with only few allowing the inclusion of NI-CS (as exclusive or possible aetiologies). Indeed, according to a recent analysis of 1305 unselected CS referred in a tertiary centre, only 31.8% met the enrolment criteria for any of the largest current and past CS trials (IABP-SHOCK II, ECLS-SHOCK, EURO-SHOCK or DanGer-SHOCK) with between studies differences (eligibility from11.9% to 26.9%). The NI cause of the CS explained 50% to 70% reason for ineligibility questioning the representativity and potential translation in clinical practice regarding NI-CS.3

Importantly, ischaemic and NI-CS differ substantially in terms of both pathophysiology including neurohormonal activation and clinical features from haemodynamic initial presentation until baseline characteristics such as age, gender, comorbidities or history. Especially in case of preexisting acutely decompensated chronic heart failure (ADCHF), NI-CS usually have a substrate of chronic compensatory mechanisms including ventricular dilatation, increase in circulating catecholamines and upregulation of the renin-angiotensinaldosterone axis.4 They present then with less profound haemodynamic instability (stage B and C Society for Cardiovascular Angiography and Interventions [SCAI]), higher systemic overload, lower left ventricular ejection fraction and cardiac index, and more renal and liver dysfunction. By contrast, AMI-CS are characterized by an acute loss of myocardial tissue and ventricular stunning with early and sometimes profound systemic inflammatory response syndrome (SIRS) due to the lack of compensatory mechanism with systemic vascular resistance rapidly crumbled which participate in the complexity of its management.5However, these data should be tempered because they differentiate well AMI-CS from ADCHF-CS but probably not from acute de novo NI-CS (fulminant myocarditis or toxics for example).

To improve patients' diagnosis, management and outcomes, studies dedicated to NI-CS are urgently needed. Trials on CS are difficult to conduct due to consent in timely acute severe condition but also due to heterogeneity of patients and practice.2To date only the ECMO-CS randomized trial reported the place of extracorporeal membrane oxygenation (ECMO) in miscellaneous CS but results were non-significant probably limited by its underpowered size (117 patients) and the large cross-over in the no early ECMO group (39% received an ECMO).6Moreover, the EURO-SHOCK trial designed to include ischaemic and NI-CS patients randomized between ECMO and medical management was recently prematurely stopped for slow inclusion rate (only 22 inclusions on 428 planned).7

In this issue of the Journal, Schrage et al.8used a large, retrospective, international, multicentre database to evaluate mechanical circulatory support (MCS) use in NI-CS. This is the first large-scale analysis of the relationship between MCS use and mortality in these severe patients. First, beyond the largest series on NI-CS published to date (n=890), they showed that MCS was used in 43% of the cases despite absence of evidence for this indication. Second, based on a propensity-score matching analysis to consider baseline differences between MCS and non-MCS patients, they showed a significant 30-day mortality reduction in favour of the MCS use in all tested subgroups (hazard ratio 0.76, 95% confidence interval0.59-0.97) at the price of higher early complications (e.g. severe bleeding [16.5% vs. 6.4%] and access site-related ischaemia [6.7%vs. 0%]).

Nevertheless, several limitations should be highlighted and dis-cussed. First, the definition of NI-CS as a CS due to severe de novo or ADCHF with no need for urgent revascularization, exclude AMI-CS but without precision on the CS cause or trigger, and no comparison between severe de novo and ADCHF-CS was reported preventing further analysis of this heterogeneous population. Indeed, to differentiate patients and management according to the probability of heart recovery seems legit and important from high probability (myocarditis, de novo dilated cardiomyopathy) versus low probability (previous advanced heart failure). Beyond myocardial recovery, the project is important to consider the possibility of heart replacement therapies such as heart transplantation or ventricular assist device (VAD). Second, they reported only early outcomes at 30 days and did not differentiate the type of support used (ECMO vs percutaneous VAD [pVAD]), although we need to consider also the risk of heart failure with the paradox of 'saving the patient but no the heart'. Long-term data in this severe population are mandatory even more with the change of several heart allocation systems (US, France) which prioritizes access for patients under acute MCS. Notably, the prioritization of patients under ECMO is questionable based on a large analysis on US heart transplantation showing a 2.4 increase of1-and 5-year mortalities for patients under ECMO at transplantation compared with patients without support or with other support (pVAD or left VAD).9 Third, as acknowledged by the authors many biases, including the retrospective design, prevent from strong conclusions. The size of the effect is suggested important (until 50% less events)but the confidence interval is wide, near to statistical significance. So, we cannot but support prospective larger trials on this topic, as proposed by the authors, in spite of the difficulties mentioned above.

Cardiogenic shock is a complex disease and haemodynamics is just one problem cascade which ultimately leads to death. Indeed, the simplistic view of CS as a unique problem of output is definitely false. Associated SIRS notably may not be reversed by device alone. So, the possible benefit of the device in terms of output may be counterbalanced by their associated adverse effects (haemolysis, bleedings, vascular complications). Patient selection and timing of implantation seem crucial when dealing with CS patients' outcomes since it can be estimated that up to 40-50%will die in absence of specific management and MCS may be helpful for them. But more 50% to 60% will survive without MCS and may be exposed to potential undeserved side effects from bleedings or acute limb ischaemia to death.10 More, through the overall survival, myocardial recovery and survival without heart replacement therapy remains crucial. This benefit/risk ratio depends on multiple factors from the type of device used, the timing of implantation, associated management, to centre experience and expertise.

The support needed and the type of device may probably be different between AMI-CS and NI-CS. For the latter, output need is probably less important thanks to chronic tissue adaptation to low perfusion but also to the lower associated SIRS. So, less efficient and less invasive MCS may be used for NI-CS support and intra-aortic balloon pump (IABP) with less associated adverse events (<3% vs. up to 40-50% for pVAD and ECMO) could be a good candidate. In a recent small randomized controlled trial (n=32) IABP was associated with a significant improvement in cardiac power output, decongestion and organ perfusion(mixed venous oxygen saturation) without serious adverse events, compared to inotropes (enoximone or dobutamine) in ADCHFpatients.11

Although in AMI-CS the earlier implantation seems the better with pVAD or extracorporeal life support in observational studies, Schrage et al. showed only a significant improvement in 30-daymortality for NI-CS stage D when the SCAI classification was considered. They confirmed that too late is probably futile12but also that the benefit/risk ratio is not in favour of too early implantation (stage C or earlier) in NI-CS. First step associating aetiological management, inotropes and vasopressors should stay mandatory in this specific population.

Based on (i) pre-clinical studies showing an increase in left ventricular pressure and alteration in coronary perfusion on ECMO, and (ii) the association of left overload with higher mortality or need for heart replacement therapies in patients with left ventricular predominant failure,13the concept of left ventricular unloading has emerged. Recently, the association between early(until 2 h after ECMO) left ventricular active unloading by pVAD was reported to be associated with a lower 30-day mortality in an adjusted and propensity-matched population of NI-CS. These results question the best type of device to choose and its timing, since this strategy was associated with a higher complication rate.14

Perspectives

Future trials dedicated to NI-CS are needed to address unmet goals of NI-CS management and improve its outcomes. To this purpose, the ongoing Altshock-2 study compares early IABP versus medical management in NI-CS.15Beyond the device itself, specific dedicated algorithms to optimize medical management and pathways are urgently needed through network organization with multidisciplinary CS centres. The ideal study should be multicentre, international, with four groups and ban cross-over. It should com-pare medical management alone versus its association with IABP, or pVAD or ECMO (with and without systematic left ventricular unloading) with composite primary endpoint of survival free from heart replacement therapy and a specific medico-economic analysis, but its feasibility seems very unlikely.

In a large international registry, Schrage et al. showed for the first time a significant association between MCS use and lower early mortality in NI-CS. Despite attractive results, this study is only hypothesis-generating and future dedicated studies are necessary to address unanswered questions in a prospective way.

Conflict of interest:

none declared.