Vincent Brault
email: vincent.brault@univ-grenoble-alpes.fr

Antoine Channarond
email: antoine.channarond@univ-rouen.fr

Fast and Consistent Algorithm for the Latent Block Model

Keywords: Latent Block Model, Largest Gaps Algorithm, Model Selection, Data analysis Mathematics Subject Classication

The latent block model is used to simultaneously rank the rows and columns of a matrix to reveal a block structure. The algorithms used for estimation are often time consuming. However, recent work shows that the log-likelihood ratios are equivalent under the complete and observed (with unknown labels) models and the groups posterior distribution to converge as the size of the data increases to a Dirac mass located at the actual groups conguration. Based on these observations, the algorithm Largest Gaps is proposed in this paper to perform clustering using only the marginals of the matrix, when the number of blocks is very small with respect to the size of the whole matrix in the case of binary data. In addition, a model selection method is incorporated with a proof of its consistency. Thus, this paper shows that studying simplistic congurations (few blocks compared to the size of the matrix or very contrasting blocks) with complex algorithms is useless since the marginals already give very good parameter and classication estimates.

Introduction

Block clustering methods aim at clustering rows and columns of a matrix simultaneously to form homogeneous blocks. There are a lot of applications of this method: genomics [START_REF] Hedenfalk | Gene-expression proles in hereditary breast cancer[END_REF][START_REF] Jagalur | Analyzing in situ gene expression in the mouse brain with image registration, feature extraction and block clustering[END_REF], recommender systems [START_REF] Bennett | The netix prize[END_REF][START_REF] Shan | Bayesian co-clustering[END_REF], archeology [START_REF] Govaert | Classication croisée[END_REF], sociology [START_REF] Hartigan | Clustering Algorithms[END_REF][START_REF] Keribin | Estimation and selection for the latent block model on categorical data[END_REF][START_REF] Wyse | Block clustering with collapsed latent block models[END_REF] or network [START_REF] Barbillon | Stochastic block models for multiplex networks: an application to a multilevel network of researchers[END_REF][START_REF] Bouveyron | The stochastic topic block model for the clustering of vertices in networks with textual edges[END_REF] for example. Among the methods proposed to solve this question, the Latent Block Model (LBM) [START_REF] Govaert | Clustering with block mixture models[END_REF] provides a chessboard structure induced by the classication of the rows and the classication of the columns. In this model, it is assumed that a sample of n individuals is collected, which contains the observation of d binary variables of the same nature. Saying that the binary variables are of the same nature means that it is possible to encode them in the same (and natural) way. This assumption is needed to ensure that decomposing the dataset in a block structure makes sense and [START_REF] Laclau | Co-clustering through optimal transport[END_REF] shows the equivalence between searching for a classication in the LBM and solving an optimal transport problem. In the case where the rows and columns represent the same individuals, we generally speak instead of Stochastic Block Model [START_REF] Barbillon | Stochastic block models for multiplex networks: an application to a multilevel network of researchers[END_REF][START_REF] Bouveyron | The stochastic topic block model for the clustering of vertices in networks with textual edges[END_REF][START_REF] Tabouy | Variational inference for stochastic block models from sampled data[END_REF]. However, in the case where the matrices are not symmetric, [START_REF] Keribin | Cluster or co-cluster the nodes of oriented graphs?[END_REF] shows that it may be more interesting to use the LBM model rather than the SBM model. Although this model is very similar, we will not discuss it in this article.

Given the number of blocks and in order to estimate the parameters, [START_REF] Govaert | Clustering with block mixture models[END_REF] suggest using a variational algorithm, [START_REF] Keribin | Model selection for the binary latent block model[END_REF] propose an adaptation of the Stochastic Expectation Maximisation introduced by [START_REF] Celeux | On Stochastic Versions of the EM Algorithm[END_REF] in the mixture case, [START_REF] Keribin | Estimation and selection for the latent block model on categorical data[END_REF] studied a Bayesian version of these two algorithms and [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF] propose a Bayesian algorithm including the estimation of the number of blocks. However, each of these iterative algorithms has a complexity at least O (ndN Block N Iter ) where N Block is the number of blocks and N Iter is the number of iterations necessary for the convergence of the algorithm. Moreover, the procedures are often associated with a model selection criterion requiring the computation of maximum likelihood estimators for all combinations of the expected number of blocks. However, the theoretical results obtained show that the distributions of the estimators (see [START_REF] Brault | Consistency and asymptotic normality of latent block model estimators[END_REF][START_REF] Celisse | Consistency of maximum-likelihood and variational estimators in the stochastic block model[END_REF][START_REF] Mariadassou | Convergence of the groups posterior distribution in latent or stochastic block models[END_REF]) are asymptotically trivial: the log-likelihood ratios are equivalent under the complete and observed (with unknown labels) models and the groups posterior distribution to converge as the size of the data increases to a Dirac mass located at the actual groups conguration.

In this article, we propose an adaptation of the Largest Gaps (LG) algorithm introduced by [8] in the Stochastic Block Model with a low complexity O(nd) (Section 3) and present the conditions to obtain asymptotically a good estimation (Section 4). We then prove that it provides a consistent procedure for all inference tasks inherent in LBM: unsupervised classication and estimation of the parameters and a selection procedure for the number of blocks (the last unknown theoretical point) is also proposed and shown to be consistent (Section 5). For ease of reading, the proofs of the results are postponed to the appendices (Sections A, B and C). These theoretical results are also illustrated on simulated data (Section 6).

By proving the consistency of the LG algorithm, some features of the asymptotic regime of the LBM will be highlighted, in particular, the concentration of some marginal distributions of the model. The secondary objective of the paper is to discuss as a conclusion (Section 7) the consequence of this when tting the LBM to large matrices; in particular, this leads to take a step back regarding the relevance of the model and the estimates when the number of blocks is small with respect to the number of cells of the matrix.

Notations and model

The binary Latent Block Model (LBM) is as follows. Let x = (x ij ) i=1,...,n;j=1,...,d be the data matrix where x ij ∈ {0, 1}; observation of a random variable X. It is assumed that there exists a partition into g row clusters z = (z ik ) i=1,...,n;k=1,...,g and a partition into m column clusters w = (w jℓ ) j=1,...,d;ℓ=1,...,m . The z ik s (resp. w jℓ s ) are binary indicators of row i (resp. column j) belonging to row cluster k (resp. column cluster ℓ), such that the random variables X ij are independent conditionally on z and w with parametric density φ(•; α kℓ ) z ik w jℓ , where α kℓ is the parameter of the conditional density of the data given Z ik = 1 and W jℓ = 1. As the binary case is studied, the density is dened for each x ∈ {0, 1} and α ∈ [0; 1] by

φ(x; α) = α x (1 -α) 1-x .
Thus, the density of X conditionally on z and w is dened for each x by

f (x|z, w; α) = n i=1 d j=1 g k=1 m ℓ=1 φ(x ij ; α kℓ ) z ik w jℓ =: i,j,k,ℓ φ(x ij ; α kℓ ) z ik w jℓ
where α = (α kℓ ) k=1,...,g;ℓ=1,...,m .

Moreover, it is assumed that the row and column labels z and w are the observations of two independent variables Z and W: p(z, w; π, ρ) = p(z; π)p(w; ρ) with p(z;

π) = i,k π z ik k = k π z +k k and p(w; ρ) = j,ℓ ρ w jℓ ℓ = ℓ ρ w +ℓ ℓ
, where (π k = P(Z ik = 1), k = 1, . . . , g) and (ρ ℓ = P(W jℓ = 1), ℓ = 1, . . . , m) are the mixing proportions and z +k = n i=1 z ik (resp. w +ℓ = d j=1 w jℓ ) represents the number of rows (resp. columns) in the class k (resp. ℓ). Hence, the density of X is dened for every x by

f (x; θ) = (z,w)∈Z×W p(z; π)p(w; ρ)f (x|z, w; α),
where Z and W denoting the sets of all possible row labels z and column labels w, and θ = (π, ρ, α), with π = (π 1 , . . . , π g ) and ρ = (ρ 1 , . . . , ρ m ). The density of X can be written for every x as

f (x; θ) = (z,w)∈Z×W k π z +k k ℓ ρ w +ℓ ℓ i,j,k,ℓ φ(x ij ; α kℓ ) z ik w jℓ .
To estimate both the classication and the parameters, many algorithms exist (for example [START_REF] Govaert | Clustering with block mixture models[END_REF], [START_REF] Keribin | Estimation and selection for the latent block model on categorical data[END_REF] or [START_REF] Wyse | Block clustering with collapsed latent block models[END_REF]) but each of these algorithms has a complexity larger than O (ndgmN Iter ) where N Iter is the number of iterations of the algorithm. This makes their use on large matrices dicult.

In the Stochastic Block Model (SBM), rows and columns are associated with the same individuals, which allows to represent a graph, whereas LBM allows to represent bipartite graphs. [8] suggested a fast algorithm, called LG for Lagest Gaps, based on a marginal of the matrix x, the degrees.

Algorithm Largest Gaps

Before the introduction of the classication algorithm Largest Gaps (LG), let us recall the main idea, inspired by [8].

Main ideas

Conditionally on Z ik = 1, meaning that row i is in class k, the probability that the variable X ij equals 1 for any j ∈ {1, . . . , d} is

P (X ij = 1|Z ik = 1) = m ℓ=1 P (X ij = 1|Z ik = 1, W jℓ = 1) P (W jℓ = 1|Z ik = 1) = m ℓ=1 P (X ij = 1|Z ik = 1, W jℓ = 1) P (W jℓ = 1) = m ℓ=1 α kℓ ρ ℓ =: τ k . (1)
In particular, conditionally on Z ik = 1, variables of row i are independent and identically distributed Bernoulli variables with parameter τ k and the sum of the cells of any row i, denoted by X i+ , is hence a binomial distributed variable Bin (d, τ k ).

As a consequence of the subgaussian concentration property of binomial distributions, variables X i = X i+ d fastly concentrates around the mean associated with its own class when d tend to innity. The point is here: if moreover these means τ 1 , . . . , τ g are pairwise distinct, the set of the variables {X i ; i = 1, . . . , n} asymptotically splits into clusters, separated by large gaps, and which exactly correspond to the clusters of the model, for d large enough. In the whole paper, τ 1 , . . . , τ g will be assumed to be pairwise distinct (Assumption (I), see discussion in Section 4).

The middle right picture of Figure 1 shows the histogram of the variable set {X i ; i = 1, . . . , n} for a matrix simulated under LBM with ve clusters. The ve clusters of rows can be seen, as well as the four large gaps which separate them. The middle left picture of Figure 1 is a representation of the vector (X i sorted in ascending order) and the bottom left picture of Figure 1, the size of the gaps between two consecutive sorted values.

To classify the rows, the idea is to identify the gaps between the clusters, which are expected to be asymptotically larger. Indeed, the internal gaps (those between two rows of the same cluster) vanish when d tends to innity due to concentration, while the external gaps (those between two rows of distinct clusters) do not, since the means are assumed to be distinct. Symmetrically the same holds for the columns.

There are several strategies to identify the large gaps. In their article, [8] assume that the number Q of clusters (the same for the rows and the columns) is known and partition the population into Q clusters by nding the Q -1 largest gaps. In order to choose Q, a model selection procedure can be rstly done separately, before the classication. The strategy chosen in this article consists in thresholding the gaps, to distinguish outer gaps from inner ones. It advantageously yields both the clusters and the numbers of clusters in only one pass.

The choice of thresholds is critical. Let's comment on the beautiful conguration of Figure 1, where the asymptotic regime has been reached, to present the key issues. If the thresholds are too large (greater than 0.06 on that example), the thresholding step will select only some of the four outer deviations and will not distinguish the ve clusters from each other. Conversely, if they are too small (less than 0.01), the thresholding step will select all four outer gaps, but also some undesirable inner gaps, and the algorithm will wrongly separate some clusters.

