Avner Bar-Hen

P Barbillon
email: pierre.barbillon@agroparistech.fr

S Donnet

Block models for multipartite networks Applications in ecology and ethnobiology

Keywords: Networks, Latent Block Models, Stochastic Block Models, Variational EM, Model Selection, Ecology, Ethnobiology

Generalized multipartite networks consist in the joint observation of several networks implying some common individuals . Such complex networks arise commonly in social sciences, biology, ecology, etc. We propose a generative probabilistic model named Multipartite Blocks Model (MBM) able to unravel the topology of multipartite networks by identify clusters (blocks) of nodes sharing the same patterns of connectivity across the collection of networks they are involved in. The model parameters are estimated through a variational version of the Expectation-Maximization algorithm. The numbers of blocks are chosen with an Integrated Completed Likelihood criterion specifically designed for our model. A simulation study illustrates the robustness of the inference strategy. Finally, two datasets respectively issued from ecology and ethnobiology are analyzed with the MBM in order to illustrate its flexibility and its relevance for the analysis of real datasets.

The inference procedure is implemented in an R-package GREMLIN, available on Github (https://github.com/Demiperimetre/GREMLIN). Simulations are implemented as vignette of the package in order to ensure reproducibility of the results.

Introduction

Networks have become fundamental tools in various fields, such as biology, ecological theory or sociology to name but a few. Statistical analysis aims to study the structure of these networks for instance by unraveling the clusters or communities of individuals (nodes) shaping the observed interactions (see [START_REF] Kolaczyk | Statistical Analysis of Network Data: Methods and Models[END_REF], for a review). The recent years have witnessed a growing interest for multi-layer networks which is a classical but wide terminology for a collection of related networks (see [START_REF] Pilosof | The multilayer nature of ecological networks[END_REF], for a review of multi-layer networks in ecology for instance). Among others, multi-layer networks encompass dynamic networks, i.e. networks evolving with time (see [START_REF] Kim | A review of dynamic network models with latent variables[END_REF], for a review), multiplex networks when several types of relations are simultaneously studied on a common set of individuals (see [START_REF] Kéfi | How structured is the entangled bank? the surprisingly simple organization of multiplex ecological networks leads to increased persistence and resilience[END_REF][START_REF] Barbillon | Stochastic block models for multiplex networks: An application to a multilevel network of researchers[END_REF], for instance) and multipartite networks. Multipartite networks are a generalization of bipartite networks. In a bipartite network, the nodes (representing the interacting entities) are partitioned into two disjoint sets and an edge links a node from one set to a node from the other set (see Figure 1, left). In a multipartite network, the nodes are divided into more than two sets and edges link entities from different sets (see Figure 1, middle). In what follows, these pre-specified sets of nodes will be referred to as functional groups. Such multipartite networks arise in ecology when studying the interactions between several groups of species such as the interactions plant/pollinator, plants/ants, etc [START_REF] Pocock | The robustness and restoration of a network of ecological networks[END_REF][START_REF] Dáttilo | Unravelling Darwin's entangled bank: architecture and robustness of mutualistic networks with multiple interaction types[END_REF] or in biology when analyzing networks issued from multi-omics datasets involving proteins, etc. (see Pavlopoulos et al., 2018, for instance). Generalized multipartite networks are an extension of multipartite networks: the nodes are still partitioned into functional groups but the interactions may occur not only between different functional groups but also inside some of the functional groups. Ethnobiology (which is the scientific study of the relations between environment and people) provides such networks. One of the problematic of ethnobiology is to understand how social relations between individuals (seed circulation in our example) may structure and guarantee biodiversity in the cultivated crop species (see Thomas and Caillon, 2016, for example). As a consequence, the network of interest has two types of nodes -namely farmers and crop species-and we observe interactions between farmers and crop species (who grows what) and inside the group of farmers (seed exchange). Marketing also provides generalized multipartite networks when individuals are connected through social networks and also interact with goods through their on-line purchases.

Many statistical tools have been proposed to analyze simple or bipartite networks. However, when individuals or biological entities are involved in several networks, there is a strong need to propose statistical tools handling these networks jointly. These tools are an essential step towards a better understanding of the global interaction systems at stake. The aim of this paper is to propose A central issue in network analysis is to be able to cluster nodes sharing the same connectivity patterns. To do so, two approaches are possible, either relying on classical metrics detecting nestedness [START_REF] Almeida-Neto | A consistent metric for nestedness analysis in ecological systems: reconciling concept and measurement[END_REF], modularity [START_REF] Barber | Modularity and community detection in bipartite networks[END_REF], etc, or adopting a probabilistic mixture strategy. Since the original paper by [START_REF] Snijders | Estimation and prediction for stochastic blockmodels for graphs with latent block structure[END_REF], probabilistic mixture models -classically referred as Stochastic Block models (SBM)-have proven their efficiency when aiming to cluster similar nodes based on their connectivity patterns without any a priori hypothesis about the patterns to be found (e.g. modularity, centrality, hierarchy). Latent Block Models (LBM) [START_REF] Govaert | Clustering with block mixture models[END_REF] are the extension of SBM for bipartite networks, resulting into a bi-clustering of the nodes.

In this paper, we propose a new probabilistic mixture model adapted to generalized multipartite networks. This model assumes that each functional group is partitioned into clusters gathering nodes sharing the same connection behavior on the basis of all the networks they are involved in. The clusterings are introduced through latent variables resulting into a mixture model on edges.

Inferring the parameters of such a model with latent variables is a complex task, since the likelihood can not be computed in a close form. We resort to a Variational Expectation-Maximization (VEM) algorithm to estimate the parameters and supplying the multi-clustering of the nodes. The crucial task of estimating the numbers of clusters in each functional group is tackled though a model selection procedure based on an adapted penalized likelihood criterion namely Integrated Classification Likelihood (ICL).

Related works Some extensions of standard descriptive tools -such as community detection-to more complex networks have been proposed in the literature (see for instance [START_REF] Yang | Community-affiliation graph model for overlapping network community detection[END_REF][START_REF] Gaskó | Community structure detection in multipartite networks: A new fitness measure[END_REF]. These metrics or descriptive tools have the drawback to only look for specific types of structures, contrary to clustering methods based on probabilistic mixture models. To the best of our knowledge, no such a probabilistic model has been proposed for generalized multipartite network. However, latent block models do exist for other kind of multi-layer networks. When the network is multiplex (i.e. when several types of interactions are studied on the same set of nodes), the SBM has been extended in [START_REF] Kéfi | How structured is the entangled bank? the surprisingly simple organization of multiplex ecological networks leads to increased persistence and resilience[END_REF]; [START_REF] Barbillon | Stochastic block models for multiplex networks: An application to a multilevel network of researchers[END_REF] by assuming a multivariate Bernoulli distribution to model the set of multiple edges between two individuals. Extending the SBM to a dynamic context where several snapshots of the networks at different time points are available, was proposed in numerous papers. A main difference between these papers is whether the evolution is considered as been discrete [START_REF] Xu | Dynamic stochastic blockmodels for timeevolving social networks[END_REF][START_REF] Matias | Statistical clustering of temporal networks through a dynamic stochastic block model[END_REF] or continuous [START_REF] Dubois | Stochastic blockmodeling of relational event dynamics[END_REF][START_REF] Xin | A continuous-time stochastic block model for basketball networks[END_REF] in time.

Outline of the paper Section 2 is dedicated to the introduction of notations supplying a flexible tool to describe generalized multipartite networks. We also provide a description of the two datasets of interest and illustrate the notations on these specific cases. The block model for multipartite networks is described in Section 3. In Section 4, the variational inference the model selection procedure are presented including the practical implementation of the algorithm. Numerical illustrations of the robustness of the inference method are provided in Section 5. Finally, the statistical analyses of the two datasets with discussion are presented in Section 6. Perspectives are discussed in the last section. Note that this document is accompanied with a appendix sections divided into three parts: in the two first sections, the derivations of the update formulas for the VEM algorithm of the ICL criterion are provided; finally, detailed estimates on the two datasets are given.

Notations and data

A generalized multipartite network can be seen as a collection of networks involving Q functional groups: each network may be simple, i.e. describing the relations inside a given functional groups or bipartite, i.e. describing the relations between individuals belonging to two different functional groups. Let n q be the number of individuals in the q-th functional group (q = 1, . . . , Q). We index the collection of networks by pairs of functional groups (q, q ) (q and q in {1, . . . , Q}). The set E denotes the list of pairs of functional groups for which we observe an interaction network. For any (q, q ) ∈ E, the interaction network is encoded in a matrix X qq such that X qq ii = 0 if there is an edge from unit i of functional group q to unit i of functional group q , 0 otherwise. Each network may be binary (X qq ii ∈ {0, 1}) or valued (X qq ii ∈ N or R). X qq may be symmetric if the relation is non-oriented, non-symmetric otherwise. X = X qq (q,q )∈E encodes the generalized multipartite network. For each network, S qq is an additional notation which refers to the list of all the possible interactions.

Illustration • The dataset 1 is issued from [START_REF] Dáttilo | Unravelling Darwin's entangled bank: architecture and robustness of mutualistic networks with multiple interaction types[END_REF]. This ecological network gathers mutualistic relations between plants and pollinators, plants and ants, and plants and frugivorous birds, resulting into Q = 4 func-tional groups, namely plants (q = 1), pollinators (q = 2), ants (q = 3) and birds (q = 4 and E = {(1, 2), (1, 3), (1, 4)}. X 1q ii = 1 if the plant species i has been observed at least once in a mutualistic interaction with the animal species i of functional group q during the observation period, 0 otherwise. Data set is available at https://figshare.com/articles/Interaction_matrix_A_in_ which_elements_aij_1_represent_the_presence_of_an_interaction_between_ the_plant_species_i_and_the_animal_species_j_and_zero_for_no_observed_ interaction/3443210/2.

• The dataset 2 comes from [START_REF] Thomas | Effects of farmer social status and plant biocultural value on seed circulation networks in Vanuatu[END_REF] and [START_REF] Thomas | A networkbased method to detect patterns of local crop biodiversity: validation at the species and infra-species levels[END_REF]. They collected the oriented network of seed circulation between farmers -resulting in a non-symmetric adjacency matrix -and the crop species grown by the farmers, resulting in an incidence matrix. Noting q = 1 for the farmers and q = 2 for crop species we get E = {(1, 1), (1, 2)}. X 11 ii = 1 if farmer i gives seeds to farmer i (oriented relation), 0 otherwise and X 12 ij = 1 if farmer i cultivates crop species j, 0 otherwise.

A block model for generalized multipartite networks

In order to account for heterogeneity among individuals, we propose a mixture model explicitly describing the way edges connect nodes in the various networks. We assume that each functional group q is divided into K q blocks (or equivalently clusters). ∀q ∈ {1, . . . , Q} and ∀i ∈ {1, . . . , n q }, let Z q i be the latent random variable such that Z q i = k if individual i of functional group q belongs to cluster k. The random variables Z q i 's are assumed to be independent and such that: ∀k ∈ {1, . . . , K q }, ∀q ∈ {1, . . . , Q}, ∀i ∈ {1, . . . , n q }:

P(Z q i = k) = π q k , (1)
with π q k ∈ [0, 1], ∀k = 1, . . . , K q and Kq k=1 π q k = 1, ∀q ∈ {1, . . . , Q}. We set Z q = (Z q i ) i∈{1,...,nq} , Z = (Z q ) q∈{1,...,Q} and π = (π q k ) k∈{1,...,Kq},q∈{1,...,Q} . Then, the nodes connect as follows: for any (q, q ) ∈ E, ∀(i, i ) ∈ S qq ,

X qq ii |{Z q i = k, Z q i = k } ∼ ind F qq (α qq kk ). (2)
F qq depends on the relation represented in X qq . We assume that for any (q, q ), F qq is either the Bernoulli distribution if the relation in X qq is binary, or the Poisson distribution if X qq ij is a counting or the Gaussian distribution if X qq encodes a continuous strenght of interaction. Equations (1) and (2) define the so-called Multipartite Block Model (MBM).

Remark 1 Our model is a generalization of the SBM and the LBM. Indeed, Equations (1)-(2) reduce to the SBM if E = {(1, 1)} and to the LBM if E = {(1, 2)}. Our extension assumes that the latent structures Z are shared among the X qq i.e. if a functional group q is at stake in several X qq , the same Z q impacts the distributions of the corresponding interaction matrices. In other words, the clusters gather individuals sharing the same properties of connection in the collection of networks. Obviously, if each functional group appears in only one element of E, the MBM reduces to independent SBMs or LBMs. In terms of probabilistic dependence, conditionally on Z, the X qq ii are independent. However, Z being latent, their marginalization introduces dependence not only between the X qq ii but also between the matrices X qq (q,q )∈E . As a consequence, the clustering variables (Z q i ) q∈{1,...,Q},i ∈{1,...,nq} are dependent once conditioned by X.

For a given vector K = (K 1 , . . . , K Q ), let θ K = (α, π) be the unknown parameters of the MBM where α = α qq kk (k,k )∈A qq ,(q,q )∈E are the connection parameters, with α qq kk ∈ Γ qq ⊂ R d qq . We have to perform two inference tasks: first, for a given vector K, estimating θ K and Z, second, selecting the right K.

Parameter inference and model selection

Let (X, Z; θ K ) denote the complete likelihood of the observations X and the latent variables Z for parameter θ K . Equations ( 1) and ( 2) lead to:

log c (X, Z; θ K ) = Q q=1 nq i=1 Kq k=1 1 {Z q i =k} log(π q k ) + (q,q )∈E (i,i )∈S qq (k,k )∈A qq 1 {Z q i =k,Z q i =k } f qq (X qq
ii , α qq kk )

(3) where f qq is the log-density of F qq . If q = q , A qq = {1, . . . , K q }×{1, . . . , K q }. If the interaction are oriented, then A qq = {1, . . . , K q } 2 ; otherwise

A qq = {(k, k ) ∈ {1, . . . , K q }| k ≤ k } 2 . Z being latent variables, the likelihood (X; θ K ) is obtained by integrating c (X, Z; θ K ) over all the possible values of Z denoted Z = {(z q i ) i∈{1,...,nq},q∈{1,...,Q} |z q i ∈ {1, . . . , K q }}: log (X; θ K ) = log Z∈Z c (X, Z; θ K ). ( 4
)
The summation over Z in Equation ( 4) becomes quickly computationally intractable when K q 's and/or the n q 's increase. Moreover, we are interested in inferring the clusterings Z.

Variational EM algorithm

In such models with latent variables, the EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] is a standard tool to maximise the likelihood, taking advantage of the simple form of log c (X, Z; θ K ). However, when conditioned by the observations X, the (Z q i ) are not independent (see Remark 1). As a consequence, the E step of the EM algorithm -which consists in the integration of log c (X, Z; θ K ) against P(Z|X; θ K )-has no explicit expression. In the context of the SBM and the LBM, the variational version of the EM (VEM) algorithm has proved to be a powerful tool for maximum likelihood inference (see [START_REF] Govaert | Block clustering with bernoulli mixture models: Comparison of different approaches[END_REF][START_REF] Daudin | A mixture model for random graphs[END_REF]. In the VEM, the problem of the E-step is tackled by replacing P(Z|X; θ K ) by an approximation R X,τ sought among distributions enforcing independence between the Z q i 's:

R X,τ (Z) = Q q=1 nq i=1 (τ q ik ) 1 Z q i =k , where τ q ik = P R X,τ (Z q i = k). At the VE-step, τ is chosen such that the Kullback-Leibler divergence KL[R τ ,X , P(•|X; θ K )] is minimized. The M-step is : θ K = arg max θ E R τ ,X [log c (X, Z; θ)] .
Iterating steps (VE) and (M) leads to the maximisation of a lower boud of the likelihood:

I θ (R τ ,X ) = log (X; θ) -KL[R τ ,X , P(•|X; θ)].
For the MBM, at iteration (t), the algorithm VEM is as follows:

VEM for MBM

• VE Step. Find (τ q ik ) i,k,q solving τ q ik ∝ exp log π q k +   q ∈Eq n q i =1 K q k =1 f qq (X qq ii , α qq kk )τ q i k   +1 (q,q)∈E j|(i,j)∈S qq Kq k =1
f qq (X qq ij , α qq kk )τ q jk + 1 (q,q)∈E 1 (i,i)∈S qq f qq (X qq ii , α qq kk ) .

• M Step. Update of the parameters:

π q k = 1 n q nq i=1 τ q ik , α qq kk = (i,i )∈S qq X qq ii τ q ik τ q i k (i,i )∈S qq τ q ik τ q i k . ( 5
)
Remark 2 Formula (5) on α qq kk is valid for the expectations of the Bernoulli, the Poisson and the Gaussian distributions. If F qq is a Gaussian distribution of variance v qq kk , then:

v qq kk = (i,i )∈S qq (X qq ii ) 2 τ q ik τ q i k (i,i )∈S qq τ q ik τ q i k -(α qq kk ) 2 .
The details of the VEM algorithm are provided in the Section A. As for any EMtype algorithm, its initialization is critical and will be discussed later. ( θK , τ ) denotes the VEM estimates. A by-product of this algorithm is an approximation of the conditional distribution P(Z|X; θK ) by R X,τ . Z is thus estimated by:

Ẑq i = arg max k∈{1,...,Kq} τ q ik . (6)
The consistency of the VEM estimates has been established for the SBM by [START_REF] Bickel | Asymptotic normality of maximum likelihood and its variational approximation for stochastic blockmodels[END_REF] and the LBM by [START_REF] Brault | Consistency and Asymptotic Normality of Latent Blocks Model Estimators[END_REF] while [START_REF] Mariadassou | Convergence of the groups posterior distribution in latent or stochastic block models[END_REF] study the behavior of P(Z|X; θK ) for the same models.

A penalized likelihood criterion

In practice, the number of clusters K = (K 1 , . . . , K Q ) is unknown and should be estimated. We adopt a model selection strategy where a model M K refers to the MBM with K clusters. Among many classical model selection criterion such as AIC, BIC and their variants, [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF] proposed the Integrated Classification Likelihood (ICL). ICL has proven its capacity to outline the clustering structure in networks in [START_REF] Daudin | A mixture model for random graphs[END_REF]

ICL(M K ) = log c (X, Ẑ; θK ) -pen(M K ) (7)
where

pen(M K ) = 1 2    Q q=1 (K q -1) log(n q ) +   (q,q )∈E d qq |A qq |   log   (q,q )∈E |S qq |     
and Ẑ has been defined in Equation ( 6). The better model is chosen as K = arg max K ICL(M K ). The penalization term pen(M K ) is made up of two parts:

Q q=1 (K q -1) log(n q ) corresponds to the clustering distribution and involves the number of nodes, while the other term depends on the size of α qq kk and the numbers of possible edges in each network |S qq |. The derivation of the ICL criterion is detailed in the Section B.

Practical algorithm

The practical choice of the model and the estimation of its parameters are computational intensive tasks. Assume that K q ∈ {1, . . . , K q }, then, ideally, we should compare Q q=1 K q models through the ICL criterion. For each model, the VEM algorithm has to be run starting from a large number of initialization points chosen carefully (due to its sensitivity to the starting point), resulting in an unreasonable computational cost. Instead, we propose to adopt a stepwise strategy, resulting in a faster exploration of the model space combined with efficient initializations of the VEM algorithm. The procedure we suggest is the following one.

Given a current model

M (m) = M(K (m) 1 , . . . , K (m)
Q ), the m-th iteration of the procedure is written as follows.

Model selection strategy for MBM

• Split proposals. For any q such that K (m) q < K q , consider the model

M (m+1)q + = M(K (m) 1 , . . . , K (m) q + 1, . . . , K (m) Q ). • Propose K (m) q
initializations by splitting any of the K (m) q current clusters into two clusters.

• From each of the K (m) q initialization points, run the VEM algorithm and keep the better variational estimate θM (m+1)q + .

• Compute the corresponding ICL(M (m+1)q +

) from formula (7).

• Merge proposals. For any q such that K (m) q > 1, consider the model

M (m+1)q - = M(K (m) 1 , . . . , K (m) q -1, . . . , K (m) Q ). • Propose K (m) q (K (m) q
-1)/2 initializations by merging any pairs of clusters among the K (m) q clusters.

• From each initialization point, run the VEM algorithm and keep the better variational estimate θM (m+1)q -

• Compute the corresponding ICL(M (m+1)q -

).

• Set M (m+1) = arg max

M (m) ICL(M) where M (m) = {M (m) }∪ q∈{1,...Q} {M (m+1)q + }∪ {M (m+1)q - }. If M (m+1) = M (m) iterate, otherwise stop.
Initializing the VEM from the clusters obtained on a smaller or larger model is much more efficient than other strategies such as random initialization or spectral clustering. Note that the tasks at each iteration can be parallelized.

Illustration on simulated datasets

We illustrate the efficiency of our inference procedure on simulated data. Simulation parameters We generate datasets mimicking the two real datasets Simulations are implemented as vignette of the R-package available on Github (https://github.com/Demiperimetre/GREMLIN) in order to ensure reproducibility of the results.

• Scenario 1: We set Q = 4 functional groups of respective sizes n = (141, 173, 46, 30), E = {(1, 2), (1, 3), (1, 4)} and K = (7, 2, 2, 1). The simulation parameters are the following: π 1 = (0.3651, 0.1270, 0.1190, 0.1460, 0.0842, 0.0794, 0.0794), π 2 = (0.1, 0.9), π 3 = (0.1, 0.9),

α 12 =          0.
         , α 14 =          0.0013 0 0.0753 0 0.0163 0.5108 0         
.

(8) • Scenario 2: Mimicking our second dataset, we set Q = 2 functional groups such that n 1 = 30, n 2 = 37, E = {(1, 1), (1, 2)} and K = (3, 2). The simulation parameters are: π 1 = (0.31, 0.42, 0.27), π 2 = (0.65, 0.35). On each simulated dataset, we run the algorithm described in Section 4.3 with K q = 10, for all q. We start the stepwize algorithm on K q = 1 and on the clusterings obtained by inferring the LBM separately on each matrix X qq . For each dataset, it takes a few minutes for the algorithm to converge.

Results for Scenario 1 Among the 100 simulated datasets, the true numbers of clusters are exactly recovered for 73 datasets as can be seen in Table 1. Three other estimated configurations are observed as detailed in Table 1, each of them corresponding to Kq = K q -1. A scrutiny shows that cases where K 3 is estimated to 1 instead of 2 correspond to datasets where one simulated cluster (meaning the simulated Z 3 ) is reduced to 1 or 2 individuals. We measure the ability of the procedure to recover the clusters by computing the Adjusted Rand Index (ARI by [START_REF] Hubert | Comparing partitions[END_REF] between Ẑ1 and the simulated Z 1 . ARI compares two clusterings with a correction for chance. It is close to 0 when the two clusterings are independent and equals 1 when the clusterings are identical (up to label switching). The boxplots of the ARI's are plotted in Figure 2a. All the ARI are high (> 0.7), even when K1 = 6, which means that the clusters are always globally recovered and the structure of the multipartite network is well discovered. We also investigate the quality of the parameters estimation. For the 73 datasets such that K = (7, 2, 2, 1), we compare the estimated parameters α to the simulation parameters (after a relabeling of the clusters to match the true clustering, if required). The results in terms of biais and Root Mean Square Error (RMSE) are reported in Table 2. We observe that the parameters are estimated without noticeable biases and with small RMSE. In particular, the null α qq kk 's are always recoverered.

Results for Scenario 2 The same analysis is done for Scenario 2. K = K = (3, 2) for 82 simulated datasets over 100. The unique alternative observed configuration is K = (2, 2). Here again, in order to assess the ability of the method to recover the clusters, we compute the ARI for Z 1 and Z 2 . These quantities are reported in Figure 2b for the two functional groups (left au right) and for K1 = 3 (triangles) and K2 = 2 (circles). As can be noticed, the clusterings are well recovered when the number of clusters are exactly recovered (triangle), with a noticeable number of cases where the clusters are exactly recovered (ARI= 1). As expected, we observe a lower ARI when K1 = 2 (circle), but with globally high values.

We also pay attention to the estimation of the parameters of connection α qq kk . Biais and RMSE are reported in Table 3. Once again, the biais and RMSE are small. Inference We run the procedure described in Subsection 4.3 starting from several automatically chosen initial points K (0) , with numbers of clusters bounded between 1 and 10. The computational time on an Intel Xeon(R) CPU E5-1650 v3 3.50GHz x12 using 6 cores is less than 10 minutes. The ICL criterion selects 7 clusters of plants, 2 clusters of pollinators, 1 cluster of birds and 2 clusters of ants. The estimated parameters are reported in Tables 5 and6 in the supplementary material. They are really similar to the ones used for the simulation study and provided in Equation 8. The resulting mesoscopic view of the multipartite network is plotted in Figure 3. Discussion From Figure 3, we conclude that plants of Clusters 7 and 2 only interact with ants (plants of Cluster 7 attract more ant species belonging to Cluster 1). The plants of Clusters 3 and 6 are only in interaction with birds, the difference between the two clusters being due to the strength of the connection. The difference between the two clusters of pollinators derives from the existence of Cluster 1 of plants.

In order to illustrate the contribution of our method, we also analyze each bipartite network separately (using an LBM) and compare the results in terms of clustering. The clusterings are compared through the Adjusted Rand Index (ARI): if ARI = 1 then the clusterings are equal (up to a label switching transformation). The ARIs are given in Table 4.

Using standard LBM, we obtain 2 clusters of ants, 1 cluster of birds and 3 clusters of pollinators. The clusterings of ants and birds are not modified by the ecosystemic approach, their ARI being equal to 1. The clustering of the pollinators is slightly modified, going from 3 clusters to 2 clusters but the additional block only contains few individuals, thus leading to a high ARI. Since the plants functional group is involved in the three bipartite networks, we obtain 3 clusterings when analyzing them separately. These three clusterings are -as expected-very different from our MBM clustering (the ARIs being respectively equal to 0.118, 0.415 and 0.163, see Table 4). When aiming at proposing a clustering taking into account the 3 bipartite networks, one may adopt a naive strategy by combining (by intersection) these three clusterings. We then obtain 12 blocks of plants and the ARI with the MBM clustering is 0.617. However, this number of clusters ( 12) is too large with respect to the model selection criterion. Our MBM clustering is a parsimonious trade-off.

Finally, Figure 3 highlights the central role played by the pollinators in the complete network also involving ants and seed dispersal birds. Our ecosystemic approach unables us to unravel such a central position of pollinators in general with respect to ants and seed dispersal birds, which would have been impossible when performing separate analysis.

Seed circulation and crop species inventory

Dataset Seed circulation among farmers is a key process that shapes crop diversity [START_REF] Coomes | Farmer seed networks make a limited contribution to agriculture ? four common misconceptions[END_REF][START_REF] Pautasso | Seed exchange networks for agrobiodiversity conservation. a review[END_REF]. Data on seed circulation within a community of first-generation migrants (30 farmers) were collected by a field survey in the island of Vanua Lava in the South Pacific archipelago nation of Vanuatu. A farmer is considered as a giver for another farmer if he/she has given at least one crop since they arrived in the new settlement site in Vanua Lava. It results in a connected and directed network of seed circulation. Besides the circulation network, inventory data for each farmer were also collected. They consist in the list of crop landraces they grow. This was aggregated at the species level, leading to 37 crop species. These inventory data can be seen as a bipartite network. The seed circulation data were analyzed in [START_REF] Thomas | Effects of farmer social status and plant biocultural value on seed circulation networks in Vanuatu[END_REF] and the inventory data were analyzed in a meta-analysis in [START_REF] Thomas | A networkbased method to detect patterns of local crop biodiversity: validation at the species and infra-species levels[END_REF].

On the basis on the joint modeling we propose in this paper, we aim to provide a clustering on farmers and crop species on the basis of the seed circulation network and the inventory bipartite network.

Inference Three clusters of farmers and two clusters of crop species were selected. The inferred parameters are given in Table 7 and a a mesoscopic view is displayed in Figure 4. The clusters were renumbered to make them correspond to the probability of connection: the larger cluster number, the larger marginal probability of connection for the members of the cluster.

Discussion The discovered clusters are straightforwardly interpretable: Cluster 3 gathers farmers who circulate seeds within the cluster and give to the two other clusters, Cluster 2 circulates seeds within the cluster contrary to Cluster 1 who only receives from Cluster 3; the two clusters of crop species are Cluster 2 with more common crop species and Cluster 1 with other species. Clusters 3 and 2 of farmers grow crop species from Clusters 1 and 2 whereas Cluster 1 of farmers grows mainly crop species from Cluster 2. It turns out that Cluster 3 gathers mainly the first migrants and Cluster 1 the last migrants. The pattern of connection is then explained by the fact that first migrants helped the others to settle by providing seeds. Moreover, the first migrants had more time to collect more crop species to grow. In order to compare the clusterings obtained by the MBM and the ones obtained from the circulation network (clustering on farmers from SBM) and the inventory network (clusterings of farmers and plants with a LBM), we compute the ARI. The clusterings on crop species remain quite close since their ARI is equal to 0.891. However, the clusterings on farmers are quite different (ARI smaller than 0.3), indeed the MBM shall make a trade off between the circulation and the inventory for farmers. To ease the comparison between clusterings on farmers, the same renumbering rule was applied for all the different clusterings so that the larger cluster numbers correspond to larger marginal probability of connection. Figure 5 is an alluvial plot which compares the three obtained clusterings of farmers. It shows how the trade-off is made between the two stand-alone clusterings in the MBM clustering. It appears quite obvious that Cluster 1 given by the MBM gathers only farmers from Cluster 1 in the seed circulation network and from Clusters 1 and 2 from the inventory data since this cluster aggregates farmers with fewer connections and who grow less crop species than the others. The same kind of observation can be made for Cluster 3 given by the MBM which aggregates farmers who were in the cluster of the most connected farmers and in the two clusters of farmers who grow more seed.

Discussion

In this paper, we proposed an extension of LBM and SBM which can handle multipartite networks, resulting in the so-called MBM. Multipartite networks encompass a lot of various situations such as the two examples dealt with in the paper. Besides, MBM can be also useful for many other contexts with different multipartite structures. The main limiting factor of the parameter inference and model selection methods we propose, is the size of networks. Indeed, the inference algorithm is suitable for networks up to 1000 nodes in order to keep computational time reasonable. If willing at studying larger networks, one should develop adapted inference algorithms. Several extensions can be thought of with no additional significant difficulty. For instance, covariates can be taken into account by writing P(X qq ii = 1|Z q i = k, Z q i = k ) = φ α qq kk + y qq ii β where y qq ii are the covariates describing the pair of individuals (i, i ). A more challenging extension could be to include extra parameters in the model to account for degree heterogeneity in the vein of the degree-corrected SBM [START_REF] Karrer | Stochastic blockmodels and community structure in networks[END_REF]. However, this should lead to cluster individuals independently of their degrees which is not always desirable. For instance in the seed circulation example, the in-degrees and Figure 5: Alluvial plot comparing the clustering on farmers obtained from an SBM on the circulation network (clExc), an LBM on the inventory network (clInv) and the MBM (clMBM) on both networks. The cluster numbers are related with the probability of connection, the larger cluster number, the larger marginal probability of connection (between farmers for clExc, between farmers and crop species for clInv). out-degrees of farmers are key factors to explain their social structure. In the interaction networks in ecology, distinguishing generalist (high degree) from specialist (low degree) species is a part of the expected outcomes of the clusterings. More recent works [START_REF] Zhu | Oriented and degree-generated block models: generating and inferring communities with inhomogeneous degree distributions[END_REF] propose to use node degrees and edge orientations to help the clustering while tolerating heavy-tailed distribution of degrees. Incorporating these variants may be of interest in the MBM framework and could lead to further developments. In a more general perspective, the study of ecological or sociological interactions supplies a wide variety of complex networks such as multilevel networks or multi-layer networks [START_REF] Pilosof | The multilayer nature of ecological networks[END_REF][START_REF] Lazega | Multilevel Network Analysis for the Social Sciences: Theory, Methods and Applications[END_REF]. Some of them can be treated as multipartite networks and then by a MBM (possibly by incorporating some extensions discussed above). The others require the development of suited models which might also rely on a latent variable modeling. They will be the subject of future works.

A Variational EM for MBM: details

The variational version of the EM algorithm maximizes a lower bound of the observed likelihood. More precisely, let R τ ,X be any probability distribution on Z, we define I θ (R τ ,X ) as:

I θ (R τ ,X ) = log (X; θ) -KL[R τ ,X , P(•|X; θ)] (9) = E R τ ,X [log c (X, Z; θ)] + H(R τ ,X ) (10) ≤ log (X; θ) ( 11
)
where KL is the Kullback-Leibler divergence and H(R τ ,X ) is the entropy of R τ ,X . The inequality in (11) derives from the positivity of the KL divergence.

The equality I θ (R τ ,X ) = log (X; θ) holds iff R τ ,X (Z) = P(Z|X; θ) which results in the classical identity of the EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF].

The principle of the VEM algorithm is to choose R τ ,X in a family of distributions P parametrized by τ , such that the conditional expectation E R τ ,X [log c (X, Z; θ)] can be computed explicitly. Iteration (t) of VEM consists in the two following steps:

• M Step θ (t) = arg max θ I θ (R τ (t-1) ,X ) = arg max θ E R τ ,X [log c (X, Z; θ)] • VE Step τ (t) = arg max τ E R τ ,X [log c (X, Z; θ)] + H(R τ ,X ) = arg min τ KL[R τ ,X , P(•|X; θ (t) )] .
The variational EM generates a sequence (θ (t) , τ (t) ) t≥0 increasing the lower bound I θ (R τ ,X ) of the likelihood log (X; θ).

Choice of P

The key point of the procedure is the choice of P, making the calculus tractable but rich enough to be a good approximation of the true conditional distribution P(Z|X; θ). Following [START_REF] Govaert | Block clustering with bernoulli mixture models: Comparison of different approaches[END_REF] and [START_REF] Daudin | A mixture model for random graphs[END_REF], we adopt the mean-field strategy [START_REF] Jaakkola | Tutorial on variational approximation methods[END_REF] and chose P as:

P = R τ |R τ (Z) = Q q=1 nq i=1 h Kq (Z q i ; τ q i ) ,
where h Kq (•; ξ) is the density of a 1 trial -multinomial distribution of parameter ξ ∈ T Kq , i.e. R X,τ is such that the latent variables Z are independent and P R τ ,X (Z q i = k) = τ q ik with Kq k=1 τ q ik = 1, ∀q ∈ {1, . . . , Q}, ∀i ∈ {1, . . . , n q }.

From this particular choice of P, we derive the VE-step and M-step for the MBM in what follows.

Explicit expression for I θ (R τ ,X ) Using the expression of the complete loglikelihood given in Equation (3), we obtain:

I θ (R τ ,X ) = -q,i,k τ q ik log τ q ik + q,i,k τ q ik log(π q k ) + (q,q ) (i,i ) (k,k ) E 1 Z q i =k 1 Z q i =k f qq (X qq ii , α qq kk ) (12) f qq (X qq ii , α qq kk ) =      X qq ii log(α qq kk ) + (1 -X qq ii ) log(1 -α qq kk ) for binary networks -α qq kk + X qq ii log(α qq kk ) -log(X qq ii !) for weigthed Poisson networks. log φ(X qq ii , µ qq kk , σ qq kk ) for Gaussian networks. ( 13
) E 1 Z q i =k 1 Z q i =k
has to be carefully calculated, when i = i . To that purpose, let us introduce the following notations :

• ∀q, E q = {q ∈ {1, . . . , Q} | q = q and (q, q ) ∈ E} . E q is the set of incidence matrices involving the functional group q.

• ∀(q, q ) ∈ E, ∀i ∈ [1, n q } we define :

S qq i = {i ∈ {1, . . . , n q } | i = i, (i, i ) ∈ S qq } .
Using these notations we detail the expression of I θ (R τ ,X ).

(q,q ) (i,i ) (k,k ) E 1 Z q i =k 1 Z q i =k f qq (X qq ii , α qq kk ) = q q ∈Eq (i,i ) (k,k ) τ q ik τ q i k f qq (X qq ii , α qq kk ) + 1 (q,q)∈E (i,i ) (k,k ) E 1 Z q i =k 1 Z q i =k f qq (X qq ii , α qq kk ) = q q ∈Eq (i,i ) (k,k ) τ q ik τ q i k f qq (X qq ii , α qq kk ) + 1 (q,q)∈E i i ∈S qq i (k,k ) E 1 Z q i =k 1 Z q i =k f qq (X qq ii , α qq kk ) +1 (q,q)∈E i (k,k ) E 1 Z q i =k 1 Z q i =k =0 if k =k f qq (X qq ii , α qq kk ) = q q ∈Eq (i,i ) (k,k ) τ q ik τ q i k f qq (X qq ii , α qq kk ) + 1 (q,q)∈E i i ∈S qq i (k,k ) τ q ik τ q i k f qq (X qq ii , α qq kk ) +1 (q,q)∈E i | (i,i)∈S qq k E 1 2 Z q i =k =1 Z q i =k f qq (X qq ii , α qq kk ) .
As a consequence, we get:

I θ (R τ ,X ) = -q,i,k τ q ik log τ q ik + q,i,k τ q ik log(π q k ) + q q ∈Eq (i,i ) (k,k ) τ q ik τ q i k f qq (X qq ii , α qq kk ) +1 (q,q)∈E i i ∈S qq i (k,k ) τ q ik τ q i k f qq (X qq ii , α qq kk ) +1 (q,q)∈E i | (i,i)∈S qq k τ q ik f qq (X qq ii , α qq kk ) . (14)
Optimization of I θ (R τ ,X ) with respect to τ , (θ being fixed)

For a fixed θ, we need to find τ such that ∀q ∈ {1, . . . , Q}, ∀k ∈ {1, . . . , K q }, ∀i ∈ {1, . . . , n q }:

∂ ∂τ q ik   I θ (R τ ,X ) + Q q =1 n q j=1 λ q j   Kq k =1 τ jk -1     = 0 (15)
where (λ q j ) 1≤q ≤Q,1≤j≤n q are the Lagrange multipliers. Combining Equations ( 12) and ( 14), we get:

0 = -(1 + log τ q ik ) + log π q k +   q ∈Eq n q i =1 K q k =1 f qq (X qq ii , α qq kk )τ q i k   +1 (q,q)∈E j∈S qq i Kq k =1 f qq (X qq ij , α qq kk )τ q jk (16) +1 (q,q)∈E 1 (i,i)∈S qq f qq (X qq ii , α qq kk ) +λ q i .
This system has no explicit solution but can be solved numerically using a fixed point strategy as in [START_REF] Daudin | A mixture model for random graphs[END_REF].

Optimization of I θ (R τ ,X ) with respect to θ, τ being fixed.

We have to compute the derivatives of I θ (R τ ,X ) with respect to θ, the variational parameters τ being fixed. We thus obtain: ∀(q, q ) ∈ E, ∀(k, k ) ∈ {1, . . . , K q } × {1, . . . , K q }: α qq kk = (i,i )∈S qq X qq ii τ q ik τ q i k (i,i )∈S qq τ q ik τ q i k 20 and ∀q ∈ {1, . . . , Q}, ∀k = {1, . . . , K q }:

π q k = 1 n q nq i=1 τ q ik .

Clustering and initializations

We denote by θ and τ the resulting estimates. The estimated clustering is the maximum a posteriori (MAP) estimate: ∀q ∈ {1, . . . , Q}, ∀i ∈ {1, . . . , n q }, Ẑq i = arg max k∈{1,...,Kq} τ q ik .

B Derivation of the ICL criterion

As exposed in Section 4.2 of the paper, we resort to the ICL criterion to perform model selection. The ICL is an asymptotic approximation of the integrated marginal complete likelihood. We supply here the details of the calculations.

Explicit expression of the marginal complete likelihood The ICL being a based on a Bayesian model selection criterion, we first set the following prior distribution:

α qq kk ∼ B(a, a) and (π q 1 , . . . π q Kq ) ∼ Dir(b, . . . , b). The prior on θ being such that π(θ) = π(α)π(π) we obtain: log m c (X, Z; M K ) = log f (X|Z; α)π(α)dα + log f (Z; π)π(π)dπ .

Taking advantage of the conditional independences in the model defined by Equations ( 1) and (2) combined with the independence of the parameters in the prior distribution, we can decompose log m c into the following sum: log m c (X, Z; M K ) = (q,q )∈E log f (X qq |Z q , Z q ; (α qq ))π(α qq )dα qq + Q q=1 log f (Z q ; π q )π(π q )dπ q .

Using the fact that f (Z q ; π q ) = Kq k=1 (π q k ) N q k with

N q k = nq i=1 1 Z q i =k (18)
and the conjugacy of the Dirichlet prior distribution, we easily deduce the following formula: f (Z q ; π q )π(π q )dπ q = Γ(bK q ) Γ(b) Kq Kq k=1 Γ(N q k + b) Γ(n q + bK q ) where Γ is the Gamma function. Now, we can reformulate f (X qq |Z q , Z q ; α qq ) as: f (X qq |Z q , Z q ; α qq ) = (i,i ,k,k ) (α qq kk )

X qq ii 1 Z q i =k 1

Z q i
=k (1 -α qq kk )

(1-X qq ii )1 Z q i =k 1

Z q i =k = Kq,K q k,k =1
(α qq kk ) S qq kk (1 -α qq kk ) N qq kk -S qq kk with S qq kk = (i,i )∈S qq X qq ii 1 Z q i =k 1 Z q i =k N qq kk = (i,i )∈S qq 1 Z q i =k 1 Z q i =k .

(

) 19
With the beta prior distribution on each α qq kk , we get:

f (X qq |Z q , Z q ; (α qq ))π(α qq )dα qq = Kq,K q k,k =1 Γ(2a) Γ(a) 2 Γ(a + S qq kk )Γ(a + N qq kk -S qq kk ) Γ(2a + N qq kk )

. log Γ(bK q ) -K q log(b) -log Γ(n q + bK q )

+ Q q=1   Kq k=1 log Γ(N q k + b)  
where N q k has been defined in Equation ( 18) and S qq kk and N qq kk in Equation ( 19) Asymptotic approximation Using the same arguments as in [START_REF] Daudin | A mixture model for random graphs[END_REF] and [START_REF] Brault | Estimation et sélection de modèle pour le modèle des blocs latents[END_REF], we obtain the following asymptotic approximation. Assume that ∀q ∈ {1, . . . , Q}, n q → ∞, then : log m c (X, Z; M K ) = max

θ K ∈Θ K log c (X, Z; θ K ) -pen(M K )
where pen(M K ) = 1 2 Q q=1 (K q -1) log(n q )

Figure 1 :

 1 Figure 1: Illustrations of bipartite (left), multipartite (center) and generalized multipartite networks (right)

Figure 2 :

 2 Figure 2: (a) : Simulation scenario 1 : distribution of the ARI for q = 1 when K1 = 6 (under-estimation) on left and K1 = 7 (exact estimation) on right. (b) Simulation scenario 2 : Boxplots of the ARI for Ẑ1 (on left) and Ẑ2 (on right) with K1 = 2 (circle) and K1 = 3 (triangle).

Figure 3 :

 3 Figure 3: Mesoscopic view of the interconnected network. The size of the nodes are proportional to the size of the clusters and the width of the edges are proportional to the probability of connection between/within clusters. Edges corresponding to probabilities of connection lower than 0.01 are not plotted.

Figure 4 :

 4 Figure 4: Summarized network provided by the MBM. Nodes correspond to the clusters detected by the MBM: clusters of farmers are in red and clusters of crop species are in green. Size of a node is proportional to the number of farmers or crop species belonging to the considered cluster. The width of the edges are proportional to the probability of connection between/within clusters. The probability of connection below 0.2 are not plotted.

  a, a) denotes the Beta distribution and Dir(b, . . . , b) is the Dirichlet distribution. By definition, log m c (X, Z;M K ) = log c (X, Z; θ K )π(θ K ; M K )dθ K .

Finally

Table 1 :

 1 Results for Scenario 1: estimated numbers of clusters ( K1 , K2 , K3 , K4 ). The simulated number of clusters is equal to (7, 2, 2, 1). studied in Section 6, i.e. we generate 100 datasets under the following two scenarios.

		K
	K1	K2	K3	K4 Nb of simulations
	7	2	2	1	73
	7	2	1	1	6
	7	1	2	1	3
	6	2	2	1	18

Table 2 :

 2 Simulation scenario 1 : Biais and RMSE for α qq kk .

			α 12				α 13		α 13
		α 12 k1	α 12 k2		α 13 k1		α 12 k2		α 14 k1
	Biais	RMSE Biais	RMSE	Biais	RMSE Biais	RMSE	Biais	RMSE
	-9e-04	0.0102 0	0.0011	0	4e-04	1e-04	6e-04	-1e-04	9e-04
	8e-04	0.0064 1e-04	2e-04	-0.0036 0.0557 0.0013	0.0085	0	0
	1e-04	7e-04	-1e-04	4e-04	0	0	0	2e-04	0.0026	0.0155
	2e-04	0.0194 4e-04	0.0038	0.0038	0.0515 -6e-04	0.0119	0	0
	0.0024 0.0281 -0.0023 0.0099	0	0	0	0	-0.0012 0.0077
	0	0	0	0	0	0	0	0	-0.0081 0.0323
	0	0	0	0	0.002	0.0553 -0.0039 0.0212	0	0

Table 3 :

 3 Simulation scenario 2 : Biais and RMSE for α qq kk .

	6 Applications
	6.1 Ecology: interactions plants / animals
	Dataset The dataset -compiled and conducted by Dáttilo et al. (2016) at Cen-
	tro de Investigaciones Costeras La Mancha (CICOLMA), located in Mexico-in-
	volves three general types of plant-animal mutualistic interaction: pollination,
	seed dispersal by frugivorous birds, and protective mutualisms between ants
	and plants with extrafloral nectaries. The dataset -which is one of the largest
	compiled so far with respect to species richness, number of interactions and
	sampling effort-includes n 1 = 141 plant species, n 2 = 173 pollinator species,
	n 3 = 46 frugivorous bird species and n 4 = 30 ant species, inducing 753 observed
	interactions of which 55% are plant-pollinator interactions, 17% are plant-birds
	interactions and 28% are plant-ant interactions.

Table 4 :

 4 Comparison of clusterings when the networks are jointly modeled by the MBM (denoted Full) and when the networks are considered apart as bipartite networks. Inter denotes the clustering obtained by intersecting the three clusterings on plants for each bipartite network. The selected number of clusters (in parenthesis) and the ARIs are provided.

		Full/Flovis Full/Ants Full/Birds Full/Inter
	Plants	(7/3) 0.118	(7/3) 0.415	(7/3) 0.163	(7/12) 0.617
	Flovis	(2/3) 0.997
	Ants		(2/2) 1.000
	Birds			(1/1) 1.000

Table 7 :

 7 Estimated parameters for MBM on the seed circulation dataset: πq k are in the first row and column, other rows and columns contain the estimates α qq kk . α qq kk identified by * are not plotted in Figure 4.

					Farmers	Species
				1	2	3	1	2
	Farmers	1 2 3	πq k	0.31 0.025 * 0.123 * 0.053 * 0.186 0.653 0.42 0.27 0.65 0.35 0.159 * 0.300 0.070 * 0.559 0.905 0.374 0.585 0.357 0.390 0.696

Acknowledgements

We thank Sophie Caillon (CEFE) for sharing the seed circulation and inventory data, Sophie Caillon and Mathieu Thomas (CIRAD) for related discussion concerning the analysis. Support of the ongoing research collaboration through the MIRES network was provided by the Institut National de la Recherche Agronomique (INRA). P. Barbillon and S. Donnet are funded by the ANR project Econet.

The first term comes from the application of the Stirling formula to the Gamma function when approximating f (Z q ; π q ). The second term comes from a BIC approximation of the part f (X qq |Z q , Z q ; α qq ). Obviously, the parameters of the prior distribution (a, b) disappear from the formula since asymptotically the importance of the prior distribution vanishes.


F Khellaf
email: feriel.khellaf@univ-grenoble-alpes.fr

R Clackdoyle

S Rit

L Desbat

Tiny changes in tomographic system matrices can cause large changes in reconstruction quality

This paper studies the impact of tiny changes in region-of-interest (ROI) tomography system matrices on the variance of the reconstructed ROI. In small-scale and medium-scale examples, the variance in the reconstructed ROI was estimated for different system matrices. The results revealed a striking and counterintuitive phenomenon: a tiny change in the system matrix can dramatically affect the variance of the ROI estimate. In one of our examples, a decrease of 0.1% in one element out of hundreds of thousands of the system matrix resulted in a systematic reduction of the variance inside the ROI, and by a factor of 5 to 10 for some pixels. Our results agree with a recently proven theorem about the ability of additional measurements to reduce the variance in ROI tomography.

Introduction

Computed tomography (CT) consists of reconstructing an image from a set of measurements

y = Sx + (1)
with S the m × n system matrix, x the n × 1 image vector, y the m × 1 measurement vector, and the m × 1 noise vector. While in classical tomography the whole image x can be reconstructed from the measurements, in region-of-interest (ROI) tomography, only a portion of the image can be uniquely and stably reconstructed [START_REF] Clackdoyle | Tomographic reconstruction in the 21st century[END_REF].

Naturally, as the measured data are noisy, the reconstructed image is not exactly equal to the true image. For both full and ROI tomography, the noise of the solution depends on the input noise vector , and on the nature of the system matrix S. In this work, we illustrate the impact that tiny variations in the system matrix S can have on the variance of the reconstructed image. Intuitively, we expect that, as long as a pixel can be reconstructed, its variance should vary continuously as a function of the matrix elements. We will show here that this is not always the case in ROI tomography, i.e. there can be a discontinuity in the variance of reconstructible pixels when the elements of the system matrix change slowly.

We note that our work is not a study of the ill-conditioning of tomographic matrices, i.e. the fact that small changes in the data cause a large variance of the reconstructed image. This problem has already been studied, e.g. for under-sampled tomography [START_REF] Jorgensen | Quantifying admissible undersampling for sparsity-exploiting iterative image reconstruction in X-ray CT[END_REF]. But in the context of ROI tomography, even ill-conditioned matrices can be used to reconstruct unique and stable ROI pixels, and the condition number is not indicative of the stability of the solution in the ROI. In this study, we are interested in how the variance changes with the matrix elements rather than in the value of the variance itself.

We study the variance for full reconstruction and ROI reconstruction using different examples, and show that our results agree with a recently published theorem [START_REF] Defrise | Do additional data improve region-of-interest reconstruction?[END_REF] about the impact of additional measurements on the ROI.

Full reconstruction

We illustrate here the fact that as the system matrix becomes more unstable (approaches a singular matrix), the quality of the reconstruction decreases (the image variance increases). Assuming = 0, the discrete model can be rewritten as

y = Sx. (2)
We also make the assumption that cov( ) = σ 2 I with I the m × m identity matrix. If the matrix S is square and invertible, then the solution is given by

x = S -1 y (3) and cov(x) = σ 2 (S T S) -1 . ( 4
)
The variance in the reconstructed image x is found from the diagonal elements of cov(x).

We provide a small-scale example of such a system to illustrate the typical outcome when a system matrix becomes unstable. In all our examples, we use a simple tomographic model of uniform square pixels of unit side length and of constant density, and of ideal line integrals where the contribution of each pixel to the line is equal to the line length through the pixel, scaled by the pixel intensity value. We first considered a 3 × 3 image and the following 9 × 9 system matrix

S =               0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 0 1 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 1 B 0 0 0 B 0 0 0 B 0 0 B 0 B 0 B 0 0 0 A 0 0 A 0 0 0 A 0 A A 0 0 0 0 0 0 αC 0 0 C 0 0 C 0 0               (5)
where, in the above matrix, the constants A, B, C are

A = √ 5 2 , B = √ 2, C = √ 10
3 , and with the parameter α belonging to (0, 1]. The measurement lines corresponding to this system matrix are shown in Figure 1. Each value of α corresponds to a different system matrix, where the last measured line has been translated slightly, thus changing a single entry in the matrix. Geometrically, α corresponds to the intersection of this line with the left hand edge of the image, with 1 -α being the distance from the upper left corner, measured downwards (see Figure 1). When α = 0, the last row becomes a linear combination of the others (the linear relationship is y 9 = y 1 +y 2 -y 3 -y 4 +y 8 /A), and the matrix becomes singular. Otherwise the system is invertible and reconstruction of all pixels is possible. As the matrix becomes singular, we expect that the reconstruction quality will get poorer when there is noise in the measurements. We computed the theoretical variance of reconstruction pixels according to equation ( 4), and Figure 2 shows that, indeed, the reconstructed variance explodes as the matrix becomes singular.

To confirm these analytical variance calculations, the variance was also estimated numerically from 10 5 reconstructions with Gaussian noise realizations added to the measurements. The 3 × 3 test image used to compute the numerical variance was a centered 2D Gaussian with a standard deviation of 2 and an amplitude of 5, and we set σ = 1. The numerical variance calculations are also shown on Figure 2 and they agree with the theoretical behavior: increased noise propagation in the reconstruction as the tomographic system tended towards a singular matrix.

It is usually the case that there are more measurement lines than pixels, in which case m > n and the system is redundant. However, this redundancy opens the (highly likely) possibility for inconsistent data. If rank(S) = n, then a unique least squares solution can be found by using the Moore Penrose generalized inverse:

x † = S † y (6)
with covariance still given by σ 2 (S T S) -1 . (Recall that S T S is invertible as it is square and has the same rank as S.) Similar reconstructed variance behavior can be observed with overdetermined (redundant) systems as with the simple 9 × 9 matrix example we showed here. (A 12 × 9 example can be formed by adding these three rows to the matrix S: (0, 0, A, 0, 0, A, 0, A, 0); (0, 0, 2D, 0, 0, 2D, 0, D, D); (0, 0, E, 0, 0, 3E/5, 0, 0, 0) with D = √ 29/10 and E = √ 89/8. These 3 rows are each linear combinations of the first 8 rows of S, so when α → 0 the rank again drops to eight. All the pixels then become underdetermined, so the variance behavior of the system is very similar to that of Figure 2).

Region-of-interest reconstruction

If rank(S) < n, the system in equation ( 2) does not have a unique solution, and the Moore-Penrose inverse selects the least squares solution of minimum norm. Region-of-interest reconstruction corresponds to those cases where, while full reconstruction of x is impossible, some components of x can be uniquely determined. In this situation of rank(S) < n, the covariance of the Moore-Penrose solution is given by [START_REF] Lewis | A generalization of the Gauss-Markov theorem[END_REF] cov(x † ) = σ 2 (S T S) † .

(7)

We are interested in the reconstructed variance of those pixels. We will illustrate some surprising behavior below, that can occur when the system matrix is changed by one element (or by one measurement line).

Small-scale examples

We consider the following 7 × 9 ROI system

S =          
1 1 1 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 1 0 0 1 0 0 1 0

A 0 A 0 0 A 0 0 A 0 0 0 A 0 0 A 0 0 0 A 0 A 0 0 A 0 0 0 0 ᾱA 0 0 A 0 0           (8)
where ᾱ = 1 -α, and as before, A = √ 5/2, and α ∈ [0, 1]. The seven corresponding projection lines are shown in Figure 3. This underdetermined system is untypical because pixels x 1 to x 3 can be solved for any value of α (for example: 3x 1 = y 1 + y 2 + (-y 4 + y 5 -y 6 )/A, 3x 2 = y 1 -2y 2 + (2y 4 + y 5 -y 6 )/A, 3x 3 = y 1 + y 2 + (-y 4 -2y 5 + y 6 )/A). Our ROI is these three pixels. (It turns out that for certain values of α, the pixels x 4 and x 7 can be reconstructed also.) It can easily be verified that the rank of S is equal to 7 for all α except α = 0 when the rank is 6. (Two independent vectors in the nullspace of S are [0, 0, 0, 0, 0, 1, 0, 0, -1] T and [0, 0, 0, 0, 1, 0, 0, -1, 0] T , and additionally when α = 0, [0, 0, 0, 1, 0, 0, -1, 0, 0] T .)

Using equation ( 7), we computed the variance analytically for the 3 ROI pixels, and we have plotted the variance (solid lines in Figure 4). Note that for both the full system matrices of equation ( 5) and the ROI matrices of equation ( 8), the rank of the matrix decreases by one when α reaches zero. However, as α → 0, the variance increases dramatically in the full system case, whereas the variance remains constant for the ROI case. Even more surprising, the variance jumps to a lower value in the ROI case when α reaches zero. So the ROI can be reconstructed better when the rank of the system matrix drops, and in this case, with just a tiny change in one entry of the matrix. We verified the variance numerically in the ROI pixels, using the same test image. Figure 4 shows the agreement of the numerical simulations with the theoretically calculated variance of the three ROI pixels. We examined a second similar small scale example, using a matrix identical to that of equation (8) except for the seventh (the last) row, which is parameterized by β ∈ [-1, 1] :

S =           1 1 1 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 1 0 0 1 0 0 1 0 A 0 A 0 0 A 0 0 A 0 0 0 A 0 0 A 0 0 0 A 0 A 0 0 A 0 0 0 0 F 4 0 0 F 7 F 8 0           (9)
where, as before, A = √ 5/2, and where

F 4 = β 2 + 2β + 5 2 (10) F 7 =    1 - 2β 1 + β F 4 β ∈ (0, 1] F 4 β ∈ [-1, 0] (11)
F 8 =    2β 1 + β F 4 β ∈ (0, 1] 0 β ∈ [-1, 0] (12)
Here we have the same ROI of {x 1 , x 2 , x 3 }, and it turns out that no other pixels can be resolved for this system, irrespective of the value of β. (A direct calculation reveals that only the x 1 , x 2 , and x 3 axes, and no others, are perpendicular to the nullspace of S.) Note that when α = 0 and β = 0, the two matrices (of equations ( 8) and ( 9)) are identical. The measurement line corresponding to the last row is described as follows. The parameter β lies on the bottom horizontal line of the image pixels, with β = 0 corresponding to the line separating (unit) pixels x 7 and x 8 as shown in Figure 3. The measurement line corresponding to a value of β ∈ [-1, 1] starts in the upper left corner of pixel x 4 and meets the lower boundary of the image at the value β. At most three matrix entries are involved for any value of β: entries S 7,4 , S 7,7 , S 7,8 . The three corresponding line-lengths for these pixels are given by equations ( 10), ( 11) and ( 12).

The analytically computed and numerically estimated variances in the ROI pixels are shown in Figure 4 for β ∈ [-1/2, 1]. The analytic calculations showed that for β > 0, the reconstructed variance in the 3 ROI pixels remained the same, even though three of the matrix values were approaching the critical β = 0 where the rank of the matrix changed. At that point the variances abruptly decreased, and remained fairly constant (with a slight upward slope in Figure 4) as β became more negative. Two new features in this example are (i) there is no change in the number of pixels that can be reconstructed when the matrix passes the critical point (β = 0), and (ii) the change in variance is sustained, not an apparently instantaneous "blip".

To verify that there is a discontinuity in the reconstructed variances even for extremely small changes in the last measurement line, the variance modifications corresponding to changes in α of 0.001 (for the first ROI matrix) and changes in β of 0.0005 (for the second ROI matrix) are given in Table 1. For pixel 2, for example, the variance jumped by more than 20% following a tiny change in the measurement line, confirming that the variance in the ROI is not continuous as a function of the system matrix elements.

Table 1: Theoretical variances in the ROI pixels when a small variation is applied to one measurement line.

α = β = 0 α = 0.001 β = 0.

A medium-scale example

The small-scale problems present simple and easily reproducible examples. We also considered a larger 32 × 32 pixel image of the Shepp-Logan phantom [START_REF] Shepp | The Fourier reconstruction of a head section[END_REF] with a 16 × 16 square ROI in the upper left corner. A 392 × 1024 system matrix was built to reconstruct the ROI. Amongst those 392 measurements, 288 were necessary to uniquely reconstruct the ROI pixels (we built an intermediate full rank system comprised of the 16 × 16 ROI pixels plus 2×16 grouped pixels corresponding to the horizontal pixel lines of the upper right corner and the vertical pixel lines of the lower left corner). Other measurements were added to include all 32×32 pixels in the system, "break" grouped pixels (identical columns in the system matrix) and add some redundancy. The reference Shepp-Logan phantom and the corresponding pseudo-inverse reconstruction from the 392 noise-free measurements are shown in Figure 5. Similar to the first ROI example of equation ( 8), we allowed one element of the system matrix to vary: S 392,513 = ( √ 257/16)ᾱ. Note from figure 5 that pixel 513 and measurement line 392 are peripheral and do not involve the 16 × 16 ROI region. Geometrically, increasing α from 0 corresponds to translating the last measurement line downwards, just as in the first ROI example. For non-zero α, the matrix had rank 368, and for α = 0, the rank was 367. We examined the reconstructed variance for the two system matrices corresponding to α = 0 and α = 0.001. These two matrices differed in only one of the 401,408 matrix elements, and by a tiny amount: 1.00195 became 1.00095. We evaluated the impact of this small variation on the ROI analytically and numerically. The variance of the noise was again set to σ 2 = 1.

The results are shown in Figure 6. The variance of the pixels reconstructed using α = 0 was smaller than when using α = 0.001. With α = 0, the reconstructed variance was reduced in nearly every ROI pixel, about 10% on average and up to 89%. The numerical variances are shown, but the theoretical and experimental calculations agreed well (the largest differences were less than 2%).

Discussion

Our observations have shown that, in the case of ROI reconstruction, the variance of reconstructible pixels does not always depend continuously on the matrix entries. We saw that a minuscule change of just one matrix entry can result in jumps in the reconstructed variances of all the ROI pixels. Close examination of the small-scale examples can provide some explanations for this remarkable behavior. Both matrices (equations ( 8) and ( 9)) had the same first six rows, and the pixel values x 1 , x 2 , x 3 could be found from the first six measurement values y 1 , y 2 , ...y 6 . Now, for α and β both nonzero, the two matrices are both full rank (with rank 7) so in both cases the last row is not a combination of the other six rows. On the other hand, when α = β = 0 the two matrices are the same, and the rank drops to 6, because the seventh row is now a linear combination of the first six rows. This new redundancy provided by the seventh row allows y 7 to contribute to the solution for x 1 , x 2 , x 3 , thereby improving the accuracy of the solution in the presence of noise.

The ROI pixels x 1 , x 2 , x 3 can be solved using the first six rows of the matrix, and for the full rank situation of α and β both non-zero, the seventh row is not helpful. This statement is a consequence of a theorem published recently [START_REF] Defrise | Do additional data improve region-of-interest reconstruction?[END_REF]. The essence of that theorem is that if the new row (or rows) of the matrix are independent of the existing rows, then there will be no improvement in the reconstructed variance of the ROI pixels.

Conclusion

We have shown a remarkable and seemingly counterintuitive phenomenon: a tiny alteration of the system matrix of a ROI problem can drastically change the variance of the reconstructed pixel values in the ROI. We note that this effect occurs when the system matrix changes rank, and, surprisingly, improved reconstruction occurred with lower rank of the matrix. We linked this result to a recently published theorem on ROI reconstruction.

Figure 1 :

 1 Figure 1: Projection lines used in the full reconstruction system. The projection line corresponding to the last row of S is shown in red.

Figure 2 :

 2 Figure2: Variance of the 9 pixels, computed analytically (solid line) and numerically (circles), as one matrix entry, S9,1, changes linearly from C to zero (i.e. one projection line is translated). Fifty α values between 0 and 1 were used for the numerical tests. Note the logarithmic vertical scale.

Figure 3 :

 3 Figure 3: Left: Projection lines used in the ROI reconstruction system of equation (8). Right: The last row of matrix S of equation (9) corresponds to the measurement line that crosses pixel x4 and exits the bottom of pixel x7 (β < 0) or pixel x8 (β ≥ 0) at a location determined by β.

Figure 4 :

 4 Figure 4: Variance of the 3 ROI pixels, computed analytically (solid line) and numerically (circles), as the projection line is translated (left) or rotated (right).

Figure 5 :

 5 Figure 5: Reference Shepp-Logan phantom (left) and reconstruction from noise-free data (right), with the ROI pixels indicated by the red square. The 392 nd measurement line is shown in red and pixel 513 in yellow.

Figure 6 :

 6 Figure 6: Top left: variances of the 16 × 16 pixels in the ROI reconstructed with α = 0.001. Top right, bottom left and right: horizontal profiles through the variance images reconstructed using either α = 0 or α = 0.001. The shown variances were computed numerically. Note the logarithmic scale.

Acknowledgements

This work was partially supported by grant ANR-17-CE19-0006 (ROIdoré project) from the Agence Nationale de la Recherche.