Michael Levine
email: mlevins@purdue.edu

Gildas Mazo
email: gildas.mazo@inrae.fr

A smoothed semiparametric likelihood for estimation of nonparametric finite mixture models with a copula-based dependence structure

Keywords: Nonparametric finite density mixture, copula, pseudo-EM algorithm 1

In this manuscript, we consider a finite multivariate nonparametric mixture model where the dependence between the marginal densities is modeled using the copula device. Pseudo EM (Expectation-Maximization) stochastic algorithms were recently proposed to estimate all of the components of this model under a location-scale constraint on the marginals. Here, we introduce a deterministic algorithm that seeks to maximize a smoothed semiparametric likelihood. No location-scale assumption is made about the marginals. The algorithm is monotonic in one special case, and, in another, leads to "approximate monotonicity"-whereby the difference between successive values of the objective function becomes non-negative up to an additive term that becomes negligible after a sufficiently large number of iterations. The behavior of this algorithm is illustrated on several simulated and real datasets. The results suggest that, under suitable conditions, the proposed algorithm may indeed be monotonic in general. A discussion of the results and some possible future research directions round out our presentation.

1 Introduction Let g(x) = g(x 1 , . . . , x d ) = K k=1 π k f k (x 1 , . . . , x d ) (1)
be a multivariate mixture model with K components (or clusters-we shall use these two words interchangeably). We view the model (1) as a nonparametric mixture model where individual components f k are not defined as belonging to any specific parametric family. The research on selecting the number of components for non-and semiparametric density mixtures is currently at a very early stage; some developments in this area can be found in e.g. [START_REF] Kasahara | Non-parametric identification and estimation of the number of components in multivariate mixtures[END_REF] and [START_REF] Kwon | Estimation of the number of components of nonparametric multivariate finite mixture models[END_REF]. Due to this, we assume that the number of components K is fixed and known in our model. In general, most of the work on nonparametric mixture modeling so far assumed that the marginal distributions f k1 , . . . , f kd of each component are conditionally independent. Such an assumption implies that, conditional on knowing which component a particular observation has been generated from, its distribution is equal to the product of its marginals. More formally, this means that

g(x) = K k=1 π k d j=1 f kj (x j ).
This model has been introduced for the first time in [START_REF] Brezis | Nonparametric estimation of component distributions in a multivariate mixture[END_REF]. The conditions sufficient to ensure identifiability for the conditionally independent model are known [START_REF] Allman | Identifiability of parameters in latent structure models with many observed variables[END_REF]. There are also a number of approaches to estimating their parameters [START_REF] Xiang | An Overview of Semiparametric Extensions of Finite Mixture Models[END_REF], both iterative [START_REF] Benaglia | An EM-like algorithm for semi-and nonparametric estimation in multivariate mixtures[END_REF]; [START_REF] Levine | Maximum smoothed likelihood for multivariate mixtures[END_REF] and closed form solutions [START_REF] Bonhomme | Non-parametric estimation of finite mixtures from repeated measurements[END_REF]. However, the assumption of conditional independence is not always a realistic one. For example, it is unlikely to be true when dealing with RNA-seq data [START_REF] Rau | Co-expression analysis of high-throughput transcriptome sequencing data with Poisson mixture models[END_REF]. Thus, it seems desirable to relax this assumption while retaining the generality of the nonparametric approach.

To the best of our knowledge, the only known results on estimation of nonparametric mixture models with conditionally non-independent components are [START_REF] Mazo | A semiparametric and location-shift copula-based mixture model[END_REF]; [START_REF] Mazo | Constraining kernel estimators in semiparametric copula mixture models[END_REF]. A somewhat related model was also considered in Vrac et al (2012). There, however, the authors model not the distribution of observations, but rather the distribution of a number of cumulative distribution curves assumed to represent observations. Thus, they assume that, at any moment in time, one can observe an entire cumulative distribution curve F (x) for any of the K possible distributions comprising the overall mixture. This is not at all possible in our setting. [START_REF] Mazo | A semiparametric and location-shift copula-based mixture model[END_REF] and [START_REF] Mazo | Constraining kernel estimators in semiparametric copula mixture models[END_REF] consider a special case of the general nonparametric mixture model, allowing for a non-trivial dependence structure where the marginals are assumed to belong to a location-scale family. Stochastic algorithms were proposed to estimate the copula parameter and the nonparametric marginals.

The estimation algorithms, while performing well in practice, do not optimize any particular objective function. Because of this, their convergence analysis will necessarily be a difficult one. In this manuscript, our goal is to suggest a deterministic algorithm capable of estimating the components of a nonparametric mixture model with conditionally non-independent components without a location-scale assumption for the marginals, since such an assumption is far from commonly satisfied in applications.

In order to continue, we are going to fix the notation first. It is well-known that, due to Sklar's theorem [START_REF] Nelsen | An introduction to copulas[END_REF] p. 18, every d-dimensional multivariate cumulative distribution function can be represented as a copula of the corresponding marginal cumulative distribution functions. Indeed, let F k1 (x 1 ), . . . , F kd (x d ) be the marginal cumulative distribution functions of the cumulative distribution function

F k (x 1 , . . . , x d ) that corresponds to the density f k (x 1 , . . . , x d ). Then, there exists a d- copula C k , which is a function C k : [0, 1] d → [0, 1], such that F k (x 1 , . . . , x d ) = C k (F k1 (x 1 ), . . . , F kd (x d )),
see [START_REF] Nelsen | An introduction to copulas[END_REF] pp. 46. If the marginal cumulative distribution functions are continuous, then the copula is unique. The copula C k can be viewed as a d-dimensional cumulative distribution function with uniform marginal distributions. Taking the derivative of order d, one immediately obtains the representation

f k (x 1 , . . . , x d ) = c k (F k1 (x 1 ), . . . , F kd (x d )) d j=1 f kj (x j )
where c k is the density of the copula C k . We assume that each copula density c k belongs to some parametric family of copula densities indexed by a parameter θ k .

Due to this, from now on we will use the index k as a subscript for θ k only but will drop this subscript for c k . Denoting by ϕ the set of all marginal densities {f kj }, and denoting by π = (π 1 , . . . , π K ) and θ = (θ 1 , . . . , θ K ) the vectors of all weights and copula parameters, respectively, we have

f k (x; θ, ϕ) = f k (x 1 , . . . , x d ; θ, ϕ) = c(F k1 (x 1 ), . . . , F kd (x d ); θ k ) d j=1 f kj (x j ), (2)
so that (1) and (2) define a class of mixture densities that can be stated as g(•; π, θ, ϕ).

To the best of our knowledge, no identifiability results are available concerning this model.

The rest of this manuscript is structured as follows. Section 2 introduces a general algorithm that can be used to estimate finite mixtures of multivariate densities with a dependence structure defined through the use of copulas. Section 3 provides some results about the monotonicity property of two simplified versions of this algorithm.

Section 4 details the implementation of the algorithm. Section 5 analyses the performance of our algorithm with several simulation studies. Section 6 presents applications to two real datasets. Finally, the conclusion section discusses the results obtained and suggests possible directions for future research.

Algorithm

The goal of our manuscript is to estimate the components and weights of the model ( 1)-( 2). The definition of such an algorithm starts with an objective function that we are going to introduce next. First, let K(•) be a proper univariate density function that can be used for kernel density estimation and K h (•) := 1 h K • h its rescaled version where h > 0 is a bandwidth. Next, for a generic function f, we define

N h f (x) := exp K h (x -u) log f (u) du (3)
which is a nonlinear smoother of the function f. Note that, even if f is a density, N f is not, in general a density due to Jensen's inequality. Now, we define the operator The objective function we seek to maximize is the population version of the smoothed semiparametric log-likelihood, given by

O by Of k (x; θ, ϕ) = c(F k1 (x 1 ), . . . , F kd (x d ); θ k ) d j=1 N f kj (x j ).
(π, θ, ϕ) = g(x) log ǧ(x; π, θ, ϕ) g(x) dx, (4)
over all (π, θ, ϕ); here g(x) is the target density. If the marginal distributions are conditionally independent then c(u 1 , . . . , u d ; θ k ) ≡ 1 for every θ k and k, and hence (4) reduces to the smoothed semiparametric log-likelihood considered in [START_REF] Levine | Maximum smoothed likelihood for multivariate mixtures[END_REF].

Lemma 1. For any choice of parameters π, θ, ϕ, the smoothed loglikelihood difference is bounded as

(π, θ, ϕ) -( π, θ, ϕ) ≤ K k=1 -log π k π k g(x)w k (x; π, θ, ϕ)dx -g(x) K k=1 w k (x; π, θ, ϕ) log d j=1 N f kj (x j ) d j=1 N f kj (x j ) dx -g(x) K k=1 w k (x; π, θ, ϕ) log c( F k1 (x 1 ), . . . , F kd (x d ); θ k ) c(F k1 (x 1 ), . . . , F kd (x d ); θ k ) dx :=Ψ 1 ( π|π, θ, ϕ) + Ψ 2 ( ϕ|π, θ, ϕ) + Ψ 3 ( θ, ϕ|π, θ, ϕ),
where the cumulative distribution functions F kj are those associated with { f kj } = ϕ and

w k (x; π, θ, ϕ) = π k Of k (x; θ, ϕ)/ǧ(x; π, θ, ϕ), (5)
k = 1, . . . , K.

Proof of Lemma 1.. By definition, the difference of smoothed log-likelihoods can be written down as

(π, θ, ϕ) -( π, θ, ϕ) = -g(x) log K k=1 π k Of k (x; θ, ϕ) K k=1 π k Of k (x; θ, ϕ) dx = -g(x) log K k=1 w k (x; π, θ, ϕ) π k Of k (x; θ, ϕ) π k Of k (x; θ, ϕ) dx
At this point, it remains only to apply Jensen's inequality to a convex combination on the right-hand side whose coefficients are w k (x; θ, ϕ).