Algorithm

The algorithm Largest Gaps1 is given in Algorithm 1 and an illustration is provided in Figure 1. The principle is to calculate the means of the values of the cells of each row and each column, to order them and to calculate the dierences between two consecutive values. Once this step is done, the clusters are formed assuming that each dierence greater than the threshold implies a change of cluster. In the sequel, the estimators provided by the algorithm are denoted by Z, W and θ.

Estimator of θ ⋆ . For the rest of the article, θ ⋆ = (π ⋆ , ρ ⋆ , α ⋆ ) represents the parameters to be estimated (i.e. the parameters that were used to simulate the data). In the algorithm 1, the estimator θ of θ ⋆ is based on Z and W. π k (resp. ρ ℓ ) is the proportion of class k (resp. ℓ) in the partition z (resp. w) and the estimator α is for all (k, ℓ) ∈ {1, . . . , g} × {1, . . . , m}:

α kℓ = n i=1 d j=1 z ik w jℓ x ij z +k w +ℓ .
Complexity of the algorithm. On algorithm 1, the complexities of each step are added in blue at the end of the line. In the end, the LG algorithm has a complexity of O (max (nd [ g + m] , n log n, d log d)). As will be seen in the section 5.3, log n is required to be much smaller than d and log d much smaller than n. In this case, the complexity is

O (nd [ g + m]). Moreover, n i=2 G i = X (n) -X (1)
being smaller than 1 and for all k ∈ {1, . . . , g -1}, G (i k ) being greater than S g then, in the worst case, g is smaller than 1/S g + 1.

As a conclusion, the complexity is O (nd [1/S g + 1/S m ]) and, if the classication only is processed, the complexity is O (nd). where for all j ∈ {2, . . . , d}, G j = X (j) -X (j-1) . Bottom-right: reorganized matrix.

4 Identiability of the model: assumptions on the model parameters In this section, the conditions for the consistency of the estimators are explained. For the rest of the article, g ⋆ and m ⋆ represent the true number of clusters in rows and columns and, as mentioned in the section 3.2, θ ⋆ = (π ⋆ , ρ ⋆ , α ⋆ ) corresponds to the true parameters. Moreover, z ⋆ and w ⋆ are the unobserved partition matrices, resulting from the laws M (1; π ⋆ ) and M (1; ρ ⋆ ), used to simulate the data.

Computation of α = Z T x W/ π k ( ρ ℓ ) T × nd. // O (nd [1/Sg + 1/Sm])
Output: Numbers of clusters g and m, matrices Z and W and parameter θ.

Algorithm 1: Algorithm Largest Gaps. The complexity of each step is highlighted in blue.

Notations 1 Key model parameters:

Let us dene π ⋆ min and ρ ⋆ min the smallest class proportions:

π ⋆ min = min 1≤k≤g ⋆ π ⋆ k and ρ ⋆ min = min 1≤ℓ≤m ⋆ ρ ⋆ ℓ
and the smallest distance between any two conditional expectations of the normalized degrees (called model smallest gaps):

δ π ⋆ = min 1≤k̸ =k ′ ≤g ⋆ τ ⋆ k -τ ⋆ k ′ and δ ρ ⋆ = min 1≤ℓ̸ =ℓ ′ ≤m ⋆ ξ ⋆ ℓ -ξ ⋆ ℓ ′
where τ ⋆ = ρ ⋆ α ⋆T and ξ ⋆ = π ⋆ α ⋆ are the proportions of the binomial distributions dened in Equation [START_REF] Barbillon | Stochastic block models for multiplex networks: an application to a multilevel network of researchers[END_REF].

Cluster identiability by the LG algorithm is possible under the following sucient conditions:

Assumptions 1 (I) The proportion of each row class (respectively column class) is non-negative, and conditional expected degrees of row (resp. column) clusters (τ k ) 1≤k≤g ⋆ (resp. (ρ ℓ ) 1≤ℓ≤m ⋆ ) are all distinct, which respectively amounts to:

π ⋆ min ρ ⋆ min > 0 and δ π ⋆ δ ρ ⋆ > 0. (I)
These assumptions will be always made in the sequel of this article. They are equivalent to the sucient conditions of identiability of LBM given in [START_REF] Keribin | Estimation and selection for the latent block model on categorical data[END_REF] but stronger than the conditions of [START_REF] Brault | Consistency and asymptotic normality of latent block model estimators[END_REF] or the identiability conditions in the case where (g, m) = (2, 2) in the [START_REF] Keribin | Estimation and selection for the latent block model on categorical data[END_REF].

The rst assumption, π ⋆ min ρ ⋆ min > 0, is classical in mixture models (for example, see [START_REF] Brault | Consistency and asymptotic normality of latent block model estimators[END_REF][START_REF] Iannario | On the identiability of a mixture model for ordinal data[END_REF][START_REF] Maugis | A non asymptotic penalized criterion for gaussian mixture model selection[END_REF]); without this assumption, the probability of having the right number of clusters is null. Thus the mixture model would be degenerated, as clusters with proportion zero would be always empty and the number of actually present clusters would be improper.

The second one, δ π ⋆ δ ρ ⋆ > 0, ensures that not only the distribution of x is identiable, but also that of (X i+ ) 1≤i≤n , which is a marginal distribution of x. This is critical in order to recover all clusters with the LG algorithm, which actually infers the clusters only from the set of variables {X i ; i = 1, . . . , n}. More precisely, it ensures that the distribution of X i+ is a proper mixture with the same number of clusters as the distribution of X ij . For example, Assumption (I) excludes the following typical case (yet identiable according to [START_REF] Keribin | Estimation and selection for the latent block model on categorical data[END_REF]):

