Mohamed Anwar Abouabdallah
email: maabouabdallah@gmail.com

Olivier Coulaud
email: olivier.coulaud@inria.fr

Nathalie Peyrard
email: nathalie.peyrard@inrae.fr

Alain Franc
email: alain.franc@inrae.fr

Computing WSBM marginals with Tensor-Train decomposition

Keywords: Binary marginals, Weighted Stochastic Block Model, Variables separation, Tensor-Train format, low rank approximation, TT matrices. approximation

The Weighted Stochastic Block Model (WSBM) is a statistical model for unsupervised clustering of individuals based on a pairwise distance matrix. The probabilities of group membership are computed as unary marginals of the joint conditional distribution of the WSBM, whose exact evaluation with brute force is out of reach beyond a few individuals. We propose to build an exact Tensor-Train (TT) decomposition of the multivariate joint distribution, from the SVD of each binary factor of a WSBM, which leads to variables separation. We present how to exploit this decomposition to compute unary and binary marginals. They are expressed without approximation as products of matrices involved in the TT decomposition. However, the implementation of the procedure faces several numerical challenges. First, the dimensions of the matrices involved grow faster than exponentially with the number of variables. We bypass this difficulty by using the format of TT-matrices. Second, the TT-rank of the products grows exponentially. Then, we use a numerical approximation of matrices product that guarantees a low TT-rank, the rounding. We compare the TT approach with two classical inference methods, the Mean-Field approximation and the Gibbs Sampler, on the problem of binary marginal inference for WSBM with 1 various distances structures and up to fifty variables. The results lead to recommend the TT approach for its accuracy and reasonable computing time. Further researches should be devoted to the numerical difficulties for controlling the rank in rounding, to be able to deal with larger problems.

Introduction

Pattern discovery in data sets has become a challenging goal in a large diversity of application domains. It is not so easy to define what an interesting pattern is, but among the quasi infinite number of possibilities, manifold learning for continuous patterns and clustering for discrete patterns are two mainstays for building a compact representation of a data set. Here, we are interested in cluster discovery formalised by a statistical model because it quantifies the uncertainties as limits to our knowledge through probabilities for group membership instead of yes/no answer. In many situations, one knows pairwise distances or dissimilarities between a set of items (see [START_REF] Deza | Encyclopedia of Distances[END_REF]). An example is distances between DNA sequences in an environmental sample, which are distances on strings computed from local alignment score [START_REF] Gusfield | Algorithms on strings, trees, and sequences[END_REF][START_REF] Smith | Identification of common molecular subsequences[END_REF]. The question of clustering items based on their pairwise distances is also classical in time series clustering [START_REF] Ansari | Spatiotemporal clustering: a review[END_REF]. For distance-based clustering, Weighted Stochastic Block Models (WSBM) are more and more used with a large scope of application domains like social sciences [START_REF] Barbillon | Stochastic block models for multiplex networks : An application to a multilevel network of researchers[END_REF], analysis of ecological interaction networks [START_REF] Miele | Revealing the hidden structure of dynamic ecological networks[END_REF], neurology [START_REF] Faskowitz | Weighted stochastic block models of the human connectome across the life span[END_REF], or numerical taxonomy [START_REF] Abouabdallah | Does clustering of DNA barcodes agree with botanical classification directly at high taxonomic levels? Trees in french guiana as a case study[END_REF]. This statistical model formalises the groups the items belong to as latent variables. One of their main advantages is the possibility to find out groups which are not communities (the items within a groupare not necessarily close to each others), like dissortative structures which exist in bipartite graphs. This property gives them great flexibility to adapt to a variety of structures and widen the scope of patterns which can be discovered beyond community detection.

Tools exist to estimate the parameters of a WSBM, as a prerequisite before inferring the group memberships. The classical approach for parameters estimation in models with latent variables is the EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] which, at step E, requires the computation of unary and binary marginals of the conditional distribution of the group memberships. If there are n items and Q groups, computing a unary WSBM marginal requires Q n-1 sums and it is a known bottleneck. The main approach used to overpass this difficulty is the Variational EM, where an assumption of independence is made and unary marginals are approximated as fixed points of the Mean-Field equations. [START_REF] Daudin | A mixture model for random graph[END_REF][START_REF] Mariadassou | Uncovering latent structure in valued graphs: A variational approach[END_REF]. Beyond the utility of marginalisation for parameters estimation, binary marginals of a WSBM are useful for inferring group memberships as their provide a richer information than the sole unary marginals. Indeed, a group allocation based solely on the unary marginals, i.e. by assigning the individual to the group with higher probability, may be associated with a large uncertainty when the mode is not peaked. In this case, the binary marginal can inform on the source of uncertainty. Indeed, focusing on a particular individual in group g, one can compute the probabilities for it to be in the same group than each other individual of group g. They will reveal if there is equal uncertainty of co-membership or on contrary if uncertainty is due to few specific members.