Instead of minimizing (π, θ, ϕ) -( π, θ, ϕ) with respect to ( π, θ, ϕ) directly, we seek to minimize the upper bound proposed by Lemma 1. This approach is in the spirit of MM (Minimization-Majorization) algorithms; see e.g. [START_REF] Wu | The MM alternative to EM[END_REF] for the detailed discussion. To do this, our heuristic is to minimize each of the three terms Ψ 1 ( π|π, θ, ϕ), Ψ 2 ( ϕ|π, θ, ϕ), Ψ 3 ( θ, ϕ|π, θ, ϕ) separately. This is sometimes called "minimization by part". To minimize the first term Ψ 1 ( π|π, θ, ϕ), we have to choose π = π where π k = g(x)w k (x; π, θ, ϕ) dx, k = 1, . . . , K. This is the result that can be obtained using standard constrained optimization techniques. Note that the resulting minimum must be non-positive since the first term can be made zero by choosing π = π. To minimize the second term Ψ 2 ( ϕ|π, θ, ϕ), define, as a first step,

f kj (u j ) = α kj g(x)w k (x; π, θ, ϕ)K h kj (x j -u j ) dx,
for any k = 1, . . . , K and j = 1, . . . , d, where α kj is the normalizing constant ensuring that the newly defined f kj is, indeed, a proper density function. Then, we have

-g(x)w k (x; π, θ, ϕ) log N f kj (x j ) dx = -g(x)w k (x; π, θ, ϕ) K h kj (x j -u j ) log f kj (u j ) du j dx = -log f kj (u j ) f kj (u j ) du j .
The same argument as in [START_REF] Levine | Maximum smoothed likelihood for multivariate mixtures[END_REF] applies: the quantity above is minimized if we select f kj (u) = f kj (u). The resulting minimum will also be less than or equal to zero because Ψ 2 ( ϕ|π, θ, ϕ) = 0 when ϕ = ϕ.

Now, we can propose the following general algorithm for estimation of (π, θ, ϕ).

A1 Choose initial values π 0 , ϕ 0 , θ 0 A2 Compute the initial set of weights

w k (x; π 0 , θ 0 , ϕ 0 ) = π 0 k Of k (x; θ 0 , ϕ 0 )/ǧ(x; π 0 , θ 0 , ϕ 0 ).
A3 At any step of iteration t = 1, 2, . . . select

π t k = g(x)w k (x; π t-1 , θ t-1 , ϕ t-1 ) dx, k = 1, . . . , K.
A4 Select as the next value of the density function vector ϕ t = {f t kj } where

f t kj (u j ) = α kj g(x)w k (x; π t-1 , θ t-1 , ϕ t-1 )K h kj (x j -u j ) dx
where α kj is the normalizing constant ensuring that the newly defined function is, indeed, a density function. As a part of this step, also compute updated cumulative

distribution functions F t kj (u j ) = uj -∞ f t kj (y) dy.
A5 Choose the value

θ t = arg min θ Ψ 3 (θ, ϕ t |π t-1 , θ t-1 , ϕ t-1 ).
A6 Redefine weights

w k (x; π t , θ t , ϕ t ) = π t k Of k (x; θ t , ϕ t )/ǧ(x; π t , θ t , ϕ t ).
and return to step A3.

At each step of the algorithm defined above, the marginals are updated first and independently of the copula parameter. This strategy was used in [START_REF] Mazo | A semiparametric and location-shift copula-based mixture model[END_REF]; [START_REF] Mazo | Constraining kernel estimators in semiparametric copula mixture models[END_REF].

Remark 1. In practice, one implements the empirical version of the algorithm. Every integral of the form g(x)ζ(x) dx, where ζ is some arbitrary function, is replaced by

1 n n i=1 ζ(X i ), where X i = (X i1 , . . . , X id ), i = 1, . . . ,

Studying the algorithm

Whether the algorithm proposed in Section 2 is monotonic with respect to the objective functional ( 4) is an open question. In some special cases, the answer is positive. One such case that we identified is when probabilities π k and the marginal densities f kj are known beforehand. In such a case, the simplified algorithm is as follows.

B1 Choose initial value of the copula parameter θ 0 .

B2 Compute the initial set of weights

w k (x; π, θ 0 , ϕ) = π k Of k (x; θ 0 , ϕ)/ǧ(x; π, θ 0 , ϕ).
B3 For any t = 1, 2, . . . choose the value

θ t = arg min θ Ψ 3 (θ, ϕ|π, θ t-1 , ϕ).
B4 Redefine weights

w k (x; π, θ t , ϕ) = π k Of k (x; θ t , ϕ)/ǧ(x; π, θ t , ϕ).
and return to step B3.

Proposition 1. The algorithm defined in B1-B4 is monotonic with respect to θ, that is, (π, θ t-1 , ϕ) -(π, θ t , ϕ) ≤ 0 for every t = 1, 2, . . . Proof.
The smoothed likelihood difference is bounded from above as

(π, θ, ϕ) -(π, θ, ϕ) ≤ Ψ 3 ( θ, ϕ|π, θ, ϕ) = -g(x) K k=1 w k (x; π, θ, ϕ) log c(F k1 (x 1 ), . . . , F kd (x d ); θ k ) c(F k1 (x 1 ), . . . , F kd (x d ); θ k ) dx. Choosing θ * = arg min θ Ψ 3 ( θ, ϕ|π, θ, ϕ) produces (π, θ, ϕ) -(π, θ * , ϕ) ≤ Ψ 3 (θ * , ϕ|π, θ, ϕ) = min θ Ψ 3 ( θ, ϕ|π, θ, ϕ);
since there exists a value θ = θ such that Ψ 3 (θ, ϕ|π, θ, ϕ) ≡ 0, the minimal value of Ψ 3 ( θ, ϕ|π, θ, ϕ) will be less than or equal to zero.

Another interesting special case results when one assumes that both component weights π k and copula parameters θ k are known while the marginal densities f kj are unknown. In this case, the simplified algorithm will be as follows.

C1 Choose initial values ϕ 0 C2 Compute the initial set of weights

w k (x; π, θ, ϕ 0 ) = π k Of k (x; θ, ϕ 0 )/ǧ(x; π, θ, ϕ 0 ).
C3 For t = 1, 2, . . . select as the next value of the density function vector

ϕ t = {f t kj } where f t kj (u j ) = α kj g(x)w k (x; π, θ, ϕ t-1 )K h kj (x j -u j ) dx.
Here, α kj is a normalizing constant, ensuring that the newly defined function is, indeed, a density function. As a part of this step, also compute updated cumulative distribution

functions F t kj (u j ) = uj -∞ f t kj (y) dy.
C4 Redefine weights

w k (x; π, θ, ϕ t ) = π k Of k (x; θ, ϕ t )/ǧ(x; π, θ, ϕ t ).
and return to step C3.

The special case of the general algorithm defined above possesses an "approximate monotonicity" property in the following sense.

Proposition 2. We assume that the target density g(x) has a compact support Ω.

We also assume that none of the known weights π k is equal to zero. Suppose that the kernel function K(•) is a proper density function defined on [-1, 1], bounded away from zero by K * > 0, and Lipschitz continuous with a positive Lipschitz constant L. We assume that the copula density function c(u 1 , . . . , u d ; θ) is also Lipschitz continuous on [0, 1] d and bounded away from zero. Then, there exists a subsequence

ϕ t l = (f t l kj , k = 1, . . . , K, j = 1, . . . , d), l = 1, 2, .
. . , such that the the algorithm C1-C4

is "approximately monotonically ascending" along this subsequence:

(π, θ, ϕ t l-1 ) -(π, θ, ϕ t l ) ≤ o(1)
as l → ∞.

Remark 2. It follows directly from the definition that K * ≤ K(•) ≤ K * where both K * and K * are positive. The assumptions of Lipschitz continuity and boundedness away from zero for the kernel function K(•) do not represent a practical problem since they are not concerned with the actual data-rather, K(•) is a tool used to analyze the data. Our simulation results suggest that they also may not be necessary.

Remark 3. The assumption of compact support for the target density g(x) and, by extension, for all of the marginal densities f kj does not represent a problem from the practical viewpoint. From the theoretical viewpoint, a result analogous to Proposition 2 can be proved if one assumes that all of the marginal densities decay to zero sufficiently fast at infinity and using the Fréchet-Kolmogorov theorem instead of the Arzelà-Ascoli theorem Brezis (2011) p. 126.

Remark 4. As an example of copulas satisfying conditions of Proposition 2 we can point out Farlie-Gumbel-Morgenstern (FGM) copulas as well as so-called copulas with cubic sections (that are direct generalizations of FGM copulas) Nelsen ( 2007) pp.

-84.

Proof. The difference in log-likelihoods can be bounded as

(π, θ, ϕ t l-1 ) -(π, θ, ϕ t l ) ≤ Ψ 2 (ϕ t l |π, θ, ϕ t l-1 ) + Ψ 3 (θ, ϕ t l |π, θ, ϕ t l-1 ) = -g(x) K k=1 w k (x; π, θ, ϕ t l-1 ) log d j=1 N f t l kj (x j ) d j=1 N f t l-1 kj (x j ) dx -g(x) K k=1 w k (x; π, θ, ϕ t l-1 ) log c(F t l k1 (x 1 ), . . . , F t l kd (x d ); θ k ) c(F t l-1 k1 (x 1 ), . . . , F t l-1 kd (x d ); θ k ) dx. Recall that minimization of Ψ 2 (ϕ t l |π, θ, ϕ t l-1 ) always results in Ψ 2 (ϕ t l |π, θ, ϕ t l-1 ) ≤ 0 since the choice f t l kj = f t l-1 kj
for all k = 1, . . . , K and j = 1, . . . , d makes this term equal to zero. Therefore, it remains to show that Ψ 3 (θ, ϕ t l |π, θ, ϕ t l-1 ) → 0 as l → ∞.