g = 2, m = 2, π = (1/2, 1/2) and α = a b b a where a, b ∈]0, 1[. Indeed, it gives τ 1 = τ 2 = a+b
2 , and hence δ π ⋆ = 0. Conditionally on either z i1 = 1 or z i2 = 1, the distribution of X i+ is the same: Bin d, a+b 2 , and therefore the two clusters cannot be distinguished just with the vector (X i+ ) 1≤i≤n . The distribution of X i+ is not a proper mixture with two distinct clusters, it is a simple binomial distribution. Thus with our approach this model would be confused with the model

g = 1, π = 1, α = a+b 2 .
Note that the set of parameters such that δ π ⋆ δ ρ ⋆ = 0 has zero Lebesgue measure.

Consistency

This section presents the main result (Theorem 2), namely the consistency of estimators. Here, consistency means that the numbers of clusters and classications are correct and that the distance between the estimates and the model parameters is smaller at any t > 0 with a probability tending towards one, when the size of the matrix (n, d) tends towards innity. Before stating this theorem, we introduce some notations, especially related to the label switching problem.

Distance on the parameters and the label switching issue

For any two parameters θ = (π, ρ, α) with (g, m) clusters and θ ′ = (π ′ , ρ ′ , α ′ ), with (g ′ , m ′ ) clusters, we dene their distance as follows:

d ∞ θ, θ ′ = θ -θ ′ ∞ if g = g ′ , m = m ′ +∞ otherwise, where ∥•∥ ∞ denotes the norm dened for any y ∈ R g by ∥y∥ ∞ = max 1≤k≤g |y k |.
We assume that two matrices z, z ′ ∈ M n×g ({0, 1}) are equivalent, denoted z ≡ Z z ′ , if there exists a permutation s ∈ S ({1, . . . , g}) such that for all (i, k) ∈ {1, . . . , n} × {1, . . . , g}, z ′ i,s(k) = z ik . By convention, we assume that two matrices with dierent numbers of columns are not equivalent. We introduce the similar notation ≡ W for the matrix w. For all parameter θ = (π, ρ, α) with (g, m) clusters and for all permutions (s, t) ∈ S ({1, . . . , g}) × S ({1, . . . , m}), we denote θ s,t = (π s , ρ t , α s,t ), by:

π s = π s(1) , . . . , π s(g) , ρ t = ρ t(1) , . . . , ρ t(m)
and α s,t = α s(1),t(1) , α s(1),t(2) , . . . , α s(1),t(m) , α s(2),t(1) , . . . , α s(g),t(m) .

Like all mixture models, the LBM is aected by the label switching problem: clusters are dened up to a permutation. The algorithm can therefore nd the right clusters but at a permutation of the labels. This also aects the parameters of the model, because the order of their coordinates depends on the labeling. The comparison of two classications, and two parameter estimates, must then be done carefully: the distance between two parameter estimates must be calculated after permutation of their coordinates, using the permutation transforming the label allocation of the classication algorithm into the original label allocation of the model. Moreover, such a permutation exists and is unique when the classication is correct, i.e. when Z ≡ Z z ⋆ (respectively W ≡ W w ⋆ ). This permutation will thus be noted s Z (resp. t W ) on the event { Z ≡ Z z ⋆ } (resp. W ≡ W w ⋆ ). Thus, the consistency of the parameter estimators amounts to proving that the following quantity disappears in probability when (n, d) tends to innity:

∀t > 0, P d ∞ θ s Z ,t W , θ ⋆ > t Z ≡ Z z ⋆ , W ≡ W w ⋆ -→ n,d→+∞ 0.
Outside of the event { Z ≡ Z z ⋆ } (resp. W ≡ W w ⋆ ), s Z (resp. t W ) will be dened as any arbitrary permutation in S ({1, . . . , g}) (resp. S ({1, . . . , m})), the identity for instance.

A non-asymptotic upper bound

This paragraph presents the main technical tool for obtaining consistency: a nonasymptotic upper bound from which the consistency results will be derived.

Theorem 1 (Concentration inequality) Under identiability assumptions (I), and if S g ∈]0, δ π ⋆ [ and S m ∈]0, δ ρ ⋆ [ for n, d large enough, then we have for all t > 0:

P g ̸ = g ⋆ or m ̸ = m ⋆ or Z ̸ ≡ Z z ⋆ or W ̸ ≡ W w ⋆ or d ∞ θ s Z ,t W , θ ⋆ > t ≤ 2n exp - d 2 min(δ π ⋆ -S g , S g ) 2 + g ⋆ 1 -π ⋆ min n +2d exp - n 2 min(δ ρ ⋆ -S m , S m ) 2 + m ⋆ 1 -ρ ⋆ min d +2g ⋆ m ⋆ 1 -π ⋆ min ρ ⋆ min 1 -e -2t 2 nd + 2g ⋆ e -2nt 2 + 2m ⋆ e -2dt 2 .
The proof (in Appendix A) is made in two steps, emphasizing the originality of the method in comparison with EM-like algorithms: here the classication is completely done rst, and parameters are then estimated afterwards. Thus an upper bound on classications and selection of class numbers will be rst established (Proposition 1), and secondly an upper bound on the parameter estimators, given that both classications and class numbers are right (Proposition 2).

If n and d increase at the same speed, the larger the size, the smaller the bound. We nd the importance that S g and S m are small enough to detect the real jumps but not so small that there are only them. Moreover, the larger π ⋆ min and ρ ⋆ min are, the more likely it is to have at least one representative of each cluster.

In addition to showing the consistency of the LG algorithm's estimators (see section 5.3), the bound of theorem 1 quanties how (very) easy a conguration is: given a conguration (θ, n, d), it is possible to estimate the probability that the algorithm Largest Gaps nds the right conguration. In the case of a conguration obtained by another algorithm, this bound can be calculated and if it is (very) small, it means that a study of the marginals would have been sucient; if the Largest Gaps algorithm does not nd the same results at all, a doubt can be cast on the results (estimation or relevance of the use of a LBM) by precaution.

Remark 1 In the proof, the Hoeding inequality is used. In the case where the parameters are very close to 0 or 1, Bernstein's inequality allows to improve the bound. This is particularly interesting if the parameters evolve with n and d (see Section 5.4).

Consistency of the estimators

The following theorem provides sucient assumptions on the threshold sequences

(S n,d g , S n,d m ) n,d
to ensure the consistency of the inference method based on the LG algorithm when n, d tend to innity. Note that this result is therefore asymptotic only, and it does not provide any guarantee for xed (n, d). Nevertheless the rates can be used as a suggestion to choose the thresholds, as even though it must be carefully used, because it is anyway impossible to know whether the asymptotic regime is reached or not.

Theorem 2 (Consistency) Under identiability assumptions (I) and the following assumptions:

S n,d g -→ n,d→+∞ 0, S n,d m -→ n,d→+∞ 0, lim n,d→+∞ S n,d g d log n > √ 2 and lim n,d→+∞ S n,d m n log d > √ 2, ( TA
)
where lim is the lower limit; then classications, model selection and estimators are consistent, that is, for all t > 0:

P g ̸ = g ⋆ or m ̸ = m ⋆ or Z ̸ ≡ Z z ⋆ or W ̸ ≡ W w ⋆ or d ∞ θ s Z ,t W , θ ⋆ > t -→ n,d→+∞ 0.
The proof of this result is available in Appendix B.

The LG algorithm has two input parameters, the (S g , S m ) thresholds, which must be set correctly to discover all clusters. Recall that the purpose of gap thresholding is to distinguish between external and internal gaps (see comments in section 3.1). First, the thresholds must be smaller than the smallest gaps in the model,

S n,d g < δ π ⋆ and S n,d
m < δ ρ ⋆ , for n and d suciently large; otherwise, some clusters will consist of mixed clusters. Since δ π ⋆ and δ ρ ⋆ are not known a priori, we assume that the thresholds decrease with respect to n and d, to ensure that they are asymptotically small enough. On the other hand, if the threshold sequences decrease too fast, the thresholds will be asymptotically too small and at least one class will be split into several clusters by the algorithm. More technically, the convergence rate given in the theorem guarantees that the upper bound of the theorem 1 tends to 0 which implies consistency. If the sequences disappear faster, coherence is no longer guaranteed.

Consistency can thus be obtained whatever the sequences S n,d g and S n,d m taken provided that they verify the conditions (TA); nevertheless, if we couple with the condition that S n,d g < δ π ⋆ and S n,d m < δ ρ ⋆ of Theorem 1, we see that it is preferable that the thresholds decrease rapidly. For example, the threshold S n,d g = 2 log(n)/d(1+ε) with ε > 0 can be used. Moreover, once a rst conguration has been found, it is possible to estimate δ π ⋆ and δ ρ ⋆ and to re-estimate the partitions with the thresholds δ π ⋆ /2 and δ ρ ⋆ /2 in order to see if a better conguration is obtained. Remark 2 The assumption (TA) of the theorem implies that n/ log d and d/ log n tend to +∞ ; this assumption is found in the results on the consistency of the maximum likelihood estimator (see [START_REF] Brault | Consistency and asymptotic normality of latent block model estimators[END_REF]). Therefore, x is allowed to have an oblong shape. For example, d = n γ with γ > 0 satises the assumption.

Remark 3 The roles of the rows and columns are symmetric: the result remains true if the transpose of X is studied rather than the matrix X.

Consistency when model parameters vary

Finally varying model parameters are also considered in this paragraph. Indeed when n and d increase, it can be reasonable to assume that new clusters arise: in this paragraph g ⋆n,d and m ⋆n,d are hence assumed to be growing to innity when n, d tend to +∞. The consequence of this assumption is the convergence to zero of both the proportions of the smallest clusters and the model smallest gaps: π ⋆ min , ρ ⋆ min , δ π ⋆ , δ ρ ⋆ tend to 0 when n, d tend to innity. For example, since the parameters τ ⋆ = (τ ⋆ k ) 1≤k≤g ⋆ are probabilities, the following inequalities are obtained:

(g ⋆ -1)δ π ⋆ ≤ g ⋆ -1 k=1 (τ ⋆ (k+1) -τ ⋆ (k) ) = τ ⋆ (g ⋆ ) -τ ⋆ (1) ≤ 1 and (m ⋆ -1)δ ρ ⋆ ≤ m ⋆ -1 l=1 (ξ ⋆ (l+1) -ξ ⋆ (l) ) ≤ 1, where τ ⋆ (1) < • • • < τ ⋆ (g ⋆ ) (resp. ξ ⋆ (1) < • • • < ξ ⋆ (m ⋆ ) ) are the (τ ⋆ k ) 1≤k≤g ⋆ (resp.
(ξ ⋆ ℓ ) 1≤l≤m ⋆ ) sorted in increasing order. In the other side, g ⋆ is bounded by 1/π ⋆ min by the following inequality:

g ⋆ π ⋆ min ≤ g ⋆ k=1 π ⋆ k = 1.
In particular, it is possible that one of these parameters tend to zero while the number of groups does not change (for example, in the case of increasingly sparse matrices). In this framework, admissible convergence rates of the model parameters are provided, as well as the corresponding admissible convergence rates of the thresholds. It thus tells how robust the consistency is. The proof of this result is available in Appendix C. Theorem 3 shows that the estimates remain consistent even with sparse matrices. On the other hand, for an optimal result and if a number of non-zero cells is xed, the estimation will be all the easier if the latter are well distributed between the classes, for example, in a staircase shape.

+∞

where N ⋆ nd block = g ⋆n,d m ⋆n,d is the number of blocks. In particular, the number of blocks can therefore increase with n and d but not too quickly.

Simulations

In this section, simulations to test the quality of the bounds obtained in the previous theorems and to compare the computing times with classical procedures are presented2 .

Estimation of the number of clusters

We use an experimental design to illustrate the results of Theorem 2. As the number of row (resp column) clusters is the basis of the other estimations, this is the only parameter studied in this section. The experimental design is dened with g ⋆ = 5 and m ⋆ = 4 and the following parameters

α ⋆ =       ε ε ε ε 1 -ε ε ε ε 1 -ε 1 -ε ε ε 1 -ε 1 -ε 1 -ε ε 1 -ε 1 -ε 1 -ε 1 -ε      
with ε ∈ {0.05, 0.1, 0.15, 0.2, 0.25, 0.3}. For ε = 0.05, the values of the Bernoulli parameters are 0.05 (resp. 0.95) and the associated blocks contain essentially white (resp. black) cases. At the opposite, for ε = 0.3, the block are more homogeneous and more dicult to distinguish. For the class proportions, we suppose two possibilities Balanced proportions: The number of rows n and the number of columns d vary between 20 and 4000 by step 20 and for each conguration, 100 matrices were simulated. For the choice of the thresholds S g , we studied four cases: 2. In practice, there is mostly no reason why the parameter δ π ⋆ could be known.

In this case, Theorem 2 claims that we need to use a varying threshold S n,d g instead, such that S n,d g tends to 0 but not too fast (slower than 2 log n/d). If the threshold decreases too slowly, it may be larger than δ π ⋆ and the smallest gaps could be undetected. On the opposite, if the threshold decreases too fast, we may detect too many gaps. In the simulation, we studied three possibilities: (a) Faster threshold: Figures 2 and3 show the proportions of bad estimations of g ⋆ following the parameters ε (in rows) as function of the number of rows n, the numbers of columns d and the thresholds used (columns). For each gure, the number of columns (d) increases on the x-axis and the number of rows (n) increases on the y-axis. The red color corresponds to a bad estimate while the blue color corresponds to a good estimate (the more blue it is, the better).

S n,d
As expected, it appears that the best threshold is the oracle S 1 = δ π ⋆ /2 but this threshold cannot be used in practice because δ π ⋆ is unknown. For the scaled thresholds, S n,d 2 = 2 log n/d 1 + 10 -10 is the best. We can see that the larger the number of rows n is, the worse the estimation is and the larger the number of columns d is, the better the estimation is. In the case of n = d, the quality of the estimation increases with n. In particular, the model selection can be generalized for the case of [8] and the results would be similar.

π ⋆ min has a weak eect because it is rare to have an empty class but the eect of δ π ⋆ seems to be greater.

Estimation of the parameters

To illustrate the convergence of the the estimation θ to the true parameter θ ⋆ , the same experimental design is chosen with the number of rows n and the number of columns d vary between 40 and 4000 by step 40 and, for each conguration, 100 matrices were simulated. To estimate the quality of the estimation, the distance d ∞ introduced in Section 5.1 is calculated. As the distance equals +∞ when the number of clusters of θ is dierent of its of θ ⋆ , we chose to represent the mean of the nite values with the size of the point corresponding at the number of nite values (see Figure 4 for the balanced case and Figure 7 for the arithmetic case).

We observe that the error decreases with the number of observations and the results are identical when the numbers of clusters are correctly estimated. The optimization of the LG algorithm thus depends on the choice of the threshold. In particular, we nd that the oracle threshold (S 1 ) nds the right number of clusters faster than the evolutionary thresholds. In particular,the threshold S n,d does not appear on the graphic because it underestimates the number of blocks.

Comparison of computing times

To conclude the part of the simulations, the computation times are estimated and compared with a classical procedure. For that, the plan presented previously is ; the algorithm variational Bayes with the hyperparameters [START_REF] Bouveyron | The stochastic topic block model for the clustering of vertices in networks with textual edges[END_REF][START_REF] Bouveyron | The stochastic topic block model for the clustering of vertices in networks with textual edges[END_REF] as proposed by [START_REF] Keribin | Estimation and selection for the latent block model on categorical data[END_REF] knowing the true number of parameters g ⋆ and m ⋆ (named Simple VBayes for the next); as the number of blocks is usually unknown, the combination of the algorithm variational Bayes and the Integrated Complete-data Likelihood (see [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF][START_REF] Keribin | Estimation and selection for the latent block model on categorical data[END_REF]) is studied with the hyperparameters (4, 4) on a 2 × 7 square grid (named VBayes+ICL for the next).

For the implementation, the R package blockcluster (version 4.5.1; see [START_REF] Singh Bhatia | blockcluster: An R package for model-based coclustering[END_REF]) with the function coclusterBinary is used. For each conguration, 20 matrices are simulated and procedures are evaluated on two criteria: the mean computation times with 10 runs; for this, the package microbenchmark (version 1.4.9; see [START_REF] Mersmann | microbenchmark: Accurate Timing Functions[END_REF]) with the function microbenchmark is used and plans are launched on a cluster 3 . the estimation quality of the number of clusters per row after an estimation of each matrix (except for Simple VBayes as the number is assumed be known); results are averaged over the 20 matrices.

The results are displayed in the Figures 5 (balanced case) and 6 (arithmetic case) where the computation times (in seconds with logarithmic scale; on the top) and the quality of the estimations (averaged for each ε on over the 20 matrices; on the bottom) have been grouped by matrix size (regardless of the values of ε). A detailed version for each ε is available on the Figures 8,9, 10 and 11 in supple- mentary material. Moreover, the means and standard deviations (in milliseconds) are represented in the Tables 1 and2.

On the quality, the procedure VBayes+ICL estimates better than the other the number of clusters in rows and the dierences between the four threshold are the same than the section 6.1.

For the computing times, the order of Largest Gaps is that of the millisecond while that of VBayes+ICL is of the second (and the minute after n = 300). Moreover, the computing time for the Largest Gaps seems to increase linearly with n (the side length of the square matrix) as stated in the section 3.2 while VBayes+ICL increases much faster than linearly. Moreover, Largest Gaps does not seem to be impacted by the conguration while the convergence time of VBayes+ICL (and even Simple VBayes who knows the right number of clusters) is larger for the case of arithmetic proportions. And nally, the computation time of VBayes+ICL depends on the maximum choice of the number of clusters (xed here quite close to the true result) but unknown in practice while Largest Gaps does not need this kind of parameter and the time is almost independent of the number of blocks.

Conclusion and discussion

First of all, Largest Gaps is a co-clustering algorithm, which has nice theoretical properties: its computational cost is much lower than most known algorithms. And it provides a consistent procedure under the Latent Block Model, for all inference tasks: model selection, classication and estimation of the model parameters. Since the algorithm is simple, the consistency is rather easy to obtain. Note that in this article, only binary matrices have been studied, but the model as well as the method and the proofs can be directly extended for distributions which have the same concentration properties, for example compactly-supported distributions, where the support is known.

As a consistent algorithm, the advantage of the LG algorithm is the simplicity: it is a simple and original way to exploit the concentration of the marginal distributions of the matrix x under LBM. But it lacks robustness: a single outlier in the marginal distribution can aect the whole classication, which is not desirable by the use on real-world datasets. Moreover, the procedure is essentially based on the choice of threshold and simulations have shown that there is still room for improvement in this choice. Many other algorithms, also based on this nice feature, could be used instead and might be better in practice, but harder to analyze from a theoretical point of view. For example, tting a binomial mixture model on the variables X (1) , . . . X (n) with an EM-algorithm, could be fruitful as well.

The contribution of this article actually goes beyond the LG algorithm. It shows special features of the latent structure of the LBM, and their consequences. In particular, when the asymptotic regime of the model is reached, the latent structure is almost obvious, and moreover it can be pick either one but not the two from a summary of the data. Indeed even basic algorithms like LG can retrieve the latent clusters from variables X (1) , . . . X (n) (sums of the rows and of the columns of the matrix x), whereas most known classication algorithms are unusable, because of their complexity.

A consequence of this remark is that the LBM should be used sparingly on very large real data sets: if the number of clusters requested is too small and if the blocks are suciently contrasted then the marginals should highlight the clusters.

Table 1 Computing times for the dierent procedures in columns (LG with the four threshold, VBayes with the true number of clusters and VBayes coupled with the ICL criterion on a 2 × 7 square grid) following the number of rows and columns (in rows) for the balanced case. Each cell represents the average (in milliseconds) and the standard deviation (in parenthesis) over the 1200 simulations (200 per ε). If this is not the case, the use of the LBM should be questioned. The bound of the theorem 1 would give an estimate of the quality of the use of this model.

S 1 S n,d 2 S n,d 3 S n,
Finally, it appears in the simulations that the estimate of the number of clusters is underestimated. Moreover, the distribution of marginals from real-world data rarely has such obvious deviations as assumed by the model asymptotic. To overcome this problem, it would be interesting to estimate the row clusters with a mixture model on the variables (X (1) , . . . , X (n) ); this will be the subject of future work. A Proof of Theorem 1

In this appendix, we present demonstrations of the concentration inequality. We rst estimate the probability of having the right number of clusters and then the right classication knowing that we have the right number of clusters. Finally, we evaluate the quality of the parameter estimates.

First of all, note that { Z ≡ Z z ⋆ } ⊂ { g = g ⋆ } and { W ≡ W w ⋆ } ⊂ { m = m ⋆ }, hence :

P g ̸ = g ⋆ or m ̸ = m ⋆ or Z ̸ ≡ Z z ⋆ or W ̸ ≡ W w ⋆ or d ∞ θ s Z ,t W , θ ⋆ > t = P Z ̸ ≡ Z z ⋆ or W ̸ ≡ W w ⋆ or d ∞ θ s Z ,t W , θ ⋆ > t = P Z ̸ ≡ Z z ⋆ or W ̸ ≡ W w ⋆ + P d ∞ θ s Z ,t W , θ ⋆ > t Z ̸ ≡ Z z ⋆ or W ̸ ≡ W w ⋆ = P Z ̸ ≡ Z z ⋆ or W ̸ ≡ W w ⋆ + P d ∞ θ s Z ,t W , θ ⋆ > t, Z ≡ Z z ⋆ , W ≡ W w ⋆ ≤ P Z ̸ ≡ Z z ⋆ + P W ̸ ≡ W w ⋆ +P d ∞ θ s Z ,t W , θ ⋆ > t Z ≡ Z z ⋆ , W ≡ W w ⋆ P Z ≡ Z z ⋆ , W ≡ W w ⋆ ≤ P Z ̸ ≡ Z z ⋆ + P W ̸ ≡ W w ⋆ + P d ∞ θ s Z ,t W , θ ⋆ > t Z ≡ Z z ⋆ , W ≡ W w ⋆
To complete the proof, we then need to bound from above the terms of this inequality. The two rst terms are bounded using Proposition 1 in Appendix A.1, and the last term is bounded with Proposition 2 in Appendix A.2.

A.1 Concentration inequality on Z

Let us rst dene the following events.

There is at least one individual in each row class, denoted by

A g ⋆ = g ⋆ k=1 Z ⋆ +k ̸ = 0 .
Denoting D the random variable equal to the maximal distance between X i and the center of the class of row i:

D = max 1≤k≤g ⋆ sup 1≤i≤n with z ⋆ i,k =1 X i -τ k ,
we also dene:

A Sg = {2D < Sg < δ π ⋆ -2D} and A id = A g ⋆ ∩ A Sg .
Lemma 1 (Interesting event)

A id ⊂ { g = g ⋆ } ∩ { Z ≡ Z z ⋆ }
Proof On the event A Sg , for any two rows i ̸ = i ′ ∈ {1, . . . , n}, we have two possibilities:

Either the rows i and i ′ are in the same class k, and then on A Sg , we have:

X i -X i ′ ≤ X i -τ k + X i ′ -τ k ≤ 2D < Sg.
Or row i is in the class k and row i ′ in the class k ′ ̸ = k, and on the event A Sg , we have:

X i -X i ′ = X i -τ k ′ -X i ′ -τ k ′ ≥ X i -τ k ′ -X i ′ -τ k ′ ≥ X i -τ k ′ -D ≥ |τ k -τ k ′ | -X i -τ k -D ≥ δ π ⋆ -2D > Sg.
Therefore, G i = X (i) -X (i-1) is less than Sg if and only if both rows (i -1) and (i) are in the same class. On A Sg , the algorithm hence nds the true classication. Moreover, on A g ⋆ , there is at least one row in each class, then the algorithm nds the true number of classes. As a conclusion, on A id , both g = g ⋆ and Z ≡ Z z ⋆ are satised.

Lemma 2 (Concentration inequality for

A id ) Under Assumption (I) and if Sg ∈]0, δ π ⋆ [ and Sm ∈ ]0, δ ρ ⋆ [ for n, d large enough: P A id ≤ 2n exp - d 2 min(δ π ⋆ -Sg, Sg) 2 + g ⋆ (1 -π ⋆ min ) n
where A id is the complementary of the event A id .

Proof Using an union bound, we rst obtain:

P A id ≤ P A g ⋆ + P A Sg
Now we bound from above each of these terms. Again with an union bound:

P A g ⋆ = P   g ⋆ k=1 Z ⋆ +k ̸ = 0   ≤ g ⋆ k=1 P Z ⋆ +k ̸ = 0 ≤ g ⋆ k=1 P Z ⋆ +k = 0 ≤ g ⋆ k=1 n i=1 P Z ⋆ i,k = 0 ≤ g ⋆ k=1 n i=1 (1 -π ⋆ k ) ≤ g ⋆ k=1 n i=1 (1 -π ⋆ min ) ≤ g ⋆ (1 -π ⋆ min ) n ,
which gives the upper bound of the rst term. Secondly:

A Sg = {2D < Sg < δ π ⋆ -2D} = {2D < Sg, 2D < δ π ⋆ -Sg} = D < 1 2 min(δ π ⋆ -Sg, Sg) .
Denoting t = min(δ π ⋆ -Sg, Sg),

P A Sg = P D ≥ t 2 = P    g ⋆ k=1 i|z ⋆ i,k =1 X i -τ k ≥ t 2    ≤ g ⋆ k=1 i|z ⋆ i,k =1 P X i -τ k ≥ t 2 .
Moreover for all i ∈ {1, . . . , n}, given z ⋆ i,k = 1, X i+ has a binomial distribution Bin (d, τ k ). The concentration properties of this distribution are then exploited through the Hoeding inequality:

P X i -τ k ≥ t 2 = P |X i+ -dτ k | ≥ dt 2 ≤ 2e -1 2 dt 2 .
And as a conclusion, the bound of the second term is:

P A Sg ≤ g ⋆ k=1 i|z ⋆ i,k =1 2e -1 2 dt 2 = 2ne -1 2 dt 2 .
With these two lemmas, the following proposition is obtained

P Z ̸ ≡ Z z ⋆ ≤ 2n exp - d 2 min(δ π ⋆ -Sg, Sg) 2 + g ⋆ (1 -π ⋆ min ) n . P W ̸ ≡ W w ⋆ ≤ 2d exp - n 2 min(δ ρ ⋆ -Sm, Sm) 2 + m ⋆ (1 -ρ ⋆ min ) d .
Proof Lemma 1 tells that whenever the event A id is satised, then both true number of row classes and their true classication are obtained. Lemma 2 provides an upper bound of P A id .

From these lemmas, it is directly deduced that:

P Z ̸ ≡ Z z ⋆ ≤ P A id ≤ 2n exp - d 2 min(δ π ⋆ -Sg, Sg) 2 + g ⋆ (1 -π ⋆ min ) n ,
which is Proposition 1.

A.2 Concentration inequality on

d ∞ θ s Z ,t W , θ ⋆ > t
In this part, the proof of the inequality on

d ∞ θ s Z ,t W , θ ⋆ > t is detailled.
Proposition 2 For all t > 0, we have:

P d ∞ θ s Z ,t W , θ ⋆ > t Z ≡ Z z ⋆ , W ≡ W w ⋆ ≤ 2g ⋆ m ⋆ 1 -π ⋆ min ρ ⋆ min 1 -e -2t 2 nd + 2g ⋆ e -2nt 2 + 2m ⋆ e -2dt 2
The proof consists in obtaining three bounds: one for each parameter. The inequalities on π and ρ are an application of the Hoeding inequality and are similar to [8] for the row class proportions. To obtain the inequality for α, it is necessary to study the conditional probability, given the true partition (z ⋆ , w ⋆ ). Apart from the problem of two asymptotic behaviors, the proof is similar to [8].

In the sequel, and for ease of reading, we remove the superscripts s Z and t W . Therefore, for all t > 0:

P d ∞ θ, θ ⋆ > t Z ≡ Z z ⋆ , W ≡ W w ⋆ = P max ∥ π -π ⋆ ∥ ∞ , ∥ ρ -ρ ⋆ ∥ ∞ , ∥ α -α ⋆ ∥ ∞ > t Z ≡ Z z ⋆ , W ≡ W w ⋆ ≤ P ∥ π -π ⋆ ∥ ∞ > t Z ≡ Z z ⋆ , W ≡ W w ⋆ + P ∥ ρ -ρ ⋆ ∥ ∞ > t Z ≡ Z z ⋆ , W ≡ W w ⋆ +P ∥ α -α ⋆ ∥ ∞ > t Z ≡ Z z ⋆ , W ≡ W w ⋆ ≤ g ⋆ k=1 P | π k -π ⋆ k | > t Z ≡ Z z ⋆ , W ≡ W w ⋆ + m ⋆ ℓ=1 P | ρ ℓ -ρ ⋆ ℓ | > t Z ≡ Z z ⋆ , W ≡ W w ⋆ + g ⋆ k=1 m ⋆ ℓ=1 P | α kℓ -α ⋆ kℓ | > t Z ≡ Z z ⋆ , W ≡ W w ⋆ .
The upper bounds of the rst and second terms are the same as [8]; only the last term is dierent. For α kℓ , rst note that when Z ≡ Z z ⋆ and W ≡ W w ⋆ α kℓ = α kℓ = 1 z ⋆ +k w ⋆ +ℓ (i,j) z ⋆ i,k w ⋆ j,ℓ =1

X ij and given (z ⋆ , w ⋆ ), the Hoeding inequality gives for all t > 0:

P | α kℓ -α ⋆ kℓ | > t Z ≡ Z z ⋆ , W ≡ W w ⋆ = P | α kℓ -α ⋆ kℓ | > t Z ≡ Z z ⋆ , W ≡ W w ⋆ ≤ P (| α kℓ -α ⋆ kℓ | > t) ≤ E Z ⋆ ,W ⋆ [ P (| α kℓ -α ⋆ kℓ | > t| Z ⋆ , W ⋆ )] ≤ E Z ⋆ ,W ⋆ 2e -2Z ⋆ +k W ⋆ +ℓ t 2 .
But, as

Z ⋆ +k = n i=1 Z ⋆ i,k and W ⋆ +ℓ = d j=1 W ⋆
jℓ and the variables are independents, the expectation is:

E Z ⋆ ,W ⋆ 2e -2Z ⋆ +k W ⋆ +ℓ t 2 = 2E Z ⋆ ,W ⋆   exp   -2 n i=1 Z ⋆ i,k d j=1 W ⋆ jℓ t 2     = 2E Z ⋆ ,W ⋆   n i=1 exp   -2Z ⋆ i,k d j=1 W ⋆ jℓ t 2     = 2 n i=1 E Z ⋆ ,W ⋆   exp   -2Z ⋆ i,k d j=1 W ⋆ jℓ t 2     = 2 n i=1 d j=1 E Z ⋆ ,W ⋆ exp -2Z ⋆ i,k W ⋆ jℓ t 2
Since, for all i ∈ {1, . . . , n}, j ∈ {1, . . . , d}, k ∈ {1, . . . , g ⋆ } and ℓ ∈ {1, . . . , m ⋆ }, the variable

Z ⋆ i,k W ⋆
jℓ is a Bernoulli of parameter π ⋆ k ρ ⋆ ℓ , then:

E Z ⋆ ,W ⋆ exp -2z ⋆ i,k w ⋆ j,ℓ t 2 = π ⋆ k ρ ⋆ ℓ exp -2 × 1 × t 2 + (1 -π ⋆ k ρ ⋆ ℓ ) exp -2 × 0 × t 2 = π ⋆ k ρ ⋆ ℓ exp -2t 2 + 1 -π ⋆ k ρ ⋆ ℓ = 1 -π ⋆ k ρ ⋆ ℓ 1 -e -2t 2 .
Finally, the inequality is:

P | α kℓ -α ⋆ kℓ | > t Z ≡ Z z ⋆ , W ≡ W w ⋆ ≤ 1 -π ⋆ k ρ ⋆ ℓ 1 -e -2t 2 nd ≤ 1 -π ⋆ min ρ ⋆ min 1 -e -2t 2 nd
and the result is obtained by the sum on each cluster.

B Proof of Theorem 2: consistency

The proof is based on Theorem 1, as n → +∞ and d → +∞ and by the assumption (I), the following limits are obtained for every t > 0: For the same reasons, as soon as t is positive, 0 < π ⋆ min ρ ⋆ min 1 -e -2t 2 < 1 and

g ⋆ (1 -π ⋆ min ) n + m ⋆ (1 -ρ ⋆ min ) d -→ n,d→+∞
g ⋆ m ⋆ 1 -π ⋆ min ρ ⋆ min 1 -e -2t 2 nd -→ n,d→+∞ 0.
For the last terms, the assumption (TA) gives that Finally, the assumption (MA.1) implies that there exists a positive constant C > 2 such that for n and d large enough

D Supplementary Material

In this supplementary material, additional gures from the experiments in the section 6 are presented.

Fig. 1

 1 Fig. 1 Top-left: Initial matrix. Top-right: Example of a vector X (1) , . . . , X (d) . Middleleft: representation of the vector X (1) , . . . , X (d) sorted in increasing order. Middle-right: histograms of X (1) , . . . , X (d) . Bottom-left: representation of the vector of gaps (G 2 , . . . , G d )where for all j ∈ {2, . . . , d}, G j = X (j) -X (j-1) . Bottom-right: reorganized matrix.

Theorem 3 ( 1 )

 31 Consistency in sparse case) Under the assumptions of the previous theorem ( (I) and (TA)) and the following additional assumptions : Assumptions on δ n,d π ⋆ and S n,d g (resp. δ n,d ρ ⋆ and S n,d m Assumptions on π ⋆ n,d min and ρ ⋆ n,d min : classications, model selection and estimators are also consistent.

π

  ⋆ = 0.2 0.2 0.2 0.2 0.2 and ρ ⋆ = 0.25 0.25 0.25 0.25 with the following parameters π ⋆ min = 0.2 and δ π ⋆ = 0.25 -0.5ε. Arithmetic proportions: π ⋆ = 0.1 0.15 0.2 0.25 0.3 and ρ ⋆ = 0.1 0.2 0.3 0.4 with the following parameters π ⋆ min = 0.1 and δ π ⋆ = 0.1 -0.2ε.

2 = 2 3 = 2 2 4 =

 22324 log(n)/d 1 + 10 -10 . (b) Middle threshold: S n,d log(n)/d. (c) Slower threshold: S n,d (log(n)/d) 1/4 .

Fig. 2

 2 Fig. 2 Proportions of bad estimations of g ⋆ for the parameter ε ∈ {0.05, . . . , 0.3} (in rows)following the thresholds (in columns) used for the balanced case: for each graphic, the number of rows n (y-axis) and the number of columns d (x-axis) varies between 20 and 4000.

Fig. 3

 3 Fig. 3 Proportions of bad estimations of g ⋆ for the parameter ε ∈ {0.05, . . . , 0.3} (in rows)following the thresholds (in columns) used for the arithmetic case: for each graphic, the number of rows n (y-axis) and the number of columns d (x-axis) varies between 20 and 4000.

Fig. 4

 4 Fig. 4 Average estimate of the distance d ∞ θ ⋆ , θ between the true parameters and the estimated parameters following ε (rows) in the balanced case: for each graphic, the threshold is represented by dierent symbols (• for S 1 , ▲ for S n,d 2 , ■ for S n,d 3 and + for S n,d 3 ) and the size for the number of nite values used; the number of rows n varies between 40 and 4000 and d is supposed equal n.

Fig. 5

 5 Fig. 5 On the top, boxplots of the computing times (in seconds; logarithmic scale) for each procedure (colour) in function of the numbers of rows and columns (n = d) over the 1200 simulations (200 per ε) for the balanced case. On the bottom, boxplots of quality of estimations averaged over the 20 matrices for each ε for each procedure (colour) in function of the numbers of rows and columns (n = d) for the balanced case.

Fig. 6

 6 Fig. 6 On the top, boxplots of the computing times (in seconds; logarithmic scale) for each procedure (colour) in function of the numbers of rows and columns (n = d) over the 1200 simulations (200 per ε) for the arithmetic case. On the bottom, boxplots of quality of estimations averaged over the 20 matrices for each ε for each procedure (colour) in function of the numbers of rows and columns (n = d) for the balanced case.

0 and g ⋆ e -2nt 2 + m ⋆ e -2dt 2

 22

Cexp n 1 n log 1 π0 2 ≤ exp -n 1 n log π ⋆ n,d min + 2t 2 ≤

 122 Thanks the same Assumption (TA), there exists a positive constant C > √ 2 such that for n and d large enough S Proof of Theorem 3: consistency in the sparse case First, if π ⋆ n,d min tends to zero (and asg ⋆n,d ≤ 1/π ⋆ n,d min ), the series expansion of t → log(1 -t) is used g ⋆n,d 1 -π ⋆ n,d min n = exp log g ⋆n,d + n log 1 -π ⋆ n,d min ≤ ⋆ n,d min + log 1 -π ⋆ n,d min ≤ exp n -1 n log π ⋆ n,d min -π ⋆ n,d min + o π ⋆ n,d min ≤ exp n -π ⋆ n,d min + o π ⋆ n,d min ≤ exp -nπ ⋆ n,d min + o nπ ⋆ n,dmin and, by the assumption (MA.2), nπ ⋆ n,d min tends to innity andg ⋆n,d 1 -π ⋆ n,dFor the same reasons, for all t > 0m ⋆n,d 1 -ρ ⋆ n,d and g ⋆ m ⋆ 1 -π ⋆ min ρ ⋆ min 1 -e -2t 2 nd -→ n,d→+∞ 0.Moreover, as g ⋆n,d e -2nt 2 = exp log g ⋆n,d -2nt exp -n2t 2 + o (n) .

>

  δ n,d π ⋆ > S n,d g C ⇔ δ n,d π ⋆ > S n,d g (C -1) + S n,d g ⇔ δ n,d π ⋆ -S n,d g n and d large enough, min δ n,d π ⋆ -S n,d g , S n,d g is S n,dg and the assumptions (I) and (TA) allow to conclude.

Fig. 11

 11 Fig. 11 Quality of estimations averaged over the 20 matrices for each ε (in rows) for each procedure (colour) in function of the numbers of rows and columns (n = d) for the arithmetic case.

  Input: data matrix x, threshold for row Sg and for column Sm. Selection of i 1 < . . . < i g -1 such that (G i 1 , . . . , G i g -1 ) are every greater than Sg.

	// Computation of gaps for i ∈ {1, . . . , n} do Computation of X i =	x i+ d .	// O(nd)
	Ascending sort of X (1) , . . . , X (n) .	// O(n log n)
	for i ∈ {2, . . . , n} do
	Computation of the gaps G i = X (i) -X (i-1) .
				// O(n)
	// Computation of g
	// O(n)
	// Computation of Z for i ∈ {(1), . . . , (n)} do Denition of z (i)k = 1 if and only if (i k-1 ) < (i) ≤ (i k ) with i 0 = 0 and i g = n.
				// O(n)
	// Computation of m and W
	Do the same on the columns.	// O(max(dn, d log d))
	// Computation of θ
	for k ∈ {1, . . . , g} do
	Computation of π k =		z +k n .	// O(n/Sg)
	for ℓ ∈ {1, . . . , m} do
	Computation of ρ ℓ =	w +ℓ d .	// O(d/Sm)

Table 2

 2 Computing times for the dierent procedures in columns (LG with the four threshold, VBayes with the true number of clusters and VBayes coupled with the ICL criterion on a 2 × 7 square grid) following the number of rows and columns (in rows) for the arithmetic case. Each cell represents the average (in milliseconds) and the standard deviation (in parenthesis) over the 1200 simulations (200 per ε).

	n=100	S 1 2 (0.9)	S n,d 2 2 (0.7)	S n,d 3 2 (0.7)	S n,d 4 2 (1.2)	Simple VBayes VBayes+ICL 530 20819 (302.3) (7768.3)
	n=200	4 (1.4)	4 (1.4)	4 (1.4)	4 (1.4)	1677 (693.7)	58494 (18294.9)
	n=300	7 (2.4)	7 (2.4)	7 (2.4)	7 (2.4)	3183 (1145.8)	117450 (34865.2)
	n=400	10 (3.2)	10 (3.7)	10 (3.2)	9 (3.3)	5637 (2119)	203479 (55356.2)
	n=500	14 (3.9)	14 (4)	13 (3.9)	13 (3.9)	8503 (2862.3)	304668 (76486.8)
	n=600	18 (5.2)	18 (5.5)	18 (5.9)	18 (5.3)	12423 (3889.1)	435626 (100281.1)
	n=700	25 (6.5)	24 (6.6)	24 (7)	24 (6.8)	17794 (7596.6)	622631 (180216.9)
	n=800	27 (2.1)	26 (2.9)	25 (2.6)	25 (2.5)	18475 (3680)	684012 (94177.3)
	n=900	36 (3.2)	35 (3.3)	35 (3.2)	35 (3.8)	22635 (4697.3)	828716 (106707.4)
	n=1000	41 (2.8)	41 (3)	41 (3.4)	40 (2.2)	31799 (13840.7)	1130163 (371089.7)
	n=1100	47 (4.4)	47 (5.5)	47 (5.7)	46 (4.6)	32203 (11270.2)	1226172 (290324.2)
	n=1200	47 (3.3)	46 (3.6)	47 (3.6)	47 (4.2)	34399 (6683.5)	1295896 (167982.4)
	n=1300	72 (21.2) (21.3) (21.8) 72 72	71 (21)	45497 (11230.9)	1673698 (353671.7)

  Proposition 1 Under Assumption (I) and if Sg ∈]0, δ π ⋆ [ and Sm ∈]0, δ ρ ⋆ [ for n, d large enough:

An implementation in the language R is available on the following Gitlab: https:// gricad-gitlab.univ-grenoble-alpes.fr/braultv/largest-gaps

For the sake of repeatability, all codes are available on the following Gitlab: https:// gricad-gitlab.univ-grenoble-alpes.fr/braultv/largest-gaps

Cluster Luke44, 28 cores, 128Go RAM, GPU 2xK40m, 2xIntel Xeon E5-2680 2.40 GHz; more informations on the url https://scalde.gricad-pages.univ-grenoble-alpes.fr/web/ pages/presentation.html

Thanks to Stéphane Robin and the reviewers for them suggestions. This work has been partially supported by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01) funded by the French program Investissement d'avenir and partially supported by MIAI@Grenoble Alpes (ANR-19-P3IA-0003). All the computations presented in this paper were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by Grenoble research communities.