Therefore, in this work, we consider the problem of computing binary marginals in a WSBM. We will use the fact that a WSBM is a graphical model. A graphical model is a stochastic model in which the joint distribution of the variables can be expressed as a product of factors, each factor involving only a subset of the variables. In a WSBM, the joint distribution of interest is the conditional distribution of the group memberships, and it is a pairwise graphical model (factors are binary at most). For some graphical models, the factorisation property makes marginalisation simple. It depends on the structure of the graph associated to the graphical model, a topic which has been thoroughly studied (see e.g. [START_REF] Wainwright | Graphical models, exponential families, and variational inference[END_REF][START_REF] Koller | Probabilistic Graphical Models -Principles and Techniques[END_REF][START_REF] Murphy | Machine Learning: A Probabilistic Perspective[END_REF][START_REF] Peyrard | Exact or approximate inference in graphical models: why the choice is dictated by the treewidth, and how variable elimination can be exploited[END_REF]. Very often, one needs not only to compute one marginal for one variable, but all marginals for all variables, or subsets of variables. This global calculation is performed exactly when the graph of the model is a tree, with message passing or belief propagation. On graphs with loops it can sometimes be performed exactly with message passing using a Junction Tree built from the graph, and otherwise the marginals are approximated as fixed points of loopy belief propagation. The tree width, a feature of the graph associated to the graphical model, characterises the complexity of the calculation of the marginalisation. If the tree-width is w, this complexity is in O(kQ w+1 ) if there are Q possible values for a discrete random variable. In case of large tree-width, the exact calculation by message passing in the Junction-Tree becomes impossible.

The graph associated to the conditional distribution of the group memberships in a WSBM is a clique. The Junction Tree has one node (the clique), and the tree-width is n -1. It is the most complex situation with n nodes. Therefore, the issue of marginalisation in WSBM should be addressed with other approaches. The abovementioned Mean-Field approach is a classical option, and MCMC techniques [START_REF] Robert | Monte Carlo statistical methods[END_REF], like the Gibbs sampler is another. The first one is quick but is limited to unary marginals, and the second allows calculation of binary marginals but to the price of a higher computational cost.

In this work we explore an alternative to the Mean-Field and the MCMC approaches for approximate marginals inference in a WSBM. We use the fact that the joint distribution of a graphical model with n variables and Q possible values for each random variable is a tensor of order n and dimension Q per mode (see section 3 for the definition of these terms). If the joint distribution of a set of random variables can be expressed as a tensor of low rank, then marginalisation becomes fairly easy. One possible strategy for marginalisation of the joint distribution is to find the best The Tucker decomposition is exploited in [START_REF] Lyu | Latent space model for higher-order networks and generalized tensor decomposition[END_REF]. The authors define a new family of graphical models, the latent space models, to study the organisation of individuals linked by hyperedges. To estimate the model parameters and restore the latent positions of the individuals, they propose to maximise the probability of the observed adjacency tensor knowing the latent variables, under a constraint on the Tucker ranks of the tensor of interest. The approach does not requires marginal inference, and the authors do not consider this task. Moreover marginalisation with the Tucker decomposition remains computationally complex.

On the contrary, finding the best TT rank r approximation of a given tensor is a well-posed problem, numerical algorithms for computing it are robust [START_REF] Oseledets | A new tensor decomposition[END_REF][START_REF] Oseledets | Tensor-Train decomposition[END_REF], and as we will see the TT format is well adapted for marginalisation because it leads to variables separation. [START_REF] Ducamp | An efficient low-rank tensors representation for algorithms in complex probabilistic graphical models[END_REF] use the TT format of the factors of a BN to write and implement the equations of the classical message passing algorithm for probabilistic inference (i.e. to compute some marginals of the BN). They show on several BN examples how this approach can lead to a reduction in the storage and a faster inference time compared to the exact message passing algorithm, with a reasonable error on the marginals values. This property has been exploited as well in [START_REF] Novikov | Putting MRFs on a Tensor Train[END_REF] to compute the normalising constant and the unary marginals of any graphical model by first finding a best low rank TTapproximation of each factor and, second, writing the product of the approximations of the factors as a tensor in TT-format by a clever utilisation of the mixed-product property of the Kronecker product.

Our work aims at proposing a way to compute the binary marginals of a WSBM, and more generally of a graphical model with binary factors, by using the TT format of the tensor associated to the joint distribution of the model. The article is organised as follows. We first recall some definitions and notations that we will use throughout this paper, on WSBM in section 2, and on tensors in section 3. The material in those sections is standard. In section 4, we present in detail the approach of [START_REF] Novikov | Putting MRFs on a Tensor Train[END_REF]. Our contribution is based on an extension of their work, and is developed in the next 3 sections. In section 5, we observe that the Singular Value Decomposition (SVD) of the matrices which are tables of binary factors is a variable separation between indices of rows and of columns. When plugged in step 1 in Novikov procedure, this leads to writing the joint distribution of a pairwise graphical model exactly as a tensor in TT format, without approximation. Distributivity of multiplication over addition in the ring of matrices leads to an algebraically simple calculation of the maginals. We develop a approach by recursion for the simultaneous calculation of all binary marginals, over all variables and all possible values the variables, which extends Novikov's procedure for computing all unary marginals. If such an approach is algebraically straightforward, it leads to sums and products between matrices with huge dimensions. In section 6, elaborating on the TT-matrix approach used in [START_REF] Novikov | Putting MRFs on a Tensor Train[END_REF], we show that it leads to numerical locks for which some solutions are proposed, essentially based on the rounding [START_REF] Oseledets | Tensor-Train decomposition[END_REF][START_REF] Al Daas | Parallel algorithms for tensor-train arithmetic[END_REF]. Finally, in section 7, we compare both in time and accuracy the computation of binary marginals by our TT approach to three classical methods: brute force enumeration, simulation with a Gibbs sampler, and as fixed point of a Mean-Field approximation.

Weighted Stochastic Block Model

The Stochastic Block Model (SBM, [START_REF] Holland | Stochastic blockmodels: First steps[END_REF][START_REF] Daudin | A mixture model for random graph[END_REF] was originally defined for identifying groups in a set of individuals connected by binary relationships (two individuals are either connected or not). In this work, we consider the extension to dissimilarity matrices, where to each pair of individuals is attached a weight which is a measure of dissimilarity between the two individuals.

The Weighted SBM (WSBM) is a model from the family of statistical models with latent variables. Let us consider n individuals. The observed variable is the dissimilarity matrix D of size n × n, with element D[i, j]. The latent (unobserved) variables are the group memberships of each individual: Z i ∈ {1, . . . , Q} is the group of individual i. The model relies on two assumptions. First, the Z i 's are independent and their distribution is described by the vector of probabilities α = (α 1 , . . . , α Q ), such that P (Z i = q) = α q . Second, the distribution of the dissimilarity between i and j depends only on the groups of i and j (an not on the individuals i and j). We assume that the distribution of D[i, j] conditionally to (Z i , Z j ) depends on a parameter Λ and we denote θ = (α, Λ) the model's parameters. For instance, P(D[i, j] | Z i = q, Z j = q ′ ) can be the Poisson distribution with parameter λ qq ′ . The connectivity matrix of the model is the Q by Q matrix Λ such that Λ[q, q ′ ] = λ qq ′ . Let Z = (Z 1 , . . . , Z n ). Then, the joint probability distribution of Z and D is:

P θ (Z, D) = P Λ (D | Z)P α (Z) = n i=1 i>j P Λ (D[i, j] | Z i , Z j ) n i=1 P α (Z i ) (1)
In practice, when modelling a dissimilarity matrix with a WSBM the main goal is to recover each individual's memberships, i.e. to restore the Z i s. This is performed by maximising the conditional distribution P θ (Z | D). Beyond this 'definitive' affectation of an individual to a group, one may be interested by the vector of probabilities that an individual belongs to each group, namely {P θ (Z i = q | D)} 1≤q≤Q , to quantify the uncertainty on the membership. Another probability of interest is the joint conditional probability that two individuals i and j belong to the same group: P θ (Z i = q, Z j = q | D). Therefore, the distributions of interest are not only the conditional distribution P θ (Z | D) but also its unary and binary marginals. The probability

P θ (Z | D) is proportional to P θ (Z, D), P θ (Z | D) = P θ (Z, D)/W (2)
and the normalising constant W is equal to the likelihood of the observed dissimilarity matrix D:

W = z∈{1,...,Q} n P θ (Z = z, D) = P θ (D)
The binary marginal for two individuals, for instance 1 and 2, is defined by

P θ (Z 1 = q, Z 2 = q ′ | D) = Q z3=1 . . . Q zn=1 P θ (Z 1 = q, Z 2 = q ′ , Z 3 = z 3 , . . . , Z n = z n | D) = 1 W Q z3=1 . . . Q zn=1 P θ (Z 1 = q, Z 2 = q ′ , Z 3 = z 3 , . . . , Z n = z n , D) (3)
The unary marginal P θ (Z i = q | D) is obtained by marginalising any binary marginal involving Z i :

P θ (Z i = q | D) = Q q ′ =1 P θ (Z i = q, Z j = q ′ | D)
In practice, un-normalised binary marginals, Pθ (Z i = q, Z j = q ′ | D) = P θ (Z i = q, Z j = z j , D) are computed, by marginalising P θ (Z, D) instead of P θ (Z | D). Then, un-normalised unary marginals Pθ (Z i = q | D) = P θ (Z i = q, D) are derived by marginalisation of the Pθ (Z i = q, Z j = q ′ | D). Finally W is obtained from any un-normalised unary marginal as

W = Q q=1 Pθ (Z i = q | D).
For the following, it will be convenient to cast the WSBM in the family of graphical models [START_REF] Koller | Probabilistic Graphical Models -Principles and Techniques[END_REF]. A vector of random variables Z = {Z 1 , . . . , Z n } is a graphical model if the joint distribution can be expressed (up to the normalising constant) as a product of functions, called factors, involving only subsets of the variables. From expression (1), we can see that in a WSBM, the distribution P θ (Z | D) is that of a graphical model with only binary and unary factors. There are n(n-1) 2 binary factors, equal to P Λ (D[i, j] | Z i , Z j ) and n unary factors equal to P α (Z i ). If computed by applying naïvely expression (3), the complexity of evaluation of one binary marginal is in O(Q n-2 ).

The Mean-Field approximation of a WSBM is its best approximation, in the sense of the Kullback-Leibler divergence, by a product of independent unary factors ϕ.

P θ (Z 1 = z 1 , . . . , Z n = z n , D) ≈ n i=1 ϕ i (z i )
In such a case, we have separation of variables, and marginalisation is simple due to the distributivity of multiplication over addition:

Q z1=1 . . . Q zn=1 n i=1 ϕ i (z i ) = n i=1 Q zi=1 ϕ i (z i ) (4)
which requires nQ sums and n -1 products. Such a distributivity exists as well in the ring of matrices. Therefore, marginalisation with separation of variables can be done as well in the case where the joint distribution of the n variables can be approximated by a product of unary factors G, each being a matrix, like, for n = 3

P θ (Z 1 = z 1 , . . . , Z n = z n , D) ≈ G 1 (z 1 ) G 2 (z 2 ) G 3 (z 3 ) with      G 1 (z 1 ) ∈ R 1×r1 G 2 (z 2 ) ∈ R r1×r2 G 3 (z 3 ) ∈ R r2×1
(5) Such a decomposition is known as Tensor-Train decomposition (TT), and has been thoroughly studied in the literature (see next section). We will see in Section 5 that the joint distribution of a WSBM (or any pairwise graphical model) can be exactly decomposed in a TT format, which paves the way for its marginalisation. Indeed, in the simple case of the probability written in TT-format in equation ( 5), we have

W = Q z1=1 Q z2=1 Q z3=1 G 1 (z 1 ) G 2 (z 2 ) G 3 (z 3 ) = Q z1=1 G 1 (z 1 ) Q z2=1 G 2 (z 2 ) Q z3=1 G 3 (z 3 ) (6)
3 Tensors and Tensor Trains

The main idea which has driven this study is to consider that the joint distribution of a WSBM defines a tensor, i.e. a multi-way array. Then, following [START_REF] Novikov | Putting MRFs on a Tensor Train[END_REF], we use possibilities given by computing in TT-format (mainly due to [START_REF] Oseledets | A new tensor decomposition[END_REF][START_REF] Oseledets | Tensor-Train decomposition[END_REF] to separate the variables for each variable z i of the WSBM in the expression of P θ (Z, D). This separation enables to compute the marginals of the joint distribution. Before positioning our method, we recall here the main definitions and results on tensors we will use in the rest of the article.

Tensors:

Let E = (E 1 , . . . , E i , . . . , E n ) be a family of n finite dimensional real vector spaces. A tensor T on E is a multilinear form on E 1 × . . . × E n . Then, n is the order of the tensor, E i for 1 ≤ i ≤ n is a mode of T (often simplified as i) and the dimension of

E i , denoted Q i , is its dimension for mode E i . The multilinear dimension of the tensor is denoted Q 1 × . . . × Q n , with T ∈ R Q1×...×Qn .
If a basis has been selected on each space E i , T can be represented in those basis as a n-dimensional array, the elements of which are denoted

T[z 1 , . . . , z i , . . . , z n ] with 1 ≤ i ≤ n, 1 ≤ z i ≤ Q i ,
or T z1...zn .

Slice of a tensor:

Let us explain what a slice of a tensor is on the example of a 3-modes tensors on

E 1 × E 2 × E 3 with indices z 1 , z 2 , z 3 . A general term of T is T[z 1 , z 2 , z 3 ] and T has dimension Q 1 × Q 2 × Q 3 with 1 ≤ z 1 ≤ Q 1 , 1 ≤ z 2 ≤ Q 2 and 1 ≤ z 3 ≤ Q 3 . The slice T E2 (z 2 )
on mode E 2 with index z 2 or, more simply, T 2 (z 2 ), is the tensor of dimension Q 1 × Q 3 obtained by letting the indices in modes E 1 and E 3 (̸ = E 2 ) run over all dimensions of respectively E 1 and E 3 and fixing the index for mode E 2 at z 2 :

T 2 (z 2 ) = (T[z 1 , z 2 , z 3 ]) z1,z3 with 1 ≤ z 1 ≤ Q 1 , 1 ≤ z 3 ≤ Q 3 .
or, more simply:

T 2 (z 2 ) = T[:, z 2 , :]
This can be extended with some technicalities in notations to any tensor, by selecting a mode E i (referred to as T i ) and an index z i (1 ≤ z i ≤ Q i ):

T i (z i ) = T[:, . . . , :, z i , :, . . . , :]

Frobenius norm of a tensor:

The set of tensors of a gi ven multilinear dimension is a vector space, which can inherit the Euclidean structure of all modes. The Frobenius norm ∥T∥ of T is defined by

∥T∥ 2 = Q1 i1=1 . . . Qn in=1 T 2 [i 1 . . . i n ].
This will be useful to define a distance between tensors, as d(T, T ′ ) = ∥T -T ′ ∥.

Tensor-Train format:

The Tensor-Train format, denoted TT-format in the rest of this article, is a format proposed in [START_REF] Oseledets | A new tensor decomposition[END_REF][START_REF] Oseledets | Tensor-Train decomposition[END_REF], well adapted to separation of variables (it can be read as a local, or componentwise, separation of variables). A tensor T = (T[z 1 , . . . , z n ]) z1,...,zn is in TT-format if there exists matrices (Oseledets, 2011, equation (1.2)). The matrices G i (z i ) are called the cores of the TTdecomposition of T. This is explicitly variable separation for the modes of the tensor. The dimensions of matrix G i (z i ) are r i-1 × r i , with r 0 = r n = 1. These are called the ranks of the TT-format. The TT-rank of T is r = (r 0 , . . . , r n ).

G i (z i ) for 1 ≤ i ≤ n such that ∀ (z 1 , . . . , z n ), with 1 ≤ z i ≤ Q i , T[z 1 , . . . , z n ] = G 1 (z 1 ) . . . G n (z n ),
The TT-decomposition of T can be developed by defining componentwise for each mode i

G i (z i ) = (G i (z i )[α i-1 , α i ]) αi-1,αi ∈ R ri-1×ri with 1 ≤ α i-1 ≤ r i-1 1 ≤ α i ≤ r i and the tensor G i of order 3 of dimension r i-1 ×Q i ×r i with coefficients G i [α i-1 , z i , α i ].
The matrix G i (z i ) is the slice z i of second mode of tensor G i . This leads to

∀(z 1 , . . . , z n ), T[z 1 , . . . , z n ] = r1 α1=1 r2 α2=1 . . . rn-1 αn-1=1 G 1 [z 1 , α 1 ]G 2 [α 1 , z 2 , α 2 ] . . . G n [α n-1 , z n ]
Vidal (2003) has shown in the framework of the MPS that such a decomposition always exists, and is not unique. However, a simple dimension analysis shows that the rank may be huge (r ∼ Q n-1 /n) in general). A classical problem is, a tensor T and a rank r being given, to find a tensor in TT format of rank r which is closest to T with Frobenius norm. [START_REF] Oseledets | Tensor-Train decomposition[END_REF] has proposed a robust algorithm, called TT-SVD, to solve this optimisation problem (see Oseledets, 2011, Algorithm 1). He has shown that many basic operations in tensor calculus can be done in the framework of the TT-format, i.e. knowing the cores only, without using explicitly the coefficients of the tensors involved: addition, contraction, Hadamard product, inner product, matrix-byvector product, and the best approximation at a given TT-rank called the rounding. Bypassing the explicitation of all coefficients leads to powerful tools for many modes tensors which cannot be stored in memory.

Suitability for marginalisation:

Marginalisation in tensors is an operation associated with slicing. It consists in (i) selecting a slice, which can be empty and (ii) summing up all terms which are not in the slice. Here, we show why TT format is suitable for marginalisation. Let us first develop it for a tensor of order 3 with equal dimension Q on all modes, as

∀1 ≤ z 1 , z 2 , z 3 ≤ Q, T[z 1 , z 2 , z 3 ] = G 1 (z 1 ) G 2 (z 2 ) G 3 (z 3 ) with      G 1 (z 1 ) ∈ R 1×r1 G 2 (z 2 ) ∈ R r1×r2 G 3 (z 3 ) ∈ R r2×1
Let us consider the slice defined by mode E 2 and index z 2 to define the marginal m 2 (z 2 ). Then, because of distributivity of multiplication on addition in the ring of matrices

m 2 (z 2 ) = Q z1,z3=1 T[z 1 , z 2 , z 3 ] = Q z1,z3=1 G 1 (z 1 )G 2 (z 2 )G 3 (z 3 ) = Q z1=1 G 1 (z 1 ) G 2 (z 2 ) Q z3=1 G 3 (z 3 ) = B 1 G 2 (z 2 ) B 3 (7)
with

B 1 = Q z1=1 G 1 (z 1 ) ∈ R 1×r1 , B 3 = Q z3=1 G 3 (z 3 ) ∈ R r2×1 .
This can be extended in a straightforward way to tensors of any order n with unequal dimensions.

TT to compute the normalising constant of a graphical model

In this section we present in detail the approach developed in [START_REF] Novikov | Putting MRFs on a Tensor Train[END_REF]. They use the TT-approximation of the factors of a graphical model to compute the normalising constant. Our work on inference for WSBM is elaborated on their results. We present in Table 1 a few notations used in the rest of the article that allow a transition between the domains of statistical modelling and tensor algebra.

symbol WSBM tensor n # of individuals # of modes i an individual a mode Q number of groups dimension of a mode z i latent group membership index in mode i for individual i z i ∈ {1, . . . , Q} ψ ij (z i , z j ) a factor in WSBM Ψ ij [z i , z j ] matrix of the table associated to ψ ij ψ(z 1 , . . . , zn) joint distribution Ψ[z 1 , . . . , zn] associated tensor m # of factors n(n -1)/2 - ℓ
index of a factor -Table 1 Notations for transition between WSBM and tensor calculus.

When a is a vector or M a matrix, we denote by a[i] the component i of a and by M [i, j] the component at row i and column j of M . If, for example, ψ ij is a binary factor of a graphical model on edge (i, j), represented by a matrix Ψ ij , if the group membership of individual i (resp. j) is z i (resp. z j ), both in state space S, we denote by ψ ij (z i , z j ) the factor as a function from S × S on R and by Ψ ij [z i , z j ] the coefficients of the matrix associated to it.

Let us have a graphical model with n variables, each representing an individual, with z i being the group of individual i. The graphical model is defined on a hypergraph G = (I, E) where I = {1, . . . , n} and E is a set of hyperedges A ℓ ⊂ I for ℓ ∈ {1, . . . , m}. Let us denote n ℓ = |A ℓ |, and by z ℓ the groups of the individuals in A ℓ i.e. z ℓ = (z i ) i∈A ℓ . Then, the un-normalised distribution of the graphical model can be written as a product of m factors indexed by ℓ:

ψ(z 1 , . . . , z n ) = m ℓ=1 ψ ℓ (z ℓ )
The objective of the work in [START_REF] Novikov | Putting MRFs on a Tensor Train[END_REF] is to compute the normalising constant

W = Q z1=1 . . . Q zn=1 m ℓ=1 ψ ℓ (z ℓ )
Their method is composed of four steps presented here as we will elaborate on them to develop our method for computing marginals in a WSBM. Let us denote by Ψ := Ψ[z 1 , . . . , z n ] = ψ(z 1 , . . . , z n ) the tensor which is the distribution table associated to the un-normalised joint distribution ψ.

Step 1: approximating each factor for a given TT-format A tensor Ψ ℓ can be built from each factor ψ ℓ . Its order is n ℓ , the number of indices in A ℓ , and its dimensions are Q for each mode. We define

Ψ ℓ [z ℓ ] = ψ ℓ (z ℓ ), with Ψ ℓ ∈ R Q×...×Q
The first step consists in selecting a rank r and computing a best TT-approximation Ψ ℓ of Ψ ℓ at rank r (more rigorously, at rank r = (1, r, . . . , r, 1)). In the TT-format, each component

Ψ ℓ [z i ℓ 1 , . . . , z i ℓ n ℓ
] can be written as a product of matrices with separation of variables (here,

z ℓ = (z i ℓ 1 , . . . , z i ℓ n ℓ
) where i ℓ µ is µ-th index in A ℓ ):

Ψ ℓ [z i ℓ 1 , . . . , z i ℓ n ℓ ] = G ℓ i ℓ 1 (z i ℓ 1 ) . . . G ℓ i ℓ n ℓ (z i ℓ n ℓ
), with

       G ℓ i ℓ 1 (z i ℓ 1 ) ∈ R 1×r G ℓ i ℓ µ (z i ℓ µ ) ∈ R r×r i ℓ 2 ≤ i ℓ µ ≤ i ℓ n ℓ -1 G ℓ i ℓ n ℓ (z i ℓ n ℓ ) ∈ R r×1 .
For sake of simplicity, we select a same rank r for each node, but it can be developed with ranks specific to each node. It may happen that, for a given factor and a given rank, the TT-decomposition at rank r is exact. For example, if n ℓ = 3, tensor Ψ ℓ has Q 3 terms, and a TT of order 3 with rank r has Qr(r + 2) terms. A simple dimension analysis shows that if r = Q, an exact TT-decomposition likely exists. In order to encompass in a same notation exact decomposition and best approximation at a given rank for any factor, we will use the same notation Ψ ℓ from now on, and drop the (.).

We then start with

Ψ ℓ [z i ℓ 1 , . . . , z i ℓ n ℓ ] = G ℓ i ℓ 1 (z i ℓ 1 ) . . . G ℓ i ℓ n ℓ (z i ℓ n ℓ
).

Step 2: Adding non essential variable:

This step is simple, but notations can be cumbersome. So, it will be presented first on the small example of [START_REF] Novikov | Putting MRFs on a Tensor Train[END_REF], section 5.1. The idea is to build for each hyperedge A ℓ a new factor ψ ℓ that depends on all variables and not only on z ℓ but whose value does not actually depend on z j for j / ∈ A ℓ :

ψ ℓ : {1, . . . , Q} n ----→ R (z 1 , . . . , z n ) ----→ ψ ℓ (z ℓ )
The tensor Ψ ℓ attached to ψ ℓ has order n. We have

Ψ[z 1 , . . . , z n ] = m ℓ=1 Ψ ℓ [z 1 , . . . , z n ] (8)
The TT-format of tensors Ψ ℓ with non essential variables can be derived from the TT-format of the tensor corresponding to the factor ψ ℓ as shown in the following toy example, with n = 7, A ℓ = {2, 4, 6}. We then have

Ψ ℓ [z 2 , z 4 , z 6 ] = G ℓ 2 (z 2 ) G ℓ 4 (z 4 ) G ℓ 6 (z 6 ) with      G ℓ 2 (z 2 ) ∈ R 1×r G ℓ 4 (z 4 ) ∈ R r×r G ℓ 2 (z 6 ) ∈ R r×1
This can be sketched as

Ψ ℓ [z 2 , z 4 , z 6 ] = G ℓ 2 (z 2 ) G ℓ 4 (z 4 ) G ℓ 6 (z 6 ) 1 × r r × r r × 1
The product of those matrices of dimensions (1, r), (r, r) and (r, 1) is a real. To build the TT-format of Ψ ℓ , we add new cores, one per non essential variable, which are

G ℓ 1 (z 1 ) = 1, G ℓ 3 (z 3 ) = G ℓ 5 (z 5 ) = I r and G ℓ 7 (z 7 ) = 1 as in the following figure: Ψ ℓ [z 1 , z 2 , z 3 , z 4 , z 5 , z 6 , z 7 ] = • • G ℓ 1 (z 1 ) = 1 G ℓ 2 (z 2 ) G ℓ 3 (z 3 ) = I r G ℓ 4 (z 4 ) G ℓ 5 (z 5 ) = I r G ℓ 6 (z 6 )G ℓ 7 (z 7 ) = 1 1 × 1 1 × r r × r r × r r × r r × 1 1 × 1 Finally, Ψ ℓ [z 1 , z 2 , z 3 , z 4 , z 5 , z 6 , z 7 ] = G ℓ 1 (z 1 ) G ℓ 2 (z 2 ) G ℓ 3 (z 3 ) G ℓ 4 (z 4 ) G ℓ 5 (z 5 ) G ℓ 6 (z 6 ) G ℓ 7 (z 7 )
This procedure is applied to each factor of the graphical model. Then, there exists a family of cores

(G ℓ k (z k )) k,ℓ with 1 ≤ ℓ ≤ m and 1 ≤ k ≤ n such that ∀ A ℓ ∈ E, Ψ ℓ [z 1 , . . . , z n ] = n k=1 G ℓ k (z k )
Step 3: Mixed-product property of Kronecker product:

Let us recall that ⊗ denotes the Kronecker product between matrices. Knowing that xy = x ⊗ y if x, y ∈ R, equation ( 8) can be written

Ψ[z 1 , . . . , z n ] = m ℓ=1 Ψ ℓ [z 1 , . . . , z n ] = m ℓ=1 n k=1 G ℓ k (z k )
Let us recall the mixed-product property of Kronecker product (see [START_REF] Horn | Matrix Analysis[END_REF]: if A, B, C, D are matrices with relevant dimensions for the products AB and CD to be possible, then

(AB) ⊗ (CD) = (A ⊗ C) (B ⊗ D).
This leads to

Ψ[z 1 , . . . , z n ] = n k=1 m ℓ=1 G ℓ k [z k ] Let us denote A k (z k ) = m ℓ=1 G ℓ k (z k ) (9)
Then, the un-normalised joint distribution ψ can be expressed as

ψ(z 1 , . . . , z n ) = Ψ[z 1 , . . . , z n ] = n k=1 A k (z k ) (10)
which is separation of variables.

Step 4: computation of the normalising constant and of unary marginals

It is now easy to compute the normalising constant W of ψ. We recall that

W = Q z1=1 . . . Q zn=1 ψ(z 1 , . . . , z n ) Then, W = Q z1=1 . . . Q zn=1 n k=1 A k (z k ) = n k=1 Q z k =1 A k (z k ) = n k=1 B k , with B k = Q z k =1 A k (z k ).
It requires nQ matrix additions and n -1 matrix multiplications only, if all matrices

A k (z k ) have been computed for all 1 ≤ k ≤ n and 1 ≤ z k ≤ Q (hence nQ matrices).
The complexity is no longer exponential with n, but depends on the dimensions of the matrices A k . Expression (10) permits to compute the unary marginals as well, with message passing, as evoked in [START_REF] Novikov | Putting MRFs on a Tensor Train[END_REF]. Let us denote z\z i = (z 1 , . . . , z i-1 , z i+1 , . . . , z n ), i.e. all z j for j ̸ = i. The un-normalised unary marginal for variable i is defined as

m i (z i ) = z\zi∈{1,...,Q} n-1 ψ(z 1 , . . . , z n ) Then m i (z i ) = z\zi n k=1 A k (z k ) =   z1,...zi-1 i-1 k=1 A k (z k )   A i [z i ]   zi+1,...zn n k=i+1 A k (z k )   = i-1 k=1 B k A i [z i ] n k=i+1 B k
For one unary marginal, this requires only n-1 matrix multiplications. If all marginals for all z k have to be computed, the calculations can be organised in order to mutualise them. Computing all marginals m i (z i ) for all i and all z i requires 2(n -2) + 2nQ matrix multiplications: n -2 to compute recursively each i-1 k=1 B k for i = 3 to n, again n -2 to compute recursively each n k=i+1 B k for i = 1 to n -2, and finally 2(n + 1)(Q -1) (because we use q m i (q) = 1) to compute each m i (z i ) for all i and all z i by multiplying three matrices (or only 2 for m 1 and m n ).

This expression does not allow to reveal the real complexity of the calculation, which depends on the dimensions of the matrices. Indeed, if A, B are two p×p matrices, the complexity of the calculation of the product AB is in O(p 3 ). Such dimensions are specific to the graphical model defined by ψ and its factors.

From now on we focus on the case of graphical models with binary factors (as for the WSBM model) and starting from the Novikov approach we establish the following results :

1. We first observe in section 5.1 that if a factor ψ ℓ is binary, the tensor Ψ ℓ attached to it is a matrix, and its TT-SVD is the SVD of this matrix 2. For pairwise graphical models (i.e. all factors are binary, like an Ising model), this leads to a rewriting of the tensor Ψ with separation of variable by equation ( 10) without any approximation (step 1 is exact). This leads to a calculation of all the binary marginals by a recursion approach inspired by Message Passing (section 5.2) 3. We develop this approach on a specific family of pairwise graphical model: the WSBM, for which all binary factors are present (section 5.3). We show that, even though the calculation is algebraically exact, it can be intractable for large n (due to the matrix × matrix product of the A k ), and we propose some low rank approximation with the rounding (section 6).

Separation of variables and computation of binary marginals for pairwise graphical models

In this section, we show that Novikov approach when implemented on a graphical model with binary factors at most leads to rewrite the un-normalised joint distribution in Tensor-Train format, hence with variable separation for computing normalising constant or marginals. This applies to WSBM models. We first develop the approach of section 4 on pairwise graphical models, and focus after on the particular case of WSBM, which are pairwise graphical models.

Separation of variables in the un-normalised joint distribution of a pairwise graphical model

The main advantage of pairwise graphical model is that the TT-format of each binary factor (i.e. separation of variables) is obtained with the SVD of the factor. We implement it as first step in Novikov approach for binary graphical models. Further steps are unchanged, and presented in appendix C.

Let us have a pairwise graphical model with un-normalised distribution ψ defined as follows:

ψ(z 1 . . . , z n ) = (i,j)∈E ψ ij (z i , z j )
where E is the set of binary factors of the graphical model. Factor ψ ij can be given as

a Q × Q matrix Ψ ij with Ψ ij [z i , z j ] = ψ ij (z i , z j )
where z i and z j label the rows and columns of Ψ ij . The SVD of Ψ ij can be written as (see e.g. [START_REF] Horn | Matrix Analysis[END_REF][START_REF] Strang | Linear Algebra and Learning from Data[END_REF])

Ψ ij = U ij Σ ij V t ij , with U ij , Σ ij , V ij ∈ R Q×Q
where Σ ij is a diagonal matrix, and U ij , V ij are columnwise orthonormal. We assume that matrices are full rank, i.e. rank

Ψ ij = Q. Let us denote M ij = U ij Σ ij . Then, Ψ ij = M ij V t ij . So Ψ ij [z i , z j ] = Q q=1 M ij [z i , q] V ij [q, z j ] Let us denote by M ij [a] ∈ R Q the row a of M ij and by V ij [b] ∈ R Q the column b of V ij .
These can be considered as one row and one column matrices respectively

M ij [z i ] ∈ R 1×Q , V ij [z j ] ∈ R Q×1 Then ψ ij (z i , z j ) = Ψ ij [z i , z j ] = M ij [z i ] V ij [z j ] (Matrices product)
This is variable separation for ψ ij (z i , z j ): each matrix depends only on one variable z i or z j .

The subsequent steps (adding non essential variables, computation of the matrices A k (z k ) and B k ) are implemented as in section 4, and are described with a specific implementation for a WSBM in appendix C.

Computing the binary marginals of a pairwise graphical model

Like the unary marginals, the un-normalised binary marginals can be expressed in terms of the B k and the A k (z k ):

m ij (q, q ′ ) = Pθ (Z i = q, Z j = q ′ | D) = B 1 × . . . × B i-1 × A i (q) × B i+1 × . . . × B j-1 × A j (q ′ ) × B j+1 × . . . × B n
Let us define the following quantities:

∀ 1 ≤ i ≤ n, F i,i = B i ∀ 1 ≤ i ≤ n -1, ∀ i + 1 ≤ j ≤ n, F i,j = B i × B i+1 × . . . × B j
Depending on the values of i and j, m ij (q, q ′ ) can be expressed as follows:

i = 1 j = n ⇒ m 1n (q, q ′ ) = A 1 [q] × F 2,n-1 × A n [q ′ ] i = 1 2 < j < n ⇒ m 1j (q, q ′ ) = A 1 [q] × F 2,j-1 × A j [q ′ ] × F j+1,n 1 < i < n -1 j = n ⇒ m in (q, q ′ ) = F 1,i-1 × A i [q] × F i+1,n-1 × A n [q ′ ] 1 < i < n -2 i + 1 < j < n ⇒ m ij (q, q ′ ) = F 1,i-1 × A i [q] × F i+1,j-1 × A j [q ′ ] × F j+1,n i = 1 j = 2 ⇒ m 12 (q, q ′ ) = A 1 [q] × A 2 [q ′ ] × F 3,n i = n -1 j = n ⇒ m in (q, q ′ ) = F 1,n-2 × A n-1 [q] × A n [q ′ ] 1 < i < n -1 i + 1 = j ⇒ m ij (q, q ′ ) = F 1,i-1 × A i [q] × A i+1 [q ′ ] × F i+2,n
All the F i,j can be computed in an efficient way recursively with the following algorithm:

Algorithm 1 Recursive computation of the F ij 1: inputs: B i ∀ 1 ≤ i ≤ n 2: initialisation: ∀ 1 ≤ i ≤ n, F i,i = B i 3: for i = 1 to n -1 do 4:
for j = i + 1 to n do 5:

F i,j = F i,j-1 × B j 6: end for 7: end for 8: return F i,j ∀ 1 ≤ i < j ≤ n
Once the F ij are computed, all binary marginals can be computed with the following algorithm:

Algorithm 2 Computation of all binary marginals of a pairwise graphical model

1: inputs: F i,j ∀ 1 ≤ i < j ≤ n, A k (q) ∀ 1 ≤ k ≤ n, ∀ 1 ≤ q ≤ Q 2: for j = 1 to n -1 do 3: for q ′ = 1 to Q do 4: R j [q ′ ] = A j [q ′ ] × F j+1,n 5:
end for 6: end for 7: for i = 2 to n do 8:

for q = 1 to Q do 9: L i [q] = F 1,i-1 × A i [q] 10:
end for 11: end for 12: for i = 1 to n -1 do 13:

for j = i + 1 to n do 14:

for q = 1 to Q do 15:

for q ′ = 1 to Q do 16: if i = 1 and j = n then 17: m ij (q, q ′ ) = A 1 [q] × F 2,n-1 × A n [q ′ ] 18:
else if i = 1 and 2 < j < n then 19:

m ij (q, q ′ ) = A 1 [z 1 ] × F 2,j-1 × R j [q ′ ] 20: else if 1 < i < n -1 and j = n then 21: m ij (q, q ′ ) = L i [q] × F i+1,n-1 × A n [q ′ ] 22: else if 1 < i < n -2 and i + 1 < j < n then 23: m ij (q, q ′ ) = L i [q] × F i+1,j-1 × R j [q ′ ] 24:
else if i = n -1 and j = n then 25: end for 32: end for 33: return m ij (q, q ′ ) ∀ 1 ≤ i < j ≤ n, ∀ 1 ≤ q, q ′ ≤ Q Computing all the F i,j requires n-1 i=1 (n -i) = n(n -1)/2 matrices products. Then, computing all L i [q] for a given q requires n -1 matrices products, and there are Q possible values for q, hence (n -1)Q matrices products are needed to obtain all matrices L i [q]. There is a similar result for all matrices R j [z j ]. Knowing all F i,j , L i [q] and R j [q ′ ], the calculation of m ij (q, q ′ ) requires 2 matrices products, hence n(n -1)Q 2 products for all binary marginals because there are n(n -1)/2 pairs (i, j). In total, computing all binary marginals requires n(n -1)/2 + 2(n -1)Q + n(n -1)Q 2 ≈ n 2 Q 2 matrix multiplications. Let us recall that the real barrier to this calculation is the dimension of the matrices, hence each matrices product. These dimensions are evaluated in next section.

m ij (q, q ′ ) = L n-1 [q] × A n [q ′ ] 26: else if 1 < i < n -1 and i + 1 = j then 27: m ij (q, q ′ ) = L i [q] × R i+1 [q ′ ]
Note that the procedure is more general and is not restricted to binary marginals. It can also be used to compute marginals of order 3 or more since any marginal can be expressed as products of some B k and A k [z k ].

Dimensions of the matrices

A k (z k ) for a WSBM
If the graphical model is a WSBM, for which E = {(i, j) : 1 ≤ i < j ≤ n}, the price to pay in complexity for variable separation using TT is very high. Indeed, the dimensions of the matrices A k (z k ) are huge. This can be understood intuitively before developing the calculation, as each

A k (z k ) is a Kronecker product of m matrices G ij k [z k ] with m = n(n -1)/2. Each matrix G ij k [z k ] has dimension 1 × 1, 1 × Q, Q × 1, or Q × Q,
and their Kronecker product has dimension at most Q m × Q m . Hence, complexity comes from the number of pairs (i, j) in E. What follows is established for WSBM, but the approach is relevant for any pairwise graphical model.

We first establish that, for WSBM models with n individuals, the dimensions of the matrix

A k (z k ) is dim A k (z k ) = Q (k-1)(n-k+1) × Q k(n-k) Proof. Let us start from A k (z k ) = 1≤i<j≤n G ij k [z k ]
(see equation ( 9) with (i,

j) ≡ ℓ). Let us recall that if dim A = (a, a ′ ) and dim B = (b, b ′ ), then dim A ⊗ B = (ab, a ′ b ′ ). A k (z k ) is a Kronecker product of m = n(n -1)/2 matrices. Let us assume that m r matrices are in R, m q in R 1×Q , m q ′ in R Q×1 and m Q in R Q×Q . Then dim A k (z k ) = Q m q ′ +m Q × Q mq+m Q
Let us select a k ∈ {1, . . . , n} and express m q , m q ′ and m q in terms of k and n. The dimension of matrix G ij k [z k ] depends on the position of k when compared to i and j according to equation (C14). The value of G ij k [z k ] depending on i and j can be visualised in figure 1 for k / ∈ {i, j}. k being fixed, the counts of pairs (i, j) per condition in (C14) follow as:

→ condition k < i: it corresponds to all pairs (i, j) with k + 1 ≤ i < j ≤ n, hence (n -k)(n -k -1)/2 pairs → condition k = i: it corresponds to all pairs (i, j) with k = i < j ≤ n, hence n -k pairs → condition i < k < j: it corresponds to all pairs (i, j) with i < k < j, hence k -1 values for i and n -k for j, hence (k -1)(n -k) pairs → condition k = j: it corresponds to all pairs (i, j) with 1 ≤ i < k = j, hence k -1 pairs → condition k > j: it corresponds to all pairs (i, j) with 1 ≤ i < j < k, hence (k -1)(k -2)/2 pairs 1 n 1 n k k i < k i > k j < k j > k 1 1 I q
Fig. 1 This matrix displays all matrices G ij k [z k ] for k being fixed and all pairs (i, j) ∈ {1, . . . , n} × {1, . . . , n}. The gray zone represents those for which 1 ≤ i < j ≤ n. Indices i for which k < i or k > i, and indices j for which k < j or k > j are separated by dashed lines. This enables to show the blocks of indices (i, j) for which G

k number of pairs (i, j) size of ψ ij k [z k ] k < i (n -k)(n -k -1)/2 1 × 1 k = i n -k 1 × Q i < k < j (k -1)(n -k) Q × Q k = j k -1 Q × 1 k > j k(k -1)/2 1 × 1 Hence, m q ′ = k -1, m Q = (k -1)(n -k) and m q = n -k.
Let us note that the counts for k < i and k > j are not necessary since they correspond to m r , which is not involved in the size of A k (z k ).

The highest values are obtained for k ≃ n/2, which yields dim A k (z k ) ≃ Q (n/2) 2 × Q (n/2) 2 , which is worse than Q n (the number of terms in the tensor of the joint distribution). For example, for n = 20 and Q = 3, this yields dim A n/2 [z] ≃ 3 10 2 = 3 100 ≃ 5.15×10 47 . The growth of the dimension is due to the multiplicity of Kronecker products. Hence, algebra yields exact formula for computing the normalising constant, but leads to inextricable numerical difficulties for the calculation. We address in next section such difficulties. Let us mention that the dimensions are dominated by the pairs (i, j) with i < k < j: the highest contribution to the dimension of A k (z k ) is a Kronecker product of (k -1)(n -k) matrices I Q . This is the identity matrix of size k) , the number of terms in this identity matrix is x 2 , among which, in the diagonal, x are non zero and equal to 1). To the best of our knowledge, taking this sparsity into account to simplify the calculations remains an open question.

Q (k-1)(n-k) , which is sparse (if x = Q (k-1)(n-

Numerical difficulties beyond small sized problems

We have seen in section 5.3 that the dimensions of the matrices involved in matrices products for exact computation of the normalising constant or some marginals (like

W = i B i with B i = zi A i (z i )) have dimensions which can reach ≃ Q (n/2) 2 × Q (n/2) 2 .
It is not possible to perform these calculations using conventional matrix calculation algorithms, in particular the BLAS library. [START_REF] Novikov | Putting MRFs on a Tensor Train[END_REF] proposed to use a TT-format specific for matrix calculation which has been defined in [START_REF] Oseledets | Tensor-Train Ranks for Matrices and Their Inverses[END_REF] and [START_REF] Oseledets | Solution of linear systems and matrix inversion in the tt-format[END_REF]. We call here this format TTmatrix format, and not TT-format, because it is not the TT-format of a matrix read as a tensor of order 2 (this TT-format is given by the SVD of the matrix, see section 5.1). To avoid any confusion, we specify "matrix expressed as TT-matrix", or more simply "TT-matrix", instead of matrix in TT format. Calculation can be done exactly in TT-matrix format, because matrices A i (z i ) are expressed as TT-matrices of rank one, hence matrices B i are expressed as TT-matrices too, as well as their products. Definition of TT-matrices is given in Section 6.1, after the presentation of multiplexing / demultiplexing of indices. Even when using TT-matrices, numerical difficulties remain and we present in Section 6.2 how we propose to bypass them.

Matrices expressed as TT-matrices

Multiplexing / demultiplexing of indices:

Let z = (z 1 , . . . , z n ) ∈ {1, . . . , Q} n be a multi-index. An index a in {1, . . . , Q n } can be associated to z using the following construction: a = n-1 i=1 Q n-i (z i -1) + z n , which is a bijection between {1, . . . , Q} n and {1, . . . , Q n }. The transformation from a to z is called multiplexing, and the reverse operation is called demultiplexing.

Definition of a matrix expressed as a TT-matrix:

Let A be a Q n × Q n matrix. Let a ≡ (z 1 , . . . , z i , . . . , z n ) and b ≡ (t 1 , . . . , t i , . . . , t n ) be two demultiplexed indices and their multiplexed equivalent. A is said to be expressed in TT-matrix format of rank r if there exists n families of r × r matrices, where each family is composed of Q 2 matrices and is indexed by i, and where each matrix in a given family is indexed by (z i , t i ):

∀ i, (z i , t i ) ----→ M i (z i , t j ), with      i ∈ {1, . . . , n} z i , t j ∈ {1, . . . , Q} M i (z i , t j ) ∈ R r×r such that A[a, b] := A[(z 1 , . . . , z n ) ; (t 1 , . . . , t n )] = M 1 (z 1 , t 1 ) . . . M n (z n , t n ). ( 11
)
(M 1 is 1 × r and M n is r × 1 for the r.h.s. to be a scalar). Matrices M i (z i , t i ) are called the "cores" of the TT-matrix. Each matrix M i (z i , t i ) has r 2 elements (or r for M 1 (z 1 , z 1 ) and M n (z n , t n )), and there are nQ 2 of them. Hence, storage of A in TTmatrix format (i.e. storing the matrices M i (z i , t i ) instead of storing A) requires nQ 2 r 2 elements instead of Q n .

Elementary operations of matrix calculus in TT-matrix format:

In addition to the impressive storage savings achieved by using the TT-matrix format when possible, basic matrix calculus operations can be carried out directly in TTmatrix format, without the need for massive and complete storage of matrices in memory. It is possible to compute the cores of the sum A + B of two matrices, or of their product AB (see appendix B), knowing the cores of A and B only. It is possible to show that TT rank (A + B) ≤ TT rank A + TT rank B and TT rank (AB) ≤ (TT rank A)(TT rank B).

Application to the computation of W : controlling the TT-rank:

So, all calculations of B k = z k A k (z k ), W = k B k
and marginals can be computed in TT-matrix format. We have TT rank B ≤ Q, and TT rank W ≤ Q m . Of course, W is a scalar, and TT rank W = 1. However, the rank of intermediate products can be exponential with the number of terms in the product. For controlling the exponential growth of the TT-rank when computing the product i B i , algorithm 1 in [START_REF] Novikov | Putting MRFs on a Tensor Train[END_REF] consists in approximating the product by a matrix of low TT-rank, with an operation called the "rounding", which can also be performed in TT-matrix format.

Let us recall that each matrix B i in W = n i=1 B i is the sum of Q matrices A i (q) for q ∈ {1, . . . , Q}. Hence, storage of a matrix B i requires a memory space of size mQ 4 (r = Q). Using the TT-toolbox, we have carried out the control the TT-rank of the matrices after each product with rounding with a prescribed accuracy of ϵ T T = 5.10 -2 , starting with n = 8 variables and increasing n progressively. At given values of n, some new numerical difficulties appeared. We present them next in the order in which they have appeared, and the choices we made for overcoming each of them.

Numerical difficulties encountered and options for overcoming them

1. for n ≥ 24, numerical instabilities occurred because of small values for W or marginals. So we have developed a simple action to fix it (see Calculating with very small values); 2. for n ≥ 30, the rounding at prescribed accuracy yields matrices with too large TT-rank (the TT-rank is no longer controlled). Therefore, a new strategy has been carried out for the control of the TT-rank (see Strategy for controlling the TT-rank in the rounding); 3. for n ≥ 45, the number of cores became an issue. Indeed, each matrix B i has m cores. The rounding requires ∼ m 2 factorisations LQ and SVD which becomes prohibitive when m is large, and the TT-toolbox is limited to matrices with 1024 cores or less. As the number m of cores is n(n -1)/2, this corresponds to n = 45 nodes or less. Therefore, we set up core fusion to work with fewer but larger cores (see Fusion of cores).

We therefore have implementation choices that depend on the problem size, as detailed below. This allowed us to calculate the different marginals up to n = 65 for a WSBM with an assortative connectivity matrix (see Section 7 for a description of several WSBM connectivity structures).

Calculating with very small values:

Both the un-normalised marginals and the normalising constant, computed separately, are very small (see table 2): numerical simulations have shown that W ∼ exp -m . Because of the finite representation of real numbers in IEEE format and the small number of points to represent very small values, calculating marginals leads to errors when n ≥ 24. This can be solved easily by scaling the WSBM factors with a factor denoted α. Let us denote with an exponent s the scaled factors:

ψ s i,j (z i , z j ) = α ψ i,j (z i , z j ). We have P s (z, D) = α m P(z, D) because it is a product of m fac- tors. So, W s = z P s (z, D) = α m W . If m i,j (z i , z j )
is the un-normalised binary marginal of variables (i, j) for states (z i , z j ), we have m i,j (z i , z j ) = z\zi,zj P(z, D), so m s i,j (z i , z j ) = α m m i,j (z i , z j ). So, normalised marginals are unchanged by scaling, because m s i,j (z i , z j )/W s = α m m i,j (z i , z j )/α m W = m i,j (z i , z j )/W . We have chosen α = 10, which is not the optimal parameter, but it is sufficient in our case.

Strategy for the controlling the TT-rank in the rounding:

There is indeed a trade off between TT-rank and accuracy: as with Singular Value Decomposition of matrices which provides a best low rank approximation of a matrix, the higher the rank, the better the accuracy. We wish at the same time a low TT-rank for memory issues, and a high accuracy. To choose the best balance between low TTrank and high precision, we carried out a numerical experiment to compare the quality of the calculation of W with prescribed precision and several choices of prescribed TTrank. For this experiment, we considered a WSBM with Poisson distributed distances. It means that P Λ (D[i, j] | Z i = q, Z j = q ′ ) is a Poisson distribution with parameter λ qq ′ , and the connectivity matrix Λ is defined by Λ[q, q ′ ] = λ qq ′ . We used Λ =   2 10 10 10 3 8 10 10 4

 

The result is given in table 2 for different values of n (from n = 8 to n = 32). To this end, in the TT rounding algorithm we have combined the two approaches of prescribed accuracy and prescribed rank, by setting the accuracy to ϵ T T = 10 -2 , but limiting the rank growth to r T T max = 27. If the final rank reaches 27, that means we have lost some of our prescribed accuracy.

Merging of cores:

To have fewer but larger cores in TT-matrix format, with m ≤ 1024, the cores of the matrices A k (z k ) can be merged together thanks to the associativity of the Kronecker product, as seen on this very simple example: 26 1.80 × 10 -353 (16) 8.13 × 10 -375 1.80 × 10 -353 1.80 × 10 -353 1.80 × 10 -353 28 1.34 × 10 -393 (21) 1.09 × 10 -393 1.34 × 10 -393 1.34 × 10 -393 1 Table 2 Normalising constant calculated according to different accuracy and maximal rank settings in rounding: R1: ϵ T T = 5.10 -2 ; R2: r T T max = 3; R3: r T T max = 9; R4: r T T max = 27; R5: r T T max = 81. For ϵ T T = 5.10 -2 , the maximum of observed rank is given in brackets. Calculations have been done for a WSBM with Poisson distributed distances and Q = 3. smaller cores. This has been implemented for each matrix A k (z k ). Let us recall that (see formula (C15)):

M 1 ⊗ M 2 ⊗ M 3 ⊗ M 4 = (M 1 ⊗ M 2 ) ⊗ (M 3 ⊗ M 4 ):
A k (z k ) = 1≤i<j≤n G ij k (z k )
The fusion here is carried out at the level of the matrices

G ij k [z k ].
Merging can be implemented either by building groups with the same number of initial cores (uniform numbers), or by creating merged cores of the same dimensions (uniform volumes). We have tested the impact of each method on the calculation of W for different values of n for the same WSBM model than above. (see table 3). As the results are very similar, and observing that the method with uniform numbers is much simpler to implement, we have chosen to merge cores with uniform numbers. 3 Normalising constant without merging cores, or when merging them with uniform numbers or uniform volumes, for increasing number of variables (n). Calculations have been done for a WSBM with Poisson distributed distances and Q = 3.

Remark:

The solutions proposed above to overcome the numerical obstacles in carrying out the computation of marginal for large n with separation of variables in TT format have led to the possibility of computing them on data sets with n ≤ 65. Even though the algebraic calculation is exact, the rounding steps introduce approximations in the numerical implementation.

7 Comparison of the TT-based method and state-of-the-art approaches on simulated data

In this section we compare the behaviour of the TT-based method (refered to as TT in the following) for computing binary marginals in a WSBM model, to that of three classical inference methods: the exact method by complete enumeration, the Gibbs Sampler and the Mean-Field approximation. For these experiments, we considered a WSBM with Poisson distributed distances.

Different methods for computing the marginals

The first method tested for computing the binary marginals (3) is the complete enumeration of all the terms in the sum. It means Q n-2 terms for computing a single binary marginal. So this method is only available for small n. Beyond n = 12 it is not possible to store the tensor corresponding to P θ (Z, D) in memory. To circumvent this complexity two strategies are classically used. The exact computation can approximated using simulations of the model, it leads for instance to the Gibbs Sampler [START_REF] Robert | Monte Carlo statistical methods[END_REF]. The other option consists in computing the marginals on a simpler model close enough to the original one in a certain sense. This is the principle of the Mean-Field approximation used in the E step of the EM algorithm for WSBM in [START_REF] Daudin | A mixture model for random graph[END_REF]. The Gibbs Sampler can be very precise but to the price of a large computing time due to the number of simulations required. On contrary, the Mean-Field method is faster to solve, but the quality of the approximated marginals is lower. Let us detail these two methods.

Gibbs Sampler.

The Gibbs Sampler (GS, [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF]) is a Markov Chain Monte-Carlo method that exploits conditional probabilities associated to the complex propability of interest. In our case, GS consists in simulating iteratively a Markov chain, using only P θ (Z i = q|{Z j = z j } j̸ =i , D), and whose stationary distribution converges towards P θ (Z | D). In practice, we simulate a single GS run (or path) for a number N GS burnin iterations, then we sample a simulated value of Z every N GS thining new iterations until we collect L realisations. Binary marginals of P θ (Z|D) can be approximated by empirical frequencies.

Mean-Field.

With the Mean-Field approximation (MF), the conditional distribution

P θ (Z | D) is approximated by a distribution Q θ (Z | D) for which marginalisation is less costly. Q θ (Z | D) is chosen among the family Q of distributions that satisfy the hypothesis of mutual independence between the Z i s. It means that Q θ (Z | D) = n i=1 q i θ (Z i | D)
, where the q i θ (Z i | D) are the unary marginals of Q θ . The binary marginals are easily obtained as Q θ

(Z i = q, Z j = q ′ | D) = q i θ (Z i = q | D)q j θ (Z j = q ′ | D). The distribution Q θ (Z | D)
is chosen as the one that minimises the Kullback-Leibler divergence with the true distribution P θ (Z | D):

Q θ (. | D) = arg max Q∈Q KL(Q(.) | P θ (. | D)) with KL(Q(.)|P(.)) = z Q(z) ln Q(z) P(z)
. The solution is on the form of a fixed point equation which is solved by an iterative scheme (see Appendix A for the demonstration). If we denote by τ t i,q the current value of Q(Z i = q | D), we have

∀1 ≤ i ≤ n, ∀1 ≤ q ≤ Q,        τ t+1 i,q = µ t+1/2 i,q Q l=1 µ t+1/2 i,l µ t+1/2 i,q = α q j̸ =i Q q ′ =1 λ D(i,j) q,q ′ D(i,j)! exp -λ q,q ′ τ t j,q ′
The scheme is initialised by τ 0 i,q = α q . We stop the scheme when

∥∆∥ ∞ < ϵ M F with ∆ = {∆ i,q , 1 ≤ i ≤ n, 1 ≤ q ≤ Q} and ∆ i,q = |τ t+1 i,q -τ t i,q | |τ t i,q | .

The 6 WSBM structures

In the literature, different structures of the connectivity matrix Λ has been considered [START_REF] Kalmbach | Generating Synthetic Internet-and IP-Topologies using the Stochastic-Block-Model[END_REF][START_REF] Funke | Stochastic block models: A comparison of variants and inference methods[END_REF]. They correspond to different organisations of the distances between the individuals. We have selected the following ones: assortative (two examples, easy and difficult), dissortative, core periphery, ordered, and hierarchical structures. The assortative structure corresponds to the situations where the individuals are organised into well separated communities : intra group dissimilarities λ qq are small and inter group dissimilarities λ qq ′ are large. In the dissortative structure, it is the opposite. In a core periphery structure, groups are organised with a core plus groups at increasing distance from this core. Individuals in the core are close to each other while those in the periphery are more and more scattered. In an ordered structure, all the intra group dissimilarities are small then groups have a chain-like organisation where a group is close only to its direct neighbours. In the hierarchical structure, all inter dissimilarities are large and the intra dissimilarities range from small to large. For each structure, the Λ matrix that we used in the experiments, as well as a realisation of the dissimilarity matrix D for n = 100 are presented in Figure 2.

Protocol

For n = 12, exact computation is still possible and so we are able to compare the binary marginals computed by TT, GS and MF to the true ones. For each of the 6 connectivity matrix of Figure 2 we generated 10 dissimilarity matrices D and we ran inference of each of them with the 4 methods. We also ran the comparison for 10 random Λ matrices. Then, for n = 25 and n = 40, we compared TT, MF and GS behaviours, only on the two assortative structures and the core-periphery structure.

We chose the two assortative structures because, they represent typical easy and difficult inference problems. Indeed, even though they belong to the same family of assortative structures, one can see from Figure 2 than the dissimilarity matrix generated using the assortative difficult structure provides much less information on the underlying groups than for the assortative easy structure. This is due to the fact that in the connectivity matrix Λ of the assortative difficult structure, two lines are almost identical, meaning that two groups share almost the same pattern of connections.

We also considered the core-periphery structure because for n = 12 we observed that the associated unary marginals have a different pattern than for the two assortative structures (see Section 7.4.1 and Table 4). Therefore we can expect that the 3 selected structures correspond to different difficulties regarding marginal inference.

Generation of the dissimilarity matrices. We set Q = 3 and α q = 1/3 in all the experiments. Then, for a given number of individuals n and a given connectivity matrix Λ, we generated a realisation z of P α (Z) and then we generated a dissimilarity matrix D according to P Λ (D | Z = z).

Settings of the inference methods. The parameters for each method are set to the following values:

-exact method : with this method, the tensor of the joint distribution is stored in memory, and the binary and unary marginals are computed by using expression (3). It does not require any parameter setting; -TT : parameters are the accuracy ϵ T T and the maximal TT-rank r T T max . We set ϵ T T = 10 -2 and r T T max = 27 for all experiments except n = 40 and the assortative easy structure where we lowered r T T max to 9 since it was enough to reach a good agreement with the other methods.

-Gibbs Sampler : parameters are the number of iterations L GS = 5 × 10 5 , the warming phase N GS burnin = 10 3 and the thinning frequency N GS thinning = 5. There exist no universal operational rules to set these parameters [START_REF] Casella | Explaining the Gibbs Sampler[END_REF]. We choose these values empirically by checking that for n = 12 the GS estimates of unary marginals converge to the true values in most cases, and that for n = 25 or 40 GS reached an agreement with the two other approximate methods; -Mean-Field : accuracy for stopping the iteration: ϵ M F = 10 -1 . Note that these choices may be specific to Q = 3. Output. Computation of binary marginals by Gibbs Sampler can be subject to label switching (Murphy, 2012, chapter 24). It means that there should be a phase of relabeling before comparison with the binary marginals obtained with the other methods. To avoid this step which can be tricky, we compared the values of the probability that two individuals i and j belong to the same class. This probability is obtained as P SC i,j = Q q=1 P θ (Z i = q, Z j = q | D). This quantity is not sensible to label switching. Furthermore, in most applications, this quantity is more relevant than the knowledge of each probabilities P θ (Z i = q, Z j = q ′ | D).

Results

Numerical experiments for n = 12

In a preliminary step, we compared the value of the exact and approximated unary marginals. The list of all marginals for a dissimilarity matrix is a list 12 × 3 values, so we were able to compared the exact and approximated values one by one instead or summarising the difference by some statistics. We observed that the shape of the unary marginals varies a lot between the structures, and it can also vary among the 10 repetitions for a given Λ matrix due to the variance of the Poisson distribution. For assortative easy the true unary marginals are all very close to a Dirac, i.e. with one group of high probability (> 0.9) and the two others of very low probability. For assortative difficult and ordered structures, some marginals are close to Dirac and other are close to a uniform distribution between 2 groups (the third one having a probability of zero). Finally, for core-periphery, hierachic and random, the true unary marginals are more variables. The first observation is that the unary marginals computed using TT always follow the same pattern than the true ones, which is not the case for GS and MF (see Table 4). In addition quite often MF did not converge. Beyond this qualitative observation, quantitatively the TT unary marginals are always almost identical to the true one (up to the 3rd decimal). GS can be good to very good, however the quality varies with the structure, and sometimes between different runs of a same structure. When MF reaches convergence, the unary are not always of same quality than TT or GS.

Regarding now binary marginals, Table 5 provides, for each Λ structure the median of the absolute difference between the value of P SC i,j computed with the exact method and each of the approximate ones. We chose the median rather than the mean since histograms of these absolute difference are often bimodal (see Figure 1 in the SI). The value reported is the mean of this median over the 10 repetitions. Similarly Table 6 displays the mean computational time, over the 10 repetitions, for computing all unary and binary marginals, for each Λ structure.

Having in mind that the true value of P SC i,j is between 0 and 1, the global conclusion from Table 5 is that the three approximate methods provide good to very good estimates, and TT almost always leads to the lowest error. The MF estimator is also precise but less than TT. Furthermore, as above-mentioned, it can happen that the iterative scheme for solving the fixed point equation does not converge. This concerned between 0 and 10 repetitions depending on the Λ structure (see Table 5). We also observed that the MF approximation of P SC i,j is of very good quality when the unary marginals are close to 0 or 1. This occurs typically with the 'assortative easy' structure. In this case, the MF approximates very well the unary marginals. Since it can be shown (see Appendix C) that when unary marginals are equal to 0 or 1 then the binary marginals are equal to the product of the unary ones, mechanically MF leads to a good approximation of P SC i,j in this case. If the unary marginals are not close to 0 or 1, the quality if the MF approximation can decrease, this is the case for instance for the random structure.

As for the unary marginals, we observed that the quality of the approximation can vary with the structure of Λ. For instance, whatever the method, the approximation is of better quality for the "assortative easy" structure than the "assortative difficult" structure. For the latter, group 1 and group 2 have very similar connectivity pattern, therefore they can be difficult to distinguish from the dissimilarity matrix D.

Note that a method can provide a good approximation of P SC while leading to a poor approximation of the unary marginals. This is due to the fact that the computation of P SC(i, j) does not use all the binary marginals P θ (Z i = q, Z j = q ′ | D) but only those where q = q ′ . On contrary, all binary marginals are required to compute the unary marginals. It happends for instance for GS with the ordered structure or MF with the dissortative structure, so it means that in this cases GS and MF did not approximate well the binary marginals P θ (Z i = q, Z j = q ′ | D) when q ̸ = q ′ . Regarding computational time, the four methods rank as expected, with from the faster to the slowest MF, TT, GS. For TT and MF, the computational time varies with the structure of Λ.

Numerical experiments for n > 12

For n = 25 and n = 40, we compared the different methods on three structures which present three different patterns of unary marginals on the experiments for n = 12 : the two assortative structures and the core-periphery one. Table 7 provides the median of the absolute difference between the value of P SC i,j computed with two pairs of approximate methods. The median is computed over all pairs of variables i, j, and the value reported in Table 7 is the mean of this median value, over the 10 experiments. Table 8 displays the mean computational time (over the 10 experiments) for computing all unary and binary marginals.

The agreement between the TT, GS and MF estimates is very good for the assortative easy structure. In this case, the unary marginals are close to 0/1, so, as explained previously, it is a situation where MF performs well. The agreement between MF and the two other methods increases when n = 40, probably because, heuristically, Mean-Field approximation originating from statistical physics is more appropriate for systems with a large number of interacting sites.

Agreement between the three methods is weaker for the assortative difficult structure than for the assortative easy one. Still, it remains of good quality. At first sight, TT and MF are closer to each other than GS is with them. However, when looking at the biplots of all P SC i,j for all pairs (i, j) we observe that MF quite systematically returns values equal to 0, 1 or 1/2, while TT and GS provide more variable values (see Figure 3).

For the core-periphery structure, the agreement between the three method on the value of P SC, while lower than for thee assortative easy structure, is still very good.

Computing time obviously increases with n. MF remains very fast (less than a second). However, as for n = 12, the MF does not always converge (see Table 8). So even though it is the fastest method, it has to be used with care. TT remains faster than GS, however the running time of TT can vary a lot with the nature of the connectivity matrix Λ, as opposed to GS and MF. Indeed for the assortative difficult structure we had to change r T T max to 27 for n = 40 (instead of 9 for assortative easy). This was necessary to reach a good agreement with GS and MF, however rounding is the demanding task in the TT algorithm and increasing r T T max has a high impact on computing time. Note that we were able to run TT of the assortative easy structure for n = 50, but not for the assortative difficult one. The problem was probably due to a saturation of memory in the underlying fortran layer.

Fig. 3 Biplots of the approximated values of the probability to be in the same group P SC i,j per pair of methods (Q = 3, n = 40). A dot corresponds to a pair of variables (i, j) with i < j. Example for the assortative difficult structure of the WSBM connectivity matrix.

and unary marginals of any graphical model through a clever combination of TT low rank approximations of the factors. Our algebraic calculation leads to the separation of variables in the joint distribution which, by simple distributivity, allows the exact calculation of the marginals through some products of matrices. However, the size of matrices involved grows exponentially with the number m of factors in the graphical model: we have exchanged a number of sums exponential in n with products of sparse matrices of sizes exponential in m. Calculations of such products cannot be done exactly numerically. We then considered numerical approximations in order to practically implement them in the case of the Weighted Stochastic Block Model (WSBM) where the underlying graph is a clique, therefore corresponding to the most complex situation. Such an implementation requires the use of specific techniques, the main ones being working with TT-matrices, and rounding for controlling the TT-rank in matrix × matrices products.

One of the key issue for controlling the complexity of the numerical approximations is controlling the TT-rank of product of TT-matrices with rounding. What are its quality and efficiency ? It appeared in our numerical simulations that for reasonable sizes of n, like n = 12, the rounding works well, i.e. the approximate marginals are of good quality when compared to exact values, and computation times are short, whatever the structure of the connectivity matrix Λ of the WSBM. Further numerical simulations for n = 25 and 40 on three structures (called assortative easy, assortative difficult and core-periphery) showed that two elements have an impact on the computation time of the TT approach : the number n of the nodes and the easiness or difficulty with which some groups can be distinguished. The structure of Λ has also an impact on the precision (from our experiments we cannot determine if n as well because the true marginals are only available for n = 12). In structure assortative difficult, two groups are hard to distinguish since they have almost the same connection pattern, whereas in structure assortative easy, all groups are easily distinguished (see Figure 2). TT performs better on assortative easy than on assortative difficult, in the sense that it obtains results much closer to the GS and MF results. Then, in the current version of the implementation, TT cannot be ran on problems with n = 50, due to numerical limits in the rounding procedure.

We have compared our TT-based approach with classical approximation approaches like Gibbs Sampler (GS) or Mean-Field approximation (MF). GS has a sound statistical background which guarantees convergence to an unbiased estimator of the marginals [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF][START_REF] Häggström | Finite Markov Chains and Algorithmic Applications[END_REF], but there are to our knowledge nor theoretical results to tune the parameters for securing a prescribed accuracy, neither generally adopted ad-hoc technique to decide for the number of iterations. Furthermore it may lead to long calculation times, and we observe in practice that it is longer than TT. MF is often quick and accurate for computing unary marginals [START_REF] Daudin | A mixture model for random graph[END_REF], but it is not designed to compute binary marginals (apart from an coarse estimate as a product of independent unary marginals). It relies on a fixed point solution of some nonlinear system, which sometimes does not converge. TT has an intermediate running time and showed good accuracy on our experiments. So the alternative solution we propose here, based on TT approximation, appears as a useful solution for computing WSBM marginals for intermediate values of n (below 50).

The algebraic exactness of the TT approach (separation of variables) is a very good basis for developing efficient and accurate numerical approaches for calculating the marginals of pairwise graphical models. One challenge for further studies is to go beyond the current numerical limits, moving towards higher sizes, and addressing structures where certain groups are more difficult to distinguish. This should be achieved through investing efforts for repelling the numerical limits in the "rounding" step, possibly by incorporating recent developments in this field, be it with implementations with distributed memory [START_REF] Al Daas | Parallel algorithms for tensor-train arithmetic[END_REF] or, even more recently, with randomised algorithms [START_REF] Al Daas | Randomized algorithms for rounding in the Tensor-Train format[END_REF]. 4 Pattern of unary marginals for different structures of the connectivity matrix (n = 12 and Q = 3). For each structure, 10 dissimilarity matrices have been generated (corresponding to the 10 lines below each structure name) and for each of them the unary marginals of each variable have been computed by 4 different methods (corresponding to the columns in each structure). Then for a dissimilarity matrix and a method we have classified the pattern of the n unary marginals : orange cells correspond to situations where all marginals are Dirac, yellow cells to situations where marginals are either Dirac or uniform on two groups, and green cells to situations where marginals are variable. We also indicate when Label Switching occured (LS) and when MF did not converge (NC). For instance, for the first dissimilarity matrix of structure dissortative, unary marginals computed exactly or with TT were variables (e.g. (0.6, 0.3, 0.1)), while the one computed with GS were all Dirac and MF did not converge.

Λ TT Gibbs MF Assortative difficult 2.126641e-03 1.137171e-01 2.096444e-02 (10) Assortative easy 7.487330e-10 8.200865e-09 1.177558e-08 (4) Core periphery 1.430457e-03 1.981476e-04 2.106457e-02 (8) Dissortative 2.037286e-07 2.723249e-07 7.601169e-07 (8) Hierarchic (II) 3.293547e-03 6.799929e-03 2.340016e-02(5) Ordered 3.891483e-08 2.265697e-06 -(0) Random 2.067789e-02 1.037811e-03 1.925229e-01 (5) Table 5 Absolute difference between the probability to be in the same group computed by the exact method and by three different approximate methods, for n = 12 and Q = 3. The value reported is the mean (over 10 runs) of the median of the absolute errors for the n(n -1)/2 pairs of variables. For MF, we provide under parentheses the number of runs that converged. 7 Absolute difference between the probability to be in the same group computed by two approximate methods, for n = 25 and n = 40, and Q = 3. The value reported is the mean (over 10 runs) of the median of the absolute errors for the n(n -1)/2 pairs of variables. For TT, experiments on the assortative easy structure were run with r T T max = 9 and those on the assortative difficult and core-periphery structures were run with r 8 Computing time (in seconds) for computing all unary and binary marginals, for n = 25 and n = 40, and Q = 3. The value reported is the mean (over 10 runs). For MF, we provide under parentheses the number of runs that converged. For TT, experiments on the assortative easy structure were run with r T T max = 9 and those on the assortative difficult and core -periphery structures were run with r T T max = 27.

Λ

Fig. 1 Example of an histogram of the absolute difference between the value of P SC i,j computed with the exact method and with GS, for each pair (i, j) of variables.

  This leads to the following table:

  .34 × 10 -393 30 4.03 × 10 -459 (21) 4.03 × 10 -459 4.03 × 10 -459 4.03 × 10 -459 4.03 × 10 -459 32 3.43 × 10 -511 (34) 3.19 × 10 -511 3.43 × 10 -511 3.43 × 10 -511 3.43 × 10 -511

Fig. 2

 2 Fig. 2 The 6 WSBM structures used for the numerical experiments: Λ connectivity matrix and a realisation of the dissimilarity matrix D for n = 100 individuals, for each structure.

  we have two larger cores as far as dimensions are concerned instead of four 10 -101 5.38 × 10 -101 5.38 × 10 -101 5.38 × 10 -101 16 3.30 × 10 -135 (10) 3.05 × 10 -135 3.30 × 10 -135 3.30 × 10 -135 3.30 × 10 -135

	n	R 1	R 2	R 3	R 4	R 5
	8	6.21 × 10 -28 (9)	2.68 × 10 -28	6.21 × 10 -28	6.21 × 10 -28	6.21 × 10 -28
	10	3.2 × 10 -45 (12)	2.24 × 10 -46	3.20 × 10 -45	3.2 × 10 -45	3.2 × 10 -45
	12	1.0 × 10 -75 (31)	3.32 × 10 -77	6.73 × 10 -76	1.0 × 10 -75	1.0 × 10 -75
	14 1.11 × 18 5.38 × 10 -101 (8) 2.05 × 10 -161 (15) 1.14 × 10 -161 2.05 × 10 -161 2.05 × 10 -161 2.05 × 10 -161
	20	1.18 × 10 -199 (18)	2.65 × 10 -205 1.18 × 10 -199 1.18 × 10 -199 1.18 × 10 -199
	22	1.59 × 10 -247 (15)	1.80 × 10 -247 1.59 × 10 -247 1.59 × 10 -247 1.59 × 10 -247
	24	9.57 × 10 -294 (16)	1.12 × 10 -297 9.57 × 10 -294 9.57 × 10 -294 9.57 × 10 -294

  Computing time (in seconds) for computing all unary and binary marginals, for n = 12 and Q = 3. The value reported is the mean (over 10 runs). For MF, we provide under parentheses the number of runs that converged.

			Exact	TT	Gibbs	MF
	Assortative difficult 0.1748605 7.366531 354.3655 0.009956175 (10)
	Assortative easy 0.1761076 3.626953 352.9239	0.046608290 (4)
	Core periphery 0.1770441 5.178768 368.2358	0.062990053 (8)
		Dissortative 0.1755904 4.698245 357.3765	0.059073172 (8)
	Hierarchic (II) 0.1756495 7.124235 364.6317	0.021913800 (5)
		Ordered 0.1800913 3.962546 359.0718	-(0)
		Random 0.1784855 6.781812 361.7040	0.020079037 (5)
	Table 6 Λ		n = 25			n = 40
		TT vs GS	TT vs MF	GS vs MF	TT vs GS	TT vs MF	GS vs TT
	Assort. easy	0	3.8e-16	4.1e-16	9.07e-61	5.03e-49	6.16e-56
	Assort. difficult	1.02e-01	2.16e-02	1.53e-01	1.15e-01	1.93e-02	1.86e-01
	Core periphery	1.26e-06	4.91e-07	1.54e-06	3.76e-09	1.88e-11	5.73e-10
	Table

Conclusions and perspectivesIn this paper, we propose an algebraic approach to compute exactly the marginals (normalising constant, unaries, binaries, ternaries, ...) of a pairwise graphical model (the size of the factors is at most two). For this purpose, we have adapted a previous calculation of[START_REF] Novikov | Putting MRFs on a Tensor Train[END_REF] developed for approximating normalising constant

Acknowledgment

Some experiments presented in this paper were carried out using the PlaFRIM experimental testbed, supported by Inria, CNRS (LABRI and IMB), Université de Bordeaux, Bordeaux INP and Conseil Régional d'Aquitaine (see https://www.plafrim.fr), and on cluster "Hawai" of BioGeCo Research Unit at Pierroton, supported by European Grid Infrastructure programme EOSC-Pillar. M.A. Abouabdallah PhD was funded by INRAE and by INRIA.

9 Author contributions A.F., N.P., O.C. and M.A.A developed the methodology. M.A.A developed the code. A.F., N.P., O.C. and M.A.A designed the numerical experiments. M.A.A. and A.F. ran the numerical experiments. A.F. and N.P. wrote the manuscript. All authors commented on and approved the final manuscript.

Disclosure statement

The authors report that there are no competing interests to declare.

Appendix A Mean-Field fixed point equation A.1 Principle

The Mean-Field approximation consists in approximating P θ (Z | D) by another joint distribution Q θ (Z | D) for which computation of marginals is easier. Under the distribution Q θ the variables Z i are independent, i.e. Q θ (Z | D) = n i=1 q i θ (Z i | D) where q i θ (Z i | D) is the unary marginal of Z i under Q θ . The distribution Q θ is the one that miminises the Kullback-Leibler divergence between a distribution of independent variables and P θ (Z | D). We recall the definition of the Kullback-Leibler divergence between two distributions Q and P:

A.2 Mean-Field as solution of a fixed point equation

Let us define the family Q of distributions Q such that Q(Z) = n i=1 q i (Z i ), and let us denote τ i,q = q i (Z

The minimisation of KL(Q(.) | P θ (. | D)) is under the constraint that for all i, Q q=1 τ i,q = 1. So, using Lagrangian multipliers, we minimise

The variables are the τ l,q for l ∈ {1, . . . , n} and q ∈ {1, . . . Q}, and the σ l for l ∈ {1, . . . , n}. Let us compute the derivatives of G.

∂G ∂σ

Then using the fact that D and Λ are symmetric matrices, we obtain:

So ∂g(σ l ,τ l,q ) ∂τ l,q = 0 leads to

We define

Then τ l,q = µ l,q e -(σ l +1) (A3) Finally, by combining equation (A1) and equation (A3) we obtain

So the τ l,q are solution of a fixed point equation, since µ l,q is a function of the τ j,q ′ for j ̸ = l.

The associated fixed point iteration scheme is

In the case of a Poisson distribution of the dissimilarities, we obtain

Let A and B be two matrices with dimensions such that the product C = AB exists. Here, we show how to compute the cores of C expressed in TT-matrix format knowing only the cores of A and B expressed in TT-matrix format. Although it is classical, it is seldom presented in articles. We develop the calculation in detail on a simple example with 3 cores, and evaluate its complexity in a more general case.

Let A, B be TT-matrices with

Would the matrices product be commutative, reordering the terms would lead to simplifications. But it is not the case. However, reordering can be done in R between coefficients. Therefore, let us introduce the following notations:

Let us reorder the terms after the by grouping (α 1 , β 1 ) and (α 2 , β 2 ) as

Let us note that, as indices α i , β i run each over {1, . . . , r}, indices γ i run over {1, . . . , r 2 }. Then

which is a the expression of C as a TT-matrix of rank r 2 .

This calculation can be extended to TT-matrices with more than 3 cores. Let us recall that, in general, A, B have m cores, labelled by µ with 1 ≤ µ ≤ m, of size r × r (except for the first and last one), and that

with 1 ≤ i µ , j µ ≤ n. Then, there are n 2 matrices Q µ (i µ , k µ ) indexed by (i µ , k µ ) for a given µ. A matrix Q µ (i µ , k µ ) has r 2 terms, each being a sum of n terms of a matrix P µ (i µ , j µ , k µ ) (see equation ( B12)), each term of a P µ being a product of two real coefficients. So, there are n products for a term of a Q µ (i µ , k µ ), nr 2 products for a matrix Q µ (i µ , j µ ), and n 3 r 2 products for all matrices Q µ (i µ , k µ ) for a given µ, and mr 2 n 3 products for all matrices Q µ (i µ , k µ ) for all cores µ. Hence, the complexity of the matrix product in TT-format is in O(mr 2 n 3 ). This is for the product of 2 matrices. The complexity of the product of p matrices is in O(mr p n 3 ), and grows exponentially with p. Each matrix A, B has dimensions n m × n m , hence their product requires in general (n m ) 3 = n 3m products, versus mr 2 n 3 in TT-matrix format.

Appendix C Computation of matrices A i (z i ) for a WSBM

The expression of matrices A k (z k ) is given here specifically for a WSBM, with two steps, as for general graphical model (see section 4): (i) adding non essential variables, and (ii) computation of matrices A k (z k ) with mixed product property of Kronecker product.

Adding non essential variables.

Let us show it first on a toy example, with n = 4, i = 2 and j = 4 and E being {(i, j) : 1 ≤ i < j ≤ n}. We have

Completing with non essential variables (see step 2 in section 4) leads to the definition of factor ψ 24 with non essential variables z 1 , z 3 which can be written in TT format as follows (× is the matrices product)

The product of those matrices (1, 1), (1, Q), (Q, Q) and (Q, 1) is a real. This is variable separation for the factor ψ 24 (z 1 , z 2 , z 3 , z 4 ).

This can be generalized to any n by computing all

Such a calculation is done for any pair (i, j) ∈ E. We then have

where the product is the matrices product. This is separation of variables for factor ψ ij .

Computation of the matrices

Here, we can borrow the same paths as in [START_REF] Novikov | Putting MRFs on a Tensor Train[END_REF]. In our toy example, it leads to

Here, A k (z k ) is the Kronecker product of matrices which depend each on z k . Then it is a matrix which depends on z k . This can be generalised as

We then have

where the terms A k (z k ) are matrices and is the matrices product. This is variable separation for the non normalised joint distribution ψ. Let us note that the TT-rank of TT-matrix A k (z k ) is one, which is a useful property for deriving the fusion of cores in the calculation. To show this, let us write A k (z k ) in TT-matrix format. Therefore, let

Then, from the definition of the Kronecker product,

which can be written

by reordering the upper and lower indices, and shows that the TT-rank is one as a product of scalars. It can be seen on a simpler example: Let A = G 1 ⊗ G 2 be a matrix, with G 1 := G 1 [i 1 , j 1 ] and G 2 := G 2 [i 2 , j 2 ] with coefficients. Then, by definition of Kronecker product, A can be expressed as a TT-matrix of TT-rank one, by

Appendix D A situation where binary marginals are product of unary ones.

Here we show that, given two variables Z i and Z j , taking value in {1, . . . , Q}, such that there exists q with P(Z i = q) = 1, then, the binary marginals P(Z i , Z j ) are the products of the unary marginals. Let us have Q = 3. Binary marginals can be organized as elements in a 3 × 3 matrix. The unary marginals are the row-wise or column-wise sums. It can be displayed as

and all terms and marginals in this matrix being in [0, 1]. Let us assume without loss of generality that x = 1. So, x ′ = x ′′ = 0. From a ′ , b ′ , c ′ ≥ 0 and a ′ + b ′ + c ′ = x ′ = 0, it follows that a ′ = b ′ = c ′ = 0, and the same for a ′′ , b ′′ , c ′′ . So, the matrix is a b c 1 0 0 0 0 0 0 0 0 y y ′ y ′′ 1 We can observe that each term in the matrix is the product of the marginals of its row and column. For example, a = y × 1 as y = a + 0 + 0 = a. Note that is it not required that y, or y ′ or y ′′ is 1 or 0. This demonstration can be extended to Q > 3.

Computing SBM marginals with Tensor-Train decomposition: supplementary material Mohamed Anwar Abouabdallah *Corresponding author(s). E-mail(s): nathalie.peyrard@inrae.fr; Contributing authors: maabouabdallah@gmail.com; olivier.coulaud@inria.fr; alain.franc@inrae.fr;