To do this, let us introduce a lemma.

Lemma 2. For each k = 1, . . . , K and j = 1, . . . , d, the sequence f t kj , t = 1, 2, . . . has a uniformly converging subsequence f t l kj , l = 1, 2, . . . .

The proof of Lemma 2 is similar to the proof of Lemma A2 in [START_REF] Levine | Maximum smoothed likelihood for multivariate mixtures[END_REF] and is not given. Denote by f * kj the limit of f t l kj as l → ∞. Denote by ϕ * the collection of all such limits. Since Ω is compact, it follows in a straightforward manner from Lemma 2 that each subsequence F t l kj (u) converges uniformly to

F * kj (u) := u -∞ f * kj (x) dx. To show that Ψ 3 (θ, ϕ t l |π, θ, ϕ t l-1
) goes to zero as l goes to infinity, we proceed as follows. We have

|Ψ 3 (θ, ϕ t l |π, θ, ϕ t l-1 )| ≤ K k=1 g(x)w k (x; π, θ, ϕ t l-1 ) log c(F t l k1 (x 1 ), . . . , F t l kd (x d ); θ k ) c(F t l-1 k1 (x 1 ), . . . , F t l-1 kd (x d ); θ k ) dx .
Each summand is bounded as

g(x)w k (x; π, θ, ϕ * ) log c(F t l k1 (x 1 ), . . . , F t l kd (x d ); θ k ) c(F t l-1 k1 (x 1 ), . . . , F t l-1 kd (x d ); θ k ) dx + g(x)(w k (x; π, θ, ϕ t l-1 ) -w k (x; π, θ, ϕ t l )) log c(F t l k1 (x 1 ), . . . , F t l kd (x d ); θ k ) c(F t l-1 k1 (x 1 ), . . . , F t l-1 kd (x d ); θ k ) dx . (6)
Since the copula density is bounded from above and below, the second term is less than or equal to a constant times g(x)|w k (x; π, θ, ϕ t l-1 )-w k (x; π, θ, ϕ t l )| dx. But, by the dominated convergence theorem, this integral vanishes because the kernel K and the copula density are bounded from above and below, the copula density is Lipschitz continuous and, from [START_REF] Levine | Maximum smoothed likelihood for multivariate mixtures[END_REF], N f t l kj converges uniformly to N f * kj as l → ∞.

The first term in ( 6) is bounded by

g(x)w k (x; π, θ, ϕ * ) log c(F t l k1 (x 1 ), . . . , F t l kd (x d ); θ k ) c(F * k1 (x 1 ), . . . , F * kd (x d ); θ k ) dx + g(x)w k (x; π, θ, ϕ * ) log c(F * k1 (x 1 ), . . . , F * kd (x d ); θ k ) c(F t l-1 k1 (x 1 ), . . . , F t l-1 kd (x d ); θ k ) dx .
But again this bound goes to zero by similar arguments. This finishes the proof.

Implementation

Details about the implementation of the algorithm of Section 2 are given below.

Initialization. For initialization, the data are partitioned into K groups by a k-means algorithm. The initial weights π 0 are set equal to proportions of observations belonging to each group. The marginal densities ϕ 0 are initialized by standard kernel density estimation methods. The marginal density of the kth group in the jth dimension is set to its kernel density estimate calculated from the projection of the data belonging to the kth group into the jth dimension. The bandwidths are specified by standard bandwidth selection methods (Silverman, 1998, p. 47-48). A value for the bandwidth of the marginal of the kth group in the jth dimension h kj is selected by applying a bandwidth selection method to the projection of the data belonging to the kth group into the jth dimension. The bandwidth selection method used consists of taking

h kj = 1.06A kj n -1/5 kj
, where A kj is the minimum between the standard deviation of the data and, the interquartile range divided by 1.34 [START_REF] Scott | Multivariate density estimation: theory, practice, and visualization[END_REF]. The initial copula parameters θ 0 are set to the value corresponding to the independence copula.

Choice of the kernel. It is well known in kernel density estimation that the choice of the kernel has little impact on the estimates [START_REF] Wiley | Density estimation for statistics and data analysis[END_REF]. Therefore, the Gaussian kernel was chosen for convenience.

Bandwidth selection. Once the bandwidths have been initialized, they can be kept fixed or be updated from one iteration to another. In the latter, each observation x i is assigned the cluster that maximizes the current value of w k (x i ; π, θ, ϕ) over k = 1, . . . , K and the same bandwidth selection method as in the initialization step is applied.

Numerical evaluation of the integral (3). A bottleneck of the algorithm is the numerical evaluation of the integral in (3). Indeed, the quantity log f (u) might be close or even equal to -∞ in some regions of the integration domain. Moreover, if the values of K h (x -u) are zero or close to zero, this may create numerical issues of the kind "0 × ∞". To avoid those issues, two remedies are implemented. First, we substitute max{f (u), ε} for f (u) where ε is some tolerance threshold. We arbitrarily set ε = 10 -5 . Second, we truncate the domain of integration. After thresholding and a change of variables, the integral to evaluate becomes ∞ -∞ K h (u) log[max{f (x -u), ε}] du. We evaluate the integral on (0 ± 1.96h) instead of the whole real line, retaining about 95% of the mass of the kernel.

Stopping criterion. To terminate the algorithm, we may let the algorithm run an arbitrary number of steps and, in retrospect, visually check the convergence of the sequence of the objective function values, or we may stop the algorithm once some criterion has been reached. One possible stopping criterion is the relative increase of the objective function. That is, if ¯ t = g(x) log ǧ(x; π t , θ t , ϕ t ) dx denotes the objective function to be maximized at step t of the algorithm, then the algorithm may be stopped as soon as the inequality | ¯ t+1 -¯ t | < ε| ¯ t | occurs k times in a row. In practice, we arbitrarily set ε = 10 -2 and k = 3.

Choice of the number of clusters. To estimate the number of mixture components K in the mixture model (1-2), the algorithm of Section 2 is run with with several values of K. To select the "best" model, we use the pseudo-AIC criterion introduced in Mazo (2017), namely "maximum smoothed semiparametric log-likelihood times sample size minus number of copula parameters". Note that in the definition above we need to multiply by the sample size because the smoothed semiparametric log-likelihood in ( 4) is defined as an expectation and hence, contrarily to [START_REF] Mazo | A semiparametric and location-shift copula-based mixture model[END_REF], the sample version is a sample average.

Simulation studies 5.1 A first study

Five hundred replications of four independent artificial datasets of sizes n = 300, 500, 700, 900 were generated from the mixture model ( 1)-( 2) with K = 3 clusters of equal proportions, FGM copulas with parameters -0.5, 0.5, 0 and marginals as in Table 1, where N(µ, σ 2 ) and L(µ, σ 2 ) refer to the normal and Laplace distributions with mean µ and standard deviation σ, respectively. (The density of a L(µ, σ 2 ) distribution is then given by f (x) = e - √ 2|x-µ|/σ /( √ 2σ) for any real x.) The algorithm of Section 2 was run with K = 3 to estimate the cluster proportions, the copula parameters and the marginal densities. Initialization was carried out as described in Section 4.

The bandwidths were kept fixed after initialization. The algorithm was stopped after 50 iterations.

cluster 1 cluster 2 cluster 3 dim 1 N(-3, 2 2 ) N(0, 0.7 2 ) N(3, 1.4 2 ) dim 2 L(0, 0.7 2 ) L(3, 1.4 2 ) L(0, 2.8 2 )
Table 1: Marginals used for the numerical experiment of Section 5.1.

Figure 1 shows the values of the empirical smoothed log-likelihood (4) at each step of the algorithm for the first ten replications in the case n = 300 and n = 900.

All of the trajectories look monotonic. It was numerically calculated that, out of the N = 500 trajectories, only 17 were non-monotonic for n = 300 at the 10 -5 precision.

This number goes down to 1 for n = 500, and zero for n = 700 and n = 900. This suggests that the algorithm of Section 2 may indeed be monotonic for the copula and marginal families chosen above.

Figure 2 shows the sum of the estimated squared biases and variances for the copula parameter vector. The variance is at least 10 times higher than the squared bias for all values of n = 300. The variance decreases with n at a rate about that of the "parametric" rate 1/n: the variance at n = 900 is between 2.18 and 3.02 times smaller than the variance at n = 300.

Figure 3 shows the marginal density estimates at the last step of the algorithm for n = 900, for the last replication. The estimates agree well with the true marginal densities. We noticed, however, that they were similar to the initial estimates.

Sensitivity to initialization

To assess the impact of initialization on the results of the algorithm of Section 2, algorithm-compare the bottom row of Figure 4 with the top row of Figure 3. This is in sharp contrast with the bottom row of Figure 3, where the estimates were good.

It seems that initialization plays a key role in the final performance of the algorithm. This is confirmed by comparing the values of the three components of the estimated copula parameter vector across iterations, depicted in Figure 5. We see in Figure 5(a)

that one of the sequence of estimates seems to have not converged, while the others have their values stuck at -1 and 1, which is in general not an indication that estimation was performed correctly. By contrast, Figure 5(b) depicts stable and reasonable estimates. In sum, a Gaussian mixture modeling step during initialization of the algorithm produced poor estimates.

Estimation of the number of mixture components K

A numerical experiment was carried out to see whether the pseudo-AIC criterion described in Section 4 is able to select the correct number of components. A number of 500 synthetic datasets of size 300 were generated from a mixture model with three components of equal weights. The components are bivariate normal distributions with means (0,3), (3,0), (-3,0), standard deviations (

√ 2, 1/ √ 2), ( √ 2, 1/ √ 2), ( √ 2, 1/ √
2) and correlations 0.5, 0.5, 0.5. Gaussian copulas were assumed for all components. For each dataset, the mixture model (1-2) was fitted with the algorithm of Section 2 for K = 2, . . . , 5, where K denotes the number of components of the mixture model, and the pseudo-AIC criterion of Section 4 was computed to select the number of clusters.