D Chen Legrand

C Barus

V Garçon

Square wave voltammetry measurements of low concentrations of nitrate using Au/AgNPs electrode in chloride solutions

Keywords:

The aim of this work is to highlight the potential of using a modified gold electrode with controlled quantity of silver nanoparticles as a working electrode to detect low concentrations of nitrate in chloride solutions. Optimal charge for silver deposition has been determined to obtain the highest signal for the nitrate reduction as the electrocatalytic properties of the bimetallic electrode were directly influenced by its composition. According to the Volcano plot obtained the charge chosen was -52 µC for a 3 mm diameter electrode, corresponding to 4.6 x 10 15 Ag atoms cm -2 . It has been shown that dioxygen did not participate to the nitrate reduction mechanism. In order to decrease the limit of quantification, square wave voltammetry was preferred to less sensitive cyclic voltammetry. Nitrate was quantified in chloride solutions in the concentration range found in the open ocean, i.e. 0.39-50 µmol L -1 with a good linear regression (R 2 =0.9969). The stability of the bimetallic Au-Ag systems has been evaluated and showed almost no difference on the signal recorded over a 26 days period which is suitable to consider an in situ sensor development for marine applications.

Nitration detection, Square Wave Voltammetry, Environmental chemistry.

-Introduction

The nitrogen cycle is one of the most important biogeochemical cycles in the marine environment. Nitrogen is found in five oxidation states: (-III) NH 4 + /NH 3 , (0) N 2 , (+I) N 2 O, (+III) NO 2 and (+V) NO 3 -where N 2 is the most abundant species but essentially unreactive. Nitrate (NO 3 -)

