# Computing WSBM marginals with Tensor-Train decomposition

Mohamed Anwar Abouabdallah, Olivier Coulaud, Nathalie Peyrard, Alain Franc

► To cite this version:

Mohamed Anwar Abouabdallah, Olivier Coulaud, Nathalie Peyrard, Alain Franc. Computing WSBM marginals with Tensor-Train decomposition. 2024. hal-04394024

HAL Id: hal-04394024

https://hal.inrae.fr/hal-04394024

Preprint submitted on 15 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_1.png" /></figure>

HAL
open science

# Computing WSBM conditional marginals with Tensor-Train decomposition

Mohamed Anwar Abouabdallah, Olivier Coulaud, Nathalie Peyrard, Alain Franc

<figure><img src="image_2.png" /></figure>

## To cite this version:

Mohamed Anwar Abouabdallah, Olivier Coulaud, Nathalie Peyrard, Alain Franc. Computing WSBM conditional marginals with Tensor-Train decomposition. 2025. hal-04394024v2

HAL Id: hal-04394024

https://hal.inrae.fr/hal-04394024v2

Preprint submitted on 30 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Computing WSBM marginals with Tensor-Train
decomposition

Mohamed Anwar Abouabdallah¹,², Olivier Coulaud³,
Nathalie Peyrard⁴*, Alain Franc¹,²

¹Université de Bordeaux, INRAE, UMR BIOGECO, Cestas, 33612, France.

²Pleiade, EPC INRIA-INRAE, Université de Bordeaux, Talence, 33405, France.

³Concace, EPII INRIA, Talence, 33405, France.

⁴*Université de Toulouse, INRAE, UR MIAT, Castanet-Tolosan, 31320, France.

*Corresponding author(s). E-mail(s): nathalie.peyrard@inrae.fr;

Contributing authors: maabouabdallah@gmail.com;

olivier.coulaud@inria.fr; alain.franc@inrae.fr;

Abstract

The Weighted Stochastic Block Model (WSBM) is a statistical model for unsu-
pervised 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 sev-
eral numerical challenges. First, the dimensions of the matrices involved grow
faster than exponentially with the number of variables. We bypass this diffi-
culty by using the format of TT-matrices. Second, the TT-rank of the products
grows exponentially. Then, we use a numerical approximation of matrices prod-
uct 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
# Computing WSBM conditional marginals with Tensor-Train decomposition

Mohamed Anwar Abouabdallah<fnref n="1" />,<fnref n="2" />, Olivier Coulaud<fnref n="3" />, Nathalie Peyrard<fnref n="4" />*and Alain Franc<fnref n="1" />,<fnref n="2" />

## Summary

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 distribution conditionally to the distances whose exact evaluation with exact enumeration is out of reach beyond a few individuals. This distribution is factored into binary factors and we propose to build its exact Tensor-Train (TT) decomposition from the SVD of each factor. It has the advantage to lead to an expression with variables separation. We present how to exploit this decomposition to compute unary and binary conditional marginals of group membership. They are expressed without approximation as products of matrices involved in the TT decomposition. 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 and 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, Mean-Field and the Gibbs Sampler, on the problem of binary conditional marginal inference for WSBM with 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.

**Key words:** Weighted Stochastic Block Model, Binary conditional marginals, Variables separation, Tensor-Train format, TT matrices, Low rank approximation

<footnote>*Author to whom correspondence should be addressed.</footnote>

<footnote n="1">¹ Université de Bordeaux, INRAE, UMR BIOGECO, Cestas, 33612, France</footnote>

<footnote n="2">² Pleiade, EPC INRIA-INRAE, Université de Bordeaux, Talence, 33405, France</footnote>

<footnote n="3">³ Concace, EPII INRIA, Talence, 33405, France</footnote>

<footnote n="4">⁴ Université de Toulouse, INRAE, UR MIAT, Castanet-Tolosan, 31320, France
Email: nathalie.peyrard@inrae.fr</footnote>

<footnote>**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'i 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.</footnote>
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.

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

# 1 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 Deza and Deza, 2016). An example is distances between DNA sequences in an environmental sample, which are distances on strings computed from local alignment score (Gusfield, 1997; Smith and Waterman, 1981). The question of clustering items based on their pairwise distances is also classical in time series clustering (Ansari et al., 2020). For distance-based clustering, Weighted Stochastic Block Models (WSBM) are more and more used with a large scope of application domains like social sciences (Barbillon et al., 2017), analysis of ecological interaction networks (Miele and Matias, 2017), neurology (Faskowitz et al., 2018), or numerical taxonomy (Abouabdallah et al., 2022). 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 group are not necessarily close to each others), like dissontative 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 (Dempster et al., 1977) 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. (Daudin et al., 2008; Mariadassou et al., 2010). 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
7

# 1. 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 a yes/no answer. In many situations, one knows pairwise distances or dissimilarities between a set of items (see Deza & Deza 2016). An example among many is distances between DNA sequences in an environmental sample, which are distances on strings computed from local alignment score (Gusfield 1997; Smith & Waterman 1981). The question of clustering items based on their pairwise distances is also classical in time series clustering (Ansari et al. 2020). For distance-based clustering, Weighted Stochastic Block Models (WSBMs) are more and more used with a large scope of application domains like social sciences (Barbillon et al. 2017), analysis of ecological interaction networks (Miele & Matias 2017), neurology (Faskowitz et al. 2018), or numerical taxonomy (Abouabdallah, Peyrard & Franc 2022). This statistical model formalises the groups the items belong to as latent variables. One of the main advantages is the possibility to find out groups which are not communities (the items within a group are 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 distances structures and widen the scope of patterns which can be discovered beyond community detection.

29

The unary marginals of the joint conditional distribution of a WSBM are classically computed since they inform on the individual group memberships. Binary conditional marginals are also useful. They provide a richer information that take into account correlations, and can be used to compute the probability that two individuals belong to the same group. They can also be used to build a Bethe (Yedidia, Freeman & Weiss 2005) approximation of the joint conditional distribution of a WSBM, less coarse than its Mean-field approximation that is based only on unary marginals. Therefore, in this work, we address the problem of computing binary marginals of the conditional distribution of the membership variables given the distances, 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 (up to a constant) as a product of factors, each factor involving only a subset of the variables. In a WSBM, the joint distribution of interest is the distribution of the group memberships conditionally to the distances, and
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 associ- ated 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. Wainwright and Jordan, 2008; Koller and Friedman, 2009; Murphy, 2012; Peyrard et al., 2019). 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 associ- ated 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 above- mentioned Mean-Field approach is a classical option, and MCMC techniques (Robert and Casella, 2004), 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
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. Wainwright & Jordan 2008; Koller & Friedman 2009; Murphy 2012; Peyrard et al. 2019). 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 marginals. 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 above-mentioned Mean-Field approach is a classical option, and MCMC techniques (Robert & Casella 2004), like the Gibbs sampler are another. The first one is quick but is limited to unary marginals, and the second allows calculation of binary marginals but with a much 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 approximation in the sense of Frobenius norm of this distribution by a tensor of low rank, and to perform marginalisation on the low rank approximation. There exists several notions of rank for tensors (see e.g. Coppi & Bolasco 1989; Kolda & Bader 2009), among which CP-rank (Harshman 1970; Carroll & Chang 1970), Tucker rank (Tucker 1966; Kroonenberg & de Leeuw 1980; de Lathauwer, de Moor & Vandewalle 2000), and TT-rank (Oseledets 2009, 2011; Al Daas, Ballard & Benner 2022). It is now acknowledged that Tensor-Train (TT) decomposition is equivalent to the notion of Matrix Product States (MPS) which has been widely developed in statistical physics (Fannes,
approximation in the sense of Frobenius norm of this distribution by a tensor of low
rank, and to perform marginalisation on the low rank approximation. There exists
several notions of rank for tensors (see e.g. Coppi and Bolasco, 1989; Kolda and Bader,
2009), among which CP-rank (Harshman, 1970; Carroll and Chang, 1970), Tucker
rank (Tucker, 1966; Kroonenberg and de Leeuw, 1980; de Lathauwer et al., 2000),
and TT-rank (Oseledets, 2009, 2011; Al Daas et al., 2022). It is now acknowledged
that TT decomposition is equivalent to the notion of Matrix Product States which
has been widely developed in statistical physics (Fannes et al., 1992; Orus, 2014).

Approximation of tensors which correspond to the Conditional Probability Tables (CPT) of a Bayesian Network (BN, directed graphical models) by decompositions based on CP-rank ones tensors has been exploited in Savicky and Vomlel (2007). The authors establish theoretical results on the closed-form expression of the minimal rank-one decomposition for some particular forms of CPT and propose a numerical method to compute a rank-one decomposition for the other situations. In the same spirit, Wrigley et al. (2017) propose an approximation of the factors of a graphical model by rank-one tensors decompositions to implement the Junction Tree algorithm. Marginalisation becomes easy because approximated factors are fully factored (i.e. product of functions of scope of size 1). Multiplication is not as easy and they propose to sample to approximate the product of two factors. However, finding the best low CP-rank approximation of a given tensor is not an easy task. It has been shown that such a problem may be ill-posed because the manifold of tensors of a given CP-rank may not be closed (de Silva and Lim, 2008). Close to its boundary, numerical algorithm may not converge.

The Tucker decomposition is exploited in Lyu et al. (2023). 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 (Oseledets,
2009, 2011), and as we will see the TT format is well adapted for marginalisation
because it leads to variables separation. Ducamp et al. (2020) 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 Novikov et al. (2014) to compute the normalising constant and
Nachtergaele & Werner 1992; Orus 2014).

81

Approximation of tensors which correspond to the Conditional Probability Tables
(CPT) of a Bayesian Network (BN, directed graphical model) by decompositions based on
CP-rank ones tensors has been exploited in Savicky & Vomlel (2007). The authors establish
theoretical results on the closed-form expression of the minimal rank-one decomposition
for some particular forms of CPT and propose a numerical method to compute a rank-one
decomposition for the other situations. In the same spirit, Wrigley, Lee & Ye (2017) propose
an approximation of the factors of a graphical model by rank-one tensors decompositions to
implement the Junction Tree algorithm. Marginalisation becomes easy because approximated
factors are fully factored (i.e. product of functions of scope of size 1). Multiplication is not as
easy and they propose to sample to approximate the product of two factors. However, finding
the best low CP-rank approximation of a given tensor is not an easy task. It has been shown
that such a problem may be ill-posed because the manifold of tensors of a given CP-rank
may not be closed (de Silva & Lim 2008). Close to its boundary, a numerical algorithm may
not converge.

96

The Tucker decomposition is exploited in Lyu, Xia & Zhang (2023). 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.

105

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 (Oseledets 2009, 2011),
and, as we will see, the TT format is well adapted for marginalisation because it leads to
variables separation. Ducamp et al. (2020) 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 Novikov et al. (2014) to compute the normalising
constant and the unary marginals of any graphical model by first finding a best low rank
TT-approximation of each factor and, second, writing the product of the approximations of

100

101

102

103

104

106

107

108

109

110

111

112

113

114

115

116
the unary marginals of any graphical model by first finding a best low rank TT-approximation 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 Novikov et al. (2014). 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 marginals. 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 Novikov et al. (2014), we show that it leads to numerical locks for which some solutions are proposed, essentially based on the rounding (Oseledets, 2011; Al Daas et al., 2022). 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.

## 2 Weighted Stochastic Block Model