Initialization was carried out as described in Section 4. The bandwidths were updated at each step of the algorithm. The stopping criterion described in Section 4 was used to terminate the algorithm.

The results are reported in Figure 6. Among the 500 estimates, 402 (standard error 9) were correct, and 98 were incorrect (standard errors 7 and 6 for K = 4 and K = 5, respectively). The chart suggests that the pseudo-AIC criterion is reasonable. This is consistent with the findings in [START_REF] Mazo | A semiparametric and location-shift copula-based mixture model[END_REF].

6 Real data analysis

The iris dataset

The iris dataset has n = 150 observations of d = 4 variables (sepal and petal length and width) belonging to three groups ("setosa", "versicolor", "virginica"). For simplicity and illustrative purposes, only two variables were considered (sepal and petal length). The algorithm of Section 2 was run with Gaussian copulas. Initialization was carried out as described in Section 4. The stopping criterion was used to terminate the algorithm. For bandwidth selection, the two strategies described in Section 4 were tested: the first consists of keeping the bandwidth fixed after initialization and the second consists of updating the bandwidth at each step of the algorithm. For the first strategy, the algorithm was successfully run with K = 2, . . . , 6 clusters. With 7 clusters, one of the clusters became empty and the algorithm was stopped. For the second strategy, the algorithm was successfully run with K = 2, 3, 4 clusters.

The values of the pseudo-AIC criterion are depicted in Figure 7. In Figure 7(a), the values increase and stabilize at K = 4. In Figure 7(b), the presence of a plateau is less clear. The pseudo-AIC increases as K increases, and hence a reasonable choice would also be K = 4. The obtained classification results are reported in Figure 8. We see a clear difference between the two bandwidth selection strategies. The classification results for the case K = 3 are reported in Figure 9. Here, the results of the two bandwidth selection strategies are similar (top row) and better reflect the true partitioning of the data than the results of the Gaussian mixture model fitted with 3 clusters (bottom row). For the Gaussian mixture model, the optimal number of clusters according to the BIC criterion is two.

The wine dataset

To illustrate the practical performance of our method, we will apply it to the analysis of the wine dataset that has been analyzed earlier in [START_REF] Bouveyron | Model-based clustering and classification for data science: with applications in R[END_REF] We conduct our analysis of this dataset based on these two variables. As a tool for selection of the number of clusters, we use the pseudo-AIC criterion introduced earlier in Section 4. To do so, we select a range of the possible number of clusters from 2 to 8 and compute the value of pseudo-AIC for each of these choices. Gaussian copula was used to model the dependence between the two variables. Initialization was carried out as described in Section 4. The bandwidths were updated at each iteration of the algorithm. The stopping criterion described in Section 4 was used to terminate the algorithm. The result is illustrated in Figure 10.

The result suggests the choice of either K = 5 or K = 8 as a possible number of clusters. Since the choice of five clusters produces an obviously more parsimonious model, we proceed with it. (Note that, when using BIC as a model selection criterion, Bouveyron et al ( 2019) also comes up with two possible models based on either 3 or 7 clusters with different respective covariance matrix structure.) At first sight, the choice of 5 clusters does not seem to be a very reasonable one since there are only three types of wine described by this dataset. However, we will see later that, nevertheless, this solution describes the true classification quite well. The resulting classification is illustrated in Figure 11.

Note that the red group (Grignolino) has about 9 observations separated from the main cluster. These observations, that were the source of confusion for Gaussian density mixture based solutions of [START_REF] Bouveyron | Model-based clustering and classification for data science: with applications in R[END_REF], also present some difficulties for our approach as well. Most of these observations have been separated into a separate (green) cluster. Moreover, observations with large values of the Flavonoids indicator have become, apparently, a reason for creation of yet another (blue) cluster.

At the same time, the confusion matrix of our classification that compares it with the partition into three wine types suggests that our classification is not widely off the mark. Indeed, consider the confusion matrix given in Table 2. Note that Barbero samples are split between Clusters 1 and 4 while Barolo samples are split (with the exception of just one sample) between Clusters 3 and 5. et al (2019). Thus, we believe that our approach provides an adequate clustering and classification analysis of the wine dataset.

Conclusion

An algorithm was designed and implemented to estimate the parameters of copulabased semiparametric mixture models. The model considered is a very general one since it does not impose any specific structure (such as the location-scale assumption) on marginal densities. The algorithm is deterministic, and hence always returns the same result if fed with the same initial point. Good performance was obtained in illustrative numerical examples, which suggests that the algorithm may indeed be monotonic under appropriate conditions. However, its theoretical analysis proved to be challenging and only partial results were obtained for versions of the algorithm where either the copula parameter or the marginals were fixed. A future avenue of research may consist of rejecting those updates where the smoothed log-likelihood does not increase and investigate whether convergence results of [START_REF] Meyer | Sufficient conditions for the convergence of monotonic mathematical programming algorithms[END_REF]; [START_REF] Zangwill | Nonlinear Programming-A Unified Approach[END_REF] could be applied. To simplify, the full parametric case may first be considered. To improve the numerical implementation of the algorithm, the integral (3) may be computed using other methods, such as Qiang (2010).

Fig. 4: Estimated marginal densities of the three clusters and the two dimensions for the last dataset generated in Section 5.1 with n = 900 and initialization by fitting a Gaussian mixture model. The top and bottom rows contain the results after the initialization step and at the last iteration of the algorithm, respectively. The column on the left and on the right contain the first and the second dimensions, respectively. 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0

  n, are observations from the target density g. The objective function to be maximized is then the empirical version of the smoothed log-likelihood, given by 1 n n i=1 log ǧ(X i ; π, θ, ϕ) (up to an additive constant). Here the bandwidths of the nonlinear smoothers are allowed to depend on the data.

Fig. 1 :

 1 Fig. 1: Values of the empirical smoothed log-likelihood at each step of the algorithm, for the first ten replications. Black plain lines: n = 900. Blue dotted lines: n = 300.

Fig. 2 :

 2 Figure3). Intriguingly, the estimates have deteriorated at the last iteration of the

  pp. 60-65. This dataset contains 27 physical and chemical measurements on 178 wine samples of three types-Barolo, Grignolino, and Barbera. The dataset is publicly available as a part of the pgmm R package (McNicholas, 2016). Bouveyron et al (2019) conducted a preliminary principal components analysis and selected 5 variables with the highest loadings on each of the first five principal components. Moreover, they noted that just two variables-Flavonoids and Color Intensity-seem to give a strong visual suggestion of clustering, based on the pairs plot. Their analysis is based on the use of multivariate Gaussian density mixtures with various covariance matrix structures.

Fig. 3 :

 3 Fig.3: True and estimated marginal densities of the three clusters and the two dimensions for n = 900 (last replication). The top row contains the true marginals and the column on the left contains the first dimension. The marginal estimates are those found at the last step of the algorithm.

Fig. 5 :

 5 Fig. 5: Values of the three components of the estimated copula parameter vector across iterations for the last dataset generated in Section 5.1 with n = 900 and initialization by fitting (a) a Gaussian mixture model and (b) a k-means algorithm.

Fig. 6 :

 6 Fig. 6: Frequency of the various values for the estimated number of components. The frequency for K = 2 is zero. The vertical bars correspond to the Monte Carlo asymptotic confidence intervals of level 95%.

Fig. 7 :

 7 Fig. 7: The value of pseudo-AIC model selection criterion for several values of the cluster number (iris dataset). Left: the bandwidth is kept fixed after initialization. Right: the bandwidth is updated at every step of the algorithm.

Fig. 8 :

 8 Fig. 8: Iris data: classification based on the choice of 4 clusters. Top left: algorithm of Section 2 with the bandwidths kept fixed after initialization. Top right: with the bandwidths updated at every step. Bottom left: results for the Gaussian mixture model with 4 clusters. Bottom right: true classification.

Fig. 9 :Fig. 10 :

 910 Fig. 9: Iris data: classification based on the choice of 3 clusters. Top left: algorithm of Section 2 with the the bandwidths kept fixed after initialization. Top right: with the bandwidths updated every step. Bottom left: results for the Gaussian mixture model with 3 clusters. Bottom right: true classification.

Fig. 11 :

 11 Fig. 11: Wine data: true classification (left) and classification based on the choice of 5 clusters (right). The measurements shown are Flavonoids and Color Intensity.

Table 2 :

 2 Confusion matrix of the 5clusters solutionIf Clusters 1 and 4, on one hand, and Clusters 3 and 5, on the other hand, are merged, one ends up with a 3-clusters solution whose misclassification rate is only 12 178 . For comparison purposes, when[START_REF] Bouveyron | Model-based clustering and classification for data science: with applications in R[END_REF] merge the necessary clusters of their 7-clusters solution, the resulting misclassification rate is 11 178 . Even if such a merger is not contemplated, the misclassification rate of our 5-clusters solution is 61 178 which is less then 52% misclassification rate of the 7-clusters solution of Bouveyron

		Barolo Grignolino Barbero
	Cluster 1	0	0	28
	Cluster 2	1	60	0
	Cluster 3	29	7	0
	Cluster 4	0	2	20
	Cluster 5	29	2	0

Acknowledgments

Michael Levine's research has been partially funded by the NSF-DMS grant # 2311103.

We thank two anonymous reviewers for helpful comments that led to an improved version of the manuscript.

Declaration

The authors have no relevant financial or non-financial interests to disclose.


Keywords: catecholamine release, cultured chromaffin cell, fetal rat adrenal medulla, gambierol, potassium current

teaching and research institutions in France or abroad, or from public or private research centers.