is the main form assimilated by marine organisms and as the other nutrients (silicate and phosphate) it can be used as a water mass tracer [1]. The concentration range found in the open ocean is from a few nanomolars to 50 µmol L -1 . Long term monitoring of nitrate ions concentration with in situ sensors is the key to improve our understanding of biogeochemical cycling in the ocean as well as to follow water masses and their mixing.

The classical technique to measure nitrate concentration in seawater is wet chemical colorimetry [2,3]. However, this method requires chemical reagents that need to be changed regularly and high power demand. Ultraviolet spectroscopy measurements at 220 nm are also used to detect nitrate ions in seawater [4]. In situ nitrate sensors for marine application using ultraviolet absorbance such as the ISUS (In situ Ultraviolet Spectrophotometer) or the SUNA (Submersible Ultraviolet Nitrate Sensor) are implementable on under-water vehicles such as ARGO profiling floats and they are already available on the market. However their sensitivity is still quite high, with a limit of detection around 1.5 µmol L -1 . In addition they suffer from strong chemical interferences such as with bromide ions, carbonates, or Colored Dissolved Organic Matter (CDOM) [5,6].

In order to decrease the limit of quantification, the cost and the size of sensors, electrochemistry is proposed to measure nitrate concentration in seawater. Nitrate being an electroactive species, it can be detected directly with electrochemical methods (voltammetry, amperometry…), no reagent is needed, which will facilitate miniaturisation of the detection system suitable for in situ analyses [2,7,8]. The electroreduction of nitrate has been intensively studied on various electrodes (Pt, Pd, Ni, Pb, Zn, Fe, Cu, Ti, Au, Au(Hg), Ru, Rh, Ir, Bi, Sn, alloys of Cu, Sn…) and is known to be a complex process [9,10]. Depending on the experimental conditions such as electrolyte anions [11], pH of the solution [12], nitrate concentration [13,14] and the structure of the working electrode [15], the mechanism is different and leads to nitrite (NO 2 -), ammonia (NH 4

+

), nitrogen (N 2 ) as bi-products (HNO 2 , NO + , NH 2 OH, N 2 O…) [14,16,17]. Electrochemical sensors have been already used to detect nitrate in fresh water [18]. Hafezi and Majidi developed an amperometric sensor to detect nitrate in foodstuffs and mineral waters in the concentration range between 1 and 35 µmol L -1 with a limit of detection of 0.59 µmol L -1 . However amperometry requires to stir the solution to control the convection which is not suitable for in situ sensor in seawater [19].

