The variables are the $\tau_{l,q}$ for $l \in \{1, \dots, n\}$ and $q \in \{1, \dots, Q\}$, and the $\sigma_l$ for $l \in \{1, \dots, n\}$. Let us compute the derivatives of $G$.

$$ \frac{\partial G}{\partial \sigma_l} = \sum_{q=1}^{Q} \tau_{l,q} - 1 $$

So $\frac{\partial G}{\partial \sigma_l} = 0$ leads to

$$ \sum_{q=1}^{Q} \tau_{l,q} = 1 \qquad (A1) $$

Then using the fact that $D$ and $\Lambda$ are symmetric matrices, we obtain:

$$ \begin{aligned} \frac{\partial G}{\partial \tau_{l,q}} &= \ln(\tau_{l,q}) + 1 + \sigma_l - \ln(\alpha_q) - \sum_{j=1}^{l-1} \sum_{q'=1}^{Q} \tau_{j,q'} [\ln \mathbb{P}_{\Lambda}(D[j, l] | Z_l = q, Z_j = q')] \\ &\quad - \sum_{j>l} \sum_{q'=1}^{Q} \tau_{j,q'} [\ln \mathbb{P}_{\Lambda}(D[j, j] | Z_l = q, Z_j = q')] \\ &= \ln(\tau_{l,q}) + 1 + \sigma_l - \ln(\alpha_q) - \sum_{j \neq l} \sum_{q'=1}^{Q} \tau_{j,q'} [\ln \mathbb{P}_{\Lambda}(D[j, j] | Z_l = q, Z_j = q')] \end{aligned} $$

So $\frac{\partial g(\sigma_l, \tau_{l,q})}{\partial \tau_{l,q}} = 0$ leads to

$$ \tau_{l,q} = \alpha_q \exp \left( \sum_{j \neq l} \sum_{q'=1}^{Q} \tau_{j,q'} [\ln \mathbb{P}_{\Lambda}(D[l, j] | Z_l = q, Z_j = q')] \right) \exp(-(\sigma_l + 1) \quad (A2) $$

We define

$$ \mu_{l,q} = \alpha_q \prod_{j \neq l}^{n} \prod_{q'=1}^{Q} [\mathbb{P}_{\Lambda}(D[l, j] | Z_l = q, Z_j = q')]^{\tau_{j,q'}} $$

Then

$$ \tau_{l,q} = \mu_{l,q} e^{-(\sigma_l+1)} \qquad (A3) $$

Finally, by combining equation (A1) and equation (A3) we obtain

$$ \tau_{l,q} = \frac{\mu_{l,q}}{\sum_{q=1}^{Q} \mu_{l,q}} \qquad (A4) $$

So the $\tau_{l,q}$ are solution of a fixed point equation, since $\mu_{l,q}$ is a function of the $\tau_{j,q'}$ for $j \neq l$.

The associated fixed point iteration scheme is

$$ \forall 1 \le l \le n, \forall 1 \le q \le Q, \quad \left\{ \begin{array}{l} \tau_{l,q}^{t+1} = \frac{\mu_{l,q}^{t+1/2}}{\sum_{q'=1}^{Q} \mu_{l,q'}^{t+1/2}} \\ \mu_{l,q}^{t+1/2} = \alpha_q \prod_{j \neq l} \prod_{q'=1}^{Q} [\mathbb{P}_{\Lambda}(D[l, j] | Z_l = q, Z_j = q')]^{\tau_{j,q'}} \end{array} \right. $$
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 $\widehat{(\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. The idea is to build for each hyperedge $A_\ell$ a new factor $\bar{\psi}_\ell$ that depends on all variables and not only on $z_\ell$ but whose value does not actually depend on $z_j$ for $j \notin A_\ell$ (these variables $z_j$ are called non essential variables):

$$ \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} $$

702 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 (10) $$

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 example (from Novikov et al. (2014), section 5.1), 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}{ll} 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. $$

703 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} $$
In the case of a Poisson distribution of the dissimilarities, we obtain

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

## Appendix B Matrix product in TT-matrix format

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

Let A, B be TT-matrices with

$$ A[(i_1, i_2, i_3); (j_1, j_2, j_3)] = M_1(i_1, j_1) M_2(i_2, j_2) M_3(i_3, j_3) \quad (B5) $$

and

$$ B[(i_1, i_2, i_3); (j_1, j_2, j_3)] = N_1(i_1, j_1) N_2(i_2, j_2) N_3(i_3, j_3) \quad (B6) $$

with $1 \le i_\mu, j_\mu \le n$. Then

$$ C[(i_1, i_2, i_3); (k_1, k_2, k_3)] = \sum_{j_1, j_2, j_3=1}^{n} A[(i_1, i_2, i_3); (j_1, j_2, j_3)] B[(j_1, j_2, j_3); (k_1, k_2, k_3)] \\ = \sum_{j_1, j_2, j_3=1}^{n} M_1(i_1, j_1) M_2(i_2, j_2) M_3(i_3, j_3) N_1(j_1, k_1) N_2(j_2, k_2) N_3(j_3, k_3) \quad (B7) $$

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

<table><thead><tr><th>Matrix</th><th>size</th><th>coefficients</th></tr></thead><tbody><tr><td>M<sub>1</sub>(i<sub>1</sub>, j<sub>1</sub>)</td><td>1 &times; r</td><td>M<sub>1</sub>[i<sub>1</sub>, j<sub>1</sub>, &alpha;<sub>1</sub>]</td></tr><tr><td>M<sub>2</sub>(i<sub>2</sub>, j<sub>2</sub>)</td><td>r &times; r</td><td>M<sub>2</sub>[&alpha;<sub>1</sub>, i<sub>2</sub>, j<sub>2</sub>, &alpha;<sub>2</sub>]</td></tr><tr><td>M<sub>3</sub>(i<sub>3</sub>, j<sub>3</sub>)</td><td>r &times; 1</td><td>M<sub>3</sub>[&alpha;<sub>2</sub>, i<sub>3</sub>, j<sub>3</sub>]</td></tr><tr><td>N<sub>1</sub>(j<sub>1</sub>, k<sub>1</sub>)</td><td>1 &times; r</td><td>N<sub>1</sub>[j<sub>1</sub>, k<sub>1</sub>, &beta;<sub>1</sub>]</td></tr><tr><td>N<sub>2</sub>(j<sub>2</sub>, k<sub>2</sub>)</td><td>r &times; r</td><td>N<sub>2</sub>[&beta;<sub>1</sub>, j<sub>2</sub>, k<sub>2</sub>, &beta;<sub>2</sub>]</td></tr><tr><td>N<sub>3</sub>(j<sub>3</sub>, k<sub>3</sub>)</td><td>r &times; 1</td><td>N<sub>3</sub>[&beta;<sub>2</sub>, j<sub>3</sub>, k<sub>3</sub>]</td></tr></tbody></table>

with $1 \le \alpha_i, \beta_j \le r$. Then

$$ M_1(i_1, j_1) M_2(i_2, j_2) M_3(i_3, j_3) N_1(j_1, k_1) N_2(j_2, k_2) N_3(j_3, k_3) \\ = \sum_{\alpha_1, \alpha_2=1}^{r} \sum_{\beta_1, \beta_2=1}^{r} M_1[i_1, j_1, \alpha_1] M_2[\alpha_1, i_2, j_2, \alpha_2] M_3[\alpha_2, i_3, j_3] \times $$
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,r}$ and $G_7^{\ell}(z_7) = 1$ as in the following synthetic representation of
$\bar{\Psi}_{\ell}[z_1, z_2, z_3, z_4, z_5, z_6, z_7]$:

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

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

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

**Step 3: TT decomposition of the joint distribution.** Here $\otimes$ denotes the Kronecker product between matrices. Knowing that $xy = x \otimes y$ if $x, y \in \mathbb{R}$, equation (10) can be written

$$
\begin{align*}
\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{align*}
$$

710 Let us recall the mixed-product property of Kronecker product (see Horn & Johnson 2012):

711 if A, B, C, D are matrices with relevant dimensions for the products AB and CD to be

712 possible, then

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

This leads to

$$
\mathbf{\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)
$$
$$N_1[j_1, k_1, \beta_1] N_2[\beta_1, j_2, k_2, \beta_2] N_3[\beta_2, j_3, k_3] \qquad (B8)$$

Let us reorder the terms after the $\sum$ by grouping $(\alpha_1, \beta_1)$ and $(\alpha_2, \beta_2)$ as

$$M_1[i_1, j_1, \alpha_1] N_1[j_1, k_1, \beta_1] M_2[\alpha_1, i_2, j_2, \alpha_2] N_2[\beta_1, j_2, k_2, \beta_2] M_3[\alpha_2, i_3, j_3] N_3[\beta_2, j_3, k_3] \qquad (B9)$$

and define, with $\gamma_i$ coding for $(\alpha_i, \beta_i)$ with $\gamma_i = r(\alpha_i - 1) + \beta_i$,

$$\underbrace{M_1[i_1, j_1, \alpha_1] N_1[j_1, k_1, \beta_1]}_{=P_1[i_1, j_1, k_1, \gamma_1]} \underbrace{M_2[\alpha_1, i_2, j_2, \alpha_2] N_2[\beta_1, j_2, k_2, \beta_2]}_{=P_2[\gamma_1, i_2, j_2, k_2, \gamma_2]} \underbrace{M_3[\alpha_2, i_3, j_3] N_3[\beta_2, j_3, k_3]}_{=P_3[\gamma_2, i_3, j_3, k_3]} \qquad (B10)$$

Let us note that, as indices $\alpha_i, \beta_i$ run each over $\{1, \dots, r\}$, indices $\gamma_i$ run over $\{1, \dots, r^2\}$. Then

$$C[(i_1, i_2, i_3); (k_1, k_2, k_3)] \\
= \sum_{j_1, j_2, j_3=1}^{n} \sum_{\gamma_1, \gamma_2=1}^{r^2} P_1[i_1, j_1, k_1, \gamma_1] P_2[\gamma_1, i_2, j_2, k_2, \gamma_2] P_3[\gamma_2, i_3, j_3, k_3] \\
= \sum_{\gamma_1, \gamma_2=1}^{r^2} \sum_{j_1, j_2, j_3=1}^{n} P_1[i_1, j_1, k_1, \gamma_1] P_2[\gamma_1, i_2, j_2, k_2, \gamma_2] P_3[\gamma_2, i_3, j_3, k_3] \\
= \sum_{\gamma_1, \gamma_2=1}^{r^2} \left( \sum_{j_1=1}^{n} P_1[i_1, j_1, k_1, \gamma_1] \right) \left( \sum_{j_2=1}^{n} P_2[\gamma_1, i_2, j_2, k_2, \gamma_2] \right) \left( \sum_{j_3=1}^{n} P_3[\gamma_2, i_3, j_3, k_3] \right) \qquad (B11)$$

Let us define $Q_1, Q_2, Q_3$ by

$$\left\{ \begin{array}{l} Q_1[i_1, k_1, \gamma_1] = \sum_{j_1=1}^{n} P_1[i_1, j_1, k_1, \gamma_1] \\ Q_2[\gamma_1, i_2, k_2, \gamma_2] = \sum_{j_2=1}^{n} P_2[\gamma_1, i_2, j_2, k_2, \gamma_2] \\ Q_3[\gamma_2, i_3, k_3] = \sum_{j_3=1}^{n} P_3[\gamma_2, i_3, j_3, k_3], \end{array} \right. \qquad (B12)$$

with

<table><thead><tr><th>Matrix</th><th>Size</th><th>Coefficients</th><th>as a tensor</th></tr></thead><tbody><tr><td>Q<sub>1</sub>(i<sub>1</sub>, k<sub>1</sub>)</td><td>1 × r<sup>2</sup></td><td>Q<sub>1</sub>(i<sub>1</sub>, k<sub>1</sub>)[γ<sub>1</sub>]</td><td>= Q<sub>1</sub>[i<sub>1</sub>, k<sub>1</sub>, γ<sub>1</sub>]</td></tr><tr><td>Q<sub>2</sub>(i<sub>2</sub>, k<sub>2</sub>)</td><td>r<sup>2</sup> × r<sup>2</sup></td><td>Q<sub>2</sub>(i<sub>2</sub>, k<sub>2</sub>)[γ<sub>1</sub>, γ<sub>2</sub>]</td><td>= Q<sub>2</sub>[γ<sub>1</sub>, i<sub>2</sub>, k<sub>2</sub>, γ<sub>2</sub>]</td></tr><tr><td>Q<sub>3</sub>(i<sub>3</sub>, k<sub>3</sub>)</td><td>r<sup>2</sup> × 1</td><td>Q<sub>3</sub>(i<sub>3</sub>, k<sub>3</sub>)[γ<sub>2</sub>]</td><td>= Q<sub>3</sub>[γ<sub>2</sub>, i<sub>3</sub>, k<sub>3</sub>].</td></tr></tbody></table>
713 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 (11) $$

714 which is separation of variables.

**Step 4: Computation of the normalising constant.** 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} $$

715 It requires $nQ$ matrix additions and $n-1$ matrix multiplications, provided that 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).

716 Evaluating the complexity of this calculation is not straightforward, because it involves 717 products of matrices, some of which are sparse. Such an evaluation will be done for the 718 simple case $n = 3$ in section 5.4. 719

**Computation of the unary marginals:** Expression (11) permits to compute the unary marginals as well, as quickly mentionned in Novikov et al. (2014). We provide here the expression of these marginals. The un-normalised unary marginal for variable $i$ is defined as

$$ m_i(q) = \sum_{z \in \{1, \dots, Q\}^n, \text{s.t. } z_i=q} \psi(z_1, \dots, z_n) $$
Then

$$ 
\begin{aligned}
C[(i_1, i_2, i_3); (k_1, k_2, k_3)] &= \sum_{\gamma_1, \gamma_2=1}^{r^2} Q_1[i_1, k_1, \gamma_1] Q_2[\gamma_1, i_2, j_2, \gamma_2] Q_3[\gamma_2, i_3, k_3] \quad (B13) \\
&= Q_1(i_1, k_1) Q_2(i_2, k_2) Q_3(i_3, k_3),
\end{aligned}
$$

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

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

$$ 
\left\{ 
\begin{aligned}
A[(i_1, \dots, i_m); (j_1, \dots, j_m)] &= \prod_{\mu=1}^{m} M_{\mu}(i_{\mu}, j_{\mu}) \\
B[(i_1, \dots, i_m); (j_1, \dots, j_m)] &= \prod_{\mu=1}^{m} N_{\mu}(i_{\mu}, j_{\mu}),
\end{aligned}
\right. 
$$

with $1 \le i_\mu, j_\mu \le n$. Then, there are $n^2$ matrices $Q_\mu(i_\mu, k_\mu)$ indexed by $(i_\mu, k_\mu)$ for a given $\mu$. A matrix $Q_\mu(i_\mu, k_\mu)$ has $r^2$ terms, each being a sum of n terms of a matrix $P_\mu(i_\mu, j_\mu, k_\mu)$ (see equation (B12)), each term of a $P_\mu$ being a product of two real coefficients. So, there are n products for a term of a $Q_\mu(i_\mu, k_\mu)$, $nr^2$ products for a matrix $Q_\mu(i_\mu, j_\mu)$, and $n^3 r^2$ products for all matrices $Q_\mu(i_\mu, k_\mu)$ for a given $\mu$, and $mr^2 n^3$ products for all matrices $Q_\mu(i_\mu, k_\mu)$ for all cores $\mu$. Hence, the complexity of the matrix product in TT-format is in $O(mr^2n^3)$. This is for the product of 2 matrices. The complexity of the product of p matrices is in $O(mr^p n^3)$, and grows exponentially with p. Each matrix A, B has dimensions $n^m \times n^m$, hence their product requires in general $(n^m)^3 = n^{3m}$ products, versus $mr^2 n^3$ in TT-matrix format.

# Appendix C Computation of matrices $A_i(z_i)$ for a WSBM

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

## Adding non essential variables.

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

$$ \psi_{24}(z_2, z_4) = M_{24}[z_2] V_{42}[z_4] \quad \text{with} \quad \begin{cases} M_{24}[z_2] & \in \mathbb{R}^{1 \times Q} \\ V^{24}[z_4] & \in \mathbb{R}^{Q \times 1} \end{cases} $$
Then, plugging (11) into this expression we obtain

$$
\begin{align*}
m_i(z_i) &= \sum_{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, in a message passing algorithm.

Appendix II

Let A and B be two matrices with dimensions such that the product $C = AB$ exists.

Here, we show how to compute the cores of $C$ expressed in TT-matrix format knowing only

the cores of A and B expressed in TT-matrix format. Although it is classical, it is seldom

presented in articles. We develop the calculation in detail on a simple example with 3 cores,

and evaluate its complexity in a more general case.

Let A, B be TT-matrices with

$$
A[(i_1, i_2, i_3); (j_1, j_2, j_3)] = M_1(i_1, j_1) M_2(i_2, j_2) M_3(i_3, j_3)
$$

and

$$
B[(i_1, i_2, i_3); (j_1, j_2, j_3)] = N_1(i_1, j_1) N_2(i_2, j_2) N_3(i_3, j_3)
$$

with $1 \le i_\mu, j_\mu \le n$. Then

$$
\begin{align*}
C[(i_1, i_2, i_3); (k_1, k_2, k_3)] &= \sum_{j_1, j_2, j_3=1}^{n} A[(i_1, i_2, i_3); (j_1, j_2, j_3)] B[(j_1, j_2, j_3); (k_1, k_2, k_3)] \\
&= \sum_{j_1, j_2, j_3=1}^{n} M_1(i_1, j_1) M_2(i_2, j_2) M_3(i_3, j_3) N_1(j_1, k_1) N_2(j_2, k_2) N_3(j_3, k_3)
\end{align*}
$$

Would the matrices product be commutative, reordering the terms would lead to

simplifications. But it is not the case. However, reordering can be done in $\mathbb{R}$ between

coefficients. Therefore, let us introduce the following notations:
Completing with non essential variables (see step 2 in section 4) leads to the definition of factor $\bar{\psi}_{24}$ with non essential variables $z_1, z_3$ which can be written in TT format as follows (× is the matrices product)

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

with, $\forall z_1, z_2, z_3, z_4$:

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

which can be sketched as

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

The product of those matrices $(1, 1)$, $(1, Q)$, $(Q, Q)$ and $(Q, 1)$ is a real. This is variable separation for the factor $\bar{\psi}_{24}(z_1, z_2, z_3, z_4)$.

This can be generalized to any $n$ by computing all $G_{ij}^{(k)}[z_k]$ for $1 \le k \le n$ as

$$ \left\{ \begin{array}{ll} k < i & \Rightarrow G_{ij}^{(k)}[z_k] = 1 \quad \in \mathbb{R}^{1 \times 1} \\ k = i & \Rightarrow G_{ij}^{(k)}[z_i] = M_{ij}[z_i] \quad \in \mathbb{R}^{1 \times Q} \\ i < k < j & \Rightarrow G_{ij}^{(k)}[z_k] = \mathbb{I}_{Q,Q} \quad \in \mathbb{R}^{Q \times Q} \\ k = j & \Rightarrow G_{ij}^{(k)}[z_j] = V_{ij}[z_j] \quad \in \mathbb{R}^{Q \times 1} \\ k > j & \Rightarrow G_{ij}^{(k)}[z_k] = 1 \quad \in \mathbb{R}^{1 \times 1} \end{array} \right. \qquad (C14) $$

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

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

where the product is the matrices product. This is separation of variables for factor $\bar{\psi}_{ij}$.
<table><thead><tr><th>Matrix</th><th>size</th><th>coefficients</th><th>as a tensor</th></tr></thead><tbody><tr><td>M<sub>1</sub>(i<sub>1</sub>, j<sub>1</sub>)</td><td>1 &times; r</td><td>M<sub>1</sub>(i<sub>1</sub>, j<sub>1</sub>)[&alpha;<sub>1</sub>]</td><td>M<sub>1</sub>[i<sub>1</sub>, j<sub>1</sub>, &alpha;<sub>1</sub>]</td></tr><tr><td>M<sub>2</sub>(i<sub>2</sub>, j<sub>2</sub>)</td><td>r &times; r</td><td>M<sub>2</sub>(i<sub>2</sub>, j<sub>2</sub>)[&alpha;<sub>1</sub>, &alpha;<sub>2</sub>]</td><td>M<sub>2</sub>[&alpha;<sub>1</sub>, i<sub>2</sub>, j<sub>2</sub>, &alpha;<sub>2</sub>]</td></tr><tr><td>M<sub>3</sub>(i<sub>3</sub>, j<sub>3</sub>)</td><td>r &times; 1</td><td>M<sub>3</sub>(i<sub>3</sub>, j<sub>3</sub>)[&alpha;<sub>2</sub>]</td><td>M<sub>3</sub>[&alpha;<sub>2</sub>, i<sub>3</sub>, j<sub>3</sub>]</td></tr><tr><td>N<sub>1</sub>(j<sub>1</sub>, k<sub>1</sub>)</td><td>1 &times; r</td><td>N<sub>1</sub>(j<sub>1</sub>, k<sub>1</sub>)[&beta;<sub>1</sub>]</td><td>N<sub>1</sub>[j<sub>1</sub>, k<sub>1</sub>, &beta;<sub>1</sub>]</td></tr><tr><td>N<sub>2</sub>(j<sub>2</sub>, k<sub>2</sub>)</td><td>r &times; r</td><td>N<sub>2</sub>(j<sub>2</sub>, k<sub>2</sub>)[&beta;<sub>1</sub>, &beta;<sub>2</sub>]</td><td>N<sub>2</sub>[&beta;<sub>1</sub>, j<sub>2</sub>, k<sub>2</sub>, &beta;<sub>2</sub>]</td></tr><tr><td>N<sub>3</sub>(j<sub>3</sub>, k<sub>3</sub>)</td><td>r &times; 1</td><td>N<sub>3</sub>(j<sub>3</sub>, k<sub>3</sub>)[&beta;<sub>2</sub>]</td><td>N<sub>3</sub>[&beta;<sub>2</sub>, j<sub>3</sub>, k<sub>3</sub>]</td></tr></tbody></table>

with $1 \le \alpha_i, \beta_j \le r$. Then

$$M_1(i_1, j_1) M_2(i_2, j_2) M_3(i_3, j_3) N_1(j_1, k_1) N_2(j_2, k_2) N_3(j_3, k_3) = \sum_{\alpha_1, \alpha_2=1}^{r} \sum_{\beta_1, \beta_2=1}^{r} M_1[i_1, j_1, \alpha_1] M_2[\alpha_1, i_2, j_2, \alpha_2] M_3[\alpha_2, i_3, j_3] \times N_1[j_1, k_1, \beta_1] N_2[\beta_1, j_2, k_2, \beta_2] N_3[\beta_2, j_3, k_3]$$

Let us reorder the terms after the $\sum$ by grouping $(\alpha_1, \beta_1)$ and $(\alpha_2, \beta_2)$ as

$M_1[i_1, j_1, \alpha_1] N_1[j_1, k_1, \beta_1] M_2[\alpha_1, i_2, j_2, \alpha_2] N_2[\beta_1, j_2, k_2, \beta_2] M_3[\alpha_2, i_3, j_3] N_3[\beta_2, j_3, k_3]$ and define, with $\gamma_i$ coding for $(\alpha_i, \beta_i)$ with $\gamma_i = r(\alpha_i - 1) + \beta_i$,

$$\underbrace{M_1[i_1, j_1, \alpha_1]}_{=P_1[i_1, j_1, k_1, \gamma_1]} \underbrace{N_1[j_1, k_1, \beta_1]}_{=P_2[\gamma_1, i_2, j_2, k_2, \gamma_2]} \underbrace{M_2[\alpha_1, i_2, j_2, \alpha_2]}_{=P_3[\gamma_2, i_3, j_3, k_3]} \underbrace{N_2[\beta_1, j_2, k_2, \beta_2]}_{=P_4[\gamma_2, i_3, j_3, k_3]} \underbrace{M_3[\alpha_2, i_3, j_3]}_{=P_5[\gamma_2, i_3, j_3, k_3]} \underbrace{N_3[\beta_2, j_3, k_3]}_{=P_6[\gamma_2, i_3, j_3, k_3]}$$

Let us note that, as indices $\alpha_i, \beta_i$ run each over $\{1, \dots, r\}$, indices $\gamma_i$ run over $\{1, \dots, r^2\}$. Then

$$C[(i_1, i_2, i_3); (k_1, k_2, k_3)] = \sum_{j_1, j_2, j_3=1}^{n} \sum_{\gamma_1, \gamma_2=1}^{r^2} P_1[i_1, j_1, k_1, \gamma_1] P_2[\gamma_1, i_2, j_2, k_2, \gamma_2] P_3[\gamma_2, i_3, j_3, k_3] \\ = \sum_{\gamma_1, \gamma_2=1}^{r^2} \sum_{j_1, j_2, j_3=1}^{n} P_1[i_1, j_1, k_1, \gamma_1] P_2[\gamma_1, i_2, j_2, k_2, \gamma_2] P_3[\gamma_2, i_3, j_3, k_3] \\ = \sum_{\gamma_1, \gamma_2=1}^{r^2} \left( \sum_{j_1=1}^{n} P_1[i_1, j_1, k_1, \gamma_1] \right) \left( \sum_{j_2=1}^{n} P_2[\gamma_1, i_2, j_2, k_2, \gamma_2] \right) \left( \sum_{j_3=1}^{n} P_3[\gamma_2, i_3, j_3, k_3] \right)$$
Computation of the matrices $A_k(z_k)$.

Here, we can 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).$$

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

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

We then have

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

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$. Let us note that the TT-rank of TT-matrix $A_k(z_k)$ is one, which is a useful property for deriving the fusion of cores in the calculation. To show this, let us write $A_k(z_k)$ in TT-matrix format. Therefore, let $G_{ij}^{(k)}(z_k)$ be with coefficients $G_{ij}^{(k)}(z_k) [\alpha_{ij}^{(k)}, \beta_{ij}^{(k)}]$ with $1 \le \alpha_{ij}^{(k)}, \beta_{ij}^{(k)} \le r$. Then, from the definition of the Kronecker product,

$$A_k(z_k) [ (\alpha_{12}^{(k)}, \dots, \alpha_{n-1,n}^{(k)}); (\beta_{12}^{(k)}, \dots, \beta_{n-1,n}^{(k)}) ] = \prod_{1 \le i < j \le n} G_{ij}^{(k)}(z_k) [\alpha_{ij}^{(k)}, \beta_{ij}^{(k)}]$$

which can be written

$$A_k(z_k) [ (\alpha_{12}^{(k)}, \dots, \alpha_{n-1,n}^{(k)}); (\beta_{12}^{(k)}, \dots, \beta_{n-1,n}^{(k)}) ] = \prod_{1 \le i < j \le n} G_{ij}^{(k),z_k}(\alpha_{ij}^{(k)}, \beta_{ij}^{(k)})$$

by reordering the upper and lower indices, and shows that the TT-rank is one as a product of scalars. It can be seen on a simpler example: Let $A = G_1 \otimes G_2$ be a matrix, with $G_1 := G_1[i_1, j_1]$ and $G_2 := G_2[i_2, j_2]$ with coefficients. Then, by definition of Kronecker product, $A$ can be expressed as a TT-matrix of TT-rank one, by $A[(i_1, i_2); (j_1, j_2)] = G_1[i_1, j_1] G_2[i_2, j_2] = g_1(i_1, j_1) g_2(i_2, j_2)$ where $g_1, g_2$ are scalars. As the TT-rank of $A_k(z_k)$ is one, the TT-rank of TT-matrix $B_k = \sum_{z_k} A_k(z_k)$ is $Q$.

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

Here we show that, given two variables $Z_i$ and $Z_j$, taking value in $\{1, \dots, Q\}$, such that there exists $q$ with $P(Z_i = q) = 1$, then, the binary marginals $P(Z_i, Z_j)$ are the products of the unary marginals. Let us have $Q = 3$. Binary marginals can be organized as elements in a $3 \times 3$ matrix. The unary marginals are the row-wise or
733 Let us define $Q_1, Q_2, Q_3$ by

$$ egin{cases} Q_1[i_1, k_1, \gamma_1] = \sum_{j_1=1}^{n} P_1[i_1, j_1, k_1, \gamma_1] \\ Q_2[\gamma_1, i_2, k_2, \gamma_2] = \sum_{j_2=1}^{n} P_2[\gamma_1, i_2, j_2, k_2, \gamma_2] \\ Q_3[\gamma_2, i_3, k_3] = \sum_{j_3=1}^{n} P_3[\gamma_2, i_3, j_3, k_3], \end{cases} \quad (12) $$

734 with

<table><thead><tr><th>Matrix</th><th>Size</th><th>Coefficients</th><th>as a tensor</th></tr></thead><tbody><tr><td>$Q_1(i_1, k_1)$</td><td>$1 \times r^2$</td><td>$Q_1(i_1, k_1)[\gamma_1]$</td><td>$= Q_1[i_1, k_1, \gamma_1]$</td></tr><tr><td>$Q_2(i_2, k_2)$</td><td>$r^2 \times r^2$</td><td>$Q_2(i_2, k_2)[\gamma_1, \gamma_2]$</td><td>$= Q_2[\gamma_1, i_2, k_2, \gamma_2]$</td></tr><tr><td>$Q_3(i_3, k_3)$</td><td>$r^2 \times 1$</td><td>$Q_3(i_3, k_3)[\gamma_2]$</td><td>$= Q_3[\gamma_2, i_3, k_3]$.</td></tr></tbody></table>

735

Then

$$ egin{aligned} C[(i_1, i_2, i_3); (k_1, k_2, k_3)] &= \sum_{\gamma_1, \gamma_2=1}^{r^2} Q_1[i_1, k_1, \gamma_1] Q_2[\gamma_1, i_2, j_2, \gamma_2] Q_3[\gamma_2, i_3, k_3] \\ &= Q_1(i_1, k_1) Q_2(i_2, k_2) Q_3(i_3, k_3), \end{aligned} $$

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

737

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

$$ egin{cases} A[(i_1, \dots, i_m); (j_1, \dots, j_m)] = \prod_{\mu=1}^{m} M_{\mu}(i_{\mu}, j_{\mu}) \\ B[(i_1, \dots, i_m); (j_1, \dots, j_m)] = \prod_{\mu=1}^{m} N_{\mu}(i_{\mu}, j_{\mu}), \end{cases} $$

741 with $1 \le i_\mu, j_\mu \le n$. Then, there are $n^2$ matrices $Q_\mu(i_\mu, k_\mu)$ indexed by $(i_\mu, k_\mu)$ for a given $\mu$. A matrix $Q_\mu(i_\mu, k_\mu)$ has $r^2$ terms, each being a sum of $n$ terms of a matrix $P_\mu(i_\mu, j_\mu, k_\mu)$ (see equation (12)), each term of a $P_\mu$ being a product of two real coefficients. So, there are $n$ products for a term of a $Q_\mu(i_\mu, k_\mu)$, $nr^2$ products for a matrix $Q_\mu(i_\mu, j_\mu)$, and $n^3 r^2$ products for all matrices $Q_\mu(i_\mu, k_\mu)$ for a given $\mu$, and $mr^2 n^3$ products for all matrices $Q_\mu(i_\mu, k_\mu)$ for all cores $\mu$. Hence, the complexity of the matrix product in TT-format is in $O(mr^2n^3)$. This is for the product of 2 matrices. The complexity of the product of $p$ matrices is in $O(mrp^n)$,
column-wise sums. It can be displayed as

<table><tr><td>a</td><td>b</td><td>c</td><td>x</td></tr><tr><td>a'</td><td>b'</td><td>c'</td><td>x'</td></tr><tr><td>a''</td><td>b''</td><td>c''</td><td>x''</td></tr><tr><td>y</td><td>y'</td><td>y''</td><td>1</td></tr></table>

with

$$ \left\{ \begin{array}{lcl} a+b+c & = & x \\ a'+b'+c' & = & x' \\ a''+b''+c'' & = & x'' \\ a+a'+a'' & = & y \\ b+b'+b'' & = & y' \\ c+c'+c'' & = & y'' \end{array} \right. $$

and all terms and marginals in this matrix being in [0, 1]. Let us assume without loss of generality that $x = 1$. So, $x' = x'' = 0$. From $a', b', c' \ge 0$ and $a' + b' + c' = x' = 0$, it follows that $a' = b' = c' = 0$, and the same for $a'', b'', c''$. So, the matrix is

<table><tr><td>a</td><td>b</td><td>c</td><td>1</td></tr><tr><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>y</td><td>y'</td><td>y''</td><td>1</td></tr></table>

We can observe that each term in the matrix is the product of the marginals of its row and column. For example, $a = y \times 1$ as $y = a + 0 + 0 = a$. Note that it is not required that $y$, or $y'$ or $y''$ is 1 or 0. This demonstration can be extended to $Q > 3$.

## References

* Abouabdallah, M.A., N. Peyrard, and A. Franc. 2022. Does clustering of DNA bar-codes agree with botanical classification directly at high taxonomic levels? *Trees in french guiana* as a case study. *Molecular Ecology Resources* 22(5): 1746–1761. https://doi.org/10.1111/1755-0998.13579.

* Al Daas, H., G. Ballard, and P. Benner. 2022. Parallel algorithms for tensor-train arithmetic. *SIAM Journal of Scientific Computing* 44: C25–C53. https://doi.org/10.1137/20M1387158.

* Al Daas, H., G. Ballard, P. Cazeaux, E. Hallman, A. Miedlar, M. Pasha, T.W. Reid, and A.K. Saibaba. 2023. Randomized algorithms for rounding in the Tensor-Train format. *SIAM Journal of Scientific Computing* 45: A74–195. https://doi.org/10.1137/21M1451191.

* Ansari, M., A. Ahmad, S. Khan, G. Bhushan, and M. Siddique. 2020. Spatiotemporal clustering: a review. *Artificial Intelligence Review* 53(4): 2381–2423. https://doi.org/10.1007/s10462-019-09736-1.

* Barbillon, P., S. Donnet, E. Lazega, and A. Bar-Hen. 2017. Stochastic block models for multiplex networks: An application to a multilevel network of researchers. *Journal*
748 and grows exponentially with $p$. Each matrix $A, B$ has dimensions $n^m \times n^m$, hence their
749 product requires in general $(n^m)^3 = n^{3m}$ products, versus $mr^2n^3$ in TT-matrix format.

750

Appendix III

**Mean field principle** The Mean-Field approximation consists in approximating $\mathrm{Pr}_{\theta}(Z | D)$ by another joint distribution $\mathbb{Q}_{\theta}(Z | D)$ for which computation of marginals is easier. Under the distribution $\mathbb{Q}_{\theta}$ the variables $Z_i$ are independent, i.e. $\mathbb{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 $\mathbb{Q}_{\theta}$. The distribution $\mathbb{Q}_{\theta}$ is the one that miminises the Kullback-Leibler divergence between a distribution of independent variables and $\mathrm{Pr}_{\theta}(Z | D)$. We recall the definition of the Kullback-Leibler divergence between two distributions $\mathbb{Q}$ and $\mathrm{Pr}$:

$$KL(\mathbb{Q}|\mathrm{Pr}) = \mathrm{E}_{\mathbb{Q}}\left[\ln\left(\frac{\mathbb{Q}(Z)}{\mathrm{Pr}(Z)}\right)\right]$$

**Mean-Field as solution of a fixed point equation** Let us define the family $\mathcal{Q}$ of distributions $\mathcal{Q}$ such that $\mathcal{Q}(Z) = \prod_{i=1}^{n} q^{i}(Z_i)$, and let us denote $\tau_{i,q} = q^{i}(Z_i = q)$. The quantity $\mathbb{Q}_{\theta}(\cdot | D)$ is solution of

$$\mathbb{Q}_{\theta}(\cdot|D) = \arg \min_{\mathbb{Q} \in \mathcal{Q}} KL(\mathbb{Q} | \mathrm{Pr}_{\theta}(\cdot | D))$$

We have

$$KL(\mathbb{Q}(\cdot) | \mathrm{Pr}_{\theta}(\cdot | D)) = \mathrm{E}_{\mathbb{Q}}[\ln(\mathbb{Q})] - \mathrm{E}_{\mathbb{Q}}[\ln(\mathrm{Pr}_{\theta}(Z | D))]$$
of the Royal Statistical Society: Series A Statistics in Society 180(1): 295-314. https://doi.org/10.1111/rssa.12193.

Carroll, J.D. and J.J. Chang. 1970. Analysis of individual differences in multidimensional scaling via *n*-ways generalization of eckart-young decomposition. *Psychometrika* 35: 283-319.

Casella, G. and E.I. George. 1992. Explaining the Gibbs Sampler. *The American Statistician* 46(3): 167-174.

Coppi, R. and S. Bolasco. 1989. *Multiway data analysis*. Amsterdam: Elsevier.

Daudin, J.J., F. Picard, and S. Robin. 2008. A mixture model for random graph. *Statis- tics and Computing* 18(2): 173-183. https://doi.org/10.1007/s11222-007-9046-7

de Lathauwer, L., B. de Moor, and J. Vandewalle. 2000. A Multilinear Singular Value Decomposition. *SIAM Journal on Matrix Analysis and Applications* 21(4): 1253.

de Silva, V. and L.H. Lim. 2008. Tensor rank and the ill posedness of the best low rank approximation, problem. *SIAM Journal on Matrix Analysis and Applications* 30(3): 1084-1027.

Dempster, A., N.M. Laird, and D.B. Rubin. 1977. Maximum likelihood from incomplete data via the EM algorithm. *Journal of the Royal Statistical Society* 39: 1-22.

Deza, M.M. and E. Deza. 2016. *Encyclopedia of Distances* (Fourth ed.). Springer.

Ducamp, G., P. Bonnard, A. Nouy, and P.H. Wuillemin 2020. An efficient low-rank tensors representation for algorithms in complex probabilistic graphical models. In *Proceedings of the 10th International Conference on Probabilistic Graphical Models*, Skørping, Denmark.

Fannes, M., B. Nachtergaele, and R.F. Werner. 1992. Finitely Correlated States on Quantum Spin Chains. *Communications in Mathematical Physics* 144: 443-490.

Faskowitz, J., X. Yan, X.N. Zuo, and O. Sporns. 2018. Weighted stochastic block models of the human connectome across the life span. *Scientific Reports* 8.

Funke, T. and T. Becker. 2019. Stochastic block models: A comparison of variants and inference methods. *PLOS ONE* 14(4): e0215296. https://doi.org/10.1371/journal.pone.0215296.

Geman, S. and D. Geman. 1984. Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. *IEEE Transactions on Pattern Analysis and Machine Intelligence PAMI-6* (6): 721-741. https://doi.org/10.1109/TPAMI.1984.4767596.
From equations (1) and (2), it is equal to (with $Z_{i,q} = 1$ if $Z_i = q$ and zero otherwise)

$$
\begin{align*}
KL(Q(.) | \Pr_{\theta}(., D)) &= \sum_{i=1}^{n} \sum_{q=1}^{Q} \tau_{i,q} \ln(\tau_{i,q}) \\
&\quad - E_{Q} \left[ \ln \left( \frac{1}{W} \left( \prod_{i=1}^{n} \prod_{j>i} \Pr_{\Lambda}(D[i,j] | Z_i, Z_j) \right) \prod_{i=1}^{n} \Pr_{\alpha}(Z_i) \right) \right] \\
&= \sum_{i=1}^{n} \sum_{q=1}^{Q} \tau_{i,q} \ln(\tau_{i,q}) \\
&\quad - E_{Q} \left[ \ln \left( \frac{1}{W} \times \prod_{i=1}^{n} \left[ \prod_{j>i} \Pr_{\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}) \\
&\quad - \sum_{i=1}^{n} \sum_{j>i} \sum_{q=1}^{Q} \sum_{q'=1}^{Q} \tau_{i,q} \tau_{j,q'} \ln \Pr_{\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\left(\frac{1}{W}\right)
\end{align*}
$$

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

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

The variables are the $\tau_{i,q}$ for $l \in \{1, \dots, n\}$ and $q \in \{1, \dots, Q\}$, and the $\sigma_i$ for $l \in \{1, \dots, n\}$. Let us compute the derivatives of $G$.

$$\frac{\partial G}{\partial \sigma_l} = \sum_{q=1}^{Q} \tau_{l,q} - 1$$

751 So $\frac{G}{\partial \sigma_l} = 0$ leads to

$$\sum_{q=1}^{Q} \tau_{l,q} = 1 \qquad (13)$$
Gusfield, D. 1997. *Algorithms on strings, trees, and sequences*. Cambridge, UK: Cambridge University Press.

Harshman, R.A. 1970. Foundations of the PARAFAC procedure: Models and conditions for an explanatory multimodal factor analysis. *UCLA Working Papers in Phonetics* 16: 1-84.

Holland, P.W., K.B. Laskey, and S. Leinhardt. 1983. Stochastic blockmodels: First steps. *Social Networks* 5(2): 109-137. https://doi.org/10.1016/0378-8733(83)90021-7.

Horn, R.A. and C.R. Johnson. 2012. *Matrix Analysis* (Second ed.). Cambridge.

Häggström, O. 2002. *Finite Markov Chains and Algorithmic Applications*. Cambridge University Press. (London Mathematical Society Student Texts, Series Number 52).

Kalmbach, P., A. Blenk, M. Kluegel, and W. Kellerer 2017. Generating Synthetic Internet- and IP-Topologies using the Stochastic-Block-Model. In *2017 IFIP/IEEE Symposium on Integrated Network and Service Management*, Lisbon, Portugal.

Kolda, T.G. and B.W. Bader. 2009. Tensor decomposition and applications. *SIAM Review* 51(3): 455-500.

Koller, D. and N. Friedman. 2009. *Probabilistic Graphical Models - Principles and Techniques*. MIT Press.

Kroonenberg, P.M. and J. de Leeuw. 1980. Principal Component Analysis of three modes data by means of Alternating Least Square algorithms. *Psychometrika* 45(1): 69-97.

Lyu, Z., D. Xia, and Y. Zhang. 2023. Latent space model for higher-order networks and generalized tensor decomposition. *Journal of Computational and Graphical Statistics* 32(4): 1-17. https://doi.org/10.1080/10618600.2022.2164289.

Mariadassou, M., S. Robin, and C. Vacher. 2010. Uncovering latent structure in valued graphs: A variational approach. *The Annals of Applied Statistics* 4(2): 715-742. https://doi.org/10.1214/10-AOAS361.

Miele, V. and C. Matias. 2017. Revealing the hidden structure of dynamic ecological networks. *Royal Society Open Science* 4(6).

Murphy, K.P. 2012. *Machine Learning: A Probabilistic Perspective*. MIT Press.

Novikov, A., A. Rodomanov, A. Osokin, and D. Vetrov 2014. Putting MRFs on a Tensor Train. In *Proceedings of the 31st International Conference on Machine Learning, Beijing, China, 2014.*
Then using the fact that $D$ and $\Lambda$ are symmetric matrices, we obtain:

$$ 
\begin{aligned}
\frac{\partial G}{\partial \tau_{l,q}} &= \ln(\tau_{l,q}) + 1 + \sigma_l - \ln(\alpha_q) - \sum_{j=1}^{l-1} \sum_{q'=1}^{Q} \tau_{j,q'} [\ln \operatorname{Pr}_{\Lambda}(D[j,l] \mid Z_l = q, Z_j = q')] \\
&\quad - \sum_{j>l} \sum_{q'=1}^{Q} \tau_{j,q'} [\ln \operatorname{Pr}_{\Lambda}(D[l,j] \mid Z_l = q, Z_j = q')] \\
&= \ln(\tau_{l,q}) + 1 + \sigma_l - \ln(\alpha_q) - \sum_{j \neq l} \sum_{q'=1}^{Q} \tau_{j,q'} [\ln \operatorname{Pr}_{\Lambda}(D[l,j] \mid Z_l = q, Z_j = q')]
\end{aligned}
$$

So $\frac{\partial G}{\partial \tau_{l,q}} = 0$ leads to

$$ \tau_{l,q} = \alpha_q \exp \left( \sum_{j \neq l} \sum_{q'=1}^{Q} \tau_{j,q'} [\ln \operatorname{Pr}_{\Lambda}(D[l,j] \mid Z_l = q, Z_j = q')] \right) \exp(-(\sigma_l + 1) $$

We define

$$ \mu_{l,q} = \alpha_q \prod_{j \neq l} \prod_{q'=1}^{n} [\operatorname{Pr}_{\Lambda}(D[l,j] \mid Z_l = q, Z_j = q')]^{\tau_{j,q'}} $$

752 Then

$$ \tau_{l,q} = \mu_{l,q} e^{-(\sigma_l + 1)} \quad (14) $$

753 Finally, by combining equation (13) and equation (14) we obtain

$$ \tau_{l,q} = \frac{\mu_{l,q}}{\sum_{q'=1}^{Q} \mu_{l,q'}} $$

754 So the $\tau_{l,q}$ are solution of a fixed point equation, since $\mu_{l,q}$ is a function of the $\tau_{j,q'}$ for $j \neq l$.

The associated fixed point iteration scheme is

$$ \text{for all } 1 \le l \le n, \text{ for all } 1 \le q \le Q, \quad \left\{ \begin{array}{l} \tau_{l,q}^{t+1} = \frac{\mu_{l,q}^{t+1/2}}{\sum_{q'=1}^{Q} \mu_{l,q'}^{t+1/2}} \\ \mu_{l,q}^{t+1/2} = \alpha_q \prod_{j \neq l} \prod_{q'=1}^{Q} [\operatorname{Pr}_{\Lambda}(D[l,j] \mid Z_l = q, Z_j = q')]^{\tau_{j,q'}} \end{array} \right. $$

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

$$ \text{for all } 1 \le l \le n, \text{ for all } 1 \le q \le Q, \quad \left\{ \begin{array}{l} \tau_{l,q}^{t+1} = \frac{\mu_{l,q}^{t+1/2}}{\sum_{q'=1}^{Q} \mu_{l,q'}^{t+1/2}} \\ \mu_{l,q}^{t+1/2} = \alpha_q \prod_{j \neq l} \prod_{q'=1}^{Q} \left[ \frac{\lambda_{q,q'}^{D(l,j)}}{D(l,j)!} \exp(-\lambda_{q,q'}) \right]^{\tau_{j,q'}} \end{array} \right. $$
Orus, R. 2014. A Practical Introduction to Tensor Networks: Matrix Product States and Projected Entangled Pair States. *Annals of Physics* 349: 117–158.

Oseledets, I., E. Tyrtyshnikov, and N. Zamarashkin. 2011. Tensor-Train Ranks for Matrices and Their Inverses. *Computational Methods in Applied Mathematics* 11(3): 394–403.

Oseledets, I.V. 2009. A new tensor decomposition. *Doklady Mathematics* 80(1): 495–496.

Oseledets, I.V. 2011. Tensor-Train decomposition. *SIAM Journal of Scientific Computing* 33: 2295–2317.

Oseledets, I.V. and S.V. Dolgov. 2012. Solution of linear systems and matrix inversion in the tt-format. *SIAM Journal of Scientific Computing* 34(34): A2718–A2739.

Peyrard, N., M.J. Cros, S. de Givry, A. Franc, S. Robin, R. Sabbadin, T. Schiex, and M. Vignes. 2019. Exact or approximate inference in graphical models: why the choice is dictated by the treewidth, and how variable elimination can be exploited. *Australian and New Zealand Journal of Statistics* 61(2): 89–133. https://doi.org/ doi:10.1111/anzs.12257.

Robert, C. and G. Casella. 2004. *Monte Carlo statistical methods (second edition)*. Springer Verlag.

Savicky, P. and J. Vomlel. 2007. Exploiting tensor rank-one decomposition in probabilistic inference. *Kybernetica* 43(5): 747–764.

Smith, P.D. and M.S. Waterman. 1981. Identification of common molecular sub sequences. *Journal of Molecular Biology* 147: 195–197.

Strang, G. 2019. *Linear Algebra and Learning from Data*. Wellesley Cambridge Press.

Tucker, L.R. 1966. Some mathematical notes on three-modes factor analysis. *Psychometrika* 31(3): 279–311.

Vidal, G. 2003. Efficient Classical Simulation of Slightly Entangled Quantum Com putations. *Physical Review Letters* 91. https://doi.org/DOI:10.1103/PhysRevLett. 91.147902.

Wainwright, M.J. and M.I. Jordan. 2008. Graphical models, exponential families, and variational inference. *Foundations and Trends in Machine Learning* 1: 1–305.

Wrigley, A., W.S. Lee, and N. Ye 2017, 06–11 Aug. Tensor belief propagation. In D. Precup and Y. W. Teh (Eds.), *Proceedings of the 34th International Conference on Machine Learning*, Volume 70, pp. 3771–3779. PMLR.
## Appendix IV

<figure><img src="image_9.png" /><figcaption>Figure 4. Example of an histogram of the absolute difference between the value of $PSC_{i,j}$ computed with the exact method and with GS, for each pair $(i, j)$ of variables.</figcaption></figure>

## Appendix V

Here we show that, given two variables $Z_1$ and $Z_2$, taking value in $\{1, \dots, Q\}$, if there exists $q$ such that $Pr(Z_1 = q) = 1$ then, the binary marginals $Pr(Z_1, Z_2)$ are the products of the unary marginals. Let us consider $Q = 3$. Binary marginals can be organized as elements in a $3 \times 3$ matrix. The unary marginals are the marges of this matrix, i.e. the row-wise or column-wise sums. It can be displayed as

<table><tr><td>a</td><td>b</td><td>c</td><td>x</td></tr><tr><td>a'</td><td>b'</td><td>c'</td><td>x'</td></tr><tr><td>a''</td><td>b''</td><td>c''</td><td>x''</td></tr><tr><td>y</td><td>y'</td><td>y''</td><td>1</td></tr></table>
Computing SBM marginals with Tensor-Train
decomposition: supplementary material

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;
with

$$ egin{cases} a+b+c &= x \\ a'+b'+c' &= x' \\ a''+b''+c'' &= x'' \\ a+a'+a'' &= y \\ b+b'+b'' &= y' \\ c+c'+c'' &= y'' \end{cases} $$

764 All elements in this system are in $[0, 1]$.

765 Let us assume without loss of generality that $Pr(Z_1 = 1) = x = 1$. So, $x' = x'' = 0$.

766 From $a', b', c' \ge 0$ and $a' + b' + c' = x' = 0$, it follows that $a' = b' = c' = 0$, and the same for $a'', b'', c''$. So, the matrix is

$$ \begin{array}{c|ccc|c} a & b & c & 1 \\ \hline 0 & 0 & 0 & 0 \\ \hline 0 & 0 & 0 & 0 \\ \hline \end{array} \quad \begin{array}{c|ccc|c} y & y' & y'' & 1 \end{array} $$

769 with $a = y$, $b = y'$ and $c = y''$. We can observe that each term in the matrix is the product of the marginals of its row and column. It is not required that $y$, or $y'$ or $y''$ is 1 or 0. This demonstration can be extended to $Q > 3$.

772 References

773 ABOUABDALLAH, M.A., PEYRARD, N. & FRANC, A. (2022). Does clustering of DNA barcodes agree with botanical classification directly at high taxonomic levels? Trees in French Guiana as a case study. *Molecular Ecology Resources* **22**, 1746–1761. doi:https://doi.org/10.1111/1755-0998.13579.

776 AL DAAS, H., BALLARD, G. & BENNER, P. (2022). Parallel algorithms for tensor-train arithmetic. *SIAM Journal of Scientific Computing* **44**, C25–C53. doi:10.1137/20M1387158.

778 AL DAAS, H., BALLARD, G., CAZEAUX, P., HALLMAN, E., MIEDLAR, A., PASHA, M., REID, T.W. & SAIBABA, A.K. (2023). Randomized algorithms for rounding in the Tensor-Train format. *SIAM Journal of Scientific Computing* **45**, A74–195. doi:10.1137/21M1451191.

781 ANSARI, M., AHMAD, A., KHAN, S., BHUSHAN, G. & SIDDIQUE, M. (2020). Spatiotemporal clustering: a review. *Artificial Intelligence Review* **53**, 2381–2423. doi:10.1007/s10462-019-09736-1.

783 BARBILLON, P., DONNET, S., LAZEGA, E. & BAR-HEN, A. (2017). Stochastic block models for multiplex networks: An application to a multilevel network of researchers. *Journal of the Royal Statistical Society: Series A Statistics in Society* **180**, 295–314. doi:10.1111/rssa.12193. URL https://hal.archives-ouvertes.fr/hal-01520820.

787 CARROLL, J.D. & CHANG, J.J. (1970). Analysis of individual differences in multidimensional scaling via *n*-ways generalization of eckart-young decomposition. *Psychometrika* **35**, 283–319.

789 CASELLA, G. & GEORGE, E.I. (1992). Explaining the Gibbs Sampler. *The American Statistician* **46**, 167–174. URL http://www.jstor.org/stable/2685208.

791 COPPI, R. & BOLASCO, S. (1989). *Multiway data analysis*. Amsterdam: Elsevier.

792 DAUDIN, J.J., PICARD, F. & ROBIN, S. (2008). A mixture model for random graph. *Statistics and Computing* **18**, 173–183. doi:10.1007/s11222-007-9046-7.
<figure><img src="image_6.png" /><figcaption>**Fig. 1** Example of an histogram of the absolute difference between the value of $PSC_{i,j}$ computed with the exact method and with GS, for each pair $(i, j)$ of variables.</figcaption></figure>
794 DE LATHAUWER, L., DE MOOR, B. & VANDEWALLE, J. (2000). A Multilinear Singular Value Decomposition. *SIAM Journal on Matrix Analysis and Applications* **21**, 1253.

795 DE SILVA, V. & LIM, L.H. (2008). Tensor rank and the ill posedness of the best low rank approximation, problem. *SIAM Journal on Matrix Analysis and Applications* **30**, 1084–1027.

798 DEZA, M.M. & DEZA, E. (2016). *Encyclopedia of Distances*. Springer, 4th edn.

799 DUCAMP, G., BONNARD, P., NOUY, A. & WUILLEMIN, P.H. (2020). An efficient low-rank tensors representation for algorithms in complex probabilistic graphical models. In *Proceedings of the 10th International Conference on Probabilistic Graphical Models*. Skørping, Denmark.

802 FANNES, M., NACHTERGAELE, B. & WERNER, R.F. (1992). Finitely Correlated States on Quantum Spin Chains. *Communications in Mathematical Physics* **144**, 443–490.

804 FASKOWITZ, J., YAN, X., ZUO, X.N. & SPORNS, O. (2018). Weighted stochastic block models of the human connectome across the life span. *Scientific Reports* **8**.

806 FUNKE, T. & BECKER, T. (2019). Stochastic block models: A comparison of variants and inference methods. *PLOS ONE* **14**, e0215296. doi:10.1371/journal.pone.0215296. URL https://dx.plos.org/10.1371/journal.pone.0215296.

809 GEMAN, S. & GEMAN, D. (1984). Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. *IEEE Transactions on Pattern Analysis and Machine Intelligence* **PAMI-6**, 721–741. doi: 10.1109/TPAMI.1984.4767596.

812 GUSFIELD, D. (1997). *Algorithms on strings, trees, and sequences*. Cambridge, UK: Cambridge University Press.

814 HARSHMAN, R.A. (1970). Foundations of the PARAFAC procedure: Models and conditions for an explanatory multimodal factor analysis. *UCLA Working Papers in Phonetics* **16**, 1–84.

816 HOLLAND, P.W., LASKEY, K.B. & LEINHARDT, S. (1983). Stochastic blockmodels: First steps. *Social Networks* **5**, 109–137. doi:https://doi.org/10.1016/0378-8733(83)90021-7. URL https://www.sciencedirect.com/science/article/pii/0378873383900217.

819 HORN, R.A. & JOHNSON, C.R. (2012). *Matrix Analysis*. Cambridge, 2nd edn.

820 HÄGGSTRÖM, O. (2002). *Finite Markov Chains and Algorithmic Applications*. Cambridge University Press. (London Mathematical Society Student Texts, Series Number 52).

822 KALMBACH, P., BLENK, A., KLUEGEL, M. & KELLERER, W. (2017). Generating Synthetic Internet-and IP-Topologies using the Stochastic-Block-Model. In *2017 IFIP/IEEE Symposium on Integrated Network and Service Management*. Lisbon, Portugal.

825 KOLDA, T.G. & BADER, B.W. (2009). Tensor decomposition and applications. *SIAM Review* **51**, 455–500.

826 KOLLER, D. & FRIEDMAN, N. (2009). *Probabilistic Graphical Models - Principles and Techniques*. MIT Press.

828 KROONENBERG, P.M. & DE LEEUW, J. (1980). Principal Component Analysis of three modes data by means of Alternating Least Square algorithms. *Psychometrika* **45**, 69–97.

830 LYU, Z., XIA, D. & ZHANG, Y. (2023). Latent space model for higher-order networks and generalized tensor decomposition. *Journal of Computational and Graphical Statistics* **32**, 1–17. doi:10.1080/10618600.2022.2164289.

833 MIELE, V. & MATIAS, C. (2017). Revealing the hidden structure of dynamic ecological networks. *Royal Society Open Science* **4**.

835 MURPHY, K.P. (2012). *Machine Learning: A Probabilistic Perspective*. MIT Press.

836 NOVIKOV, A., RODOMANOV, A., OSOKIN, A. & VETROV, D. (2014). Putting MRFs on a Tensor Train. In *Proceedings of the 31st International Conference on Machine Learning, Beijing, China, 2014*.

838 ORUS, R. (2014). A Practical Introduction to Tensor Networks: Matrix Product States and Projected Entangled Pair States. *Annals of Physics* **349**, 117–158.

840 OSELEDETS, I., TYRTYSHNIKOV, E. & ZAMARASHKIN, N. (2011). Tensor-Train Ranks for Matrices and Their Inverses. *Computational Methods in Applied Mathematics* **11**, 394–403.

842 OSELEDETS, I.V. (2009). A new tensor decomposition. *Doklady Mathematics* **80**, 495–496.
843 OSELEDETS, I.V. (2011). Tensor-Train decomposition. *SIAM Journal of Scientific Computing* **33**, 2295-2317.

845 OSELEDETS, I.V. & DOLGOV, S.V. (2012). Solution of linear systems and matrix inversion in the TT-format. *SIAM Journal of Scientific Computing* **34**, A2718-A2739.

847 PEYRARD, N., CROS, M.J., DE GIVRY, S., FRANC, A., ROBIN, S., SABBADIN, R., SCHIEX, T. & VIGNES, M. (2019). Exact or approximate inference in graphical models: why the choice is dictated by the treewidth, and how variable elimination can be exploited. *Australian and New Zealand Journal of Statistics* **61**, 89-133. doi:doi:10.1111/anzs.12257.

851 ROBERT, C. & CASELLA, G. (2004). *Monte Carlo statistical methods (second edition)*. Springer Verlag.

852 SAVICKY, P. & VOMLEL, J. (2007). Exploiting tensor rank-one decomposition in probabilistic inference. *Kybernetica* **43**, 747-764.

854 SMITH, P.D. & WATERMAN, M.S. (1981). Identification of common molecular subsequences. *Journal of Molecular Biology* **147**, 195-197.

856 STRANG, G. (2019). *Linear Algebra and Learning from Data*. Wellesley Cambridge Press.

857 TUCKER, L.R. (1966). Some mathematical notes on three-modes factor analysis. *Psychometrika* **31**, 279-311.

859 VIDAL, G. (2003). Efficient Classical Simulation of Slightly Entangled Quantum Computations. *Physical Review Letters* **91**, doi:DOI:10.1103/PhysRevLett.91.147902.

861 WAINWRIGHT, M.J. & JORDAN, M.I. (2008). Graphical models, exponential families, and variational inference. *Foundations and Trends in Machine Learning* **1**, 1-305.

863 WRIGLEY, A., LEE, W.S. & YE, N. (2017). Tensor belief propagation. In *Proceedings of the 34th International Conference on Machine Learning*, eds. D. Precup & Y.W. Teh, vol. 70. PMLR, pp. 3771-3779. URL https://proceedings.mlr.press/v70/wrigley17a.html.

866 YEDIDIA, J., FREEMAN, W. & WEISS, Y. (2005). Constructing free energy approximations and generalized belief propagation algorithms. *IEEE Transactions on Information Theory* **51**, 2282-2312.