Gambierol action on K(+) currents and catecholamine release in cultured single chromaffin cells from fetal rat adrenal medulla

Jordi Gambierol, characterized by a transfused octacyclic polyether core, was first isolated and chemically typified from cultured Gambierdiscus toxicus dinoflagellates collected in French Polynesia. Subsequently, distinct groups in Japan and the USA, using various strategies, achieved its total chemical synthesis. This allowed detailed studies on its mode of action. Gambierol inhibits voltage-gated K + (KV) channels in various excitable and non-excitable cells, as well as in motor nerve terminals of the skeletal neuromuscular junction. In the present study, we have investigated first, the effects of nanomolar concentrations of gambierol on K + current of cultured chromaffin cells from fetal rat adrenal medulla using perforated patch-clamp current-recordings. Our results show that gambierol only blocked a small component of the total K + current, and affected neither calcium-activated K + (KCa) nor ATP-sensitive K + (KATP) channels, as revealed using apamin and iberiotoxin (selective KCa channel blockers) and glibenclamid (KATP channel blocker). After inhibiting KATP and KCa channel activation, the gambierol concentration blocking 50% the K + current-component (IC50) was 7.6 ± 1.1 nM.

The repolarizing phase of all-or-none elicited action potentials, recorded under current-clamp conditions (triggered by 1-ms depolarizing pulses), was sensitive to the action of gambierol but insensitive to the action of apamin and iberiotoxin, indicating that KCa channels do not participate in the modulation of action potential duration triggered by short depolarizing pulses. The use of simultaneous patch-clamp and single-cell amperometry allowed controlling the membrane potential and detecting exocytosis events (with a carbon electrode polarized to +650 mV to allow the oxidation of released catecholamines). Such recordings revealed that gambierol did not modify the membrane potential following 14-second depolarizing current-steps, and did not increase significantly the number of exocytotic catecholamine release events with respect to controls. The addition of KCa channel blockers (in the continuous presence of gambierol) enhanced the membrane depolarization by about 15 mV (during the 14-second current step), and at the same time increased significantly the number of exocytotic events related to catecholamine secretion. Such enhanced depolarization induced by the KCa channel blockers probably brings the membrane potential above the activation threshold of high-voltage activated CaV channels triggering both Ca 2+ influx and subsequent catecholamine secretion. These results emphasize the diversity of KV channels in chromaffin cells from fetal rat adrenal medulla and highlight the modulatory role played by K Ca channels in the control of exocytosis in the absence of splanchnic innervation.

  Molgó, Roland Bournaud, Sébastien Schlumberger, Makoto Sasaki, Haruhiko Fuwa, Denis Servent, Evelyne Benoit To cite this version: Jordi Molgó, Roland Bournaud, Sébastien Schlumberger, Makoto Sasaki, Haruhiko Fuwa, et al.. Gambierol action on K(+) currents and catecholamine release in cultured single chromaffin cells from fetal rat adrenal medulla. RT27 -27th Meeting of the French Society of Toxinology -Toxins: Mr Hyde or Dr Jekyll?, Dec 2021, Paris (Virtual), France. hal-04470261

Poster Gambierol action on K + currents and catecholamine release in cultured single chromaffin cells from fetal rat adrenal medulla

  Mr Hyde or Dr Jekyll? 27 èmes Rencontres en Toxinologie -27 th Meeting on Toxinology -e-RT27 Université Paris-Saclay, CEA, INRAE, Institut des sciences du vivant Frédéric Joliot, Département Médicaments et Technologies pour la Santé (DMTS), Service d'Ingénierie Moléculaire pour la Santé (SIMoS), ERL CNRS 9004, F-91191 Gif-sur-Yvette, France. 2 CNRS, Laboratoire de Neurobiologie Cellulaire et Moléculaire-UPR 9040, F-91198 Gif-sur-Yvette, France. 3 Tohoku University, Graduate School of Life Sciences, Sendai, Japan.

	48
	(Rencontres virtuelles -Virtual Meeting)
	9-10 décembre 2021 -December 9-10, 2021
	Session "Plant toxins, algal toxins, mycotoxins"
	Jordi MOLGÓ 1,2,* , Roland BOURNAUD 2 , Sébastien SCHLUMBERGER 2 , Makoto SASAKI 3 ,
	Haruhiko FUWA 4 , Denis SERVENT 1 , Évelyne BENOIT 1,2

Toxines: Mr Hyde ou Dr Jekyll? -Toxins:

1 4 Chuo University, Faculty of Science and Engineering, Department of Applied Chemistry,

 Bunkyo, Japan. *Corresponding author,


Go Green: Selected Configuration Interaction as a More Sustainable Alternative for High Accuracy Pierre-François Loos, 1, a) Yann Damour, 1 Abdallah Ammar, 1 Michel Caffarel, 1 Fábris Kossoski, 1 and Anthony Scemama 1 Laboratoire de Chimie et Physique Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS, France Recently, a new distributed implementation of the full configuration interaction (FCI) method has been reported [Gao et al. J. Chem Theory Comput. (2024) 10.1021/acs.jctc.3c01190]. Thanks to a hybrid parallelization scheme, the authors were able to compute the exact energy of propane (C 3 H 8 ) in the minimal basis STO-3G. This formidable task involves handling an active space of 26 electrons in 23 orbitals or a Hilbert space of 1.3 × 10 12 determinants. This is, by far, the largest FCI calculation reported to date. Here, we illustrate how, from a general point of view, selected configuration interaction (SCI) can achieve microhartree accuracy at a fraction of the computational and memory cost, via a sparse exploration of the FCI space. The present SCI calculations are performed with the Configuration Interaction using a Perturbative Selection made Iteratively (CIPSI) algorithm, as implemented in a determinant-driven way in the quantum package software. The present study reinforces the common wisdom that among the exponentially large number of determinants in the FCI space, only a tiny fraction of them significantly contribute to the energy. More importantly, it demonstrates the feasibility of achieving comparable accuracy using more reasonable and sustainable computational resources, hence reducing the ever-growing carbon footprint of computational chemistry. The full configuration interaction (FCI) method provides the exact solution of the Schrödinger equation in a given basis set, the number of electrons and orbitals defining the size of the Hilbert space. As such, it is often considered the ultimate answer when one wants to compare the performances of approximate methods. [1][2][3] However, the computational cost and memory requirement of a FCI calculation increases exponentially fast with the size of the system making it feasible on small systems only. FCI calculations on realistic molecular systems were achieved in the late 1980's, largely owing to the groundbreaking work of Knowles and Handy. [4][5][6][7] The significant milestone of surpassing the one-billion determinant barrier was accomplished in 1990 by Olsen et al. 8 a) Electronic mail: loos@irsamc.ups-tlse.fr While there were several improvements during the first decade of the twenty-first century, [9][10][11] the pace notably slowed thereafter. 12,13 This slowdown can be ascribed partially to the emergence of alternative methodologies, notably selected configuration interaction (SCI) methods, which had lain dormant [14][15][16][17] but underwent a resurgence in the early 2010's. 18,19 The essence of these methods lies in their sparse exploration of the Hilbert space, focusing only on the most energetically relevant determinants. This strategy stems from the realization that within the vast FCI space, only a tiny fraction of determinants significantly impact the energy.

Nowadays, SCI methods alongside related techniques like density-matrix renormalization group (DMRG) approaches [20][21][22] and others, [23][24][25][26][27][28][29][30][31][32] have become central to modern electronic structure theory. [1][2][3] Their main objective is to compute reference correlation and excitation energies in small molecular systems with exceptional accuracy, often rivaling FCI results. 2,[33][34][35][36][37][38][39][40] Modern implementations of SCI encompass methodologies such as CIPSI (CI using a Perturbative Selection made Iteratively), 1,16,18,19,33,[START_REF] Caffarel | Using CIPSI Nodes in Diffusion Monte Carlo[END_REF][START_REF] Garniron | [END_REF][43][44][45][46] adaptive sampling CI (ASCI), [47][48][49][50] semistochastic heatbath CI (SHCI), 35,36,[51][52][53][54][55] and iterative CI (iCI). [56][57][58][59][60] Stochastic CI techniques, like Monte Carlo CI (MCCI) 61,62 and FCI quantum Monte Carlo (FCIQMC), [63][64][65][66][67][68] adopt a similar approach, utilizing stochastic wave function representations to identify crucial determinants.

Recently, Gao et al. 69 reported a new implementation of FCI using a hybrid parallelization scheme combining multiprocessing with MPI and multithreading with OpenMP. They were able to compute the exact energy of propane (C 3 H 8 ) in the STO-3G basis set. This calculation corresponds to a Hilbert space of size 1.3 × 10 12 determinants. This formidable feat effectively breaks the one-trillion determinant barrier and is by far the largest FCI calculation reported to date. According to Ref. 69, this calculation required running 512 processes on 256 nodes for a total wall time of 113.6 hours. Their hardware consisted of nodes equipped with two Intel Xeon Gold 6148 CPUs, boasting a total of 384 GB of memory and 40 logical cores per CPU, enabling the execution of 40 threads per process. The most memory-consuming task requires an astonishing amount of 19 TB of memory.

In this study, our objective is to demonstrate that SCI calculations can be conducted routinely with significantly more reasonable computational resources while still producing energies that closely match those obtained through FCI, thus contributing to the reduction of the ever-growing carbon footprint associated with computational chemistry.

Our calculations are performed on a single node of the CALMIP supercomputer center (Toulouse, France), which has been in production since September 2018. This node is a dual-socket Intel Skylake 6140 CPU@2.3 Ghz with 192 GB of memory for a total of 36 physical CPU cores. The numbers given for the energy consumption of a calculation represent all the energy used by the compute node for the duration of the run. These values were obtained from the database of the SLURM job scheduler, which was configured such that the acct_gather_energy plugin obtains the energy consumption from the Baseboard Management Controller (BMC) via the Intelligent Platform Management Interface (IPMI) protocol.