The Stochastic Block Model (SBM, Holland et al., 1983; Daudin et al., 2008) 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 \times n$, with element $D[i, j]$. The latent (unobserved) variables are the group memberships of each individual: $Z_i \in \{1, \dots, 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 $\alpha = (\alpha_1, \dots, \alpha_Q)$, such that $P(Z_i = q) = \alpha_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
117 the factors as a tensor in TT-format by a clever utilisation of the mixed-product property of
118 the Kronecker product.

119

120 Our work focuses on WSBMs, and more generally of a graphical model with binary
121 factors and proposes a way to compute binary marginals by deriving the exact TT format of
122 the tensor associated to the joint distribution of the graphical model. The article is organised
123 as follows. We first recall some definitions that we will use throughout this paper, on WSBM
124 in section 2, and on tensors in section 3. The material in these sections is standard. In section
125 4, we recall the main steps of the approach of Novikov et al. (2014). Our contribution is
126 based on an extension of their work, and is developed in the next 3 sections. In section
127 5, we observe that the Singular Value Decomposition (SVD) of the matrices which are
128 tables of binary factors is a variable separation between indices of rows and of columns.
129 When plugged in step 1 in Novikov procedure, this leads to writing the joint distribution
130 of a pairwise graphical model exactly as a tensor in TT format, without approximation.
131 Distributivity of multiplication over addition in the ring of matrices leads to an algebraically
132 simple calculation of the marginals. We develop a approach by recursion for the simultaneous
133 calculation of all binary marginals, over all variables and all possible values the variables,
134 which extends Novikov’s procedure that is presented only for the normalising constant and
135 unary marginals. If such an approach is algebraically straightforward, it leads to sums and
136 products between matrices with huge dimensions. In section 6, elaborating on the TT-matrix
137 approach used in Novikov et al. (2014), we show that it leads to numerical locks for which
138 some solutions are proposed, essentially based on the rounding (Oseledets 2011; Al Daas,
139 Ballard & Benner 2022). Finally, in section 7, we compare both in time and accuracy
140 the computation of binary marginals by our TT approach and by three classical methods:
141 exact enumeration, simulation with a Gibbs sampler, and as fixed point of a Mean-Field
142 approximation.

143

## 2. Weighted Stochastic Block Model

144 The Stochastic Block Model (SBM, Holland, Laskey & Leinhardt 1983; Daudin, Picard
145 & Robin 2008) was originally defined for identifying groups in a set of individuals connected
146 by binary relationships (two individuals are either connected or not). In this work, we
147 consider the extension to dissimilarity matrices, where to each pair of individuals is attached
148 a weight which is a measure of dissimilarity between the two individuals.

149

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*
the distribution of $D[i, j]$ conditionally to $(Z_i, Z_j)$ depends on a parameter $\Lambda$ and we denote $\theta = (\alpha, \Lambda)$ the model's parameters. For instance, $\mathbb{P}(D[i, j] | Z_i = q, Z_j = q')$ can be the Poisson distribution with parameter $\lambda_{qq'}$. The connectivity matrix of the model is the $Q$ by $Q$ matrix $\Lambda$ such that $\Lambda[q, q'] = \lambda_{qq'}$. Let $Z = (Z_1, \dots, Z_n)$. Then, the joint probability distribution of $Z$ and $D$ is:

$$ \mathbb{P}_{\theta}(Z, D) = \mathbb{P}_{\Lambda}(D | Z)\mathbb{P}_{\alpha}(Z) \\ = \left( \prod_{i=1}^{n} \prod_{i>j}^{n} \mathbb{P}_{\Lambda}(D[i, j] | Z_i, Z_j) \right) \prod_{i=1}^{n} \mathbb{P}_{\alpha}(Z_i) \quad (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 $\mathbb{P}_{\theta}(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 $\{\mathbb{P}_{\theta}(Z_i = q | D)\}_{1 \le q \le 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: $\mathbb{P}_{\theta}(Z_i = q, Z_j = q | D)$. Therefore, the distributions of interest are not only the conditional distribution $\mathbb{P}_{\theta}(Z | D)$ but also its unary and binary marginals. The probability $\mathbb{P}_{\theta}(Z | D)$ is proportional to $\mathbb{P}_{\theta}(Z, D)$,

$$ \mathbb{P}_{\theta}(Z | D) = \mathbb{P}_{\theta}(Z, D)/W \quad (2) $$

and the normalising constant $W$ is equal to the likelihood of the observed dissimilarity matrix $D$:

$$ W = \sum_{z \in \{1, \dots, Q\}^n} \mathbb{P}_{\theta}(Z = z, D) \\ = \mathbb{P}_{\theta}(D) $$

The binary marginal for two individuals, for instance 1 and 2, is defined by

$$ \mathbb{P}_{\theta}(Z_1 = q, Z_2 = q' | D) = \sum_{z_3=1}^{Q} \dots \sum_{z_n=1}^{Q} \mathbb{P}_{\theta}(Z_1 = q, Z_2 = q', Z_3 = z_3, \dots, Z_n = z_n | D) \\ = \frac{1}{W} \sum_{z_3=1}^{Q} \dots \sum_{z_n=1}^{Q} \mathbb{P}_{\theta}(Z_1 = q, Z_2 = q', Z_3 = z_3, \dots, Z_n = z_n, D) \quad (3) $$

The unary marginal $\mathbb{P}_{\theta}(Z_i = q | D)$ is obtained by marginalising any binary marginal involving $Z_i$:

$$ \mathbb{P}_{\theta}(Z_i = q | D) = \sum_{q'=1}^{Q} \mathbb{P}_{\theta}(Z_i = q, Z_j = q' | D) $$
of size $n \times n$, where element $D[i, j]$ is the distance between individuals $i$ and $j$. The latent (unobserved) variables are the group memberships of each individual: $Z_i \in \{1, \dots, 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 $\alpha = (\alpha_1, \dots, \alpha_Q)$, such that $\Pr(Z_i = q) = \alpha_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 $\Lambda$ and we denote $\theta = (\alpha, \Lambda)$ the model's parameters. Let $Z = (Z_1, \dots, Z_n)$. The joint probability distribution of $Z$ and $D$ is:

$$ \begin{aligned} \Pr_{\theta}(Z, D) &= \Pr_{\Lambda}(D | Z)\Pr_{\alpha}(Z) \\ &= \left( \prod_{i=1}^{n} \prod_{i>j} \Pr_{\Lambda}(D[i, j] | Z_i, Z_j) \right) \prod_{i=1}^{n} \Pr_{\alpha}(Z_i) \end{aligned} \quad (1) $$

For instance, $\Pr(D[i, j] | Z_i = q, Z_j = q')$ can be the Poisson distribution with parameter $\lambda_{qq'}$. The connectivity matrix of the model is the $Q$ by $Q$ matrix $\Lambda$ such that $\Lambda[q, q'] = \lambda_{qq'}$.

The probability $\Pr_{\theta}(Z | D)$ is proportional to $\Pr_{\theta}(Z, D)$,

$$ \Pr_{\theta}(Z | D) = \Pr_{\theta}(Z, D)/W, \quad (2) $$

where the normalising constant $W$ is equal to the likelihood of the observed dissimilarity matrix $D$:

$$ \begin{aligned} W &= \sum_{z \in \{1, \dots, Q\}^n} \Pr_{\theta}(Z = z, D) \\ &= \Pr_{\theta}(D). \end{aligned} \quad (3) $$

When modelling a dissimilarity matrix with a WSBM the main goal is to recover each individual's membership, i.e. to restore the $Z_i$s. This is performed by maximising the conditional distribution $\Pr_{\theta}(Z | D)$. Beyond this 'definitive' affectation of an individual to a group, one may be interested by the vector of conditional probabilities that an individual belongs to each group, namely $\{\Pr_{\theta}(Z_i = q | D)\}_{1 \le q \le 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: $\Pr_{\theta}(Z_i = q, Z_j = q | D)$. Therefore, the distributions of interest are not only the joint conditional distribution $\Pr_{\theta}(Z | D)$ but also its unary and binary marginals.
In practice, un-normalised binary marginals, $\tilde{P}_{\theta}(Z_i = q, Z_j = q' | D) = P_{\theta}(Z_i = q, Z_j = z_j, D)$ are computed, by marginalising $P_{\theta}(Z, D)$ instead of $P_{\theta}(Z | D)$. Then, un-normalised unary marginals $\tilde{P}_{\theta}(Z_i = q | D) = P_{\theta}(Z_i = q, D)$ are derived by marginalisation of the $\tilde{P}_{\theta}(Z_i = q, Z_j = q' | D)$. Finally $W$ is obtained from any un-normalised unary marginal as $W = \sum_{q=1}^{Q} \tilde{P}_{\theta}(Z_i = q | D)$.

For the following, it will be convenient to cast the WSBM in the family of graphical models (Koller and Friedman, 2009). A vector of random variables $Z = \{Z_1, \dots, 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_{\theta}(Z | D)$ is that of a graphical model with only binary and unary factors. There are $\frac{n(n-1)}{2}$ binary factors, equal to $P_{\Lambda}(D[i, j] | Z_i, Z_j)$ and $n$ unary factors equal to $P_{\alpha}(Z_i)$. If computed by applying naively 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 $\phi$.

$$P_{\theta}(Z_1 = z_1, \dots, Z_n = z_n, D) \approx \prod_{i=1}^{n} \phi_i(z_i)$$

In such a case, we have separation of variables, and marginalisation is simple due to the distributivity of multiplication over addition:

$$\sum_{z_1=1}^{Q} \cdots \sum_{z_n=1}^{Q} \left( \prod_{i=1}^{n} \phi_i(z_i) \right) = \prod_{i=1}^{n} \left( \sum_{z_i=1}^{Q} \phi_i(z_i) \right) \quad (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_{\theta}(Z_1 = z_1, \dots, Z_n = z_n, D) \approx G_1(z_1) G_2(z_2) G_3(z_3) \quad \text{with} \quad \begin{cases} G_1(z_1) & \in \mathbb{R}^{1 \times r_1} \\ G_2(z_2) & \in \mathbb{R}^{r_1 \times r_2} \\ G_3(z_3) & \in \mathbb{R}^{r_2 \times 1} \end{cases} \quad (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 binary conditional marginal for two individuals, for instance 1 and 2, is defined by

$$
\begin{align}
\Pr_{\theta}(Z_1 = q, Z_2 = q' | D) & \\
&= \sum_{z_3=1}^{Q} \dots \sum_{z_n=1}^{Q} \Pr_{\theta}(Z_1 = q, Z_2 = q', Z_3 = z_3, \dots, Z_n = z_n | D) \nonumber \\
&= \frac{1}{W} \sum_{z_3=1}^{Q} \dots \sum_{z_n=1}^{Q} \Pr_{\theta}(Z_1 = q, Z_2 = q', Z_3 = z_3, \dots, Z_n = z_n, D) \tag{4}
\end{align}
$$

The unary conditional marginal Pr<sub>θ</sub>(Z<sub>i</sub> = q | D) is obtained by marginalising any binary marginal involving Z<sub>i</sub>:

$$
\Pr_{\theta}(Z_i = q | D) = \sum_{q'=1}^{Q} \Pr_{\theta}(Z_i = q, Z_j = q' | D)
$$

Since W is unkown, un-normalised binary marginals, $\tilde{\Pr}_{\theta}(Z_i = q, Z_j = q' | D) =$
$\Pr_{\theta}(Z_i = q, Z_j = z_j, D)$ are computed first, by marginalising $\Pr_{\theta}(Z, D)$ instead of
$\Pr_{\theta}(Z | D)$. Then, un-normalised unary marginals $\tilde{\Pr}_{\theta}(Z_i = q | D) = \Pr_{\theta}(Z_i = q, D)$ are
derived by marginalisation of the $\tilde{\Pr}_{\theta}(Z_i = q, Z_j = q' | D)$. Finally W is obtained from any
un-normalised unary marginal as $W = \sum_{q=1}^{Q} \tilde{\Pr}_{\theta}(Z_i = q | D)$. If computed by applying
naively expression (4), referred to as exact enumeration in the following, the complexity
of evaluation of one un-normalised binary marginal is in $O(Q^{n-2})$, and that of the binary
marginal is in $O(Q^n)$.

One way to approximate the evaluation of unary and binary marginals is to use the Mean-
Field approximation. 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 $\phi$
(meaning that under the Mean-Field approximation the variables are independent).

$$
\Pr_{\theta}(Z_1 = z_1, \dots, Z_n = z_n | D) \approx p_{\theta}^{MF}(Z_1 = z_1, \dots, Z_n = z_n | D) = \prod_{i=1}^{n} \phi_i(z_i).
$$

In such a case, we have separation of variables in the expression of the joint distribution.
Furthermore, the $\phi$s are the unary marginals of $p_{\theta}^{MF}$. The binary marginals are simply the
product of the unary marginals. Marginalisation with separation of variables can be done as
well in the case where the joint distribution of the $n$ variables is approximated by a product
of unary factors $G_i$, each being a matrix, like, for $n = 3$

$$
\mathrm{Pr}_{\theta}(Z_1 = z_1, Z_2 = z_2, Z_3 = z_3, D) \approx G_1(z_1) G_2(z_2) G_3(z_3) \quad (5)
$$
the simple case of the probability written in TT-format in equation (5), we have

$$ 
\begin{aligned}
W &= \sum_{z_1=1}^{Q} \sum_{z_2=1}^{Q} \sum_{z_3=1}^{Q} G_1(z_1) G_2(z_2) G_3(z_3) \\
&= \left( \sum_{z_1=1}^{Q} G_1(z_1) \right) \left( \sum_{z_2=1}^{Q} G_2(z_2) \right) \left( \sum_{z_3=1}^{Q} G_3(z_3) \right)
\end{aligned} 
\quad (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 Novikov et al. (2014), we use possibilities given by computing in TT-format (mainly due to Oseledets, 2009, 2011) to separate the variables for each variable $z_i$ of the WSBM in the expression of $\mathbb{P}_\theta(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 $\mathbf{E} = (E_1, \dots, E_i, \dots, E_n)$ be a family of $n$ finite dimensional real vector spaces. A tensor $\mathbf{T}$ on $\mathbf{E}$ is a multilinear form on $E_1 \times \dots \times E_n$. Then, $n$ is the order of the tensor, $E_i$ for $1 \le i \le n$ is a mode of $\mathbf{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 \times \dots \times Q_n$, with $\mathbf{T} \in \mathbb{R}^{Q_1 \times \dots \times Q_n}$. If a basis has been selected on each space $E_i$, $\mathbf{T}$ can be represented in those basis as a $n$-dimensional array, the elements of which are denoted

$$ \mathbf{T}[z_1, \dots, z_i, \dots, z_n] \quad \text{with} \quad 1 \le i \le n, \ 1 \le z_i \le Q_i, $$

or $\mathbf{T}_{z_1 \dots z_n}$.

### 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 \times E_2 \times E_3$ with indices $z_1, z_2, z_3$. A general term of $\mathbf{T}$ is $\mathbf{T}[z_1, z_2, z_3]$ and $\mathbf{T}$ has dimension $Q_1 \times Q_2 \times Q_3$ with $1 \le z_1 \le Q_1$, $1 \le z_2 \le Q_2$ and $1 \le z_3 \le Q_3$. The slice $\mathbf{T}_{E_2}(z_2)$ on mode $E_2$ with index $z_2$ or, more simply, $\mathbf{T}_2(z_2)$, is the tensor of dimension $Q_1 \times Q_3$ obtained by letting the indices in modes $E_1$ and $E_3$ ($\ne E_2$) run over all dimensions of respectively $E_1$ and $E_3$ and fixing the index for mode $E_2$ at $z_2$:

$$ \mathbf{T}_2(z_2) = (\mathbf{T}[z_1, z_2, z_3])_{z_1, z_3} \quad \text{with} \quad 1 \le z_1 \le Q_1, \ 1 \le z_3 \le Q_3. $$

or, more simply:

$$ \mathbf{T}_2(z_2) = \mathbf{T}[;, z_2, ;] $$
$$ \text{with} \quad \begin{cases} G_1(z_1) & \in \mathbb{R}^{1 \times r_1} \\ G_2(z_2) & \in \mathbb{R}^{r_1 \times r_2} \\ G_3(z_3) & \in \mathbb{R}^{r_2 \times 1} \end{cases} $$

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 can be exactly decomposed in a TT format, which paves the way for its marginalisation. Indeed, in the simple case of a distribution written in TT-format like in equation (5), we have

$$ \begin{aligned} W &= \sum_{z_1=1}^{Q} \sum_{z_2=1}^{Q} \sum_{z_3=1}^{Q} G_1(z_1) G_2(z_2) G_3(z_3) \\ &= \left( \sum_{z_1=1}^{Q} G_1(z_1) \right) \left( \sum_{z_2=1}^{Q} G_2(z_2) \right) \left( \sum_{z_3=1}^{Q} G_3(z_3) \right) \end{aligned} $$

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 Novikov et al. (2014), we use possibilities given by computing in TT-format (mainly due to Oseledets 2009, 2011) to separate each variable $z_i$ of the WSBM in the expression of $\text{Pr}_{\theta}(Z, D)$. This separation enables to compute efficiently the marginals of the joint distribution.

The work by Novikov et al. (2014) applies to distributions from the family of graphical models (Koller & Friedman 2009). A vector of random variables $Z = \{Z_1, \dots, 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 $\text{Pr}_{\theta}(Z | D)$ is that of a graphical model with only binary and unary factors. There are $\frac{n(n-1)}{2}$ binary factors, corresponding to the $\text{Pr}_{\Lambda}(D[i, j] | Z_i, Z_j)$ and $n$ unary factors corresponding to the $\text{Pr}_{\alpha}(Z_i)$.

### 3. Tensors and Tensor Trains

Before introducing the work by Novikov et al. (2014) and our method for computing binary marginals in a WSBM, and more generally in any pairwise graphical model, we recall here the main definitions and results on tensors we will use in the rest of the article.

**Tensor:** Let $\mathbf{E} = (E_1, \dots, E_i, \dots, E_n)$ be a family of $n$ finite dimensional real vector spaces. A tensor $\mathbf{T}$ on $\mathbf{E}$ is a multilinear form on $E_1 \times \dots \times E_n$. The quantity $n$ is the order of the tensor, $E_i$ for $1 \le i \le n$ is a mode of $\mathbf{T}$ (often simplified as $i$) and the dimension
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 \le z_i \le Q_i$):

$$T_i(z_i) = T[:..., :; z_i, :; ..., :]$$

### Frobenius norm of a tensor:

The set of tensors of a given 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 = \sum_{i_1=1}^{Q_1} \dots \sum_{i_n=1}^{Q_n} T^2[i_1 \dots 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 Oseledets (2009, 2011), well adapted to separation of variables (it can be read as a local, or componentwise, separation of variables). A tensor $T = (T[z_1, \dots, z_n])_{z_1, \dots, z_n}$ is in TT-format if there exists matrices $G_i(z_i)$ for $1 \le i \le n$ such that

$$\forall (z_1, \dots, z_n), \quad \text{with } 1 \le z_i \le Q_i, \quad T[z_1, \dots, z_n] = G_1(z_1) \dots G_n(z_n),$$

(Oseledets, 2011, equation (1.2)). The matrices $G_i(z_i)$ are called the *cores* of the TT-decomposition 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} \times r_i$, with $r_0 = r_n = 1$. These are called the ranks of the TT-format. The TT-rank of $T$ is

$$\mathbf{r} = (r_0, \dots, r_n).$$

The TT-decomposition of $T$ can be developed by defining componentwise for each mode $i$

$$G_i(z_i) = (G_i(z_i)[\alpha_{i-1}, \alpha_i])_{\alpha_{i-1}, \alpha_i} \in \mathbb{R}^{r_{i-1} \times r_i} \quad \text{with} \quad \begin{cases} 1 \le \alpha_{i-1} \le r_{i-1} \\ 1 \le \alpha_i \le r_i \end{cases}$$

and the tensor $G_i$ of order 3 of dimension $r_{i-1} \times Q_i \times r_i$ with coefficients $G_i[\alpha_{i-1}, z_i, \alpha_i]$. The matrix $G_i(z_i)$ is the slice $z_i$ of second mode of tensor $G_i$. This leads to

$$\forall (z_1, \dots, z_n), \quad T[z_1, \dots, z_n] = \sum_{\alpha_1=1}^{r_1} \sum_{\alpha_2=1}^{r_2} \dots \sum_{\alpha_{n-1}=1}^{r_{n-1}} G_1[z_1, \alpha_1] G_2[\alpha_1, z_2, \alpha_2] \dots G_n[\alpha_{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
$Q_i$ of mode $i$ is the dimension of the vector space $E_i$. The multilinear dimension of the tensor is denoted $Q_1 \times \dots \times Q_n$, with $T \in \mathbb{R}^{Q_1 \times \dots \times Q_n}$. 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, \dots, z_i, \dots, z_n] \quad \text{with} \quad 1 \le i \le n, \quad 1 \le z_i \le Q_i.$$

**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 \times E_2 \times 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 \times Q_2 \times Q_3$ with $1 \le z_1 \le Q_1$, $1 \le z_2 \le Q_2$ and $1 \le z_3 \le Q_3$. The slice $T_{E_2}(z_2)$ on mode $E_2$ with index $z_2$ or, more simply, $T_2(z_2)$, is the tensor of dimension $Q_1 \times Q_3$ obtained by letting the indices in modes $E_1$ and $E_3$ 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])_{1 \le z_1 \le Q_1, 1 \le z_3 \le Q_3},$$

or, more compactly:

$$T_2(z_2) = T[\cdot, z_2, \cdot]$$

This can be extended with some technicalities in notations to any tensor, by selecting a mode $E_i$ and an index $z_i$ ($1 \le z_i \le Q_i$). The slide $T_i(z_i)$ of tensor $T$ is

$$T_i(z_i) = T[\cdot, \cdot, \cdot, z_i, \cdot, \cdot, \cdot]$$

**Frobenius norm of a tensor:** The set of tensors of a given 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 = \sum_{z_1=1}^{Q_1} \dots \sum_{z_n=1}^{Q_n} T^2[z_1 \dots z_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 Oseledets (2009, 2011), well adapted to separation of variables (it can be read as a local, or componentwise, separation of variables). A tensor $T = (T[z_1, \dots, z_n])_{z_1, \dots, z_n}$ is in TT-format if there exists matrices $G_i(z_i)$ for $1 \le i \le n$ such that

$$\text{for all } (z_1, \dots, z_n), \text{ with } 1 \le z_i \le Q_i, \quad T[z_1, \dots, z_n] = G_1(z_1) \dots G_n(z_n),$$
rank may be huge ($r \sim \sqrt{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. Oseledets (2011) 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-by-vector 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

$$ \forall 1 \le z_1, z_2, z_3 \le Q, \quad \mathbf{T}[z_1, z_2, z_3] = G_1(z_1) G_2(z_2) G_3(z_3) \quad \text{with} \quad \begin{cases} G_1(z_1) & \in \mathbb{R}^{1 \times r_1} \\ G_2(z_2) & \in \mathbb{R}^{r_1 \times r_2} \\ G_3(z_3) & \in \mathbb{R}^{r_2 \times 1} \end{cases} $$

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

$$ \begin{aligned} m_2(z_2) &= \sum_{z_1, z_3=1}^{Q} \mathbf{T}[z_1, z_2, z_3] \\ &= \sum_{z_1, z_3=1}^{Q} G_1(z_1) G_2(z_2) G_3(z_3) \\ &= \left( \sum_{z_1=1}^{Q} G_1(z_1) \right) G_2(z_2) \left( \sum_{z_3=1}^{Q} G_3(z_3) \right) \\ &= B_1 G_2(z_2) B_3 \end{aligned} \quad (7) $$

with

$$ B_1 = \sum_{z_1=1}^{Q} G_1(z_1) \in \mathbb{R}^{1 \times r_1}, \quad B_3 = \sum_{z_3=1}^{Q} G_3(z_3) \in \mathbb{R}^{r_2 \times 1}. $$

This can be extended in a straightforward way to tensors of any order $n$ with unequal dimensions.

## 4 TT to compute the normalising constant of a graphical model

In this section we present in detail the approach developed in Novikov et al. (2014). They use the TT-approximation of the factors of a graphical model to compute the
(Oseledets 2011, equation (1.2)). The matrices $G_i(z_i)$ are called the cores of the TT-decomposition of $T$. Each core depends only on one of the $Z_i$s: this is what is referred to as variable separation for the modes of the tensor. The dimensions of matrix $G_i(z_i)$ are $r_{i-1} \times r_i$, with $r_0 = r_n = 1$. These are called the ranks of the TT-format. The TT-rank of $T$ is

$$ \mathbf{r} = (r_0, \dots, r_n). $$

The TT-decomposition of $T$ can be developed by defining componentwise for each mode $i$

$$ G_i(z_i) = (G_i(z_i)[\alpha_{i-1}, \alpha_i])_{\alpha_{i-1}, \alpha_i} \in \mathbb{R}^{r_{i-1} \times r_i} \quad \text{with} \quad \begin{cases} 1 \le \alpha_{i-1} \le r_{i-1} \\ 1 \le \alpha_i \le r_i \end{cases} $$

and the tensor $G_i$ of order 3 of dimension $r_{i-1} \times Q_i \times r_i$ with coefficients $G_i[\alpha_{i-1}, z_i, \alpha_i]$. The matrix $G_i(z_i)$ is the slice $z_i$ of second mode of tensor $G_i$. This leads to

$$ \text{for all } (z_1, \dots, z_n), \quad \mathbf{T}[z_1, \dots, z_n] = \sum_{\alpha_1=1}^{r_1} \sum_{\alpha_2=1}^{r_2} \dots \sum_{\alpha_{n-1}=1}^{r_{n-1}} \mathbf{G}_1[z_1, \alpha_1] \mathbf{G}_2[\alpha_1, z_2, \alpha_2] \dots \mathbf{G}_n[\alpha_{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 TT-rank may be huge ($r \sim \sqrt{Q^{n-1}/n}$) in general). A tensor $T$ and a TT-rank $r$ being given, a classical problem is to find a tensor in TT format of TT-rank $r$ which is closest to $T$ with Frobenius norm. Oseledets (2011) has proposed a robust algorithm, called TT-SVD, to solve this optimisation problem (see Oseledets 2011, Algorithm 1). Moreover, 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-by-vector 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.

**Tensor marginalisation using the TT format:** Marginalisation in tensors is an operation associated with slicing. It consists in (i) selecting a slice, and (ii) summing up all terms which are not in the slice. If the slice is empty (no particular mode selected), the sum is over all the terms in the tensor. Here, we show why TT format is suitable for tensor marginalisation, for
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.

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

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, $\psi_{ij}$ is a binary factor of a graphical model on edge $(i, j)$, represented by a matrix $\Psi_{ij}$, if the group membership of individual $i$ (resp. $j$) is $z_i$ (resp. $z_j$), both in state space $S$, we denote by $\psi_{ij}(z_i, z_j)$ the factor as a function from $S \times S$ on $\mathbb{R}$ and by $\Psi_{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_\ell \subset I$ for $\ell \in \{1, ..., m\}$. Let us denote $n_\ell = |A_\ell|$, and by $z_\ell$ the groups of the individuals in $A_\ell$ i.e. $z_\ell = (z_i)_{i \in A_\ell}$. Then, the un-normalised distribution of the graphical model can be written as a product of $m$ factors indexed by $\ell$:

$$ \psi(z_1, \dots, z_n) = \prod_{\ell=1}^{m} \psi_{\ell}(z_{\ell}) $$

The objective of the work in Novikov et al. (2014) is to compute the normalising constant

$$ W = \sum_{z_1=1}^{Q} \dots \sum_{z_n=1}^{Q} \left( \prod_{\ell=1}^{m} \psi_{\ell}(z_{\ell}) \right) $$

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 $\Psi := \Psi[z_1, \dots, z_n] = \psi(z_1, \dots, z_n)$ the tensor which is the distribution table associated to the un-normalised joint distribution $\psi$.
a tensor of order 3 with equal dimension $Q$ on all modes, as

$$ \text{for all } 1 \le z_1, z_2, z_3 \le Q, \quad \mathbf{T}[z_1, z_2, z_3] = G_1(z_1) G_2(z_2) G_3(z_3) $$

$$ \text{with} \quad \left\{ \begin{array}{ll} G_1(z_1) & \in \mathbb{R}^{1 \times r_1} \\ G_2(z_2) & \in \mathbb{R}^{r_1 \times r_2} \\ G_3(z_3) & \in \mathbb{R}^{r_2 \times 1}. \end{array} \right. $$

Let us consider the slice defined by mode $E_2$ and index $z_2$. The marginal $m_2(z_2)$ of $\mathbf{T}$ is

$$ \begin{aligned} m_2(z_2) &= \sum_{z_1, z_3=1}^{Q} \mathbf{T}[z_1, z_2, z_3] \\ &= \sum_{z_1, z_3=1}^{Q} G_1(z_1) G_2(z_2) G_3(z_3) \end{aligned} $$

Then, because of distributivity of multiplication on addition in the ring of matrices

$$ \begin{aligned} m_2(z_2) &= \left( \sum_{z_1=1}^{Q} G_1(z_1) \right) G_2(z_2) \left( \sum_{z_3=1}^{Q} G_3(z_3) \right) \\ &= B_1 G_2(z_2) B_3 \end{aligned} $$

229 with

$$ B_1 = \sum_{z_1=1}^{Q} G_1(z_1) \in \mathbb{R}^{1 \times r_1}, \quad B_3 = \sum_{z_3=1}^{Q} G_3(z_3) \in \mathbb{R}^{r_2 \times 1}. $$

230 This can be extended in a straightforward way to tensors of any order $n$ with unequal dimensions.
231

232

## 4. TT to compute the normalising constant of a graphical model

233 In this section we recall the approach developed in Novikov et al. (2014). They use the 234 TT-approximation of the factors of a graphical model to compute its normalising constant. 235 This approach exploits the correspondence between a joint distribution over several variables 236 and a tensor (or a matrix if only two variables are involved).

Let us consider a graphical model with $n$ variables $z_1, \dots, z_n$, each representing the state of an individual, with $z_i \in \{1, \dots, Q\}$. The graphical model is defined on a hypergraph $\mathcal{G} = (\mathcal{I}, \mathcal{E})$ where $\mathcal{I} = \{1, \dots, n\}$ and $\mathcal{E}$ is a set of hyperedges $A_\ell \subset \mathcal{I}$ for $\ell \in \{1, \dots, m\}$, with $m = |\mathcal{E}|$. Let us denote $\mathbf{z}_\ell$ the list of the groups of the individuals in $A_\ell$ i.e. $\mathbf{z}_\ell = (z_i)_{i \in A_\ell}$. The un-normalised distribution of the graphical model can be written as a product of $m$ factors
### Step 1: approximating each factor for a given TT-format

A tensor $\Psi_{\ell}$ can be built from each factor $\psi_{\ell}$. Its order is $n_{\ell}$, the number of indices in $A_{\ell}$, and its dimensions are $Q$ for each mode. We define

$$ \Psi_{\ell}[\mathbf{z}_{\ell}] = \psi_{\ell}(\mathbf{z}_{\ell}), \quad \text{with} \quad \Psi_{\ell} \in \mathbb{R}^{Q \times \dots \times Q} $$

The first step consists in selecting a rank $r$ and computing a best TT-approximation $\tilde{\Psi}_{\ell}$ of $\Psi_{\ell}$ at rank $r$ (more rigorously, at rank $\mathbf{r} = (1, r, \dots, r, 1)$). In the TT-format, each component $\tilde{\Psi}_{\ell}[z_{i_1^{\ell}}, \dots, z_{i_{n_{\ell}}^{\ell}}]$ can be written as a product of matrices with separation of variables (here, $\mathbf{z}_{\ell} = (z_{i_1^{\ell}}, \dots, z_{i_{n_{\ell}}^{\ell}})$ where $i_{\mu}^{\ell}$ is $\mu$-th index in $A_{\ell}$):

$$ \tilde{\Psi}_{\ell}[z_{i_1^{\ell}}, \dots, z_{i_{n_{\ell}}^{\ell}}] = G_{i_1^{\ell}}^{\ell}(z_{i_1^{\ell}}) \dots G_{i_{n_{\ell}}^{\ell}}^{\ell}(z_{i_{n_{\ell}}^{\ell}}), \quad \text{with} \quad \left\{ \begin{array}{l} G_{i_1^{\ell}}^{\ell}(z_{i_1^{\ell}}) \in \mathbb{R}^{1 \times r} \\ G_{i_{\mu}^{\ell}}^{\ell}(z_{i_{\mu}^{\ell}}) \in \mathbb{R}^{r \times r} \\ G_{i_{n_{\ell}}^{\ell}}^{\ell}(z_{i_{n_{\ell}}^{\ell}}) \in \mathbb{R}^{r \times 1} \end{array} \right. \quad i_{2}^{\ell} \le i_{\mu}^{\ell} \le i_{n_{\ell}-1}^{\ell} $$

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_{\ell} = 3$, tensor $\Psi_{\ell}$ 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 $\Psi_{\ell}$ from now on, and drop the ($\cdot$). We then start with

$$ \Psi_{\ell}[z_{i_1^{\ell}}, \dots, z_{i_{n_{\ell}}^{\ell}}] = G_{i_1^{\ell}}^{\ell}(z_{i_1^{\ell}}) \dots G_{i_{n_{\ell}}^{\ell}}^{\ell}(z_{i_{n_{\ell}}^{\ell}}). $$

### 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 Novikov et al. (2014), section 5.1. The idea is to build for each hyperedge $A_{\ell}$ a new factor $\tilde{\psi}_{\ell}$ that depends on all variables and not only on $\mathbf{z}_{\ell}$ but whose value does not actually depend on $z_j$ for $j \notin A_{\ell}$:

$$ \begin{array}{l} \bar{\psi}_{\ell} : \{1, \dots, Q\}^{n} \longrightarrow \mathbb{R} \\ (z_1, \dots, z_n) \longrightarrow \psi_{\ell}(\mathbf{z}_{\ell}) \end{array} $$

The tensor $\bar{\Psi}_{\ell}$ attached to $\bar{\psi}_{\ell}$ has order $n$. We have

$$ \Psi[z_1, \dots, z_n] = \prod_{\ell=1}^{m} \bar{\Psi}_{\ell}[z_1, \dots, z_n] \qquad (8) $$

The TT-format of tensors $\bar{\Psi}_{\ell}$ with non essential variables can be derived from the TT-format of the tensor corresponding to the factor $\bar{\psi}_{\ell}$ as shown in the following toy
indexed by $\ell$:

$$ \psi(z_1, \dots, z_n) = \prod_{\ell=1}^{m} \psi_{\ell}(\mathbf{z}_{\ell}). $$

The distribution $\psi$ can be represented by a tensor $\Psi$ with $n$ modes, each of dimension $Q$:
the probability $\psi(z_1, \dots, z_n)$ is the element $\Psi[z_1, \dots, z_n]$ of the tensor. Each factor can also
be represented by a tensor. If, for example, $\psi_{ij}$ is a binary factor, it can be represented by a
matrix $\Psi_{ij}$ of dimension $Q \times Q$, with $\psi_{ij}(z_i, z_j)$ being the element $\Psi_{ij}[z_i, z_j]$ of the matrix.

The objective of the work in Novikov et al. (2014) is to compute the normalising constant

$$ W = \sum_{z_1=1}^{Q} \dots \sum_{z_n=1}^{Q} \left( \prod_{\ell=1}^{m} \psi_{\ell}(\mathbf{z}_{\ell}) \right). $$

Their method is composed of four steps. We provide in Appendix I a presentation of
these steps, described with much more details than in the original paper, and summarised
here. In step one, each factor of the graphical model is approximated by a low TT-rank
tensor in TT format. In step two, each factor is extended to form a factor that depends on
all $n$ variables, but with same value (e.g. if the factor is $\psi(z_1, z_3)$ and $n = 3$, it is replaced
by a factor $\bar{\psi}(z_1, z_2, z_3)$ with $\bar{\psi}(z_1, z_2, z_3) = \psi(z_1, z_3)$ for all $z_2$). This extension is also
expressed in TT format. In step three, the TT formats of each factor are combined to obtain
an expression of the joint probability $\psi(z_1, \dots, z_n)$ with separation of variables. Step four
consists in computing the normalizing constant $W$ based on this expression. In section 5,
we elaborate on their work to compute unary and binary marginals in a pairwise graphical
model.

## 5. Separation of variables and computation of binary marginals for pairwise graphical models

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

1. In section 5.1, we first observe that if a factor $\psi_{\ell}$ is binary, the tensor $\Psi_{\ell}$ associated to it in step 1 is a matrix, and its TT-SVD is the SVD of this matrix. This enables to express the tensor $\Psi$ in TT format with variable separation, without any approximation. Therefore step 1 is exact.
example, with $n = 7$, $A_\ell = \{2, 4, 6\}$. We then have

$$ \Psi_{\ell}[z_2, z_4, z_6] = G_2^{\ell}(z_2) G_4^{\ell}(z_4) G_6^{\ell}(z_6) \quad \text{with} \quad \left\{ \begin{array}{l} G_2^{\ell}(z_2) \in \mathbb{R}^{1 \times r} \\ G_4^{\ell}(z_4) \in \mathbb{R}^{r \times r} \\ G_2^{\ell}(z_6) \in \mathbb{R}^{r \times 1} \end{array} \right. $$

This can be sketched as

$$ \Psi_{\ell}[z_2, z_4, z_6] = \begin{array}{c} \boxed{\phantom{G_2^{\ell}(z_2)}} \\ G_2^{\ell}(z_2) \\ 1 \times r \end{array} \quad \begin{array}{c} \boxed{\phantom{G_4^{\ell}(z_4)}} \\ G_4^{\ell}(z_4) \\ r \times r \end{array} \quad \begin{array}{c} \boxed{\phantom{G_6^{\ell}(z_6)}} \\ G_6^{\ell}(z_6) \\ r \times 1 \end{array} $$

The product of those matrices of dimensions $(1, r)$, $(r, r)$ and $(r, 1)$ is a real. To build the TT-format of $\bar{\Psi}_{\ell}$, we add new cores, one per non essential variable, which are $G_1^{\ell}(z_1) = 1$, $G_3^{\ell}(z_3) = G_5^{\ell}(z_5) = \mathbb{I}_r$ and $G_7^{\ell}(z_7) = 1$ as in the following figure:

$$ \bar{\Psi}_{\ell}[z_1, z_2, z_3, z_4, z_5, z_6, z_7] = \bullet \quad \begin{array}{c} \boxed{\phantom{G_1^{\ell}(z_1)}} \\ G_1^{\ell}(z_1) = 1 \\ 1 \times 1 \end{array} \quad \begin{array}{c} \boxed{\phantom{G_2^{\ell}(z_2)}} \\ G_2^{\ell}(z_2) \\ 1 \times r \end{array} \quad \begin{array}{c} \boxed{\phantom{G_3^{\ell}(z_3)}} \\ G_3^{\ell}(z_3) = \mathbb{I}_r \\ r \times r \end{array} \quad \begin{array}{c} \boxed{\phantom{G_4^{\ell}(z_4)}} \\ G_4^{\ell}(z_4) \\ r \times r \end{array} \quad \begin{array}{c} \boxed{\phantom{G_5^{\ell}(z_5)}} \\ G_5^{\ell}(z_5) = \mathbb{I}_r G_6^{\ell}(z_6) \\ r \times r \end{array} \quad \begin{array}{c} \boxed{\phantom{G_6^{\ell}(z_6)}} \\ G_6^{\ell}(z_6) \\ r \times 1 \end{array} \quad \begin{array}{c} \bullet \\ G_7^{\ell}(z_7) = 1 \\ 1 \times 1 \end{array} $$

Finally,

$$ \bar{\Psi}_{\ell}[z_1, z_2, z_3, z_4, z_5, z_6, z_7] = G_1^{\ell}(z_1) G_2^{\ell}(z_2) G_3^{\ell}(z_3) G_4^{\ell}(z_4) G_5^{\ell}(z_5) G_6^{\ell}(z_6) G_7^{\ell}(z_7) $$

This procedure is applied to each factor of the graphical model. Then, there exists a family of cores $(G_k^{\ell}(z_k))_{k,\ell}$ with $1 \le \ell \le m$ and $1 \le k \le n$ such that

$$ \forall A_{\ell} \in E, \quad \bar{\Psi}_{\ell}[z_1, \dots, z_n] = \prod_{k=1}^{n} G_k^{\ell}(z_k) $$

### Step 3: Mixed-product property of Kronecker product:

Let us recall that $\otimes$ denotes the Kronecker product between matrices. Knowing that $xy = x \otimes y$ if $x, y \in \mathbb{R}$, equation (8) can be written

$$ \begin{aligned} \Psi[z_1, \dots, z_n] &= \bigotimes_{\ell=1}^{m} \bar{\Psi}_{\ell}[z_1, \dots, z_n] \\ &= \bigotimes_{\ell=1}^{m} \left( \prod_{k=1}^{n} G_k^{\ell}(z_k) \right) \end{aligned} $$

Let us recall the mixed-product property of Kronecker product (see Horn and Johnson, 2012): if $A, B, C, D$ are matrices with relevant dimensions for the products $AB$ and
2. We then study the computational complexity of this approach on a specific family of pairwise graphical model: the WSBM, for which there is a binary factor for each pair of individuals (section 5.3 and 5.4). We show that, even though the calculation is algebraically exact, it can be intractable for large n and we propose some low rank approximation with the rounding (section 6).

## 5.1. Adaptation of the Novikov approach

The adjustments made to the different stages of the Novikov approach for the particular case of a pairwise graphical model are presented here step by step.

**Step 1: Exact TT decomposition of each binary factor.** The main advantage of a pairwise graphical model is that the TT-format of each binary factor is an exact decomposition, not an approximation. It is obtained with the SVD of the factor. Let us have a pairwise graphical model with un-normalised distribution $\psi$ defined as follows:

$$ \psi(z_1, \dots, z_n) = \prod_{(i,j) \in \mathcal{E}} \psi_{ij}(z_i, z_j), $$

where $\mathcal{E}$ is the set of binary factors of the graphical model. Factor $\psi_{ij}$ can be given as a $Q \times Q$ matrix $\Psi_{ij}$ with $\Psi_{ij}[z_i, z_j] = \psi_{ij}(z_i, z_j)$ where $z_i$ and $z_j$ label the rows and columns of $\Psi_{ij}$. The SVD of $\Psi_{ij}$ can be written as (see e.g. Horn & Johnson 2012; Strang 2019)

$$ \Psi_{ij} = U_{ij} \Sigma_{ij} V_{ij}^{\mathrm{T}}, \quad \text{with } U_{ij}, \Sigma_{ij}, V_{ij} \in \mathbb{R}^{Q \times Q} $$

where $\Sigma_{ij}$ is a diagonal matrix, and $U_{ij}, V_{ij}$ are columnwise orthonormal. We assume that matrices are full rank, i.e. rank $\Psi_{ij} = Q$. Let us denote $M_{ij} = U_{ij} \Sigma_{ij}$. Then, $\Psi_{ij} = M_{ij} V_{ij}^{\mathrm{T}}$. So

$$ \Psi_{ij}[z_i, z_j] = \sum_{q=1}^{Q} M_{ij}[z_i, q] V_{ij}^{\mathrm{T}}[q, z_j]. $$

Let us denote by $M_{ij}[a] \in \mathbb{R}^Q$ the row a of $M_{ij}$ and by $V_{ij}^{\mathrm{T}}[b] \in \mathbb{R}^Q$ the column b of $V_{ij}^{\mathrm{T}}$. These can be considered as one row and one column matrices respectively

$$ M_{ij}[z_i] \in \mathbb{R}^{1 \times Q}, \quad V_{ij}^{\mathrm{T}}[z_j] \in \mathbb{R}^{Q \times 1}. $$

Then

$$ \psi_{ij}(z_i, z_j) = \Psi_{ij}[z_i, z_j] = M_{ij}[z_i] V_{ij}^{\mathrm{T}}[z_j]. $$
CD to be possible, then

$$(AB) \otimes (CD) = (A \otimes C) (B \otimes D).$$

This leads to

$$\Psi[z_1, \dots, z_n] = \prod_{k=1}^{n} \left( \bigotimes_{\ell=1}^{m} G_k^{\ell}[z_k] \right)$$

Let us denote

$$A_k(z_k) = \bigotimes_{\ell=1}^{m} G_k^{\ell}(z_k) \qquad (9)$$

Then, the un-normalised joint distribution $\psi$ can be expressed as

$$\psi(z_1, \dots, z_n) = \Psi[z_1, \dots, z_n] = \prod_{k=1}^{n} A_k(z_k) \qquad (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 $\psi$. We recall that

$$W = \sum_{z_1=1}^{Q} \dots \sum_{z_n=1}^{Q} \psi(z_1, \dots, z_n)$$

Then,

$$\begin{aligned} W &= \sum_{z_1=1}^{Q} \dots \sum_{z_n=1}^{Q} \left( \prod_{k=1}^{n} A_k(z_k) \right) \\ &= \prod_{k=1}^{n} \left( \sum_{z_k=1}^{Q} A_k(z_k) \right) \\ &= \prod_{k=1}^{n} B_k, \quad \text{with} \quad B_k = \sum_{z_k=1}^{Q} A_k(z_k). \end{aligned}$$

It requires $nQ$ matrix additions and $n-1$ matrix multiplications only, if all matrices $A_k(z_k)$ have been computed for all $1 \le k \le n$ and $1 \le z_k \le 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 Novikov et al. (2014). Let us denote $z_k = (z_1, \dots, z_{i-1}, z_{i+1}, \dots, z_n)$, i.e. all $z_j$ for $j \ne i$. The un-normalised unary marginal for variable $i$ is defined as

$$m_i(z_i) = \sum_{z_i \in \{1, \dots, Q\}^{n-1}} \psi(z_1, \dots, z_n)$$
This is variable separation for $\psi_{ij}(z_i, z_j)$: it is a product of 2 matrices, each eachdepending only on one variable, $z_i$ or $z_j$.

**Step 2: Adding non essential variables in factors.** Let us show first how to perform this step on a toy example with $n = 4$, $i = 2$, $j = 4$ and $\mathcal{E}$ being $\{(i, j) : 1 \le i < j \le n\}$. Factors are (1, 2), (1, 3), (1, 4), (2, 3), (2, 4), (3, 4). Let us add non essential variables on the example of factor (2, 4). We have

$$ \psi_{24}(z_2, z_4) = M_{24}[z_2] V_{24}^{\mathrm{T}}[z_4] \quad \text{with} \quad \left\{ \begin{array}{ll} M_{24}[z_2] & \in \mathbb{R}^{1 \times Q} \\ V_{24}^{\mathrm{T}}[z_4] & \in \mathbb{R}^{Q \times 1} \end{array} \right. $$

Completing with non essential variables (see step 2 in appendix I) leads to the definition of factor $\bar{\psi}_{24}$ with non essential variables $z_1, z_3$, written in TT format as follows

$$ \bar{\psi}_{24}(z_1, z_2, z_3, z_4) = G_{24}^{(1)}(z_1) G_{24}^{(2)}(z_2) G_{24}^{(3)}(z_3) G_{24}^{(4)}(z_4) $$

with, for all $z_1, z_2, z_3, z_4$:

$$ \left\{ \begin{array}{lll} G_{24}^{(1)}(z_1) & = & 1 \quad \in \mathbb{R}^{1 \times 1} \\ G_{24}^{(2)}(z_2) & = & M_{24}[z_2] \quad \in \mathbb{R}^{1 \times Q} \\ G_{24}^{(3)}(z_3) & = & \mathbb{I}_{Q,Q} \quad \in \mathbb{R}^{Q \times Q} \\ G_{24}^{(4)}(z_4) & = & V_{24}^{\mathrm{T}}[z_4] \quad \in \mathbb{R}^{Q \times 1} \end{array} \right. $$

This expression can be sketched as

$$ \bar{\psi}_{24}(z_1, z_2, z_3, z_4) = \begin{array}{lll} \bullet & \multicolumn{2}{c}{\boxed{\phantom{XX}}} \\ 1 & M_{24}[z_2] & \mathbb{I}_{Q,Q} \\ 1 \times 1 & 1 \times Q & Q \times Q \end{array} \quad V_{24}^{\mathrm{T}}[z_4] \quad Q \times 1 $$

The product of those matrices (1, 1), (1, Q), (Q, Q) and (Q, 1) is a real. Therefore, at the end of step 2, there is still an expression in TT format for each factor with variable separation, and these expressions all have the same scope.

281
Then

$$
\begin{align*}
m_i(z_i) &= \sum_{\mathbf{z} \setminus z_i} \left( \prod_{k=1}^{n} A_k(z_k) \right) \\
&= \left[ \sum_{z_1, \dots, z_{i-1}} \left( \prod_{k=1}^{i-1} A_k(z_k) \right) \right] A_i[z_i] \left[ \sum_{z_{i+1}, \dots, z_n} \left( \prod_{k=i+1}^{n} A_k(z_k) \right) \right] \\
&= \left( \prod_{k=1}^{i-1} B_k \right) A_i[z_i] \left( \prod_{k=i+1}^{n} B_k \right)
\end{align*}
$$

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 $\prod_{k=1}^{i-1} B_k$ for $i=3$ to $n$, again $n-2$ to compute recursively each $\prod_{k=i+1}^{n} B_k$ for $i=1$ to $n-2$, and finally $2(n+1)(Q-1)$ (because we use $\sum_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 \times 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 $\psi$ 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 $\psi_\ell$ is binary, the tensor $\Psi_\ell$ 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 $\Psi$ 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 $\times$ matrix product of the $A_k$), and we propose some low rank approximation with the rounding (section 6).

5 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
282 This can be generalised to any $n$ by computing all $G_{ij}^{(k)}[z_k]$ for $1 \le k \le n$ as

$$ egin{cases} k < i & \Rightarrow G_{ij}^{(k)}(z_k) = 1 & \in \mathbb{R}^{1 \times 1} \\ k = i & \Rightarrow G_{ij}^{(k)}(z_i) = M_{ij}[z_i] & \in \mathbb{R}^{1 \times Q} \\ i < k < j & \Rightarrow G_{ij}^{(k)}(z_k) = \mathbb{I}_{Q,Q} & \in \mathbb{R}^{Q \times Q} \\ k = j & \Rightarrow G_{ij}^{(k)}(z_j) = V_{ij}^{T}[z_j] & \in \mathbb{R}^{Q \times 1} \\ k > j & \Rightarrow G_{ij}^{(k)}(z_k) = 1 & \in \mathbb{R}^{1 \times 1} \end{cases} \quad (6) $$

Such a calculation is done for any pair $(i, j) \in \mathcal{E}$. We then have at the end of step 2

$$ \text{for all } 1 \le i < j \le n, \quad \text{for all } 1 \le k \le n, \quad \bar{\psi}_{ij}(z_1, \dots, z_n) = \prod_{k=1}^{n} G_{ij}^{(k)}(z_k) $$

**Step 3: Exact TT decomposition of the joint distribution.** We borrow the same paths as in Novikov et al. (2014). In our toy example, it leads to

$$ A_k(z_k) = G_{12}^{(k)}(z_k) \otimes G_{13}^{(k)}(z_k) \otimes G_{14}^{(k)}(z_k) \otimes G_{23}^{(k)}(z_k) \otimes G_{24}^{(k)}(z_k) \otimes G_{34}^{(k)}(z_k), $$

283 where the matrices $G_{ij}^{(k)}(z_k)$ are given by equation (6). This can be generalised as

$$ A_k(z_k) = \bigotimes_{1 \le i < j \le n} G_{ij}^{(k)}(z_k). \quad (7) $$

284 We then have

$$ \psi(z_1, \dots, z_n) = \prod_{k=1}^{n} A_k(z_k), \quad (8) $$

285 where the terms $A_k(z_k)$ are matrices and $\prod$ is the matrices product. This is variable separation for the non normalised joint distribution $\psi$.

**Step 4: Computing the normalizing constant and the unary marginals.** Our goal is to compute the binary marginals of the pairwise graphical model. The unary marginals and the normalising constant will be derived by summing on the appropriate variables in binary marginals, as described in section 2. We do not need step 4 of Novikov approach.
of section 4 on pairwise graphical models, and focus after on the particular case of WSBM, which are pairwise graphical models.

## 5.1 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 $\psi$ defined as follows:

$$ \psi(z_1, \dots, z_n) = \prod_{(i,j) \in E} \psi_{ij}(z_i, z_j) $$

where $E$ is the set of binary factors of the graphical model. Factor $\psi_{ij}$ can be given as a $Q \times Q$ matrix $\Psi_{ij}$ with $\Psi_{ij}[z_i, z_j] = \psi_{ij}(z_i, z_j)$ where $z_i$ and $z_j$ label the rows and columns of $\Psi_{ij}$. The SVD of $\Psi_{ij}$ can be written as (see e.g. Horn and Johnson, 2012; Strang, 2019)

$$ \Psi_{ij} = U_{ij} \Sigma_{ij} V_{ij}^{T}, \quad \text{with } U_{ij}, \Sigma_{ij}, V_{ij} \in \mathbb{R}^{Q \times Q} $$

where $\Sigma_{ij}$ is a diagonal matrix, and $U_{ij}, V_{ij}$ are columnwise orthonormal. We assume that matrices are full rank, i.e. rank $\Psi_{ij} = Q$. Let us denote $M_{ij} = U_{ij} \Sigma_{ij}$. Then, $\Psi_{ij} = M_{ij} V_{ij}^T$. So

$$ \Psi_{ij}[z_i, z_j] = \sum_{q=1}^{Q} M_{ij}[z_i, q] V_{ij}[q, z_j] $$

Let us denote by $M_{ij}[a] \in \mathbb{R}^Q$ the row $a$ of $M_{ij}$ and by $V_{ij}[b] \in \mathbb{R}^Q$ the column $b$ of $V_{ij}$. These can be considered as one row and one column matrices respectively

$$ M_{ij}[z_i] \in \mathbb{R}^{1 \times Q}, \quad V_{ij}[z_j] \in \mathbb{R}^{Q \times 1} $$

Then

$$ \psi_{ij}(z_i, z_j) = \Psi_{ij}[z_i, z_j] = M_{ij}[z_i] V_{ij}[z_j] \quad (\text{Matrices product}) $$

This is variable separation for $\psi_{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.
## 5.2. 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)$. Indeed at the end of step 3, starting from expression (8) and using variable separation when marginalising, we obtain

$$m_{ij}(q, q') = \tilde{\mathrm{Pr}}_{\theta}(Z_i = q, Z_j = q' | D) \\ = B_1 \times \dots \times B_{i-1} \times A_i(q) \times B_{i+1} \times \dots \times B_{j-1} \times A_j(q') \times B_{j+1} \times \dots \times B_n$$

with $B_k = \sum_{z_k=1}^Q A_k(z_k)$. Let us define the following quantities:

$$\begin{array}{ll} \text{for all } 1 \le i \le n, & F_{i,i} = B_i \\ \text{for all } 1 \le i \le n-1, \text{ for all } i+1 \le j \le n, & F_{i,j} = B_i \times B_{i+1} \times \dots \times B_j \end{array}$$

Depending on the values of i and j, $m_{ij}(q, q')$ have different form:

$$\begin{array}{lllcl} i = 1 & j = n & \Rightarrow & m_{1n}(q, q') = A_1(q) \times F_{2,n-1} \times A_n(q') \\ i = 1 & 2 < j < n & \Rightarrow & m_{1j}(q, q') = A_1(q) \times F_{2,j-1} \times A_j(q') \times F_{j+1,n} \\ 1 < i < n-1 & j = n & \Rightarrow & m_{in}(q, q') = F_{1,i-1} \times A_i(q) \times F_{i+1,n-1} \times A_n(q') \\ 1 < i < n-2 & i+1 < j < n & \Rightarrow & m_{ij}(q, q') = F_{1,i-1} \times A_i(q) \times F_{i+1,j-1} \times A_j(q') \times F_{j+1,n} \\ i = 1 & j = 2 & \Rightarrow & m_{12}(q, q') = A_1(q) \times A_2(q') \times F_{3,n} \\ i = n-1 & j = n & \Rightarrow & m_{in}(q, q') = F_{1,n-2} \times A_{n-1}(q) \times A_n(q') \\ 1 < i < n-1 & i+1 = j & \Rightarrow & m_{ij}(q, q') = F_{1,i-1} \times A_i(q) \times A_{i+1}(q') \times F_{i+2,n} \end{array}$$

All the $F_{i,j}$ can be computed in an efficient way recursively. Once the $F_{ij}$ are computed, all binary marginals can be computed with the following algorithm:

298
## 5.2 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)$:

$$ 
\begin{aligned}
m_{ij}(q, q') &= \tilde{P}_{\theta}(Z_i = q, Z_j = q' | D) \\
&= B_1 \times \dots \times B_{i-1} \times A_i(q) \times B_{i+1} \times \dots \times B_{j-1} \times A_j(q') \times B_{j+1} \times \dots \times B_n
\end{aligned}
$$

Let us define the following quantities:

$$ 
\begin{array}{l}
\forall 1 \le i \le n, \quad F_{i,i} = B_i \\
\forall 1 \le i \le n-1, \quad \forall i+1 \le j \le n, \quad F_{i,j} = B_i \times B_{i+1} \times \dots \times B_j
\end{array}
$$

Depending on the values of $i$ and $j$, $m_{ij}(q, q')$ can be expressed as follows:

$$ 
\begin{array}{l@{\qquad}l@{\qquad}l}
| i = 1 & j = n & \Rightarrow m_{1n}(q, q') = A_1[q] \times F_{2,n-1} \times A_n[q'] \\
| i = 1 & 2 < j < n & \Rightarrow m_{1j}(q, q') = A_1[q] \times F_{2,j-1} \times A_j[q'] \times F_{j+1,n} \\
| 1 < i < n-1 & j = n & \Rightarrow m_{in}(q, q') = F_{1,i-1} \times A_i[q] \times F_{i+1,n-1} \times A_n[q'] \\
| 1 < i < n-2 & i+1 < j < n & \Rightarrow m_{ij}(q, q') = F_{1,i-1} \times A_i[q] \times F_{i+1,j-1} \times A_j[q'] \times F_{j+1,n} \\
| i = 1 & j = 2 & \Rightarrow m_{12}(q, q') = A_1[q] \times A_2[q'] \times F_{3,n} \\
| i = n-1 & j = n & \Rightarrow m_{in}(q, q') = F_{1,n-2} \times A_{n-1}[q] \times A_n[q'] \\
| 1 < i < n-1 & i+1 = j & \Rightarrow m_{ij}(q, q') = F_{1,i-1} \times A_i[q] \times A_{i+1}[q'] \times F_{i+2,n}
\end{array}
$$

All the $F_{i,j}$ can be computed in an efficient way recursively with the following algorithm:

<caption>**Algorithm 1** Recursive computation of the $F_{ij}$</caption>

1: **inputs:** $B_i \quad \forall 1 \le i \le n$
2: **initialisation:** $\forall 1 \le i \le n, \quad F_{i,i} = B_i$
3: **for** $i = 1$ **to** $n-1$ **do**
4: &nbsp;&nbsp;&nbsp;&nbsp;**for** $j = i+1$ **to** $n$ **do**
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$F_{i,j} = F_{i,j-1} \times B_j$
6: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
7: **end for**
8: **return** $F_{i,j} \quad \forall 1 \le i < j \le n$

Once the $F_{ij}$ are computed, all binary marginals can be computed with the following algorithm:
<caption>**Algorithm 1** Computation of all binary marginals of a pairwise graphical model</caption>

1: **inputs:** $F_{i,j}$ for all $1 \le i < j \le n$, $A_k(q)$ for all $1 \le k \le n$, for all $1 \le q \le Q$
2: **for** $j = 1$ **to** $n-1$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;**for** $q' = 1$ **to** $Q$ **do**
4: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$R_j(q') = A_j(q') \times F_{j+1,n}$
5: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
6: **end for**
7: **for** $i = 2$ **to** $n$ **do**
8: &nbsp;&nbsp;&nbsp;&nbsp;**for** $q = 1$ **to** $Q$ **do**
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$L_i(q) = F_{1,i-1} \times A_i(q)$
10: **end for**
11: **end for**
12: **for** $i = 1$ **to** $n-1$ **do**
13: &nbsp;&nbsp;&nbsp;&nbsp;**for** $j = i+1$ **to** $n$ **do**
14: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for** $q = 1$ **to** $Q$ **do**
15: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for** $q' = 1$ **to** $Q$ **do**
16: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $i = 1$ and $j = n$ **then**
17: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$m_{ij}(q, q') = A_1(q) \times F_{2,n-1} \times A_n(q')$
18: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**else if** $i = 1$ and $2 < j < n$ **then**
19: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$m_{ij}(q, q') = A_1(z_1) \times F_{2,j-1} \times R_j(q')$
20: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**else if** $1 < i < n-1$ and $j = n$ **then**
21: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$m_{ij}(q, q') = L_i(q) \times F_{i+1,n-1} \times A_n(q')$
22: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**else if** $1 < i < n-2$ and $i+1 < j < n$ **then**
23: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$m_{ij}(q, q') = L_i(q) \times F_{i+1,j-1} \times R_j(q')$
24: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**else if** $i = n-1$ and $j = n$ **then**
25: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$m_{ij}(q, q') = L_{n-1}(q) \times A_n(q')$
26: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**else if** $1 < i < n-1$ and $i+1 = j$ **then**
27: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$m_{ij}(q, q') = L_i(q) \times R_{i+1}(q')$
28: **end if**
29: **end for**
30: **end for**
31: **end for**
32: **end for**
33: **return** $m_{ij}(q, q')$ for all $1 \le i < j \le n$, for all $1 \le q, q' \le Q$

Computing all the $F_{i,j}$ requires $\sum_{i=1}^{n-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
<caption>**Algorithm 2** Computation of all binary marginals of a pairwise graphical model</caption>

1: **inputs:** $F_{i,j} \quad \forall 1 \le i < j \le n$, $A_k(q) \quad \forall 1 \le k \le n, \forall 1 \le q \le Q$
2: **for** $j = 1$ **to** $n-1$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;**for** $q' = 1$ **to** $Q$ **do**
4: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$R_j[q'] = A_j[q'] \times F_{j+1,n}$
5: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
6: **end for**
7: **for** $i = 2$ **to** $n$ **do**
8: &nbsp;&nbsp;&nbsp;&nbsp;**for** $q = 1$ **to** $Q$ **do**
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$L_i[q] = F_{1,i-1} \times A_i[q]$
10: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
11: **end for**
12: **for** $i = 1$ **to** $n-1$ **do**
13: &nbsp;&nbsp;&nbsp;&nbsp;**for** $j = i+1$ **to** $n$ **do**
14: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for** $q = 1$ **to** $Q$ **do**
15: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**for** $q' = 1$ **to** $Q$ **do**
16: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $i = 1$ and $j = n$ **then**
17: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$m_{ij}(q, q') = A_1[q] \times F_{2,n-1} \times A_n[q']$
18: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**else if** $i = 1$ and $2 < j < n$ **then**
19: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$m_{ij}(q, q') = A_1[z_1] \times F_{2,j-1} \times R_j[q']$
20: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**else if** $1 < i < n-1$ and $j = n$ **then**
21: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$m_{ij}(q, q') = L_i[q] \times F_{i+1,n-1} \times A_n[q']$
22: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**else if** $1 < i < n-2$ and $i+1 < j < n$ **then**
23: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$m_{ij}(q, q') = L_i[q] \times F_{i+1,j-1} \times R_j[q']$
24: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**else if** $i = n-1$ and $j = n$ **then**
25: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$m_{ij}(q, q') = L_{n-1}[q] \times A_n[q']$
26: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**else if** $1 < i < n-1$ and $i+1 = j$ **then**
27: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$m_{ij}(q, q') = L_i[q] \times R_{i+1}[q']$
28: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end if**
29: **end for**
30: **end for**
31: **end for**
32: **end for**
33: **return** $m_{ij}(q, q') \quad \forall 1 \le i < j \le n, \forall 1 \le q, q' \le Q$

Computing all the $F_{i,j}$ requires $\sum_{i=1}^{n-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 \approx n^2Q^2$ matrix multiplications. Let us recall that the real barrier to this calculation is the
301 values for $q$, hence $(n-1)Q$ matrices products are needed to obtain all matrices $L_i(q)$. 302 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 303 calculation of $m_{ij}(q, q')$ requires 2 matrices products, hence $n(n-1)Q^2$ products for all 304 binary marginals because there are $n(n-1)/2$ pairs $(i, j)$. In total, computing all binary 305 marginals requires $n(n-1)/2 + 2(n-1)Q + n(n-1)Q^2 \approx n^2Q^2$ matrices products. 306 The real barrier to this calculation is the dimension of the matrices. These dimensions are 307 evaluated in next section. 308

309 Note that the procedure is more general and is not restricted to binary marginals. It can 310 also be used to compute marginals of order 3 or more since any marginal can be expressed as 311 products of some $B_k$ and $A_k(z_k)$. 312

## 5.3. Dimensions of the matrices $A_k(z_k)$ for a WSBM

313 If the graphical model is a WSBM, for which $\mathcal{E} = \{(i, j) : 1 \le i < j \le n\}$, the price 314 to pay in complexity for variable separation using TT is very high. Indeed, the dimensions 315 of the matrices $A_k(z_k)$ are huge. This can be understood intuitively before developing 316 the calculation, as each $A_k(z_k)$ is a Kronecker product of $m$ matrices $G_{ij}^{(k)}(z_k)$ with 317 $m = n(n-1)/2$ being the number of possible pairs of variables. Each matrix $G_{ij}^{(k)}(z_k)$ has 318 dimension $1 \times 1$, $1 \times Q$, $Q \times 1$, or $Q \times Q$, and their Kronecker product has dimension at 319 most $Q^m \times Q^m$. Hence, complexity comes from the number of pairs $(i, j)$ in $\mathcal{E}$. What follows 320 is established for a WSBM where the associated graphical representation is a clique. We first 321 establish the dimensions of the matrix $A_k(z_k)$. 322

**Property** For WSBM models with $n$ individuals, the dimensions of the matrix $A_k(z_k)$ for $1 \le k \le n$ is

$$ \dim A_k(z_k) = Q^{(k-1)(n-k+1)} \times Q^{k(n-k)}. $$

**Proof.** Let us start from equation (7):

$$ A_k(z_k) = \bigotimes_{1 \le i < j \le n} G_{ij}^{(k)}(z_k). $$

Let us recall that if $\dim A = (a, a')$ and $\dim B = (b, b')$, then $\dim A \otimes B = (ab, a'b')$. $A_k(z_k)$ is a Kronecker product of $m$ matrices. Let us assume that $m_r$ matrices are in $\mathbb{R}$, $m_q$ in $\mathbb{R}^{1 \times Q}$, $m_{q'}$ in $\mathbb{R}^{Q \times 1}$ and $m_Q$ in $\mathbb{R}^{Q \times Q}$. We have

$$ \dim A_k(z_k) = Q^{m_{q'} + m_Q} \times Q^{m_q + m_Q}. $$
dimension of the matrices, hence each matrices product. These dimensions are evaluated in next section.

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]$.

## 5.3 Dimensions of the matrices $A_k(z_k)$ for a WSBM

If the graphical model is a WSBM, for which $E = \{(i, j) : 1 \le i < j \le 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_k^{ij}[z_k]$ with $m = n(n-1)/2$. Each matrix $G_k^{ij}[z_k]$ has dimension $1 \times 1$, $1 \times Q$, $Q \times 1$, or $Q \times Q$, and their Kronecker product has dimension at most $Q^m \times 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)} \times Q^{k(n-k)} $$

*Proof.* Let us start from

$$ A_k(z_k) = \bigotimes_{1 \le i < j \le n} G_k^{ij}[z_k] $$

(see equation (9) with $(i, j) \equiv \ell$). Let us recall that if $\dim A = (a, a')$ and $\dim B = (b, b')$, then $\dim A \otimes 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 $\mathbb{R}$, $m_q$ in $\mathbb{R}^{1 \times Q}$, $m_{q'}$ in $\mathbb{R}^{Q \times 1}$ and $m_Q$ in $\mathbb{R}^{Q \times Q}$. Then

$$ \dim A_k(z_k) = Q^{m_{q'}+m_Q} \times Q^{m_q+m_Q} $$

Let us select a $k \in \{1, \dots, n\}$ and express $m_q, m_{q'}$ and $m_Q$ in terms of $k$ and $n$. The dimension of matrix $G_k^{ij}[z_k]$ depends on the position of $k$ when compared to $i$ and $j$ according to equation (C14). The value of $G_k^{ij}[z_k]$ depending on $i$ and $j$ can be visualised in figure 1 for $k \notin \{i, j\}$.

$k$ being fixed, the counts of pairs $(i, j)$ per condition in (C14) follow as:

* $\rightarrow$ condition $k < i$: it corresponds to all pairs $(i, j)$ with $k+1 \le i < j \le n$, hence $(n-k)(n-k-1)/2$ pairs

* $\rightarrow$ condition $k = i$: it corresponds to all pairs $(i, j)$ with $k = i < j \le n$, hence $n-k$ pairs

* $\rightarrow$ 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

* $\rightarrow$ condition $k = j$: it corresponds to all pairs $(i, j)$ with $1 \le i < k = j$, hence $k-1$ pairs

* $\rightarrow$ condition $k > j$: it corresponds to all pairs $(i, j)$ with $1 \le i < j < k$, hence $(k-1)(k-2)/2$ pairs
Let us select a $k \in \{1, \dots, 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 equations (6). The value of $G_{ij}^{(k)}(z_k)$ depending on $i$ and $j$ can be visualised in Figure 1 for $k \notin \{i, j\}$, i.e. when $k$ corresponds to a non essential variable.

326

<figure><img src="image_4.png" /><figcaption>Figure 1. This matrix displays the values of matrix $G_{ij}^{(k)}(z_k)$ for $k$ being fixed and all pairs $(i,j) \in \{1, \dots, n\} \times \{1, \dots, n\}$ such that $k \neq i, j$. The gray zone represents those for which $1 \le i < j \le 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_{ij}^{(k)}(z_k) = 1$ and those for which $G_{ij}^{(k)}(z_k) = \mathbb{I}_{Q,Q}$.</figcaption></figure>

327

The value of $k$ being fixed, the counts of pairs $(i, j)$ per condition in (6) follow as:

$328 \rightarrow$ condition $k < i$: it corresponds to all pairs $(i, j)$ with $k+1 \le i < j \le n$, hence $(n-k)(n-k-1)/2$ pairs

329

$330 \rightarrow$ condition $k=i$: it corresponds to all pairs $(i, j)$ with $k=i < j \le n$, hence $n-k$ pairs

$331 \rightarrow$ 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

332

$333 \rightarrow$ condition $k=j$: it corresponds to all pairs $(i, j)$ with $1 \le i < k = j$, hence $k-1$ pairs

$334 \rightarrow$ condition $k > j$: it corresponds to all pairs $(i, j)$ with $1 \le i < j < k$, hence $(k-1)(k-2)/2$ pairs

335

This leads to counts per condition summarised in Table 1

336

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

337

338
<figure><img src="image_2.png" /><figcaption>**Fig. 1** This matrix displays all matrices $G_k^{ij}[z_k]$ for $k$ being fixed and all pairs $(i,j) \in \{1, \dots, n\} \times \{1, \dots, n\}$. The gray zone represents those for which $1 \le i < j \le 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_{ij}^{(k)} = 1$ and those for which $G_{ij}^{(k)} = \mathbb{I}_Q$</figcaption></figure>

This leads to the following table:

<table><thead><tr><th>k</th><th>number of pairs (i,j)</th><th>size of ψ<sup>ij</sup><sub>k</sub>[z<sub>k</sub>]</th></tr></thead><tbody><tr><td>k &lt; i</td><td>(n - k)(n - k - 1)/2</td><td>1 &times; 1</td></tr><tr><td>k = i</td><td>n - k</td><td>1 &times; Q</td></tr><tr><td>i &lt; k &lt; j</td><td>(k - 1)(n - k)</td><td>Q &times; Q</td></tr><tr><td>k = j</td><td>k - 1</td><td>Q &times; 1</td></tr><tr><td>k &gt; j</td><td>k(k - 1)/2</td><td>1 &times; 1</td></tr></tbody></table>

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)$. $\square$

The highest values are obtained for $k \simeq n/2$, which yields $\dim A_k(z_k) \simeq Q^{(n/2)^2} \times 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] \simeq 3^{102} = 3^{100} \simeq 5.15 \times 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 $\mathbb{I}_Q$. This is the identity matrix of size $Q^{(k-1)(n-k)}$, which is sparse (if $x = Q^{(k-1)(n-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.
<table><caption>Table 1. Number of pairs $(i, j)$ with $1 \le i < j \le n$ satisfying the different conditions enumerated in (6).</caption><thead><tr><th>k</th><th>number of pairs (i, j)</th><th>size of $\overline{\psi_k^{ij}}[z_k]$</th></tr></thead><tbody><tr><td>$k < i$</td><td>$(n-k)(n-k-1)/2$</td><td>$1 \times 1$</td></tr><tr><td>$k = i$</td><td>$n-k$</td><td>$1 \times Q$</td></tr><tr><td>$i < k < j$</td><td>$(k-1)(n-k)$</td><td>$Q \times Q$</td></tr><tr><td>$k = j$</td><td>$k-1$</td><td>$Q \times 1$</td></tr><tr><td>$k > j$</td><td>$k(k-1)/2$</td><td>$1 \times 1$</td></tr></tbody></table>

339 involved in the size of $A_k(z_k)$.

340 □

341 The highest count values are obtained for $k \approx n/2$, which yields $\dim A_k(z_k) \approx Q^{(n/2)^2} \times Q^{(n/2)^2}$, which is worse than $Q^n$ (the number of sums defining $W$). For example, 342 for $n=20$ and $Q=3$, this yields $\dim A_{n/2}(z) \approx 3^{102} = 3^{100} \approx 5.15 \times 10^{47}$. The growth 343 of the dimension is due to the multiplicity of Kronecker products. Hence, algebra yields 344 exact formula for computing binary marginals, but leads to inextricable numerical difficulties 345 for the calculation. This is true for the calculation of normalising constant $W = \prod_{k=1}^n B_k$ 346 as well. We address such difficulties in Section 6. Let us mention that the dimensions are 347 dominated by the pairs $(i, j)$ with $i < k < j$: the highest contribution to the dimension of 348 $A_k(z_k)$ is a Kronecker product of $(k-1)(n-k)$ matrices $\mathbb{I}_{Q,Q}$. This is the identity matrix 349 of size $Q^{(k-1)(n-k)}$, which is sparse (if $x = Q^{(k-1)(n-k)}$, the number of terms in this identity 350 matrix is $x^2$, among which, in the diagonal, $x$ are non zero and equal to 1). 351

5.4. Complexity of the exact calculation of $W$ for a WSBM

Here, we show that the complexity of the calculation of $W$ using the TT approach 353 by $W = \prod_{k=1}^n B_k$ leads to the same complexity as the brute force, in $O(Q^n)$, despite the 354 sparsity of the matrices $A_k(z_k)$. Indeed, we have shown that 355

$$A_1(z_1) = \bigotimes_{j=2}^{n} M_{1j}[z_1] \quad \text{and} \quad B_1 = \sum_{z=1}^{Q} A_1(z_1).$$

356 $A_1(z_1)$ is a Kronecker product of $n-1$ vectors, each in $\mathbb{R}^Q$. So, computing $A_1(z_1)$ requires 357 $O(Q^{n-1})$ multiplications. $B_1$ is the sum of $Q$ vectors $A_1(z_1)$, for $Q$ values of $z_1$. So, 358 $Q \times Q^{n-1}$ multiplications are required in total to compute $B_1$. Then, without looking at 359 how to take into account the sparsity of matrices $B_i$ for $2 \le i \le n-1$, the computation 360 of $B_1$ itself (and of $B_n$) is in $O(Q^n)$ (those two matrices are not sparse). Then, as for
# 6 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 = \prod_i B_i$ with $B_i = \sum_{z_i} A_i(z_i)$) have dimensions which can reach $\approx Q^{(n/2)^2} \times Q^{(n/2)^2}$. It is not possible to perform these calculations using conventional matrix calculation algorithms, in particular the BLAS library. Novikov et al. (2014) proposed to use a TT-format specific for matrix calculation which has been defined in Oseledets et al. (2011) and Oseledets and Dolgov (2012). We call here this format TT-matrix 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.

## 6.1 Matrices expressed as TT-matrices

**Multiplexing / demultiplexing of indices:**

Let $z = (z_1, \dots, z_n) \in \{1, \dots, Q\}^n$ be a multi-index. An index $a$ in $\{1, \dots, Q^n\}$ can be associated to $z$ using the following construction: $a = \sum_{i=1}^{n-1} Q^{n-i}(z_i - 1) + z_n$, which is a bijection between $\{1, \dots, Q\}^n$ and $\{1, \dots, 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 \times Q^n$ matrix. Let $a = (z_1, \dots, z_i, \dots, z_n)$ and $b = (t_1, \dots, t_i, \dots, 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 \times 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)$:

$$ \forall i, \quad (z_i, t_i) \longrightarrow M_i(z_i, t_j), \quad \text{with} \quad \begin{cases} i & \in \{1, \dots, n\} \\ z_i, t_j & \in \{1, \dots, Q\} \\ M_i(z_i, t_j) & \in \mathbb{R}^{r \times r} \end{cases} $$

such that

$$ A[a, b] := A[(z_1, \dots, z_n); (t_1, \dots, t_n)] = M_1(z_1, t_1) \dots M_n(z_n, t_n). \quad (11) $$

$(M_1$ is $1 \times r$ and $M_n$ is $r \times 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
361 computing $W$, one must compute $B_1$ and $B_n$ beforehand, the complexity of the calculation
362 of $W$ is in $O(Q^n)$ at least, which is the complexity of the calculation with brute force
363 ($W = \sum_{z_1,...,z_n=1}^{Q} [\prod_{1<i<j\le n} \psi_{ij}(z_i, z_j)]$).

364
365 So, the advantage of separation of variables is not in lowering the complexity of the
366 calculation of W or marginals. It is in the possibility to write matrices in TT format, because
367 each matrix $A_k(z)$ is of TT-rank 1 as a Kronecker product of matrices. Additions and
368 multiplications of matrices can be performed in this format, without necessarily making
369 the coefficients of the matrices explicit. This TT matrix approach makes the calculations
370 accessible at the price of some low rank approximations, and is presented in the next section.

371
6. From exact algebraic TT calculation to a numerically tractable code

372
373 We have seen in section 5.3 that the dimensions of the matrices involved in matrices
374 products in the algebraic TT-based expression of the normalising constant or some marginals
375 (like $W = \prod_i B_i$ with $B_i = \sum_{z_i} A_i(z_i)$) have dimensions which can reach values close to
376 $Q^{(n/2)^2} \times Q^{(n/2)^2}$. It is not possible to perform these calculations using conventional matrix
377 calculation algorithms, because of memory space. Novikov et al. (2014) proposed to use a
378 TT-format specific for matrix calculation which has been defined in Oseledets, Tyrtyshnikov
379 & Zamarashkin (2011) and Oseledets & Dolgov (2012), and is referred to as *matrix in TT-*
380 format (Oseledets & Dolgov (2012), equation (2.2)), and enables storage saving.

381 In order to avoid confusion with the TT-format of a matrix read as a tensor of order
382 two, we call here this format *TT-matrix format*. So the TT-format of a matrix is its SVD
383 (see Section 5.1), while a TT-matrix or a matrix expressed in TT-matrix format is the
384 decomposition of a matrix as proposed in Oseledets & Dolgov (2012), and presented in
385 section 6.1.

385 Calculation can be done exactly in TT-matrix format, because matrices $A_i(z_i)$ are
386 expressed as TT-matrices of rank one, hence matrices $B_i$ are expressed as TT-matrices too, as
387 well as their products. Definition of TT-matrices is given in Section 6.1, after the presentation
388 of multiplexing / demultiplexing of indices. However, numerical difficulties remain when
389 using TT-matrices, and we present in Sections 6.2 and 6.3 how we propose to bypass them.
390 The solutions proposed have led to the possibility of computing them on data sets with
391 $n \le 50$ for the structure called 'assortative easy', sometimes less for other structures, as for
392 the structure 'assortative difficult'.
$M_1(z_1, z_2)$ and $M_n(z_n, t_n)$), and there are $nQ^2$ of them. Hence, storage of $A$ in TT-matrix format (i.e. storing the matrices $M_i(z_i, t_i)$ instead of storing $A$) requires $nQ^2r^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 TT-matrix 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 $\text{TT\_rank}(A+B) \le \text{TT\_rank} A + \text{TT\_rank} B$ and $\text{TT\_rank}(AB) \le (\text{TT\_rank} A)(\text{TT\_rank} B)$.

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

So, all calculations of $B_k = \sum_{z_k} A_k(z_k)$, $W = \prod_k B_k$ and marginals can be computed in TT-matrix format. We have $\text{TT\_rank} B \le Q$, and $\text{TT\_rank} W \le Q^m$. Of course, $W$ is a scalar, and $\text{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 $\prod_i B_i$, algorithm 1 in Novikov et al. (2014) 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 = \prod_{i=1}^n B_i$ is the sum of $Q$ matrices $A_i(q)$ for $q \in \{1, \dots, 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 $\epsilon^{TT} = 5 \cdot 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.

## 6.2 Numerical difficulties encountered and options for overcoming them

1. for $n \ge 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 \ge 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 \ge 45$, the number of cores became an issue. Indeed, each matrix $B_i$ has $m$ cores. The rounding requires $\sim 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).
## 6.1. TT-matrices to reduce memory space

**Multiplexing / demultiplexing of indices:** Let $z = (z_1, \dots, z_n) \in \{1, \dots, Q\}^n$ be a multi-index. An index $a$ in $\{1, \dots, Q^n\}$ can be associated to $z$ using the following construction: $a = \sum_{i=1}^{n-1} Q^{n-i}(z_i - 1) + z_n$, which is a bijection between $\{1, \dots, Q\}^n$ and $\{1, \dots, 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 \times Q^n$ matrix. TT-matrix format of $A$ at rank $r$ is a decomposition of the coefficients of $A$ relying on the multiplexed indices as

$$A[a, b] = A[(z_1, \dots, z_n); (t_1, \dots, t_n)] = M_1(z_1, t_1) \dots M_n(z_n, t_n), \quad (9)$$

where $M_i(z_i, t_i)$ is a $r \times r$ matrix (the variables $z_i, t_i$ are not the coefficients of the matrix). $M_1(z_1, t_1)$ is $1 \times r$ and $M_n(z_n, t_n)$ is $r \times 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 TT-matrix format (i.e. storing the matrices $M_i(z_i, t_i)$ instead of storing $A$) requires $nQ^2r^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 TT-matrix 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 II), knowing the cores of $A$ and $B$ only. It is possible to show that $\text{TT\_rank}(A+B) \le \text{TT\_rank} A + \text{TT\_rank} B$ and $\text{TT\_rank}(AB) \le (\text{TT\_rank} A)(\text{TT\_rank} B)$.

## 6.2. Controlling the rank while computing the binary marginals

So, all calculations of $A_k(z_k)$ or $B_k = \sum_{z_k=1}^Q A_k(z_k)$ necessary to compute the binary marginals can be computed in TT-matrix format. We have $\text{TT\_rank} B_k \le Q$ for all $k$ and $\text{TT\_rank} A_k(z_k) = 1$. The computation of a binary marginal is a product of $n-2$ matrices $B_k$ and two matrices $A_k(z_k)$. Hence, its rank is bounded by $Q^{n-2}$. For controlling the exponential growth of the TT-rank when computing the products like $\prod_k B_k$, algorithm 1 in Novikov et al. (2014) consists in approximating the product of two TT-matrices by a matrix of low TT-rank, with an operation called the 'rounding', which can also be performed in TT-matrix format.
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 \sim \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 \ge 24$. This can be solved easily by scaling the WSBM factors with a factor denoted $\alpha$. Let us denote with an exponent $s$ the scaled factors: $\psi_{i,j}^s(z_i, z_j) = \alpha \psi_{i,j}(z_i, z_j)$. We have $\mathbb{P}^s(z, D) = \alpha^m \mathbb{P}(z, D)$ because it is a product of $m$ factors. So, $W^s = \sum_z \mathbb{P}^s(z, D) = \alpha^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) = \sum_{z \setminus z_i, z_j} \mathbb{P}(z, D)$, so $m_{i,j}^s(z_i, z_j) = \alpha^m m_{i,j}(z_i, z_j)$. So, normalised marginals are unchanged by scaling, because $m_{i,j}^s(z_i, z_j)/W^s = \alpha^m m_{i,j}(z_i, z_j)/\alpha^m W = m_{i,j}(z_i, z_j)/W$. We have chosen $\alpha = 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 TT-rank 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 TT-rank. For this experiment, we considered a WSBM with Poisson distributed distances. It means that $\mathbb{P}_{\Lambda}(D[i, j] | Z_i = q, Z_j = q')$ is a Poisson distribution with parameter $\lambda_{qq'}$, and the connectivity matrix $\Lambda$ is defined by $\Lambda[q, q'] = \lambda_{qq'}$. We used

$$ \Lambda = \begin{pmatrix} 2 & 10 & 10 \\ 10 & 3 & 8 \\ 10 & 10 & 4 \end{pmatrix} $$

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 $\epsilon^{TT} = 10^{-2}$, but limiting the rank growth to $r_{max}^{TT} = 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 \le 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: $M_1 \otimes M_2 \otimes M_3 \otimes M_4 = (M_1 \otimes M_2) \otimes (M_3 \otimes M_4)$: we have two larger cores as far as dimensions are concerned instead of four
However, as in SVD of a matrix, there is a trade-off to tune between low rank and high accuracy. 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, we carried out a numerical experiment to compare the quality of the calculation of W with prescribed accuracy and several choices of prescribed TT-rank. For this experiment, we considered a WSBM with Poisson distributed distances. It means that $Pr_{\Lambda}(D[i, j] | Z_i = q, Z_j = q')$ is a Poisson distribution with parameter $\lambda_{qq'}$, and the connectivity matrix $\Lambda$ is defined by $\Lambda[q, q'] = \lambda_{qq'}$. We used

$$ \Lambda = \begin{pmatrix} 2 & 10 & 10 \\ 10 & 3 & 8 \\ 10 & 10 & 4 \end{pmatrix} $$

In the TT rounding algorithm we have combined the two levers of prescribed accuracy and prescribed rank, by setting the accuracy to $\epsilon^{TT} = 5.10^{-2}$, but limiting the rank growth to $r_{max}^{TT}$ for different values of $r_{max}^{TT}$. If the rank reaches the ceiling $r_{max}^{TT}$ along the course of the calculation, it means that the desired accuracy cannot be fulfilled. The results are given in table 2 for different values of $n$ (from $n = 8$ to $n = 32$).

<table><caption>Table 2. Normalising constant W calculated according to different prescribed accuracy and maximal rank settings in rounding. In all experiments, $\epsilon^{TT} = 5.10^{-2}$, then R1: no $r_{max}^{TT}$; R2: $r_{max}^{TT} = 3$; R3: $r_{max}^{TT} = 9$; R4: $r_{max}^{TT} = 27$; R5: $r_{max}^{TT} = 81$. For R1, the maximum of observed rank is given in brackets. Calculations have been done for a WSBM with Poisson distributed distances and $Q = 3$.</caption><thead><tr><th>n</th><th>R<sub>1</sub></th><th>R<sub>2</sub></th><th>R<sub>3</sub></th><th>R<sub>4</sub></th><th>R<sub>5</sub></th></tr></thead><tbody><tr><td>8</td><td>6.21 × 10<sup>-28</sup> (9)</td><td>2.68 × 10<sup>-28</sup></td><td>6.21 × 10<sup>-28</sup></td><td>6.21 × 10<sup>-28</sup></td><td>6.21 × 10<sup>-28</sup></td></tr><tr><td>10</td><td>3.2 × 10<sup>-45</sup> (12)</td><td>2.24 × 10<sup>-46</sup></td><td>3.20 × 10<sup>-45</sup></td><td>3.2 × 10<sup>-45</sup></td><td>3.2 × 10<sup>-45</sup></td></tr><tr><td>12</td><td>1.0 × 10<sup>-75</sup> (31)</td><td>3.32 × 10<sup>-77</sup></td><td>6.73 × 10<sup>-76</sup></td><td>1.0 × 10<sup>-75</sup></td><td>1.0 × 10<sup>-75</sup></td></tr><tr><td>14</td><td>5.38 × 10<sup>-101</sup> (8)</td><td>1.11 × 10<sup>-101</sup></td><td>5.38 × 10<sup>-101</sup></td><td>5.38 × 10<sup>-101</sup></td><td>5.38 × 10<sup>-101</sup></td></tr><tr><td>16</td><td>3.30 × 10<sup>-135</sup> (10)</td><td>3.05 × 10<sup>-135</sup></td><td>3.30 × 10<sup>-135</sup></td><td>3.30 × 10<sup>-135</sup></td><td>3.30 × 10<sup>-135</sup></td></tr><tr><td>18</td><td>2.05 × 10<sup>-161</sup> (15)</td><td>1.14 × 10<sup>-161</sup></td><td>2.05 × 10<sup>-161</sup></td><td>2.05 × 10<sup>-161</sup></td><td>2.05 × 10<sup>-161</sup></td></tr><tr><td>20</td><td>1.18 × 10<sup>-199</sup> (18)</td><td>2.65 × 10<sup>-205</sup></td><td>1.18 × 10<sup>-199</sup></td><td>1.18 × 10<sup>-199</sup></td><td>1.18 × 10<sup>-199</sup></td></tr><tr><td>22</td><td>1.59 × 10<sup>-247</sup> (15)</td><td>1.80 × 10<sup>-247</sup></td><td>1.59 × 10<sup>-247</sup></td><td>1.59 × 10<sup>-247</sup></td><td>1.59 × 10<sup>-247</sup></td></tr><tr><td>24</td><td>9.57 × 10<sup>-294</sup> (16)</td><td>1.12 × 10<sup>-297</sup></td><td>9.57 × 10<sup>-294</sup></td><td>9.57 × 10<sup>-294</sup></td><td>9.57 × 10<sup>-294</sup></td></tr><tr><td>26</td><td>1.80 × 10<sup>-353</sup> (16)</td><td>8.13 × 10<sup>-375</sup></td><td>1.80 × 10<sup>-353</sup></td><td>1.80 × 10<sup>-353</sup></td><td>1.80 × 10<sup>-353</sup></td></tr><tr><td>28</td><td>1.34 × 10<sup>-393</sup> (21)</td><td>1.09 × 10<sup>-393</sup></td><td>1.34 × 10<sup>-393</sup></td><td>1.34 × 10<sup>-393</sup></td><td>1.34 × 10<sup>-393</sup></td></tr><tr><td>30</td><td>4.03 × 10<sup>-459</sup> (21)</td><td>4.03 × 10<sup>-459</sup></td><td>4.03 × 10<sup>-459</sup></td><td>4.03 × 10<sup>-459</sup></td><td>4.03 × 10<sup>-459</sup></td></tr><tr><td>32</td><td>3.43 × 10<sup>-511</sup> (34)</td><td>3.19 × 10<sup>-511</sup></td><td>3.43 × 10<sup>-511</sup></td><td>3.43 × 10<sup>-511</sup></td><td>3.43 × 10<sup>-511</sup></td></tr></tbody></table>

## 6.3. Merging of cores

When $n$ is large, the number of cores becomes an issue. Indeed, each matrix $B_i$ has $m$ cores. The rounding requires about $m^2$ factorisations LQ and SVD which becomes prohibitive when $m$ is large. The TT-toolbox we used to ran our numerical experiments is limited to matrices with 1024 cores or less. As the number $m$ of cores is $n(n-1)/2$, this
<table><caption>**Table 2** Normalising constant calculated according to different accuracy and maximal rank settings in rounding: R1: $\epsilon^{TT} = 5.10^{-2}$; R2: $r_{max}^{TT} = 3$; R3: $r_{max}^{TT} = 9$; R4: $r_{max}^{TT} = 27$; R5: $r_{max}^{TT} = 81$. For $\epsilon^{TT} = 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$.</caption><thead><tr><th>n</th><th>R<sub>1</sub></th><th>R<sub>2</sub></th><th>R<sub>3</sub></th><th>R<sub>4</sub></th><th>R<sub>5</sub></th></tr></thead><tbody><tr><td>8</td><td>6.21 × 10<sup>-28</sup> (9)</td><td>2.68 × 10<sup>-28</sup></td><td>6.21 × 10<sup>-28</sup></td><td>6.21 × 10<sup>-28</sup></td><td>6.21 × 10<sup>-28</sup></td></tr><tr><td>10</td><td>3.2 × 10<sup>-45</sup>(12)</td><td>2.24 × 10<sup>-46</sup></td><td>3.20 × 10<sup>-45</sup></td><td>3.2 × 10<sup>-45</sup></td><td>3.2 × 10<sup>-45</sup></td></tr><tr><td>12</td><td>1.0 × 10<sup>-75</sup>(31)</td><td>3.32 × 10<sup>-77</sup></td><td>6.73 × 10<sup>-76</sup></td><td>1.0 × 10<sup>-75</sup></td><td>1.0 × 10<sup>-75</sup></td></tr><tr><td>14</td><td>5.38 × 10<sup>-101</sup>(8)</td><td>1.11 × 10<sup>-101</sup></td><td>5.38 × 10<sup>-101</sup></td><td>5.38 × 10<sup>-101</sup></td><td>5.38 × 10<sup>-101</sup></td></tr><tr><td>16</td><td>3.30 × 10<sup>-135</sup>(10)</td><td>3.05 × 10<sup>-135</sup></td><td>3.30 × 10<sup>-135</sup></td><td>3.30 × 10<sup>-135</sup></td><td>3.30 × 10<sup>-135</sup></td></tr><tr><td>18</td><td>2.05 × 10<sup>-161</sup>(15)</td><td>1.14 × 10<sup>-161</sup></td><td>2.05 × 10<sup>-161</sup></td><td>2.05 × 10<sup>-161</sup></td><td>2.05 × 10<sup>-161</sup></td></tr><tr><td>20</td><td>1.18 × 10<sup>-199</sup>(18)</td><td>2.65 × 10<sup>-205</sup></td><td>1.18 × 10<sup>-199</sup></td><td>1.18 × 10<sup>-199</sup></td><td>1.18 × 10<sup>-199</sup></td></tr><tr><td>22</td><td>1.59 × 10<sup>-247</sup>(15)</td><td>1.80 × 10<sup>-247</sup></td><td>1.59 × 10<sup>-247</sup></td><td>1.59 × 10<sup>-247</sup></td><td>1.59 × 10<sup>-247</sup></td></tr><tr><td>24</td><td>9.57 × 10<sup>-294</sup>(16)</td><td>1.12 × 10<sup>-297</sup></td><td>9.57 × 10<sup>-294</sup></td><td>9.57 × 10<sup>-294</sup></td><td>9.57 × 10<sup>-294</sup></td></tr><tr><td>26</td><td>1.80 × 10<sup>-353</sup>(16)</td><td>8.13 × 10<sup>-375</sup></td><td>1.80 × 10<sup>-353</sup></td><td>1.80 × 10<sup>-353</sup></td><td>1.80 × 10<sup>-353</sup></td></tr><tr><td>28</td><td>1.34 × 10<sup>-393</sup>(21)</td><td>1.09 × 10<sup>-393</sup></td><td>1.34 × 10<sup>-393</sup></td><td>1.34 × 10<sup>-393</sup></td><td>1.34 × 10<sup>-393</sup></td></tr><tr><td>30</td><td>4.03 × 10<sup>-459</sup>(21)</td><td>4.03 × 10<sup>-459</sup></td><td>4.03 × 10<sup>-459</sup></td><td>4.03 × 10<sup>-459</sup></td><td>4.03 × 10<sup>-459</sup></td></tr><tr><td>32</td><td>3.43 × 10<sup>-511</sup>(34)</td><td>3.19 × 10<sup>-511</sup></td><td>3.43 × 10<sup>-511</sup></td><td>3.43 × 10<sup>-511</sup></td><td>3.43 × 10<sup>-511</sup></td></tr></tbody></table>

smaller cores. This has been implemented for each matrix $A_k(z_k)$. Let us recall that (see formula (C15)): $A_k(z_k) = \bigotimes_{1 \le i < j \le n} G_k^{ij}(z_k)$ The fusion here is carried out at the level of the matrices $G_k^{ij}[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.

<table><caption>**Table 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$.</caption><thead><tr><th>n</th><th>Without merging</th><th>Merging with uniform numbers</th><th>Merging with uniform volumes</th></tr></thead><tbody><tr><td>12</td><td>9.35 × 10<sup>-71</sup></td><td>9.35 × 10<sup>-71</sup></td><td>9.35 × 10<sup>-71</sup></td></tr><tr><td>15</td><td>8.04 × 10<sup>-112</sup></td><td>7.95 × 10<sup>-112</sup></td><td>7.95 × 10<sup>-112</sup></td></tr><tr><td>18</td><td>3.77 × 10<sup>-163</sup></td><td>3.77 × 10<sup>-163</sup></td><td>3.77 × 10<sup>-163</sup></td></tr><tr><td>21</td><td>2.92 × 10<sup>-219</sup></td><td>2.92 × 10<sup>-219</sup></td><td>2.92 × 10<sup>-219</sup></td></tr><tr><td>24</td><td>1.53 × 10<sup>-281</sup></td><td>1.53 × 10<sup>-281</sup></td><td>1.53 × 10<sup>-281</sup></td></tr><tr><td>27</td><td>3.49 × 10<sup>-359</sup></td><td>3.76 × 10<sup>-359</sup></td><td>3.73 × 10<sup>-359</sup></td></tr><tr><td>30</td><td>4.40 × 10<sup>-459</sup></td><td>4.40 × 10<sup>-459</sup></td><td>4.40 × 10<sup>-459</sup></td></tr><tr><td>33</td><td>7.81 × 10<sup>-565</sup></td><td>7.81 × 10<sup>-565</sup></td><td>7.81 × 10<sup>-565</sup></td></tr><tr><td>36</td><td>1.55 × 10<sup>-673</sup></td><td>1.54 × 10<sup>-673</sup></td><td>1.55 × 10<sup>-673</sup></td></tr><tr><td>39</td><td>1.08 × 10<sup>-775</sup></td><td>1.08 × 10<sup>-775</sup></td><td>1.08 × 10<sup>-775</sup></td></tr><tr><td>42</td><td>1.06 × 10<sup>-915</sup></td><td>1.06 × 10<sup>-915</sup></td><td>1.06 × 10<sup>-915</sup></td></tr></tbody></table>
corresponds to $n = 45$ nodes or less. Therefore, we set up core fusion to work with fewer but larger cores. 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: $M_1 \otimes M_2 \otimes M_3 \otimes M_4 = (M_1 \otimes M_2) \otimes (M_3 \otimes M_4)$: in the r.h.s. we have two cores with dimensions larger than that of the four cores in the l.h.s.. This has been implemented for each matrix $A_k(z_k)$. Let us recall that (see formula (7)): $A_k(z_k) = \bigotimes_{1 \le i < j \le n} G_k^{ij}(z_k)$. The fusion here is carried out at the level of the matrices $G_k^{(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 solution 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 since the method with uniform numbers is much simpler to implement, we have chosen to merge cores with uniform numbers.

<table><caption>Table 3. Computed values of the normalising constant $W$ 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$.</caption><thead><tr><th>n</th><th>Without merging</th><th>Merging with uniform numbers</th><th>Merging with uniform volumes</th></tr></thead><tbody><tr><td>12</td><td>9.35 × 10<sup>-71</sup></td><td>9.35 × 10<sup>-71</sup></td><td>9.35 × 10<sup>-71</sup></td></tr><tr><td>15</td><td>8.04 × 10<sup>-112</sup></td><td>7.95 × 10<sup>-112</sup></td><td>7.95 × 10<sup>-112</sup></td></tr><tr><td>18</td><td>3.77 × 10<sup>-163</sup></td><td>3.77 × 10<sup>-163</sup></td><td>3.77 × 10<sup>-163</sup></td></tr><tr><td>21</td><td>2.92 × 10<sup>-219</sup></td><td>2.92 × 10<sup>-219</sup></td><td>2.92 × 10<sup>-219</sup></td></tr><tr><td>24</td><td>1.53 × 10<sup>-281</sup></td><td>1.53 × 10<sup>-281</sup></td><td>1.53 × 10<sup>-281</sup></td></tr><tr><td>27</td><td>3.49 × 10<sup>-359</sup></td><td>3.76 × 10<sup>-359</sup></td><td>3.73 × 10<sup>-359</sup></td></tr><tr><td>30</td><td>4.40 × 10<sup>-459</sup></td><td>4.40 × 10<sup>-459</sup></td><td>4.40 × 10<sup>-459</sup></td></tr><tr><td>33</td><td>7.81 × 10<sup>-565</sup></td><td>7.81 × 10<sup>-565</sup></td><td>7.81 × 10<sup>-565</sup></td></tr><tr><td>36</td><td>1.55 × 10<sup>-673</sup></td><td>1.54 × 10<sup>-673</sup></td><td>1.55 × 10<sup>-673</sup></td></tr><tr><td>39</td><td>1.08 × 10<sup>-775</sup></td><td>1.08 × 10<sup>-775</sup></td><td>1.08 × 10<sup>-775</sup></td></tr><tr><td>42</td><td>1.06 × 10<sup>-915</sup></td><td>1.06 × 10<sup>-915</sup></td><td>1.06 × 10<sup>-915</sup></td></tr></tbody></table>

## 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 (combining the algorithm of Section 5.2, rounding and merging, and 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 exact enumeration, the Gibbs Sampler and the Mean-Field approximation. For these experiments, we considered a WSBM with Poisson distributed
### 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 (referred 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.

## 7.1 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 (Robert and Casella, 2004). 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 Daudin et al. (2008). 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, Geman and Geman, 1984) is a Markov Chain Monte-Carlo method that exploits conditional probabilities associated to the complex probability of interest. In our case, GS consists in simulating iteratively a Markov chain, using only $P_θ(Z_i = q | \{Z_j = z_j\}_{j \neq 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_{burnin}^{GS}$ iterations, then we sample a simulated value of Z every $N_{thinning}^{GS}$ 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.
distances. In order to study the influence of the connectivity matrix $\Lambda$ on the quality of the
results, we perform the comparison on six constrained WSBM structures.

**7.1. Different methods for computing the marginals**

The first method tested for computing the binary marginals (4) is the exact enumeration
of all the terms in the sum. There are $Q^{n-2}$ terms for computing a single binary marginal.
So this method is only available for small $n$. Beyond $n = 12$ and $Q = 3$ it is not possible to
store the tensor corresponding to $\text{Pr}_{\theta}(Z, D)$ in memory. To circumvent this complexity two
strategies are classically used. The exact computation can be approximated using simulations
of the model, for instance with the Gibbs Sampler (Robert & Casella 2004). 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 Daudin, Picard & Robin
(2008). 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 has no theoretical guarantee. Let us detail these
two methods.

**Gibbs Sampler.** The Gibbs Sampler (GS, Geman & Geman 1984) is a Markov Chain Monte-Carlo method that exploits conditional probabilities associated to the complex probability of interest. In our case, GS consists in simulating iteratively a Markov chain, using only $\text{Pr}_{\theta}(Z_i = q|\{Z_j = z_j\}_{j \neq i}, D)$, and whose stationary distribution converges towards $\text{Pr}_{\theta}(Z|D)$. In practice, we simulate a single GS run (or path) for a number $N_{burnin}^{GS}$ iterations, then we sample a simulated value of $Z$ every $N_{thinning}^{GS}$ new iterations until we collect $L$ realisations. Binary marginals of $\text{Pr}_{\theta}(Z|D)$ are approximated by empirical frequencies.

**Mean-Field.** With the Mean-Field approximation (MF), the conditional distribution $\text{Pr}_{\theta}(Z|D)$ is approximated by a distribution $\mathbb{Q}_{\theta}(Z|D)$ for which marginalisation is less costly. $\mathbb{Q}_{\theta}(Z|D)$ is chosen among the family $\mathcal{Q}$ of distributions that satisfy the hypothesis of mutual independence between the $Z_i$s. It means that $\mathbb{Q}_{\theta}(Z|D) = \prod_{i=1}^{n} q_{\theta}^{i}(Z_i|D)$, where the $q_{\theta}^{i}(Z_i|D)$ are the unary marginals of $\mathbb{Q}_{\theta}$. The binary marginals are easily obtained as $\mathbb{Q}_{\theta}(Z_i = q, Z_j = q'|D) = q_{\theta}^{i}(Z_i = q|D)q_{\theta}^{j}(Z_j = q'|D)$. However, this expression ignores the dependencies between the $Z_i$s that exist in the true model.
$\mathbb{Q}_{\theta}(Z | D)$ is chosen among the family $\mathbb{Q}$ of distributions that satisfy the hypothesis of mutual independence between the $Z_i$s. It means that $\mathbb{Q}_{\theta}(Z | D) = \prod_{i=1}^{n} q_{\theta}^{i}(Z_i | D)$, where the $q_{\theta}^{i}(Z_i | D)$ are the unary marginals of $\mathbb{Q}_{\theta}$. The binary marginals are easily obtained as $\mathbb{Q}_{\theta}(Z_i = q, Z_j = q' | D) = q_{\theta}^{i}(Z_i = q | D)q_{\theta}^{j}(Z_j = q' | D)$.

The distribution $\mathbb{Q}_{\theta}(Z | D)$ is chosen as the one that minimises the Kullback-Leibler divergence with the true distribution $\mathbb{P}_{\theta}(Z | D)$:

$$ \mathbb{Q}_{\theta}(\cdot | D) = \arg \max_{\mathbb{Q} \in \mathbb{Q}} KL(\mathbb{Q}(\cdot) | \mathbb{P}_{\theta}(\cdot | D)) $$

with $KL(\mathbb{Q}(\cdot)|\mathbb{P}(\cdot)) = \sum_{z} \mathbb{Q}(z) \ln(\frac{\mathbb{Q}(z)}{\mathbb{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 $\tau_{i,q}^{t}$ the current value of $\mathbb{Q}(Z_i = q | D)$, we have

$$ \forall 1 \le i \le n, \forall 1 \le q \le Q, \quad \left\{ \begin{array}{l} \tau_{i,q}^{t+1} = \frac{\mu_{i,q}^{t+1/2}}{\sum_{l=1}^{Q} \mu_{i,l}^{t+1/2}} \\ \mu_{i,q}^{t+1/2} = \alpha_q \prod_{j \ne i} \prod_{q'=1}^{Q} \left[ \frac{\lambda_{q,q'}^{D(i,j)}}{D(i,j)!} \exp(-\lambda_{q,q'}) \right]^{\tau_{j,q'}^{t}} \end{array} \right. $$

The scheme is initialised by $\tau_{i,q}^{0} = \alpha_q$. We stop the scheme when $\|\Delta\|_{\infty} < \epsilon^{MF}$ with $\Delta = \{\Delta_{i,q}, 1 \le i \le n, 1 \le q \le Q\}$ and

$$ \Delta_{i,q} = \frac{|\tau_{i,q}^{t+1} - \tau_{i,q}^{t}|}{|\tau_{i,q}^{t}|} $$

## 7.2 The 6 WSBM structures

In the literature, different structures of the connectivity matrix $\Lambda$ has been considered (Kalmbach et al., 2017; Funke and Becker, 2019). 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 $\lambda_{qq}$ are small and inter group dissimilarities $\lambda_{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 $\Lambda$ 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.
The distribution $Q_\theta(Z | D)$ is chosen as the one that minimises the Kullback-Leibler divergence with the true distribution $\text{Pr}_\theta(Z | D)$:

$$Q_\theta(\cdot | D) = \arg \max_{Q \in \mathcal{Q}} KL(Q(\cdot) | \text{Pr}_\theta(\cdot | D))$$

with $KL(Q(\cdot)|\text{Pr}(\cdot)) = \sum_z Q(z) \ln \left( \frac{Q(z)}{\text{Pr}(z)} \right)$. The solution is on the form of a fixed point equation which is solved by an iterative scheme (see Appendix III for the demonstration). If we denote by $\tau_{i,q}^t$ the value of $Q(Z_i = q | D)$ at time $t$, we have

$$\text{for all } 1 \le i \le n, \text{ for all } 1 \le q \le Q, \quad \left\{ \begin{array}{ll} \tau_{i,q}^{t+1} = \frac{\mu_{i,q}^{t+1/2}}{\sum_{l=1}^{Q} \mu_{i,l}^{t+1/2}} & \text{with} \\ \mu_{i,q}^{t+1/2} = \alpha_q \prod_{j \ne i} \prod_{q'=1}^{Q} \left[ \frac{\lambda_{q,q'}^{D(i,j)}}{D(i,j)!} \exp(-\lambda_{q,q'}) \right]^{\tau_{j,q'}^t} & \end{array} \right. $$

The scheme is initialised by $\tau_{i,q}^0 = \alpha_q$. We stop the scheme as soon as there is a value of $t$ for which $\|\Delta_t\|_\infty < \epsilon^{MF}$ with $\Delta_t = \{\Delta_{i,q}^t, 1 \le i \le n, 1 \le q \le Q\}$ and

$$\Delta_{i,q}^{t} = \frac{|\tau_{i,q}^{t+1} - \tau_{i,q}^{t}|}{|\tau_{i,q}^{t}|}.$$

## 7.2. The 6 WSBM structures

In the literature, different structures of the connectivity matrix $\Lambda$ have been considered (Kalmbach et al. 2017; Funke & Becker 2019). They correspond to different patterns 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 : mean intra group dissimilarities $\lambda_{qq}$ are small and mean inter group dissimilarities $\lambda_{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 mean 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 mean inter dissimilarities are large and the mean intra dissimilarities range from small to large. For each structure, the $\Lambda$ 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.
<figure><img src="image_3.png" /><figcaption>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.</figcaption></figure>

7.3 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 gen-
erated 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
<figure><img src="image_5.png" /><figcaption>Figure 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.</figcaption></figure>

**7.3. Protocol**

For *n* = 12, exact enumeration 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 marginal 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,
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 $\alpha_q = 1/3$ in all the experiments. Then, for a given number of individuals $n$ and a given connectivity matrix $\Lambda$, we generated a realisation $z$ of $\mathbb{P}_\alpha(Z)$ and then we generated a dissimilarity matrix $D$ according to $\mathbb{P}_\Lambda(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 $\epsilon^{TT}$ and the maximal TT-rank $r_{max}^{TT}$. We set $\epsilon^{TT} = 10^{-2}$ and $r_{max}^{TT} = 27$ for all experiments except $n = 40$ and the assortative easy structure where we lowered $r_{max}^{TT}$ 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 \times 10^5$, the warming phase $N_{burnin}^{GS} = 10^3$ and the thinning frequency $N_{thinning}^{GS} = 5$. There exist no universal operational rules to set these parameters (Casella and George, 1992). 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: $\epsilon^{MF} = 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 rela- beling 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 $PSC_{i,j} = \sum_{q=1}^Q \mathbb{P}_\theta(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 $\mathbb{P}_\theta(Z_i = q, Z_j = q' | D)$.

## 7.4 Results

### 7.4.1 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 \times 3$ values, so we were able to compare 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 $\Lambda$ 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
512 meaning that two groups share almost the same pattern of connections. We also considered
513 the core-periphery structure because for $n = 12$ we observed that the associated unary
514 marginals have a different pattern than for the two assortative structures (see Section 7.4.1
515 and Table 4). Therefore we can expect that the 3 selected structures correspond to different
516 difficulties regarding marginal inference.

517
518 **Generation of the dissimilarity matrices.** We set $Q = 3$ and $\alpha_q = 1/3$ in all the
519 experiments. Then, for a given number of individuals $n$ and a given connectivity matrix
520 $\Lambda$, we generated a realisation $z$ of $\text{Pr}_\alpha(Z)$ and then we generated a dissimilarity matrix $D$
521 according to $\text{Pr}_\Lambda(D | Z = z)$.

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

525
526
527 - *exact enumeration*: with this method, the tensor of the joint distribution is stored in
memory, and the binary and unary marginals are computed by using expression (4). It
does not require any parameter setting;

528
529
530
531 - **TT**: parameters are the accuracy $\epsilon^{TT}$ and the maximal TT-rank $r_{max}^{TT}$. We set $\epsilon^{TT} = 10^{-2}$ and $r_{max}^{TT} = 27$ for all experiments except $n = 40$ and the assortative easy structure where we lowered $r_{max}^{TT}$ to 9 since it was enough to reach a good agreement with the other methods.

532
533
534
535
536
537 - *Gibbs Sampler*: parameters are the number of iterations $L^{GS} = 5 \times 10^5$, the warming phase $N_{burnin}^{GS} = 10^3$ and the thinning frequency $N_{thinning}^{GS} = 5$. There exist no universal operational rules to set these parameters (Casella & George 1992). We chose 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;

538
- *Mean-Field*: accuracy for stopping the iteration: $\epsilon^{MF} = 10^{-1}$.

539
Note that these choices may be specific to $Q = 3$.

540
541
**Output.** Computation of binary marginals by Gibbs Sampler can be subject to label
542 switching (Murphy 2012, chapter 24). It means that there should be a phase of relabeling
543 before comparison with the binary marginals obtained with the other methods. To
544 avoid this step which can be tricky, we compared the values of the probability that
545 two individuals *i* and *j* belong to the same class. This probability is obtained as
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 $PSC_{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 $PSC_{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 $PSC_{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 $PSC_{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 $PSC$ while leading to a poor approximation of the unary marginals. This is due to the fact that the computation of $PSC(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 happens 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 \neq q'$.
546

$PSC_{i,j} = \sum_{q=1}^{Q} \Pr_{\theta}(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 $\Pr_{\theta}(Z_i = q, Z_j = q' | D)$.

547

548

549

550

The numerical experiments have been produced with the code available at https://forgemia.inra.fr/anwar_phd/tt-wsbm.

551

552

## 7.4. Results

553

### 7.4.1. Numerical experiments for $n = 12$

554

555

556

557

558

559

560

561

562

563

564

565

566

567

568

569

570

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 compare the exact and approximated values one by one instead of 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 $\Lambda$ 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.

571

572

573

574

575

576

577

Regarding now binary marginals, Table 5 provides, for each $\Lambda$ structure, the median of the absolute difference between the value of $PSC_{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 4 in Appendix IV). The median is computed over all pairs of variables $i, j$. 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 $\Lambda$ structure.

578

579

580

Having in mind that the true value of $PSC_{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.
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 Λ.

### 7.4.2 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 $PSC_{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 $PSC_{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 $PSC$, while lower than for the 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 $\Lambda$, as opposed to GS and MF. Indeed for the assortative difficult structure we had to change $r_{max}^{TT}$ 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_{max}^{TT}$ 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.

## 8 Conclusions and perspectives

In 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 Novikov et al. (2014) developed for approximating normalising constant
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 $\Lambda$ structure (see Table 5). We also observed that the MF approximation of $PSC_{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 V) 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 $PSC_{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 $\Lambda$. For instance, whatever the method, the approximation is of better quality for the assortative easy structure than the assortative difficult structure.

Note that a method can provide a good approximation of $PSC_{i,j}$ while leading to a poor approximation of the unary marginals. This is due to the fact that the computation of $PSC_{i,j}$ does not use all the binary marginals $Pr_{\theta}(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 happens 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 $Pr_{\theta}(Z_i = q, Z_j = q' | D)$ when $q \neq q'$.

Regarding computational time, the four methods rank as expected, with from the faster to the slowest MF, TT, exact enumeration and GS. For TT and MF, the computational time varies with the structure of $\Lambda$.

### 7.4.2. 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 $PSC_{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 by construction. The agreement between MF and the two other methods increases when $n = 40$, probably because, heuristically, Mean-Field
<figure><img src="image_4.png" /><figcaption>**Fig. 3** Biplots of the approximated values of the probability to be in the same group $PSC_{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.</figcaption></figure>

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 $\times$ matrices products.

One of the key issues 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 $\Lambda$ 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
<figure><img src="image_6.png" /><figcaption>Figure 3. Biplots of the approximated values of the probability to be in the same group $PSC_{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.</figcaption></figure>

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 $PSC_{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 methods on the value of $PSC_{i,j}$, while lower than for the 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 $\Lambda$, as opposed to GS and MF. Indeed for the assortative difficult structure we had to change $r_{max}^{TT}$ 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_{max}^{TT}$ 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.
difficulty with which some groups can be distinguished. The structure of $\Lambda$ 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 (Geman and Geman, 1984; Häggström, 2002), 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 (Daudin et al., 2008), 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 (Al Daas et al., 2022) or, even more recently, with randomised algorithms (Al Daas et al., 2023).

## 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.
In this paper, we propose an algebraic approach to compute exactly the marginals (normalising constant, unaries, binaries, ternaries, ...) of a pairwise graphical model. For this purpose, we have adapted a previous method from Novikov et al. (2014), developed for approximating normalising constant 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 of a pairwise graphical model 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 the case of the Weighted Stochastic Block Model (WSBM) where the underlying graph is a clique, therefore corresponding to the most complex situation. The exact algebraic approach is not tractable in this case. We propose to combine specific techniques: working with TT-matrices, using a numerical approximation based on rounding for controlling the TT-rank in matrix × matrices products, and merging cores for controlling dimensions.

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 accuracy (from our experiments we cannot determine if *n* has such an impact 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. More globally, our numerical experiments enabled to explore the variations in accuracy of different inference methods for various connectivity patterns in the WSBM.
## 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.
670

671

672

673

674

675

676

677

678

679

680

681

682

683

We have compared our TT-based approach with 2 classical approximation approaches: Gibbs Sampler (GS) and Mean-Field approximation (MF). GS has a sound statistical background which guarantees convergence to an unbiased estimator of the marginals (Geman & Geman 1984; Häggström 2002), but there are to our knowledge no 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 (Daudin, Picard & Robin 2008), but it is not designed to compute binary marginals (apart from a 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 binary marginals for intermediate values of *n* (below 50).

684

685

686

687

688

689

690

691

692

693

The algebraic exactness of the TT-based approach (before rounding) is a very good basis for developing efficient and accurate numerical approaches for calculating the marginals of pairwise graphical models. In the current version of the implementation (with rounding), TT cannot be ran on problems with *n* = 50, due to numerical limits in the rounding procedure. 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 (Al Daas, Ballard & Benner 2022) or, even more recently, with randomised algorithms (Al Daas et al. 2023).
<figure><img src="image_5.png" /><figcaption>**Table 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.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Table 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 (dark grey) cells correspond to situations where all marginals are Dirac, yellow (light grey) cells to situations where marginals are either Dirac or uniform on two groups, and green (grey) 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 dissontative, 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.</figcaption></figure>
<table><caption>**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.</caption><thead><tr><th>Λ</th><th>TT</th><th>Gibbs</th><th>MF</th></tr></thead><tbody><tr><td>Assortative difficult</td><td>2.126641e-03</td><td>1.137171e-01</td><td>2.096444e-02 (10)</td></tr><tr><td>Assortative easy</td><td>7.487330e-10</td><td>8.200865e-09</td><td>1.177558e-08 (4)</td></tr><tr><td>Core periphery</td><td>1.430457e-03</td><td>1.981476e-04</td><td>2.106457e-02 (8)</td></tr><tr><td>Dissortative</td><td>2.037286e-07</td><td>2.723249e-07</td><td>7.601169e-07 (8)</td></tr><tr><td>Hierarchic (II)</td><td>3.293547e-03</td><td>6.799929e-03</td><td>2.340016e-02 (5)</td></tr><tr><td>Ordered</td><td>3.891483e-08</td><td>2.265697e-06</td><td>- (0)</td></tr><tr><td>Random</td><td>2.067789e-02</td><td>1.037811e-03</td><td>1.925229e-01 (5)</td></tr></tbody></table>

<table><caption>**Table 6** 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.</caption><thead><tr><th>Λ</th><th>Exact</th><th>TT</th><th>Gibbs</th><th>MF</th></tr></thead><tbody><tr><td>Assortative difficult</td><td>0.1748605</td><td>7.366531</td><td>354.3655</td><td>0.009956175 (10)</td></tr><tr><td>Assortative easy</td><td>0.1761076</td><td>3.626953</td><td>352.9239</td><td>0.046608290 (4)</td></tr><tr><td>Core periphery</td><td>0.1770441</td><td>5.178768</td><td>368.2358</td><td>0.062990053 (8)</td></tr><tr><td>Dissortative</td><td>0.1755904</td><td>4.698245</td><td>357.3765</td><td>0.059073172 (8)</td></tr><tr><td>Hierarchic (II)</td><td>0.1756495</td><td>7.124235</td><td>364.6317</td><td>0.021913800 (5)</td></tr><tr><td>Ordered</td><td>0.1800913</td><td>3.962546</td><td>359.0718</td><td>- (0)</td></tr><tr><td>Random</td><td>0.1784855</td><td>6.781812</td><td>361.7040</td><td>0.020079037 (5)</td></tr></tbody></table>

<table><caption>**Table 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*<sub>max</sub><sup>TT</sup> = 9 and those on the assortative difficult and core-periphery structures were run with *r*<sub>max</sub><sup>TT</sup> = 27.</caption><thead><tr><th rowspan="2">Λ</th><th colspan="3">n = 25</th><th colspan="3">n = 40</th></tr><tr><th>TT vs GS</th><th>TT vs MF</th><th>GS vs MF</th><th>TT vs GS</th><th>TT vs MF</th><th>GS vs TT</th></tr></thead><tbody><tr><td>Assort. easy</td><td>0</td><td>3.8e-16</td><td>4.1e-16</td><td>9.07e-61</td><td>5.03e-49</td><td>6.16e-56</td></tr><tr><td>Assort. difficult</td><td>1.02e-01</td><td>2.16e-02</td><td>1.53e-01</td><td>1.15e-01</td><td>1.93e-02</td><td>1.86e-01</td></tr><tr><td>Core periphery</td><td>1.26e-06</td><td>4.91e-07</td><td>1.54e-06</td><td>3.76e-09</td><td>1.88e-11</td><td>5.73e-10</td></tr></tbody></table>

<table><caption>**Table 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*<sub>max</sub><sup>TT</sup> = 9 and those on the assortative difficult and core-periphery structures were run with *r*<sub>max</sub><sup>TT</sup> = 27.</caption><thead><tr><th rowspan="2">Λ</th><th colspan="3">n = 25</th><th colspan="3">n = 40</th></tr><tr><th>TT</th><th>GS</th><th>MF</th><th>TT</th><th>GS</th><th>MF</th></tr></thead><tbody><tr><td>Assort. easy</td><td>61.90</td><td>1536.32</td><td>0.15 (5)</td><td>202.97</td><td>4670.47</td><td>0.42 (2)</td></tr><tr><td>Assort. difficult</td><td>278.48</td><td>1503.69</td><td>0.03 (10)</td><td>2790.38</td><td>4670.25</td><td>0.09 (10)</td></tr><tr><td>Core-periphery</td><td>86.39</td><td>1505.00</td><td>0.13 (6)</td><td>393.33</td><td>4526.13</td><td>0.36 (7)</td></tr></tbody></table>
<table><caption>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.</caption><thead><tr><th>Λ</th><th>TT</th><th>Gibbs</th><th>MF</th></tr></thead><tbody><tr><td>Assortative difficult</td><td>2.126641e-03</td><td>1.137171e-01</td><td>2.096444e-02 (10)</td></tr><tr><td>Assortative easy</td><td>7.487330e-10</td><td>8.200865e-09</td><td>1.177558e-08 (4)</td></tr><tr><td>Core periphery</td><td>1.430457e-03</td><td>1.981476e-04</td><td>2.106457e-02 (8)</td></tr><tr><td>Dissortative</td><td>2.037286e-07</td><td>2.723249e-07</td><td>7.601169e-07 (8)</td></tr><tr><td>Hierarchical</td><td>3.293547e-03</td><td>6.799929e-03</td><td>2.340016e-02 (5)</td></tr><tr><td>Ordered</td><td>3.891483e-08</td><td>2.265697e-06</td><td>- (0)</td></tr><tr><td>Random</td><td>2.067789e-02</td><td>1.037811e-03</td><td>1.925229e-01 (5)</td></tr></tbody></table>

<table><caption>Table 6. 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.</caption><thead><tr><th>Λ</th><th>Exact</th><th>TT</th><th>Gibbs</th><th>MF</th></tr></thead><tbody><tr><td>Assortative difficult</td><td>0.1748605</td><td>7.366531</td><td>354.3655</td><td>0.009956175 (10)</td></tr><tr><td>Assortative easy</td><td>0.1761076</td><td>3.626953</td><td>352.9239</td><td>0.046608290 (4)</td></tr><tr><td>Core periphery</td><td>0.1770441</td><td>5.178768</td><td>368.2358</td><td>0.062990053 (8)</td></tr><tr><td>Dissortative</td><td>0.1755904</td><td>4.698245</td><td>357.3765</td><td>0.059073172 (8)</td></tr><tr><td>Hierarchical</td><td>0.1756495</td><td>7.124235</td><td>364.6317</td><td>0.021913800 (5)</td></tr><tr><td>Ordered</td><td>0.1800913</td><td>3.962546</td><td>359.0718</td><td>- (0)</td></tr><tr><td>Random</td><td>0.1784855</td><td>6.781812</td><td>361.7040</td><td>0.020079037 (5)</td></tr></tbody></table>

<table><caption>Table 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_{max}^{TT} = 9$ and those on the assortative difficult and core-periphery structures were run with $r_{max}^{TT} = 27$.</caption><thead><tr><th rowspan="2">Λ</th><th colspan="3">n = 25</th><th colspan="3">n = 40</th></tr><tr><th>TT vs GS</th><th>TT vs MF</th><th>GS vs MF</th><th>TT vs GS</th><th>TT vs MF</th><th>GS vs TT</th></tr></thead><tbody><tr><td>Assort. easy</td><td>0</td><td>3.8e-16</td><td>4.1e-16</td><td>9.07e-61</td><td>5.03e-49</td><td>6.16e-56</td></tr><tr><td>Assort. difficult</td><td>1.02e-01</td><td>2.16e-02</td><td>1.53e-01</td><td>1.15e-01</td><td>1.93e-02</td><td>1.86e-01</td></tr><tr><td>Core periphery</td><td>1.26e-06</td><td>4.91e-07</td><td>1.54e-06</td><td>3.76e-09</td><td>1.88e-11</td><td>5.73e-10</td></tr></tbody></table>

<table><caption>Table 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_{max}^{TT} = 9$ and those on the assortative difficult and core-periphery structures were run with $r_{max}^{TT} = 27$.</caption><thead><tr><th rowspan="2">Λ</th><th colspan="3">n = 25</th><th colspan="3">n = 40</th></tr><tr><th>TT</th><th>GS</th><th>MF</th><th>TT</th><th>GS</th><th>MF</th></tr></thead><tbody><tr><td>Assort. easy</td><td>61.90</td><td>1536.32</td><td>0.15 (5)</td><td>202.97</td><td>4670.47</td><td>0.42 (2)</td></tr><tr><td>Assort. difficult</td><td>278.48</td><td>1503.69</td><td>0.03 (10)</td><td>2790.38</td><td>4670.25</td><td>0.09 (10)</td></tr><tr><td>Core-periphery</td><td>86.39</td><td>1505.00</td><td>0.13 (6)</td><td>393.33</td><td>4526.13</td><td>0.36 (7)</td></tr></tbody></table>
# Appendix A Mean-Field fixed point equation

## A.1 Principle

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

$$KL(Q|P) = E_{Q}\left[\ln\left(\frac{Q(Z)}{P(Z)}\right)\right]$$

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

Let us define the family $Q$ of distributions $Q$ such that $Q(Z) = \prod_{i=1}^{n} q^{i}(Z_i)$, and let us denote $\tau_{i,q} = q^{i}(Z_i = q)$. Then $Q_{\theta}(\cdot | D)$ is solution of

$$Q_{\theta(\cdot|D)} = \arg\min_{Q \in Q} KL(Q | P_{\theta}(\cdot | D))$$

We have

$$KL(Q(\cdot) | P_{\theta}(\cdot | D)) = E_{Q}[\ln(Q)] - E_{Q}[\ln(P_{\theta}(Z | D))]$$

From equations (1) and (2), it is equal to

$$\begin{aligned}
KL(Q(\cdot) | P_{\theta}(\cdot | D)) &= \sum_{i=1}^{n} \sum_{q=1}^{Q} \tau_{i,q} \ln(\tau_{i,q}) - E_{Q} \left[ \ln \left( \frac{1}{W} \left( \prod_{i=1}^{n} \prod_{j>i} P_{\Lambda}(D[i, j] | Z_i, Z_j) \right) \prod_{i=1}^{n} P_{\alpha}(Z_i) \right) \right] \\
&= \sum_{i=1}^{n} \sum_{q=1}^{Q} \tau_{i,q} \ln(\tau_{i,q}) - E_{Q} \left[ \ln \left( \frac{1}{W} \times \prod_{i=1}^{n} \left[ \prod_{j>i} P_{\Lambda}(D[i, j] | Z_i, Z_j) \prod_{q=1}^{Q} \alpha_{q}^{Z_{i,q}} \right] \right) \right] \\
&= \sum_{i=1}^{n} \sum_{q=1}^{Q} \tau_{i,q} \ln(\tau_{i,q}) - \sum_{i=1}^{n} \sum_{j>i} \sum_{q=1}^{Q} \sum_{q'=1}^{Q} \tau_{i,q} \tau_{j,q'} \ln P_{\Lambda}(D[i, j] | Z_i = q, Z_j = q') \\
&\quad - \sum_{i=1}^{n} \sum_{q=1}^{Q} (\tau_{i,q} \ln(\alpha_{q})) - \ln(\frac{1}{W})
\end{aligned}$$

The minimisation of $KL(Q(\cdot) | P_{\theta}(\cdot | D))$ is under the constraint that for all $i$, $\sum_{q=1}^{Q} \tau_{i,q} = 1$. So, using Lagrangian multipliers, we minimise

$$G = KL(Q | P(\cdot | D)) + \sum_{i=1}^{n} \sigma_i \left( \sum_{q=1}^{Q} \tau_{i,q} - 1 \right)$$
694

# Appendix I

695

We present here in detail the four steps in Novikov approach to compute the normalising contant of a graphical model. We detail as well how it can be prolonged to compute the unary marginals, as quickly suggested in Novikov et al. (2014).

696

697

Let us consider a graphical model with $n$ variables $z_1, \dots, z_n$, each representing the state of an individual, with $z_i \in \{1, \dots, Q\}$. The graphical model is defined on a hypergraph $G = (\mathcal{I}, \mathcal{E})$ where $\mathcal{I} = \{1, \dots, n\}$ and $\mathcal{E}$ is a set of hyperedges $A_\ell \subset \mathcal{I}$ for $\ell \in \{1, \dots, m\}$, with $m = |\mathcal{E}|$. Let us denote $n_\ell = |A_\ell|$, and $\mathbf{z}_\ell$ the list of the groups of the individuals in $A_\ell$ i.e. $\mathbf{z}_\ell = (z_i)_{i \in A_\ell}$. The un-normalised distribution of the graphical model can be written as a product of $m$ factors indexed by $\ell$:

$$ \psi(z_1, \dots, z_n) = \prod_{\ell=1}^{m} \psi_{\ell}(\mathbf{z}_{\ell}). $$

698

699

700

701

The distribution $\psi$ can be represented by a tensor $\mathbf{\Psi}$ with $n$ modes, each of dimension $Q$: the probability $\psi(z_1, \dots, z_n)$ is the element $\mathbf{\Psi}[z_1, \dots, z_n]$ of the tensor. Each factor can also be represented by a tensor. If, for example, $\psi_{ij}$ is a binary factor, it can be represented by a matrix $\mathbf{\Psi}_{ij}$ of dimension $Q \times Q$, with $\psi_{ij}(z_i, z_j)$ being the element $\mathbf{\Psi}_{ij}[z_i, z_j]$ of the matrix.

**Step 1: Approximating each factor with A TT-format for a given rank.** A tensor $\mathbf{\Psi}_\ell$ can be built from each factor $\psi_\ell$. Its order is $n_\ell$, the number of indices in $A_\ell$, and its dimensions are $Q$ for each mode. We define

$$ \mathbf{\Psi}_{\ell}[\mathbf{z}_{\ell}] = \psi_{\ell}(\mathbf{z}_{\ell}), \quad \text{with} \quad \mathbf{\Psi}_{\ell} \in \mathbb{R}^{Q \times \dots \times Q} $$

The first step consists in selecting a rank $r$ and computing a best TT-approximation $\tilde{\mathbf{\Psi}}_\ell$ of $\mathbf{\Psi}_\ell$ at rank $r$ (more rigorously, at rank $\mathbf{r} = (1, r, \dots, r, 1)$). In the TT-format, each component $\tilde{\mathbf{\Psi}}_\ell[z_{i_1}^\ell, \dots, z_{i_{n_\ell}^\ell}]$ can be written as a product of matrices with separation of variables (here, $\mathbf{z}_\ell = (z_{i_1}^\ell, \dots, z_{i_{n_\ell}^\ell})$ where $i_\ell^\ell$ is $\mu$-th index in $A_\ell$):

$$ \tilde{\mathbf{\Psi}}_{\ell}[z_{i_1}^{\ell}, \dots, z_{i_{n_\ell}^{\ell}}] = G_{i_1}^{\ell}(z_{i_1}^{\ell}) \dots G_{i_{n_\ell}^{\ell}}^{\ell}(z_{i_{n_\ell}^{\ell}}), \\ \text{with} \quad \begin{cases} G_{i_1}^{\ell}(z_{i_1}^{\ell}) \in \mathbb{R}^{1 \times r} \\ G_{i_\mu}^{\ell}(z_{i_\mu}^{\ell}) \in \mathbb{R}^{r \times r} & i_2^{\ell} \le i_\mu^{\ell} \le i_{n_\ell-1}^{\ell} \\ G_{i_{n_\ell}^{\ell}}^{\ell}(z_{i_{n_\ell}^{\ell}}) \in \mathbb{R}^{r \times 1} \end{cases} $$