Among the electrodes, bimetallic electrodes showed a better sensitivity towards the detection of nitrate as the signal is enhanced due to higher specific electrode surface area [20][21][22]. The choice of electrode materials as well as the experimental conditions for nanoparticles deposition to form the bimetallic electrodes are crucial and will have a large influence on the nitrate detection [23]. Previous work on gold modified with silver nanoparticles showed a good potential to detect nitrate in seawater [24]. However, the limit of quantification was still too high (i.e. 10 µmol L -1 ) to consider in situ application in marine environments [24]. The aim of this work is to focus on the detection of low concentrations of nitrate in chloride solutions. To do so, an optimisation of the parameters to deposit silver nanoparticles on gold electrode is essential. Square wave voltammetry measurements are proposed to go further in terms of the sensitivity of our analytical method. To develop an in situ electrochemical sensor competitive with those already existing in oceanography, the limit of detection has to be the smallest possible (< 1.5 µmol L -1 [5]) and the modified electrode has to show a good stability with time.

-Experimental section 2.1 -Chemicals

All solutions were prepared in Milli-Q water (Millipore Milli-Q water system). Standard nitrate solutions were prepared with potassium nitrate (KNO 3 , PURATRONIC®, 99.995 % purity) from Alfa Aesar in sodium chloride solution (NaCl, EMSURE®, supplied by Merck) at 34.5 g L -1 (0.6 mol L -1 ).

Sulphuric acid solution at 0.5 mol L -1 used to clean electrode surfaces was prepared from a 98% H 2 SO 4, EMSURE®, solution supplied by Merck.

Silver nanoparticles electrodepositions on gold electrode surfaces were made with silver nitrate AgNO 3 (AnalaR, NORMAPUR®, VWR Chemicals) dissolved in KNO 3 (Alfa Aesar) electrolyte solution.

-Material and methods

Chronoamperometry (I = f(t), fixed E), Cyclic Voltammetry (CV) (I = f(E)) and Square Wave Voltammetry (SWV) (I = f(E step ), with I = I pulse forward -I pulse backward ) were performed in sodium chloride solutions ([NaCl] = 34.5 g L -1 ) using a PGSTAT 128N potentiostat (Metrohm) controlled by NOVA software at room temperature under atmospheric conditions. A classical three electrodes system was used for all experiments. The counter electrode is a platinum grid purchased at GoodFellow. The reference electrode was a double junction Ag/AgCl/KCl (3 mol L -1 ) electrode; all the potentials were referred to this electrode. The gold working electrode was a commercial gold disc electrode (φ = 3 mm) from Metrohm.

-Preparation of modified electrode Au/AgNPs

The gold disc surface was polished firstly with alumina (0.3 µm slurry diameter). Secondly, it was electrochemically cleaned in sulphuric acid solution (H 2 SO 4 , 0.5 mol L -1 ) by polarising the electrode 10 seconds at +2 V and 10 seconds at -2 V to form O 2 and H 2 bubbles, respectively, at the electrode surface. Then, cyclic voltammograms were recorded between E 1 = 0 V and E 2 = 1.5 V at 100 mV s -1 until reproducible cycles were obtained.

The electrodepositions of silver nanoparticles onto gold substrates were performed by chronoamperometry at E = -0.2 V optimised potential in a 0.1 mol L -1 KNO 3 solution containing silver nitrate at [AgNO 3 ] = 0.17 mmol L -1 . The end of electrolyses was either controlled by the time of electrolysis (in seconds) or the electric charge (Q in Coulombs (C)). Several charges were studied in order to determine the optimised silver deposition. The resulting Au/AgNPs modified electrodes were then rinsed with MilliQ water. Cyclic voltammograms recorded with Au/AgNPs modified electrode in seawater highlighted the presence of silver on gold in the oxidation domain at +0.6 V (results not shown).

-Results and discussion

Fig. 1 shows cyclic voltammograms recorded on a gold commercial electrode (Au, Metrohm, φ = 3 mm) modified with silver nanoparticles (AgNPs) deposited at E = -0.2 V during 7 seconds, in sodium chloride solution (34.5 g L -1 ) containing nitrate ions ([NO 3 -] = 9 mmol L -1 ) at 20 and 500 mV s -1 scan rates [24].

Fig. 1