The present SCI calculations have been carried out with the CIPSI algorithm. 16,18,19,[START_REF] Garniron | [END_REF][43][44] One of the main differences between conventional FCI implementations and the CIPSI implementation reported in Ref. 44 is that the former is integral-driven while the latter is determinantdriven. (Note that the SHCI method, due to its different selection criteria, is also integral-driven. 35,51,52 ) In a nutshell, CIPSI iteratively increases the size of the so-called variational space characterized by its (zeroth-order) variational energy E var (see Fig. 1). As discussed further

|0 |1 |2 |3 |4 FIG.
1. Schematic presentation of the growth of the variational space as a function of the CIPSI iterations. Starting from a reference Slater determinant |0⟩ built with a given set of orbitals, the variational space is enlarged systematically via the inclusion of the most energetically important determinants of higher excitation degrees (filled circles) to build larger multideterminant wave functions (|1⟩, |2⟩, |3⟩, etc). In such a way, the Hilbert space is explored sparsely by leaving behind the determinants that do not contribute significantly to the energy (empty circles).

below, the variational energy (or SCI energy) converges slowly with respect to the number of determinants in the variational space N det . To improve this convergence, E var is corrected by its second-order perturbative energy, E PT2 , computed within Epstein-Nesbet perturbation theory. [START_REF] Garniron | [END_REF] The sum E var + E PT2 defines the SCI+PT2 energy. Because there exists an approximate linear relationship between E var and E PT2 when E PT2 is small enough, a linear extrapolation of E var as E PT2 → 0 is performed to produce the final FCI estimate. 35,[START_REF] Burton | Rationale for the extrapolation procedure in selected configuration interaction[END_REF] In the following, this extrapolated value is named exFCI.

As a first step, we perform a preliminary run where, in addition to iteratively increasing the size of the variational space, we further optimize the orbitals via the minimization of the variational energy. 45,54 This orbital optimization procedure is performed at each CIPSI iteration up to 1.8 × 10 6 determinants. This run takes 2h56 of wall time, 4.9 GB of memory, and 758 W h of energy consumption. The most memory-consuming task is the Davidson diagonalization of the largest variational space which requires 2.6 GB, while the most expensive PT2 calculation only requires 170 MB of memory. (The orbital optimization step could easily be made globally less expensive without altering the overall accuracy by stopping the calculation earlier.) Then, a new CIPSI run is performed with these optimized orbitals. In Table I, we report the evolution of the SCI (E var ), SCI+PT2 (E var + E PT2 ), and exFCI (E exFCI ) energies as functions of N det . The errors with respect to the FCI energy, E FCI , TABLE I. Convergence of the SCI (Evar), SCI+PT2 (Evar + EPT2) and exFCI (EexFCI) energies of propane (C3H8) computed in the STO-3G basis with respect to the number of determinants included in the variational space (N det ). These calculations are performed with energetically-optimized orbitals. The error with respect to the FCI value of the SCI (∆SCI), SCI+PT2 (∆SCI+PT2), and exFCI (∆exFCI) energies are also reported. The standard errors associated with the 3-point linear fitting procedure performed to obtain the exFCI energies are reported in parenthesis. of Gao et al. are also reported. 69 The exFCI values are obtained via a 3-point linear extrapolation of the variational energy using the three smallest E PT2 available at a given stage. The extrapolation procedure of E var as a function of E PT2 is illustrated in Fig. 2 for the largest variational space (N det = 127 × 10 6 ). For the sake of comparison, the CCSD, CCSD(T), CCSDT, and CCSDTQ energies, computed with cfour, [START_REF] Matthews | [END_REF] are also reported in Fig. 2 and Table II. While the CCSD, CCSD(T), and CCSDT calculations require only a few seconds of wall time and a minimal amount of memory, the CCSDTQ calculation is more expensive as it requires 42 minutes of wall time III presents the wall time, maximum memory usage, and total energy consumption for increasingly larger CIPSI calculations. For large variational space, the Davidson diagonalization performed to compute the variational energy 44 is the most memory-consuming task while the semistochastic calculation of the second-order perturbative correction [START_REF] Garniron | [END_REF] I.

N

accuracy of the exFCI energy (with respect to the FCI value) is also reported. Achieving microhartree accuracy necessitates approximately 16 × 10 6 determinants, 25 GB of RAM, and slightly over 3 hours of wall time, totaling approximately 1.1 kW h. Comparatively, a CIPSI calculation employing N det = 32 × 10 6 demands twice the memory, four times the energy consumption, and an additional 10 hours of wall time to achieve an overall accuracy of 0.2 µE h .

We believe that these numbers illustrate nicely that SCI calculations are a more eco-friendly alternative compared to FCI calculations, especially for achieving high accuracy and generating reference values for benchmarking other computational methods. As the variational space expands (see Table III), SCI calculations become more resourceintensive both in terms of computational resources and energy consumption. For example, to perform a calculation up to 126 × 10 6 determinants, almost 6 days of wall time are required for an energy consumption of 42 kW h. Hence, there remains ample scope for further optimization and improvement.

FIG. 2 .

 2 FIG.2. Variational energy as a function of the second-order perturbative correction for the SCI calculations performed on propane (C3H8) in the STO-3G basis and with energeticallyoptimized orbitals. The 3-point linear fit is based on the three largest variational wave functions reported in TableI.

TABLE II .

 II Energy of propane (C3H8) computed at different levels of theory with the STO-3G basis. The error with respect to the FCI value is also reported. Perturbatively-corrected variational energy obtained with N det = 32 × 10 6 . d Extrapolated FCI value obtained via a 3-point linear fit using N det = 32 × 10 6 as the largest variational space.

	det	SCI energy	SCI+PT2 energy	exFCI energy
		Evar (E h )	∆SCI (mE h ) Evar + EPT2 (E h ) ∆SCI+PT2 (mE h ) EexFCI (E h )	∆exFCI (mE h )
	1907 -117.082 104 21	18.018	-117.097 761 32	2.361
	3877 -117.087 335 85	12.787	-117.098 483 22	1.639
	7756 -117.090 941 51	9.181	-117.099 003 85	1.119
	15 517 -117.093 617 12	6.506	-117.099 369 72	0.753
	31 044 -117.095 587 36	4.535	-117.099 595 81	0.527
	62 113 -117.097 029 68	3.093	-117.099 781 79	0.341
	124 570 -117.097 993 09	2.130	-117.099 893 33	0.229
	249 144 -117.098 757 92	1.365	-117.099 982 64	0.140	-117.100 143 5(7)	-0.0208(7)
	498 290 -117.099 265 39	0.857	-117.100 044 65	0.078	-117.100 149(3)	-0.026(3)
	996 650 -117.099 599 16	0.524	-117.100 076 04	0.047	-117.100 139(9)	-0.016(9)
	1 993 314 -117.099 811 92	0.311	-117.100 096 21	0.026	-117.100 125 8(1)	-0.0031(1)
	3 986 707 -117.099 939 42	0.183	-117.100 107 96	0.015	-117.100 125 5(3)	-0.0028(3)
	7 973 418 -117.100 017 98	0.105	-117.100 114 98	0.0077	-117.100 124 8(2)	-0.0021(2)
	15 946 880 -117.100 065 98	0.057	-117.100 118 64	0.0040	-117.100 123 7(5)	-0.0010(5)
	31 893 835 -117.100 093 52	0.029	-117.100 120 66	0.0020	-117.100 122 89(6) -0.000 21(6)
	63 788 022 -117.100 108 27	0.014	-117.100 121 70	0.000 98	-117.100 122 76(3) -0.000 08(3)
	126 541 040 -117.100 115 87	0.0068	-117.100 122 19	0.000 49	-117.100 122 67(3) +0.000 01(3)
	Method	Energy (E h )		Error wrt FCI
	FCI a	-117.100 122 681 461
	CCSD	-117.098 767		1.355 mE h
	CCSD(T)	-117.099 708		0.414 mE h
	CCSDT	-117.099 942 158	0.181 mE h
	CCSDTQ	-117.100 120 230	2.451 µE h
	SCI b	-117.100 093 52	0.029 mE h
	SCI+PT2 c	-117.100 120 66	2.021 µE h
	exFCI d	-117.100 122 89(6)	-0.21(6) µE h

a Reference 69. b Variational energy obtained with N det = 32 × 10 6 . c

TABLE III .

 III Wall time, maximum memory consumption, total energy consumption, and error with respect to FCI of the exFCI energy for increasingly large CIPSI calculations performed on propane (C3H8) in the STO-3G basis. The energy consumption represents the total consumption of the compute node measured by the BMC for the entire job duration (see main text).

	N det	Wall time	Memory	Energy	Error
		(hh:mm) consumption consumption wrt FCI
	2 × 10 6	00:14	5.3 GB	74 W h	3 µE h
	4 × 10 6	00:33	8.1 GB	176 W h	3 µE h
	8 × 10 6	01:19	15 GB	438 W h	2 µE h
	16 × 10 6	03:12	25 GB	1.1 kW h	1 µE h
	32 × 10 6	13:16	47 GB	4.1 kW h	0.2 µE h
	64 × 10 6	43:54	83 GB	13 kW h	0.08 µE h
	127 × 10 6	138:44	138 GB	42 kW h	0.01 µE h
	and 508 MB of memory for an energy consumption of
	192 W h. The CCSDTQ energy and the SCI+PT2 energy
	obtained with N det = 32 × 10 6 have similar accuracy and almost reach microhartree accuracy. Sub-microhartree
	accuracy can be obtained via extrapolation, the exFCI
	energy having an error of only -0.21(6) µE h compared to FCI.
	Table

ACKNOWLEDGMENTS

The authors would like to thank Nicolas Renon (CALMIP) for his technical support. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 863481). Additionally, it was supported by the European Centre of Excellence in Exascale Computing (TREX), and has received funding from the European Union's Horizon 2020 -Research and Innovation program -under grant agreement no. 952165. This work was granted access to the HPC resources of CALMIP supercomputing center under the allocation 2024-18005.

DATA AVAILABILITY STATEMENT

The data that supports the findings of this study are available within the article and its Supporting Information.

SUPPORTING INFORMATION

See the Supporting Information for the geometry of propane and output files associated with the CC calculations, orbital optimization, and CIPSI calculations. Files gathering information about wall time, memory usage, and energy consumption are also provided.


Oana Stan

Vincent Thouvenot

Aymen Boudguiga

Katarzyna Kapusta

Martin Zuber

Renaud Sirdey

A Secure Federated Learning: analysis of different cryptographic tools

Keywords: federated learning, homomorphic encryption, multi-party computation, differential privacy

Federated Learning is established as one of the most efficient collaborative learning approaches aiming at training different client models using private datasets. By private, we mean that clients' datasets are never disclosed as they serve to train clients' models locally. Then, a central server is in charge of aggregating the different models' weights. The central server is generally a honest-but-curious entity that may be interested in collecting information about clients datasets by using model inversion or membership inference. In this paper, we discuss different cryptographic options for providing a secure Federated Learning framework. We investigate the use of Differential Privacy, Homomorphic Encryption and Multi-Party Computation (MPC) for confidential data aggregation while considering different threat models. In our homomorphic encryption approach, we compare results obtained with an optimized version of the Paillier cryptosystem to those obtained with BFV and CKKS. As for MPC technique, different general protocols are tested under various security assumptions. Overall we have found HE to have better performance, for a lower bandwidth usage.

Introduction

Federated Learning (FL), a recent ML distributed paradigm allowing to train a common model without sharing sensitive local data, seems an appealing option to build high-quality and robust models with large amounts of training sets from various sources. In FL a set of participating clients collaboratively learn a global model by uploading their local models updates to a central server, which coordinates the training. This central server updates the global model by averaging the local model parameters, and sends it back to the data owners. Even if federated learning has the advantage of not sharing local data, recent research showed that the local data of distributed devices can be leaked through the local model parameters. Sharing intermediate models with the coordinator server, or among the participants, can lead to various privacy attacks, e.g., extracting participants' inputs or membership inference (e.g. [START_REF] Hitaj | Deep models under the gan: Information leakage from collaborative deep learning[END_REF], [START_REF] Melis | Exploiting unintended feature leakage in collaborative learning[END_REF]).

To address these problems, several works employ different cryptographic techniques such as Differential Privacy (DP), Homomorphic Encryption (HE) or Multi-Party Computation (MPC) to propose more privacy-preserving and secure federated learning ap-proaches. However, none of these solutions is yet completely satisfactory (in terms of performances of AI models and/or security constraints) and moreover they are often tested for different applications and on different datasets.

In this paper, we propose a first work on secure federated learning in which for the same use-case and datasets, different cryptographic solutions are conceived, implemented and compared. More precisely, our contributions are as follows:

• We propose different architecture configurations for a Secure Federated Learning taking into account threats coming from the aggregation server (by means of Homomorphic Encryption or Multi-Party Computation) and/or the other clients (by means of Differential Privacy) and we analyse the impact of these countermeasures.

• We present a comparison of three HE schemes (Paillier, BFV, CKKS) when applied to federated learning averaging.

• We evaluate the performance of FL aggregation implemented using general MPC protocols for four different security configurations. We show the cost of increasing the protection from a semihonest adversary to a dishonest one.

• We implement global differential privacy on a Federated Learning framework and evaluate the privacy of the resulting model.

• We provide a comparative analysis between MPC and Homomorphic Encryption solutions for our Federated Learning model in terms of accuracy of the AI model, the additional training overhead but also from a security point of view.

Overview of Federated Learning

The main idea of Three settings of FL exist: horizontal FL, vertical FL and Transfer FL. In this paper we focus on horizontal FL, or sample-based federated learning, which corresponds to the scenarios in which datasets share the same features. We assume honest participants and security against different adversarial settings for the server. For more details about FL, interested readers can refer to these recent surveys [Li et al., 2021, Kairouz andet al., 2021].

3 Architectural options for a more secure Federated Learning

HE-based variant

One of the security assumptions we made with this variant of architecture is an honest-but-curious aggregation server, i.e. a server following honestly the protocol and trying to infer as much information as possible from the data it has access to. Homomorphic Encryption can prevent this from happening, by providing a secure countermeasure to protect the models parameters during FL training. Before the training actually begins, the central server shares with the M clients the global architecture of the neural network that will be trained. For now, we assume a default setting, with a simple key distribution, in which each of the K clients selected randomly to participate in a round has the same pair of secret and public homomorphic keys (sk, pk). sk is the private key while pk is the public key. The preliminary step of key setup is independent of this work. For example, we can assume that the key generation is made locally by one random client. The latter can be randomly selected by the server or be the result of a leadership election protocol. All the clients will then share the same pair of (sk, pk). The central server holds only the pk required for the homomorphic evaluation of the global model. At each round of the training which is an iterative process, the server randomly selects K out of the M clients. Therefore, the only function to be computed encrypted is the federated averaging of the encrypted weights provided by the K clients participating to the FL round t. That is, the central server has to compute per round t the function:

w t+1 ←-∑ K k=1 n k n •w k t+1
, where n k is the size of the training set of client k, and n = ∑ K k=1 n k is the total size of the training datasets used in round t.

MPC-based variant

As discussed in the previous section, the homomorphic encryption is a good choice for securing the averaging by the aggregation server in the semi-honest setting. However, MPC can provide a solution in a situation where such a server is not available due to lack of trust between the participants. It can replace the semi-honest aggregator server by two or more semihonest or untrusted aggregator servers that would perform the aggregation in a distributed way. MPCbased distributed aggregators could solve two problems that a single semi-honest aggregator server cannot address. First is that the secure aggregator server cannot be also a participant of the computation: as a participant has the decryption key that the aggregator server cannot possess. For the same reason, an aggregation server is not supposed to collude with one or more participants. In MPC, the aggregators can be independent or not from the participants to the learning. Therefore, a participant to the learning can also have the role of an aggregator. Second, homomorphicencryption prevents the server from decrypting the data but does not come with mechanisms addressing the dishonest threat model, in which the aggregation server could deviate from the agreed protocol. For instance, an aggregation server corrupted by an active adversary could omit some of the inputs provided by the participants (and therefore modify the aggregation results, and in consequence the final model). Homomorphic encryption could be enriched with additional verification mechanisms that would detect dishonest behavior of the aggregator. However, implementing those mechanisms is not straightforward. In contrast, MPC protocols addressing the active adversary embed already mechanisms that enable to detect deviations of the participants from the agreed protocol. The advantages of MPC come at the cost of increased communications and expensive computations (especially in the active threat model). Therefore, although MPC is an interesting alternative to HE, it can be judged impractical for use-cases where the communications costs are an important factor.

In the example of architecture from Figure 2, the aggregation servers are different from the clients and a secret-sharing MPC algorithm is used to protect the data (computation is realized on the shares of inputs and its secure unless all shares of a given data are gathered).

Adding DP guarantees to FL training

In order to address threats coming from other honestbut-curious clients, the above frameworks can be completed with Differential Privacy (DP) guarantees.

In the case of a malicious server, one may use local DP with noises generated by each client. Otherwise, if one can trust the aggregation server for the noise generation, global DP can be applied. Let us also note that for the secure variant of FL using local DP and HE, specific quantization techniques are necessary for the noise sampling to ensure high DP guarantees (refer to [START_REF] Amiri | Compressive differentially-private federated learning through universal vector quantization[END_REF] for more details).

4 Experimental results

FL implementation with DP

We have implemented a federated learning with global differential privacy based on Tensorflowfederated. To demonstrate the results, we use the Turbofan dataset [START_REF] Saxena | Damage propagation modeling for aircraft engine runto-failure simulation[END_REF]. This learning process is similar to [START_REF] Mcmahan | Learning Differentially Private Recurrent Language Models[END_REF]. To simulate the federated learning between three participants, we split the data into three using the column "unit", this way each participant gets the full data. Two participants have 33 engines and one has 34 engines. The standard deviation of the noise that we add is equal to z/3, where z is a hyperparameter that we make varying and 3 comes from the fact that we have 3 participants. In Figure 3, we vary z and consider 20 rounds of training. We collect for each z the maximum validation accuracy after 20 rounds. We denote that the process does not converge after 20 rounds for the 6 strongest z. Moreover, for the strongest z values, we obtain the maximum validation accuracy for less than 8 rounds. The results shown in Figure 3 are quite intuitive: the more noise we add the less accurate the model is. The evolution of the accuracy during the training process gives a more precise understanding of what happens. In Figure 4, we display the categorical accuracy obtained at each round on the training dataset and on the test dataset. It displays on the same plot, the curves for the three different learning settings:

• Federated Learning with global differential privacy

• Federated Learning without differential privacy

Local training on the global dataset

In Figure 4 the differential privacy was computed with z=0.0.1, a small amount of noise. Its affects the convergence rate, the accuracy increases slowly but the performances are good. We also compute the differential privacy for z=0.1, a larger amount of noise, but still quite small. The learning is affected. There is some randomness in the evaluation accuracy. The model is learnt on the training dataset but its perfor- mances on the test dataset are not predictable. In our first results of this sandbox dataset, it seems that z=0.1 is the larger amount noise that allow us to learn a model, while ensuring some privacy. We evaluate the privacy guarantees obtained after 20 rounds in order to find a trade-off between model performances and data privacy. In Figure 5 we show the validation accuracy, and the vertical line at z=0.1 highlights the limit between a learning for which we achieve a convergence and a learning that does not ends with convergence. Here, we see that the convergence is achieved only with a good accuracy and that the resulting privacy metrics are high, and thus the privacy guarantees are low.