The reduction starting around -0.2 V corresponds to the reduction of dioxygen. The reduction of nitrate occurs between -0.8 V and -1 V on Au/AgNPs surface. Only one peak is observed at -1.0 V at high scan rate (500 mV s -1 ) whereas two peaks are visible around -0.9 V and -1.0 V, respectively, at lower scan rate (20 mV s -1 ) showing the complexity of the mechanism. Aouina et al. also observed several reductions processes for nitrate reduction at 20 mV s -1 . They identified the formation of nitrite at the less negative potential and ammonia at the most negative potential [12]. Other authors also managed to highlight the formation of co-products depending on the experimental conditions (pH, electrolyte, nitrate concentration) used [10,[12][13][14][15][16][17]. According to the bibliography, nitrate reduction requires first the adsorption of NO 3 -on the electrode surface (reaction 1) [19,25,26]. The step determining the rate of the reaction on most electrodes is the reduction of NO 3 -(ads) to NO 2 -(ads) (reaction 2) [10]. The main final product is ammonia (reaction 3) [12,16,27].

(

(

(3)

NH 4 + + 3H 2 O ← NO 3 -+ 10H + + 8 e -
The scan rate of v= 500 mV s -1 is chosen to detect nitrate as the signal is much higher (Ip (v=500 mV s -1 ) = 5.8 Ip (v=20 mV s -1

) ) than for the lower scan rate. It is also more convenient to treat only one, well defined, reduction peak than two tangled reduction waves.

As it was mentioned earlier, the structure of the bimetallic electrode has a strong influence on the signal [28]. Therefore, in order to improve the electrocatalytic properties of the electrode it is crucial to determine the optimum conditions for silver deposition on the gold working electrode.

In order to control the exact quantity of silver nanoparticles deposited on the gold surface (Metrohm, φ = 3 mm), the electrolyses at fixed potential (i.e. E = -0.2 V) realised to form Au/AgNPs electrodes are controlled by the charge in Coulomb (C) and not by the time. Depending on the initial state of the bare gold electrode, electrolyses controlled with time will give different Au/AgNPs electrodes as the AgNPs quantities deposited (charges) are not necessary the same from one electrolysis to another. The cut-off with the electric charge is the only way to ensure the deposition of the same quantity of silver nanoparticles.

The influence of charges (between -14 and -75 µC) used for silver deposition is investigated by analysing the cyclic voltammograms recorded with these different Au/AgNPs in sodium chloride solution (34.5 g L -1 ) with 0.1 mmol L -1 of nitrate. The forward cycles obtained for each charge are presented on Fig. 2.

Fig. 2

The reduction of dioxygen appears on Fig. 2 in two peaks at -0.4 V and -0.6 V corresponding to the reduction of O 2 in two steps first into H 2 O 2 then H 2 O observed on the gold electrode. For all the charges applied, the gold electrodes are therefore not totally covered with silver. Even if the experiments were realised in atmospheric conditions, with saturated O 2 solutions, the signal observed on Fig. 2 for O 2 varies from one experiment to the other and might have an influence on NO 3 -measurements. The influence of O 2 on nitrate measurement will be checked later on.

As expected the intensity corresponding to nitrate reduction increases when the silver electrodeposition charge increases from -14 µC to -38 µC as the effective surface area of the electrode increased (Fig. 2). However, increasing too much the charge leads to a decrease of nitrate signal for -71 µC and -75 µC (Fig. 2). This result is more visible on the corresponding Volcano plot (I peak =f(|Q|)) presented on Fig. 3 showing two characteristic linear slopes. This behaviour has been observed in the bibliography and can be explained by the Ostwald ripening widely known in the field of crystallization and corresponding to the dissolution of smaller particles that will grow in larger ones [29][30][31]. The effective surface area of Au/AgNPs therefore decreases and the electrode tends to behave as a bare electrode. The mechanism of Ostwald ripening includes a diffusion process at the surface and grain boundary, and an electrochemical mechanism. Larger and smaller particles behave as different metals with different standard electrode potentials explaining that it is easier to reduce NO 3 -on electrode with smaller particles than on a bulk material or on larger particles electrode [30].

Adsorption is also promoted on small nanoparticles as the effective surface is bigger and so makes the reduction on this electrode easier. This is why pure bare gold electrode is inactive toward nitrate reduction (weak adsorption of the reactants) [26]. The maximum catalytic activity should be for a surface partially covered with AgNPs which is the case for all the charges tested as the characteristic reduction of O 2 in two steps on gold is always visible on Fig. 2. Therefore, in order to obtain the highest NO 3 -signal the charge to deposit AgNPs has to be the highest possible but before the Ostwald ripening process starts. According to the Sabatier principle, the best condition is determined on Volcano plot (Fig. 3), corresponding to the intersection of the 2 straight lines (I peak = f (|Q|)) [26,32]. The optimum charge found to form the modified Au/AgNPs electrode and to detect NO 3 -is -56.3 µC. Taking into account the experimental errors to build the Volcano plot and determine the optimum condition, the charge chosen for silver nanoparticles deposition is -52 µC corresponding to 4.6 x 10 15 Ag atoms cm -2 .

Fig. 3

In the open ocean, the concentrations of dissolved dioxygen vary from 300 µmol L -1 in saturated surface waters down to around 50 µmol L -1 for deeper waters. Some areas called Oxygen Minimum Zones (OMZ) are very depleted in dioxygen in sub-surface where the concentration can be close to zero [33]. On the potential window used, the dioxygen reduces just before the potential of nitrate reduction and variations of its signal have been observed even in oxygenated solutions (Fig. 2). It might have an influence on nitrate reduction even more if its concentration changes from an experiment to another. Therefore in order to study the influence of dioxygen reduction on the nitrate detection, cyclic voltammograms were recorded with the same modified electrode (Au/AgNPs: Q = -68 µC) on the same nitrate solution at 1 mmol L -1 with and without dioxygen. The solution has been deaerated during 15 minutes with nitrogen bubbling to remove all the dioxygen and a constant flow of nitrogen was maintained over the surface during the experiment to prevent dioxygen to diffuse into the solution. The results are presented on Fig. 4. [24]. They proposed a catalytic reaction between the intermediate oxygen peroxide formed from the incomplete reduction of dioxygen on gold electrode with nitrite ions produced from the incomplete nitrate reduction. However, the increase of NO 3 -peak intensity they observed on Au/AgNPs is probably due to an increase of electroactive surface area of the modified electrode Au/AgNPs from one deposition to the other as the amount of silver nanoparticles deposited was not controlled with the electric charge but with the time of electrolysis. In conclusion, the mechanism of nitrate reduction however complicated it may be, is clearly not influenced by dioxygen which is compatible with our aim to develop a new electrochemical sensor in the open ocean where O 2 concentration varies significantly throughout the water column. Now in order to go further on the detection of low concentrations of nitrate using the optimized Au/AgNPs (Q = -52 µC) electrode, square wave voltammetry method is proposed to detect nitrate. This electrochemical method is more sensitive than cyclic voltammetry thanks to short-term square shape potential pulses combined with staircase potential ramp and current sampling allowing to decrease the contribution of capacitive current. This pulse technique is attractive for quantitative analysis since it leads to better resolved peaks (Gaussian shape), less distorted by the capacitive current. Indeed, in order to obtain the highest sensitivity (corresponding to the electron(s) transfer), it is important to reduce the effect of the capacitive current and only measure the contribution of the faradic current. The current is therefore sampled at the end of each potential pulse (forward pulse and backward pulse) as the faradaic current ( ) decreases slower than the capacitive current ( ) with time [START_REF] Mirčeski | Square-Wave Voltammetry : Theory and Application[END_REF].

The SWV parameters have been optimized to obtain the least noisy signal possible using classical frequency of f = 100 Hz, an amplitude of E SW = 60 mV and a staircase potential E step = 5 mV, corresponding to the same scan rate than previously used with cyclic voltammetry (i.e. 500 mV s -1 ). Common values for amplitude and step potential are respectively and , with n the number of electron(s) exchanged [START_REF] Bard | Electrochemical Methods: Fundamentals and Applications Chapter 7[END_REF]. So considering the reduction of NO 3 -into NO 2 -, 2 electrons are exchanged and E SW should be around 25 mV and E step around 5 mV. However experimental results showed better resolution with higher amplitude. Indeed, amplitude affects the electrode kinetics of the electrode reaction. Higher amplitude gives higher net response in current, however too large amplitude may lead to interferences by the charging current and the ohmic drop resistance [START_REF] Mirceski | [END_REF]. The result obtained for 100 mV amplitude (result not shown) was indeed not satisfying and much noisy, so 60 mV amplitude has been fully validated.

The square wave voltammograms obtained with Au/AgNPs (φ = 3 mm, Q = -52 µC) electrode for nitrate concentrations between 0.39 and 50 µmol L -1 are presented on Fig. A shift of nitrate peak potential towards more anodic potentials (from -1.1 V to -0.9 V) (less negative potentials) is observed on Fig. 5 when NO 3 -concentration increases indicating an easier electronic transfer at high concentration than at low NO 3 -concentration. This result has also been observed by Bui et al., on different pulse voltammetry measurements on carbon paper electrode functionalised with gold nanoparticles [18].

As it has been mentioned before, the nitrate electroreduction mechanism is influenced by the nitrate concentration and requires as well the adsorption of nitrate on the electrode surface before its reduction [14,26]. Considering the adsorption reaction (1), the equilibrium of the reaction might shift towards the formation of NO 3 -(ads) if the NO 3 -(sol) increases, facilitating its reduction, therefore explaining the anodic peak potential shifting observed. Also the nitrate adsorption is quite weak and strongly influenced by the adsorption of anions and hydrogen (particularly on platinum electrodes) [13,15]. There is therefore a competition on the electrode surface between NO 3 -(ads) and adsorbed anions (Cl -) making harder the adsorption of nitrate ions especially at really low concentrations.

The peak intensities for each concentration have been measured using NOVA software and the corresponding calibration curve is presented on Fig. 6. A good linear behaviour is obtained between 0.39 µmol L -1 up to 50 µmol L -1 which covers the whole concentration range found in the open ocean. The limit of quantification is 0.39 µmol L -1 . Now, in order to use this modified electrode in a sensor for marine applications, its stability with time has to be evaluated particularly the stability of the silver nanoparticles. The same modified electrode Au/AgNPs (φ = 3mm, Q = -52 µC) is continuously stored in sodium chloride solution and used to detect regularly nitrate at [NO 3 -] = 25 µmol L -1 using square wave voltammetry during 625 hours (~26 days). The peak intensity corresponding to nitrate reduction remained at 95 % from the initial intensity peak (t=0) for 26 days showing quite a good stability of the system and its potential as an electrode in a marine sensor.

-Conclusions

The results presented in this work show a sensitive method for the nitrate detection in sea water using a very simple three electrodes electrochemical system. The study of the influence of the electric charge to deposit silver particles on nitrate detection allowed to determine the best conditions to enhance the electrocatalytic properties of the bimetallic electrode. The optimum charge found is Q = -52 µC corresponding to ~4.6 x 10 15 Ag atoms cm -2 . A very good calibration has been obtained using this modified Au/AgNPs electrode and square wave voltammetry measurements on the concentration range found in the open ocean (0.39 -50 µmol L -1 ) with a limit of quantification of 0.39 µmol L -1 . Finally the study of the electrode stability in time showed a good potential for in situ use.

To decrease the LOD of nitrate detection, a larger electrode surface could be used in order to increase the intensity. Also a comparison with other bimetallic couples should be tried as copper showed very good electrocatalytic properties for nitrate detection [19] or platinum substrate instead of gold for example.

Figure captions

Fig. 4

 4 Fig. 4

Fig. 4

 4 Fig.4shows indeed the disappearance of the two dioxygen peaks at -0.4 V and -0.6 V after 15 minutes of nitrogen bubbling. But in both cases, the peak intensities at E = -1.1 V corresponding to nitrate reduction are identical (~50 µA). Almost the same result has been obtained by Gross et al. in 2015 on Au/AgNPs electrode where the difference observed on the signals with and without dioxygen was extremely weak (experimental errors)[34]. The O 2 has therefore no influence on nitrate detection using Au/AgNPs electrode unlike previous work by Fajerwerg et al. suggested[24]. They proposed a catalytic reaction between the intermediate oxygen peroxide formed from the incomplete reduction of dioxygen on gold electrode with nitrite ions produced from the incomplete nitrate reduction. However, the increase of NO 3 -peak intensity they observed on Au/AgNPs is probably due to an increase of

  5. As previously observed on cyclic voltammograms, three peaks are observed at -0.25, -0.6 and -1.0 V corresponding respectively to the reduction on gold surface in two steps of O 2 into H 2 O 2 then H 2 O 2 into H 2 O and finally the reduction of NO 3 -where the peak intensities increase with NO 3 -concentrations. The fact that signals show positive current even for reduction processes is due to the amplitude chosen for the pulses as the differential intensities (I = I forward -I backward ) are plotted versus staircase potential [38].

Fig. 5

 5 Fig.5

  Fig. 6

Fig. 1 :Fig. 2 : 1 Fig. 3 :

 1213 Fig.1: Influence of the scan rate on the cyclic voltammograms obtained in NaCl solution (34.5 g L -1 ) containing 9 mmol L -1 KNO 3 using Au/AgNPs (E = -0.2 V, t = 7 s) at scan rate 20 mV s -1 (•••••) and 500 mV s -1 ()

Fig. 4 : 1 Fig. 5 :Fig. 6 :

 4156 Fig.4: Influence of dioxygen on cyclic voltammograms obtained in NaCl solution (34.5 g L -1 ) containing 1 mmol L -1 KNO 3 using Au/AgNPs electrode (Q = -68 µC) with dioxygen (straight line) and without dioxygen (dotted line); scan rate: 500 mV s-1

Acknowledgments

The present work was supported by the MIACTIS project: "Microsystèmes intégrés pour l'analyse des composés en traces in situ" (Integrated microsystems for the analyses of in situ trace compounds) as part of the overall project "Instrumentation and Environmental Sensors" funded by the RTRA-STAE: thematic network of the foundation for scientific cooperation "Sciences et Technologies pour l'Aéronautique et l'Espace" (Sciences and Technologies for Aeronautic and Space) in Toulouse, France.


S Chabert
email: chabert@laplace.univ-tlse.fr

J J Gonzalez

P Freton

PLASMA ARC CUTTING -REVERSED SWIRL RING, ELECTRODE THREAD AND CUT DIRECTION EFFECTS ON KERF GEOMETRY

Keywords: Plasma arc cutting, air plasma, cutting quality

Plasma arc cutting is used to cut any conductive material. It consists in blowing pressurized gas and feed current to an arc, leading to a thin plasma dart able to melt the material down and blow it away, creating a kerf. Its quality depends on its shape. This paper shows, through experimental measurements, how the inner geometry of the torch can affect the cut quality. It appears that one side of the kerf is much more oblique and sensitive to factors variation than the other. A theory based on a computational fluid dynamics model is proposed to investigate the causes of these phenomena.

Introduction

The torch under study is an AT-125 torch from the Toparc brand. It can withstand direct currents up to 125A. It uses compressed air as both the plasma and shielding gas. Its components can be seen in Figure 1.

The electrode, which has an helix and the swirl ring, which blows air tangentially, break the torch axisymmetrical layout. The tangential air injection is aimed at constricting the plasma jet. There is a backward flow that passes into the electrode thread. In plasma arc cutting, one side of the kerf always displays a side straighter than the other. In previous unpublished work, it has been noticed that changing the electrode rear part with a cylinder drilled with holes spread in a circular manner (instead of the usual thread) led to more symmetric kerfs. This paper pushes the investigation further with regards to torch geometry's effects on cut quality. Prototypes electrodes with reversed thread and swirl rings with opposite tangential injection have been made (Figure 2).

Cut direction has also been investigated. This gives three parameters, namely thread direction, injection angle and cut direction, which effects have been studied comparatively through a Taguchi's [START_REF] Pillet | Design of Experiments Using the Taguchi Method[END_REF] design of experiments. This method allows to study the effects of the parameters and their interactions (synergy).

Experimental Setup

The cuts were performed on 4mm thick S355 steel plates. The torch is moved by a Computer Numerical Control (CNC) Autotec cutting table. The used generator is a Gys Neocut105. All over the experiment, cutting voltage The cut kerfs are parallel lines cut either in one direction or its opposite, later referred as Direction 1 (Figure 1) and Direction 2.

After the plasma kerfs are made, the plate is cut perpendicularly to them with a band saw in order to have a profile observable with a microscope.

Measurement protocol

The measurements are performed with a Keyence VHX-6000 microscope and its embedded measurements tools. The precision is 10 µm. Measurements have been made according to standard ISO9013.

Figure 3 is taken from ISO9013 standard and shows one side of a kerf and the measured values. The ISO-9013 standard gives ∆a = 0.3mm because the plate thickness is a=4mm.

U (taken on each side of the kerf) is the value this paper focuses on and quantifies the oblicity (bevel) of the kerf. With regards to Taguchi's method, this is the studied response.

Design of Experiments

The experiments were conducted by varying the factors as seen in Table 1. For each line of this table, three kerfs were made, measured as explained previously and the results were averaged (side by side) to limit results scattering.

Experimental results

The overall results (for both sides) range from 100 to 700 µm. Standard deviation is 52 µm for the left edges and 65 µm for the right edges. In this paper, references to "right" and "left" reference for the edges of a kerf are taken with respect to the torch's direction 1 (torch direction 1 is forward and the kerf is left behind).

5 bar 105 A trials

A full 3 factors -2 levels design was employed, wherein each combination of factors at all levels was executed, resulting in a total of 8 experiments (23 = 8). These correspond to the first 8 rows of Table 1.

On the left side of the kerfs

On the left sides, nor the studied parameters nor their interaction has any measurable effect on the bevel angle. These sides remain very straight and display an average U of 151 µm for the 5 bar 105 A trials. They are also very consistent as their standard deviation is 15µm.

On the right side of the kerfs

The right sides are much more oblique and much more sensitive to the parameters variation and their interactions.

For the 5 bar 105 A trials, on the right sides, U spreads from 341 to 591 µm with an average of 481 µm and a 91 µm standard deviation. The calculated effects are given in Table 2.

As per Taguchi, this is the effect that the parameters or their interaction have on the response (U). It can be seen that a left hand thread electrode (electrode 1), a -30°swirl ring (swirl ring 2) and cutting direction 2 lead to straighter right sides of the kerfs (negative effect on U). Also, the effects of the electrode helix and of the swirl ring injection angle interacts together (but the helix-cut direction and swirl-cut direction interactions are negligible and thus not displayed in Table 2). Electrode 2 and swirl ring 1 lead to even straighter right sides when used together.

Varying pressure and currents trials

More tests were conducted with variable pressure and current. Not all configurations were studied (Table 1).

Let us recall that more vertical edges on the left were consistently observed over the experiment. Therefore, representing only the difference between the U values of the left and right edges for each case provides interesting information about these results, as is seen in Figure 4. A value close to zero indicates a rather symmetrical kerf, while a highly negative value corresponds to a strongly asymmetrical kerf.

First, the cutting direction effect is low but steady except for configuration 3: the kerfs are a bit more symmetrical with cutting direction 2. It can be seen easily on Figure 4 because the left pointing arrows are above the corresponding right pointing arrows.

About the configuration 1 and 4, for which the most data are available :

1. Decreasing current from 105 to 45 A at steady pressure makes kerfs more symmetrical the swirl ring injection (while 5 bar is experimentaly set on the cutter, because pressure drop occurs in the circuit).

2. Shield air injection (pressure inlet): 4.3 bar, 25°C

3. Ambient (pressure outlet): 1 bar and 25°C bar 4. Electrode head (wall): Uniform current profile. 2220°C (melt of the hafnium).

5. Nozzle channel (wall): Resolving Biot and Savart 6. Kerf (pressure outlet): Scalar potential value is 0V. This the porous anode. Unindexed walls temperature is set at 25°C.

Model results

Three simulations have been run. They all share the same parameters except for the current and the swirl injection pressure. These parameters and results of interest are displayed in Table 3. Figure 6 shows the swirl (Z) momentum field. Its shape remains the same in every case, only the values' range change. Z momentum written in Table 3 is taken at point A.

Decreasing the current reduces the (viscous) arc width, allowing more air to flow around the plasma jet. Increasing the injected pressure also forces more air into the plasma jet fringes. In both cases this leads to higher mass flow rates.

The swirl component of the flow could lead to assymetrical results, because the kerfs are open behind (with regards to torch movement) the torch and closed before the torch. As a result, air hitting the kerf might leave it immediately if the swirl effect bring it outwards. However, on the opposite side, the swirl effect shall bring the flow inward onto the molten metal: the air flow encounters different conditions depending on the angular position around the plasma jet where it its the plate. As a consequence, the swirl momentum might be held responsible for (part of) the assymetry. From the model : Decreasing current or increasing pressure leads to higher flow rates and noticeably swirl momentum in the plasma jet As a result, it might be suggested that:

Discussion

1. Higher flow rates make kerfs more symmetrical. This could be because of a lack of momentum axial symetry in the plasma that gets overriden at higher global flow rates (high flow rate lead to straighter kerfs).

2. Higher swirl momentum make kerfs more symmetrical. This sounds conflicting with the phenomena given as conclusion of previous section. However this could also be explained by some flow angular assymetry in the plasma flow that would get smoothed by higher swirl momenta.

Conclusion

The reasons for the kerfs assymetry in plasma arc cutting are not yet explained. This paper suggests the shape of the inner torch parts has an impact on the kerfs. Besides, the electrode helix and swirl ring angle interact with each other. However, the model presented in this paper is 2D axisymmetrical and thus cannot represent this. Cut direction also shows some noticeable yet not predominant effect on kerfs. However, every kerf in this paper showed much straighter left sides regardless of the set of parameters, so other factors that remain yet to discover may be provoking this phenomenon. Especially any assymetry in air efflux even upstream the plasma jet and its relation to swirl ring and helix should be investigated.

Figure 1 .

 1 Figure 1. AT-125 torch components

Figure 2 .

 2 Figure 2. ZY Cut plane of the original and prototype swirl rings

2 .

 2 Figure 5. 2D model domain and boundary limits -non indexed walls all share the same settings

1 .

 1 Figure 6. Azimuthal momentum field -Reference case

Table 1 .

 1 Design of experiment

	vol. 0 no. 0/2023

Table 2 .

 2 Effects of the parameters on the right side given in µm

Table 3 .

 3 Simulation cases and results

Acknowledgements

Funding, experimental setups and technical support provided by SAS Gys.

HAL Figure 4. Average difference between U on the left and right side of each kerf

These phenomena led us to a theory that is exposed in section 7 (Discussion) and that is supported by the use of a numerical model that is presented in the following sections.

The numerical model

The numerical model we used runs on Fluent represents the torch in a limited domain meshed with 17k cells and represented in Figure 5 along with the boundaries.

The model is described as follows :

1. Mathematical method: Finite Volume Method solved with user defined function to take the electromagnetic values into account and to calculate the source terms of the energy and momentum equations (Joule effect, Laplace forces).

2. Geometry: The model is 2D axisymmetrical (with swirl) and as a consequence, cannot consider the helix of the electrode nor the torch movement. The injections hole of the swirl ring have been replaced by an equivalent section full ring. The backwards flow is not represented but the effluxes in the other circuits met experimental accordance. No solid is meshed and only the fluid zone is resolved.

3. Hypothesis: The model is time-steady (no arc fluctuations), the flow is laminar, the radiation is taken into account through the Net Emission Coefficient method, no matter is ablated from the solid parts. A porous anode is used. Current is injected at the electrode head with a uniform profile over the surface of the boundary. No slip condition (fluid speed =0) on every wall.

4. Initialization: The current is fed into the simulation after the plasma zone has been patched to 10k°C to bring sufficient level of conductivity.

5. Equations: Navier Stokes equations, scalar potential (with porous anode) and vector potential (closed by Biot & Savart [START_REF] Freton | Magnetic field approaches in dc thermal plasma modelling[END_REF]) are solved.

The boundaries conditions (Figure 5) are :

1. Swirl ring injection (pressure inlet): fed air under 4.3 bar at 25°C with a 30°angle of swirl injection. Experimental and numerical accordance with regards to the flow rates is reached for 4.3 bar simulated at


Maëlys Cadel

Isabelle Cousin

Marine Lacoste

Lionel Cottenot

Study the influence of environmental drivers on N2O emissions under semi-controlled conditions using a rainfall simulator

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.