(F)HE for FL

All the experiments presented in this section are done with parameters that support a security level of λ = 128. We have trained three models of different sizes: 197, 1000 and 10000 weights. Note that 197 is the number of weights of the neural network that was trained, in the previous section, using FL with DP and the Turbofan dataset.

Implementation with Paillier

Our implementation using the Paillier cryptosystem relies on a simple packing technique as well as a faster Regev-style encryption function (encrypts by homomorphically adding a random subset of public encryptions of 0 to a pseudo-encryption of the message). Because we parameterize for exact arithmetic (as opposed to approximate arithmetic for some homomorphic schemes or depending on the parameters) we only need to concern ourselves with the time performance on the server here and the bandwidth use.

On the client side, the encryption time for a given client depends on the number of ciphertexts needed to be encrypted. That in turn is determined by the precision p (in bits) needed for the weights, the number of weights per client N w , and the number of training points used by every client n k . We assume for our tests, for simplicity sake, that all n k have the same values and that all clients have the same number of weights N w . Table 1 presents our implementation's time and bandwidth performance for several sets of parameters from small to big. Obviously the parameter sets are not exhaustively chosen and are designed to give a good idea of the scaling of our implementation. Because deciphering time is fixed (at roughly 0.05 seconds) for every set of parameters it is not included in the Table . The bandwidth size corresponds to what a client needs to send to the server. It gets back roughly the same amount of data, therefore for a total bandwidth, that needs to be doubled.

Implementation with BFV/CKKS

We did our tests for FL aggregation with BFV or CKKS schemes with Microsoft SEAL 3.7.2 library. For batching with BFV and CKKS, we set the encryption parameters with respect to the considered precision and the number of batched slots. For BFV, the plaintext modulus size is at least equal to the global precision p + ⌈log 2 (n k )⌉ + ⌈log 2 (K)⌉. Meanwhile for CKKS, the scaling factor size ⌈log 2 (δ)⌉ is at least equal to p + ⌈log 2 (n k )⌉ + ⌈log 2 (K)⌉ + ⌈log 2 (K × noise)⌉ where noise corresponds to the noise induced by a homomorphic addition. We do so to ensure that the final precision with CKKS is equal to ⌈log 2 (δ)⌉ -⌈log 2 (K × noise)⌉ = p + ⌈log 2 (n k )⌉ + ⌈log 2 (K)⌉. This precision is needed for ensuring a correct aggregation result.

Note that the FL aggregation takes at most 15 milliseconds (Tables 2 and3). It is very efficient with BFV or CKKS, as it only consists in computing additions (no multiplication needed). However, the ciphertext size reaches at most 1.8 MB with BFV.

MPC for FL

We implemented the MPC-based secure FL aggregation procedure using the recently released MP-SPDZ framework [START_REF] Keller | MP-SPDZ: A versatile framework for multi-party computation[END_REF] that allows fast translation of a python-like code into MPC operations and that implements more than 30 variants of MPC protocols. We measured the aggregation performance in terms of time and communication costs. The number of participants was set to 5 and models of 1000 weights were considered. A configuration of 2 or 3 computing parties was chosen, which is a common choice for MPC computations. Two setting were tested: honest majority and dishonest majority.

Three different MPC protocols were selected for the test: a Shamir's secret sharing-based protocol for semi-honest honest majority, the MASCOT protocol for dishonest setting, and a version of MAS-COT stripped of active security mechanisms for semihonest dishonest majority.

Performance results are presented in Table 4. They were obtained on a device with a 4-th generation i5 processor running at 1.30 GHz using 16 GiB of RAM on Linux. They confirmed that a straightforward implementation of distributed aggregation using general MPC protocols leads to a rather large overhead in terms of communication and computation costs. However, the advantage of this approach is that it enables to switch smoothly between different protocols and thus to easily adjust the security level and the number of parties. It shows that securing the computation against an active attacker is possible, although it comes at a high cost (that can be potentially acceptable for some of the use cases where securing is a much higher priority than the speed of the learning). Moreover, such implementation is robust to dynamic client or aggregation server dropouts.

Analysis

The preliminary results presented in the above sections showed interesting perspectives when securing Federated Learning with DP, FHE and MPC. As expected, there is a delicate trade-off to find between the guarantees offered by the global DP and the accuracy of the training model.

As for the results obtained with (F)HE and MPC, let us analyse and compare the concrete case of 5 clients and a model size of 1000 weights. The tests show that, under the hypothesis of a honestbut-curious server, the homomorphic encryption under various flavors (either optimized version of additive Pailier cryptosystem or classical batching of levelled BFV and CKKS) has better performance results in terms of execution times and bandwidths requirements than the general MPC protocols. This general conclusion of homomorphic encryption performing better for the aggregation in the context of federated learning remains valid for other testing parameters. Of course, the results with MPC could be ameliorated through specific protocols and it is still useful in situations in which homomorphic encryption cannot be used (see Section 3.2 for details).

  FL is highly similar to Distributed Learning (DL) where the data is spread among workers, while a central server is in charge of computing the updates. A global model is computed by averaging the local model trained on-device using the dataset owned locally by each worker (data-owner). The FL training is described below: 1. Central server sends the latest model parameters to the nodes 2. Data is collected at each node 3. Each local model is trained based on the latest pa-4. Updated model parameters are communicated back to the global model 5. Combine updates from each model and retrain the global model to get a new model 6. Restart from step 1

Figure 1 :

 1 Figure 1: Cross-silo federated learning architecture with Homomorphic Encryption

  Figure 1 depicts a round of FL training with federated averaging with homomorphic encryption. During this round, each client runs several epochs of minibatch stochastic gradient descent minimizing a local loss function. Once the clients perform this local training, they encrypt their local models and send them to the central server. The latter will perform a weighted averaging of the updated local models to obtain an encrypted updated global model. The current iteration ends by sending the global model to each of the K clients that decrypt it with the sk.

Figure 2 :

 2 Figure 2: Cross-silo federated learning architecture with MPC

Figure 3 :

 3 Figure 3: Maximum validation accuracy during 20 rounds in function of z.

Figure 4 :

 4 Figure 4: Categorical Accuracy during the training with different settings. Global differential privacy is computed with z=0.01, and the epsilon (dashed curve) is estimated for delta=0.01.

Figure 5 :

 5 Figure 5: Validation accuracy and privacy metrics after 20 rounds of training in function of parameter z.

Table 1 :

 1 Computation time (t c ) in s, encryption time (t e ) in s per client, and bandwidth overhead (B w ) in KB per client for different sets of precision (p) in bits, number of weights per client (N w ), number of training dataset size per client n k and number of clients (K).

	p	N w	n k	K	t e	t c	B w
		197	100	3 0.04 0.05 2.5
	8	1000	1000	5 0.13 0.06 14
		10000 10000 10 1.2 0.26 153
		197	100	3 0.05 0.06 3.7
	16	1000	1000	5 0.16 0.06 19
		10000 10000 10 1.6 0.32 200
		197	100	3 0.08 0.06 6.2
	32	1000	1000	5 0.24 0.07 28
		10000 10000 10 2.7 0.45 295

Table 2 :

 2 Results with BFV-Computation time (t c ) in ms, encryption time (t e ) in ms per client, and bandwidth overhead (B w ) in MB per client for different sets of precision (p) in bits, number of weights per client (N w ), number of training dataset size per client n k and number of clients (K).

	p	N w	n k	K	t e	t c	B w
		197	100	3	1.091 0.043 0.031
	8	1000	1000	5	1.397 0.056 0.031
		10000 10000 10 15.97 3.216	1.8
		197	100	3	1.349 0.024 0.031
	16	1000	1000	5	1.052	0.04 0.031
		10000 10000 10 14.264 3.422	1.8
		197	100	3	0.973 0.024 0.031
	32	1000	1000	5	3.001 0.091 0.087
		10000 10000 10 14.561 4.47	1.8

Table 3 :

 3 Results with CKKS-Computation time (t c ) in ms, encryption time (t e ) in ms per client, and bandwidth overhead (B w ) in MB per client for different sets of precision (p) in bits, number of weights per client (N w ), number of training dataset size per client n k and number of clients (K).

	p	N w	n k	K	t e	t c	B w
		197	100	3	2.015 0.051 0.081
	8	1000	1000	5	2.089 0.093 0.081
		10000 10000 10 21.412 1.733 0.973
		197	100	3	2.021 0.051 0.081
	16	1000	1000	5	2.03	0.086 0.081
		10000 10000 10 22.443 1.735 0.973
		197	100	3	2.025	0.05 0.081
	32	1000	1000	5	4.532 0.192 0.23
		10000 10000 10 23.76 1.696 0.973

Table 4 :

 4 Performance results using general MPC protocols: straightforward aggregation of 5 FL contributions containing 1000 updates each, aggregated by two or three computing parties (CP) acting as aggregator servers. Communication cost is measured for both one CP (2nd row) and all CPs (3rd row) in MB, as well as for both dishonest (DM) and honest majority (HM) configuration. There is a visible performance gap between semi-honest and dishonest model, as well as between honest and dishonest majority.

	Semi-honest, DM, 2CP Semi-honest, HM, 3CP Dishonest, DM, 2CP Dishonest, DM, 3CP
	Time [sec] 69,9	11,8	476,48	518,8
	Send, 1CP 7666,7	83,6	41817,3	83634,6
	Total send 15333,4	250,2	83613,9	250643

Acknowledgements

The research leading to these results has received funding from the European Union's Preparatory Action on Defence Research (PADR-FDDT-OPEN-03-2019). This paper reflects only the authors' views and the Commission is not liable for any use that may be made of the information contained therein.