## 6.1 Introduction

L'objectif du chapitre est d'évaluer la qualité du calcul des marginales unaires et binaires par l'approche Tenseur Train (TT) présentée dans le chapitre 5. Pour cela, nous la comparons avec trois autres méthodes présentes dans la littérature.

* Premièrement, quand $n$ est petit, nous la comparons avec la méthode de calcul exacte par énumération. Cette dernière consiste à construire le tenseur de la loi jointe puis faire le calcul des marginales par sommation.
Cette méthode est exacte mais, elle nécessite le stockage d'un tenseur dont un élément est la loi jointe de $Z$ sachant $D$ pour un état de $Z$: ce tenseur contient $Q^n$ éléments. De plus, le nombre de sommes est exponentiel avec l'ordre du tenseur.

* Ensuite, nous comparons notre approche à une méthode de type Méthode de Monte Carlo par chaîne de Markov [And01] (MCMC) qui sont des méthodes réputées pour leur forte précision mais, qui restent très lentes à mettre en œuvre. Parmi ces méthodes, nous avons choisi l'échantillonneur de Gibbs [RC04], qui est simple à mettre en œuvre. Cette méthode est utilisée dans l'algorithme MCEM.

* Enfin, nous comparons notre approche à une méthode variationnelle basée sur l'approximation en champ moyen [DPR08] de la loi de $Z$ sachant $D$: c'est une méthode basée sur une hypothèse d'indépendance. Donc, les marginales binaires sont obtenues comme le produit des unaires. Cette méthode est très rapide vu qu'elle repose sur la résolution d'une équation du point fixe mais, peut être peu précise certaines fois. Cette méthode est utilisée dans l'algorithme VEM pour l'estimation des paramètres d'un SBM.

Dans le graphique de la figure 6.1, nous présentons le comportement de chaque méthode en termes de précision et de temps de calcul.

<figure><img src="image_20.png" /><figcaption>FIGURE 6.1 – Compromis précision/temps de calcul pour les quatre méthodes testées pour l'inférence de la loi de Z sachant D</figcaption></figure>
Afin d'évaluer l'approche TT, nous avons utilisé des matrices de distances simulées à partir de six structures de distances que nous présenterons ultérieurement. Cela a pour but de nous permettre de valider l'approche TT de manière générique.

Ensuite, pour chaque structure :

1. nous avons calculé les marginales unaires et binaires pour les quatre méthodes de calculs (la méthode exacte par énumération, l'échantillonneur de Gibbs, TT et le Champ moyen) pour un nombre d'individus $n = 9$ (qui est la limite à laquelle la méthode du calcul exacte par énumération peut être développée). Puis nous les avons comparé entre eux.

2. nous avons calculé les marginales unaires et binaires pour trois méthodes de calculs (l'échantillonneur de Gibbs, TT et le Champ moyen) pour $n = 18$ : la méthode de référence est maintenant l'échantillonneur de Gibbs à laquelle nous avons comparé TT et le champ moyen.

Dans ce qui suit, nous présenterons les quatre méthodes pour le calcul des marginales dans la sous-section 6.2, les structures de distances simulées dans la sous-section 6.3, le protocole pour le calcul des marginales sur les six structures dans la sous-section 6.4, les résultats numériques dans les sous-sections 6.5 et 6.6 et enfin une discussion dans la sous-section 6.7.

## 6.2 Méthodes de calcul des marginales

Nous commençons par présenter l'expression de la loi jointe dont nous voulons calculer les marginales unaires et binaires, puis nous présentons chacune des trois méthodes autres que "Tenseur Train" pour le calcul des marginales (la méthode via les TT est présentée dans le chapitre 5).

### 6.2.1 Expression de la loi jointe

La loi jointe de $Z_1, \dots, Z_n$ conditionnellement aux données $D$ s'exprime en fonction des facteurs unaires (2.2.18) et binaires (2.2.19) présentés dans la section 2.2.1.4 ainsi :

$$ \mathbb{P}_{\theta}(Z_1, \dots, Z_n | D) = \frac{1}{W} \left[ \prod_{i=1}^{n} \prod_{j>i} \psi_{i,j}(Z_i, Z_j) \right] \left[ \prod_{i=1}^{n} \phi_i(Z_i) \right] \quad (6.2.1) $$

Avec,

$$ \psi_{ij} \in \mathbb{R}^{Q \times Q}, \quad \phi_i \in \mathbb{R}^Q \quad Z_i, Z_j \in [1, Q] \quad (6.2.2) $$

**Facteurs** Nous exprimons les facteurs binaires et unaires de la loi de $Z$ sachant $D$ :

$$ \psi_{i,j}(z_i, z_j) = \sum_{q=1}^{Q} \sum_{q'=1}^{Q} \frac{\lambda_{q,q'}^{D(i,j)}}{D(i,j)!} \exp(-\lambda_{q,q'}) z_{i,q} \times z_{j,q'} \quad (6.2.3) $$

$$ \phi_i(z_i) = \prod_{q=1}^{Q} \alpha_q^{z_{i,q}} \quad (6.2.4) $$

Pour l'implémentation des calculs, il peut être utile d'assembler les facteurs binaires et unaires dans des nouveaux facteurs binaires de cette manière :

$$ \left| \begin{array}{ll} \psi'_{i,i+1}(z_i, z_{i+1}) & = \psi_{i,i+1}(z_i, z_{i+1})\phi_i(z_i) \quad \forall 1 \le i \le n-2 \\ \psi'_{i,j}(z_i, z_j) & = \psi_{i,j}(z_i, z_j) \quad \text{si } j \ne i+1 \\ \psi'_{n-1,n}(z_{n-1}, z_n) & = \psi_{n-1,n}(z_{n-1}, z_n)\phi_{n-1}(z_{n-1})\phi_n(z_n) \end{array} \right. \quad (6.2.5) $$
## 6.2.2 Méthode exacte par énumération

Le calcul par la méthode exacte par énumération se fait de la manière suivante :

### Construction du tenseur de la loi jointe

* Nous créons un tenseur $\Psi \in (\mathbb{R}^Q)^{\otimes n}$ d'ordre $n$ où nous stockons la loi jointe.

* Puis, $\forall z_1, z_2, \dots, z_n = 1, \dots, Q$

$$ \Psi[z_1, z_2, \dots, z_n] = \frac{1}{W} \left[ \prod_{i=1}^{n} \prod_{j>i} \psi_{i,j}(z_i, z_j) \right] \left[ \prod_{i=1}^{n} \phi_i(z_i) \right] \quad (6.2.6) $$

Cette méthode consiste à stocker le tenseur complet de la loi jointe dans un ordinateur.

### Calcul des marginales
Dans ce cas, le calcul des marginales découle de la définition :

* Marginales unaires : $\forall i \in [1, n], q \in [1, Q]$

$$ p_i(Z_i = q) = \sum_{z_1, \dots, z_{i-1}, z_{i+1}, \dots, z_n \in [1, Q]^{[n-1]}} \Psi[z_1, \dots, z_{i-1}, z_i = q, z_{i+1}, \dots, z_n] \quad (6.2.7) $$

* Marginales binaires : $\forall i, j \in [1, n]^2, q, q' \in [1, Q]^2$

$$ p_{i,j}(Z_i = q, Z_j = q') = \sum_{z_l \in [1, Q], \forall l \in [1, n] \setminus \{i, j\}} \Psi[z_1, \dots, z_{i-1}, z_i = q, z_{i+1}, \dots, z_{j-1}, z_j = q', z_{j+1}, \dots, z_n] \quad (6.2.8) $$

L'avantage de cette méthode est qu'on ne fait pas d'approximation, le calcul est exact. L'inconvénient est qu'elle nécessite le stockage de $Q^n$ éléments, ce qui pose une limite de stockage. De plus, cette méthode nécessite un nombre exponentiel d'additions.

## 6.2.3 Échantillonneur de Gibbs

Dans cette section, nous commençons par expliquer le principe de l'échantillonneur de Gibbs, sa mise en œuvre et le pré-traitement nécessaire pour son implémentation.

### 6.2.3.1 Principe de l'échantillonneur de Gibbs pour le calcul des marginales

L'échantillonneur de Gibbs [RC04] est un algorithme général de type Monte-Carlo par chaînes de Markov [And01] (MCMC) qui permet de générer des séquences d'observations très proches d'une réalisation d'une loi multivariée difficile à échantillonner directement. Quand la simulation est répétée un grand nombre de fois, cela permet une estimation des marginales.

Le principe est que, étant donné une loi multivariée, il est facile d'échantillonner à partir d'une loi conditionnelle (une variable conditionnellement à toutes les autres). Afin de produire un échantillon indexé par exemple par $\ell \in [1, L]$, soit $z^{(\ell)} = (z_1^{(\ell)}, \dots, z_n^{(\ell)})$, nous partons d'un $z^{(\ell)}$ initial. Ensuite, nous mettons à jour successivement sur $i \in [1, n]$ et dans l'ordre chacune des composantes $z_i^{(\ell)}$, les autres étant fixées.

Nous répétons en boucle $N_{\text{iter}}$ fois ce processus. Lorsque $N_{\text{iter}} \to \infty$, alors le $z^{(\ell)}$ obtenu est une réalisation de la loi $P_{\theta}(Z | D)$, loi jointe de Z sachant D. L'idée pour le calcul des marginales est de produire un grand nombre $L$ d'échantillons qui sont chacun une réalisation approchée de la loi jointe à échantillonner, ici $P_{\theta}(Z | D)$. Des statistiques élémentaires sur les $L$ réalisations permettent alors de calculer les marginales souhaitées selon cette formule :
$$P_{\theta}(Z_{i,q} = 1 | D) \approx \frac{1}{L} \sum_{\ell=1}^{L} 1(z_{i,q}^{(\ell)} = 1) \qquad (6.2.9)$$

$$P_{\theta}(Z_{i,q} = 1, Z_{j,q'} = 1 | D) \approx \frac{1}{L} \sum_{\ell=1}^{L} 1(z_{i,q}^{(\ell)} = 1, z_{j,q'}^{(\ell)} = 1) \qquad (6.2.10)$$

avec,

$$\begin{cases} z_{i,q}^{(\ell)} = 1 (i \in q) \\ \sum_{q=1}^{Q} z_{i,q}^{(\ell)} = 1 \end{cases}$$

Il y a donc deux paramètres qui contrôlent la qualité du résultat :

1. $N_{\text{iter}}$ qui contrôle la qualité de chaque réalisation de la loi jointe $P_{\theta}(Z | D)$.

2. $L$ qui contrôle la qualité de l'estimation des marginales.

### 6.2.3.2 Préparation

**Loi conditionnelle** : L'approche via l'échantillonneur de Gibbs repose sur la mise à jour séquentielle de chacun des $z_i$ connaissant sa loi conditionnellement à tous les autres $z_j$. Cette loi s'écrit, pour $z_i$ :

$$P_{\theta}(Z_i = z_i | \{Z_j = z_j\}_{j \neq i}, D) = \frac{\phi_i(z_i) \prod_{j \neq i} \psi_{i,j}(z_i, z_j)}{\sum_{u=1}^{Q} \phi_i(u) \prod_{j \neq i} \psi_{i,j}(u, z_j)} \qquad (6.2.11)$$

Cela suppose de calculer et de stocker les facteurs unaires $\phi_i(z_i)$ et binaires $\psi_{i,j}(z_i, z_j)$.

**Initialisation** : Nous souhaitons simuler $L$ échantillons.

Pour chaque échantillon $l \in [1, L]$, nous prenons un point de départ $z^{(\ell)} = (z_1^{(\ell)}, \dots, z_n^{(\ell)})$ où nous générons la classe initiale $z_i^{(\ell)}$ de $i \in [1, n]$ selon une loi uniforme.

### 6.2.3.3 Mise en œuvre

L'échantillonneur de Gibbs se fait en deux étapes : une initialisation dans laquelle nous simulons $z^{(\ell)}$ initial et une étape de $N_{\text{iter}}$ mises à jour successives de chacune des coordonnées de $z^{(\ell)}$. Le but est d'aller jusqu'à $N_{\text{iter}}$ mises à jour de $z^{(\ell)}$ pour que le dernier échantillon produit soit une réalisation fiable de la loi $P_{\theta}(Z | D)$.

La procédure se déroule ainsi :

* **Initialisation** : Choisir une valeur initiale $z^{\text{init}} = (z_1 \dots z_n)$.

Nous attribuons :

$$z^{\text{old}} \leftarrow z^{\text{init}}$$

* **Mises à jour successives** : Supposons que l'on commence l'itération $t$. Nous notons $z^t = z^{\text{old}}$ et $z^{t+1} = z^{\text{new}}$.

Une itération de Gibbs consiste à prendre $z^{\text{old}}$ en entrée et mettre à jour successivement les états $z_i$ des $n$ nœuds pour retourner $z^{\text{new}}$ quand l'ensemble des nœuds ont été mis à jour. La séquence des opérations successives se déroule ainsi :

— pour $i=1$, nous simulons $z_1^{\text{new}}$ selon la loi définie par $P_{\theta}(Z_1 = z_1 | \{Z_j^{\text{old}} = z^{\text{old}}\}_{j \neq 1}, D)$

— pour $i=2$, nous simulons $z_2^{\text{new}}$ selon la loi définie par $P_{\theta}(Z_2 = z_2 | Z_1 = z_1^{\text{new}}, Z_3 = z_3^{\text{old}}, \dots, Z_n = z_n^{\text{old}}, D)$

∴
— pour $i \in [1, n]$, nous simulons $z_i^{new}$ selon la loi définie par $P_{\theta}(Z_i = z_i | Z_1 = z_1^{new}, \dots, Z_{i-1} = z_{i-1}^{new}, Z_{i+1} = z_{i+1}^{old}, \dots, Z_n = z_n^{old}, D)$

— pour $i=n$, nous simulons $z_n^{new}$ selon la loi définie par $P_{\theta}(Z_n = z_n | Z_1 = z_1^{new}, \dots, Z_{n-1} = z_{n-1}^{new}, D)$

— $z^{old} \leftarrow z^{new}$

* Nous refaisons cette simulation de manière séquentielle $N_{iter}$ fois

* Le dernier $z^{old}$ obtenu est la réalisation de la loi $P_{\theta}(Z|D)$.

Nous lançons $L$ échantillonneurs de Gibbs et nous obtenons $L$ réalisations de $Z: z^{(1)}, z^{(2)}, \dots, z^{(l)}$ à partir desquels nous calculons les marginales unaires et binaires selon les formules (6.2.9) et (6.2.10).

## 6.2.4 Champ moyen

L'approximation basée sur l'approximation en champ moyen de la loi de $Z$ sachant $D$ est une méthode variationnelle développée par Daudin & al. en 2008 [DPR08] dans le cadre de l'algorithme VEM pour l'estimation des paramètres d'un SBM.

### 6.2.4.1 Principe du champ moyen

Dans le cadre de cette approche, nous souhaitons approcher $P_{\theta}(Z | D)$ par une distribution $Q_{\theta}(Z | D)$ plus simple à calculer. Celle-ci doit vérifier les propriétés suivantes :

* $P_1: Q_{\theta}(Z | D) = \prod_{i=1}^{n} q_{\theta}^{i}(Z_i | D)$

* $P_2: Q_{\theta}$ minimise la divergence de Kullback-Leibler.

**Définition 1.** Divergence de Kullback-Leibler [AS66] entre une distribution $Q$ et une distribution $P$ s'exprime :

$$KL(Q|P) = \mathbb{E}_{Q} \left[ \log \left( \frac{Q(Z)}{P(Z)} \right) \right] \quad (6.2.12)$$

$$KL(Q|P) = \sum_{z} Q(z) \log \left( \frac{Q(z)}{P(z)} \right) = J(Q(Z)) - \mathbb{E}_{Q}[\log(P(Z))] \quad (6.2.13)$$

avec $J(Q(Z)) = \mathbb{E}_{Q}[\log(Q(Z | D))]$.

### 6.2.4.2 Mise en œuvre du champ moyen

Pour cela, nous choisisons la distribution $Q_{\theta}$ de cette manière :

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

Dans le cas optimal, sans contrainte sur la forme de $Q_{\theta}$, nous aurons $P_{\theta}(Z | D) = Q_{\theta}(Z | D)$ et donc une divergence de Kullback-Leibler nulle.

Rappelons la propriété $P_1$ :

$$Q_{\theta}(Z|D) = \prod_{i=1}^{n} q_{\theta}^{i}(Z_i | D)$$

On notera $\tau_{i,q} = q_{\theta}^{i}(q | D)$. Donc :

$$\forall i \in [1, n], \quad \sum_{q=1}^{Q} \tau_{i,q} - 1 = 0 \quad (6.2.14)$$
On développe $J(Q_\theta(Z))$ :

$$
\begin{align*}
J(Q_{\theta}(Z)) &= \mathbb{E}_{Q_{\theta}}[\log(Q_{\theta}(Z | D))] = \mathbb{E}_{Q_{\theta}}\left[\sum_{i=n}^{n} \log(q_{\theta}^{i}(z_i | D))\right] \\
&= \sum_{i=1}^{n} \mathbb{E}_{Q_{\theta}}[\log(q_{\theta}^{i}(z_i | D))] = \sum_{i=1}^{n} \sum_{q=1}^{Q} \tau_{i,q} \log(\tau_{i,q})
\end{align*}
$$

Ensuite, nous développons l'espérance de la loi jointe du modèle :

$$
\mathbb{E}_{Q_{\theta}}[\log(\mathbb{P}_{\theta}(Z | D))] = \mathbb{E}_{Q_{\theta}}\left[\log\left(\frac{1}{W} \times \prod_{i=1}^{n} \left[\prod_{j>i} \frac{\lambda_{z_i, z_j}^{D_{i,j}}}{D_{i,j}!} e^{-\lambda_{z_i, z_j}} \prod_{q=1}^{Q} \alpha_q^{z_i, q}\right]\right)\right] \quad (6.2.15)
$$

Nous appliquons la propriété suivante $\log(a \times b) = \log(a) + \log(b)$ sur la formule (6.2.15) et nous
développons :

$$
\mathbb{E}_{Q_{\theta}}[\log(\mathbb{P}_{\theta}(Z | D))] = \mathbb{E}_{Q_{\theta}}(\log(\frac{1}{W})) + \sum_{i=1}^{n} \sum_{j>i} \mathbb{E}_{Q_{\theta}}[\log(\frac{\lambda_{z_i, z_j} e^{-\lambda_{z_i, z_j} D_{ij}}}{D_{i,j}!})] + \sum_{q=1}^{Q} \sum_{i=1}^{n} \mathbb{E}_{Q_{\theta}}[z_{i,q} \log(\alpha_q)] \tag{6.2.16}
$$

Comme $W$ est une constante, alors : $\mathbb{E}_{Q_{\theta}}[\log(\frac{1}{W})] = \log(\frac{1}{W})$.
Ainsi, nous développons la formule de l'espérance pour obtenir :

$$
\mathbb{E}_{Q_{\theta}}[\log(\mathbb{P}_{\theta}(Z | D))] = \log\left(\frac{1}{W}\right) + \sum_{i=1}^{n} \sum_{j>i} \sum_{q=1}^{Q} \sum_{q'=1}^{Q} \tau_{i,q} \tau_{j,q'} \log\left(\frac{\lambda_{q,q'}^{D_{i,j}}}{D_{i,j}!} e^{-\lambda_{q,q'}}\right) + \sum_{q=1}^{Q} \sum_{i=1}^{n} [\tau_{i,q} \log(\alpha_q)] \tag{6.2.17}
$$

Ainsi,

$$
\begin{equation}
\begin{split}
KL(Q_{\theta} | \mathbb{P}) = \sum_{i=1}^{n} \sum_{q=1}^{Q} \tau_{i,q} \log(\tau_{i,q}) & - \log\left(\frac{1}{W}\right) - \sum_{q=1}^{Q} \sum_{i=1}^{n} [\tau_{i,q} \log(\alpha_q)] \\
& - \sum_{i=1}^{n} \sum_{j>i} \sum_{q=1}^{Q} \sum_{q'=1}^{Q} \tau_{i,q} \tau_{j,q'} \log\left(\frac{\lambda_{q,q'}^{D_{i,j}}}{D_{i,j}!} e^{-\lambda_{q,q'}}\right)
\end{split}
\tag{6.2.18}
\end{equation}
$$

Nous souhaitons minimiser $KL(Q_{\theta} | \mathbb{P})$ en respectant la contrainte de la formule (6.2.14). Alors, nous allons minimiser $g(\sigma_{l}, \tau_{l,q}) = KL(Q_{\theta} | \mathbb{P}) + \sum_{l=1}^{n} \sigma_{l}(\sum_{q=1}^{Q} \tau_{l,q} - 1)$ par rapport à $\sigma_{l}$ et $\tau_{l,q}$ (avec $\sigma_{l}$ le multiplicateur de Lagrange).

Afin de maximiser cette fonction, il faudrait trouver le zéro de sa dérivée par rapport à $\sigma_l$ et $\tau_{l,q}$.

Nous commençons par calculer la dérivée selon $\sigma_l$ :

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

$$
\text{Comme } \frac{\partial g(\sigma_l, \tau_{l,q})}{\partial \sigma_l} = 0, \text{ alors}
$$
$$
\sum_{q=1}^{Q} \tau_{l,q} = 1 \tag{6.2.19}
$$
Comme les matrices $D$ et $\Lambda$ sont symétriques, la dérivée de $g(\sigma_l, \tau_{l,q})$ selon $\tau_{l,q}$ s'exprime :

$$ \frac{\partial g(\sigma_l, \tau_{l,q})}{\partial \tau_{l,q}} = \log(\tau_{l,q}) + 1 + \sigma_l - \log(\alpha_q) - \sum_{j=1}^{l-1} \sum_{q'=1}^{Q} \tau_{j,q'} \left[ \log\left(\frac{\lambda_{q,q'}^{D_{l,j}}}{D_{l,j}!} e^{-\lambda_{q,q'}}\right) \right] \\ - \sum_{j>l} \sum_{q'=1}^{Q} \tau_{j,q'} \left[ \log\left(\frac{\lambda_{q,q'}^{D_{l,j}}}{D_{l,j}!} e^{-\lambda_{q,q'}}\right) \right] \quad (6.2.20) $$

Or,

$$ \sum_{j=1}^{l-1} \sum_{q'=1}^{Q} \tau_{j,q'} \left[ \log \left( \frac{\lambda_{q,q'}^{D_{l,j}}}{D_{l,j}!} e^{-\lambda_{q,q'}} \right) \right] + \sum_{j>l} \sum_{q'=1}^{Q} \tau_{j,q'} \left[ \log \left( \frac{\lambda_{q,q'}^{D_{l,j}}}{D_{l,j}!} e^{-\lambda_{q,q'}} \right) \right] = \sum_{j \neq l} \sum_{q'=1}^{Q} \tau_{j,q'} \left[ \log \left( \frac{\lambda_{q,q'}^{D_{l,j}}}{D_{l,j}!} e^{-\lambda_{q,q'}} \right) \right] $$

Ainsi,

$$ \frac{\partial g(\sigma_l, \tau_{l,q})}{\partial \tau_{l,q}} = \log(\tau_{l,q}) + 1 + \sigma_l - \log(\alpha_q) - \sum_{j \neq l} \sum_{q'=1}^{Q} \tau_{j,q'} \left[ \log \left( \frac{\lambda_{q,q'}^{D_{l,j}}}{D_{l,j}!} e^{-\lambda_{q,q'}} \right) \right] $$

On cherche :

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

Donc :

$$ \tau_{l,q} = \alpha_q e^{\sum_{j \neq l} \sum_{q'=1}^{Q} \tau_{j,q'} \left[ \log \left( \frac{\lambda_{q,q'}^{D_{l,j}}}{D_{l,j}!} e^{-\lambda_{q,q'}} \right) \right] e^{-\sigma_l - 1}} \quad (6.2.21) $$

Ensuite, nous utilisons la propriété du log :

$$ \sum_{j \neq l} \sum_{q'=1}^{Q} \tau_{j,q'} \left[ \log \left( \frac{\lambda_{q,q'}^{D_{l,j}}}{D_{l,j}!} e^{-\lambda_{q,q'}} \right) \right] = \log \left( \prod_{j \neq l} \prod_{q'=1}^{Q} \left[ \frac{\lambda_{q,q'}^{D_{l,j}}}{D_{l,j}!} e^{-\lambda_{q,q'}} \right]^{\tau_{j,q'}} \right) $$

On pose :

$$ \mu_{l,q} = \alpha_q \prod_{j \neq l} \prod_{q'=1}^{n} \left[ \frac{\lambda_{q,q'}^{D_{l,j}}}{D_{l,j}!} e^{-\lambda_{q,q'}} \right]^{\tau_{j,q'}} \quad (6.2.22) $$

Ainsi $\tau_{l,q}$ est proportionnel à $\mu_{l,q}$ et s'exprime :

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

On combine les deux formules (6.2.19) et (6.2.24) pour obtenir :

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

Une manière de trouver le zéro de cette fonction est de procéder par un schéma du point fixe.
# Schéma du point fixe et dynamique des quantités d'intérêt

L'équation de point fixe est la suivante :

$$ \tau_{i,q}^{t+1} = \frac{\mu_{i,q}^{t+1/2}}{\sum_{l=1}^{Q} \mu_{i,l}^{t+1/2}} \qquad (6.2.25) $$

La dynamique des quantités d'intérêt s'exprime :

$$ \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. \qquad (6.2.26) $$

On initialise le schéma du point fixe de cette manière :

$$ \mu_{i,q}^{0} = \alpha_{q} $$

**Critère d'arrêt :** Soit,

$$ \forall 1 \le i \le n, \forall 1 \le q \le Q, \quad \Delta_{i,q} = \frac{|\mu_{i,q}^{t+1} - \mu_{i,q}^{t}|}{|\mu_{i,q}^{t}|} $$

Nous choisissons le critère d'arrêt suivant :

$$ \|\Delta\|_{\infty} < \epsilon_{pf} $$

Après la convergence de ce schéma, nous obtenons les valeurs de $\tau_{i,q}^*$ et les marginales unaires s'expriment :

$$ q_{\theta}^{i}(z_{i} = q) = \tau_{i,q}^{*} \qquad (6.2.27) $$

Enfin, comme pour cette approche où nous faisons une approximation champ moyen, la loi jointe est approchée par le produit des marginales unaires, le calcul des marginales binaires s'obtient comme produit de marginales unaires :

$$ \forall i,j \in [1,n]^2, q_{\theta}^{i,j}(z_i, z_j) = q_{\theta}^{i}(z_i) q_{\theta}^{j}(z_j) \qquad (6.2.28) $$
## 6.3 Les six structures de réseaux à utiliser pour générer les graphes

Il existe plusieurs structures de SBM qui permettent de modéliser des systèmes assez différents comme la détection de communautés dans les réseaux sociaux et de composants fonctionnels dans les réseaux d'interaction protéine-protéine qui ont été illustrées par exemple dans l' article [KBKK]. En nous basant sur la littérature [KBKK, FYZS18, FB19, BE99], nous avons choisi les six structures suivantes pour simuler nos matrices des distances car elles correspondent à des structures de lambda bien contrastées. Nous présentons ci dessous les structures de distances utilisées et nous illustrons les matrices de distances simulée pour un exemple $D \in \mathbb{R}^{100 \times 100}$.

**La structure associative** [KBKK] les distances classes sont faibles tandis que les distances interclasses sont fortes.

**Structure $\Lambda_{asso}$**

$$ \Lambda_{asso} = \begin{pmatrix} 2 & 10 & 10 \\ 10 & 3 & 10 \\ 10 & 10 & 4 \end{pmatrix} $$

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

**Structure dissociative** [KBKK] les distances intra classes sont fortes tandis que les distances interclasses sont faibles.

**Structure $\Lambda_{diss}$**

$$ \Lambda_{diss} = \begin{pmatrix} 8 & 3 & 3 \\ 3 & 10 & 3 \\ 3 & 3 & 9 \end{pmatrix} $$

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

**Structure hiérarchique** [KBKK] où les distances intra classes sont plus fortes que les distances interclasses. De plus, les distances sont organisées de manière hiérarchique. C'est-à-dire :
Les distances intra classes vérifient :

$$ \Lambda_{\text{hier}}[1, 1] > \Lambda_{\text{hier}}[2, 2] > \Lambda_{\text{hier}}[3, 3] $$

Les distances interclasses vérifient :

$$ \Lambda_{\text{hier}}[1, 2] \geq \Lambda_{\text{hier}}[2, 3] \geq \Lambda_{\text{hier}}[1, 3] $$

**Structure $\Lambda_{\text{hier}}$**

$$ \Lambda_{\text{hier}} = \begin{pmatrix} 12 & 3 & 1 \\ 3 & 6 & 3 \\ 1 & 3 & 3 \end{pmatrix} $$

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

**Ordonné** [KBKK] Les distances intra classes sont beaucoup plus fortes que les distances interclasses.

De plus, les distances interclasses sont organisés de manière décroissante et vérifient :

$$ \Lambda_{\text{hier}}[1, 2] \geq \Lambda_{\text{hier}}[2, 3] \geq \Lambda_{\text{hier}}[1, 3] $$

**Structure $\Lambda_{\text{ord}}$**

$$ \Lambda_{\text{ord}} = \begin{pmatrix} 10 & 6 & 2 \\ 6 & 10 & 6 \\ 2 & 6 & 10 \end{pmatrix} $$

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

**Cores periphery** [BE99, FYZS18, FB19], cette structure, suppose qu'il existe deux classes de nœuds, ceux du centre et ceux de la péripherie. Les nœuds adoptent ce comportement :

* Les nœuds du centre qui sont fortement interconnectés et intra connectés. Leurs distances intra et interclasses sont faibles.

* Les nœuds périphériques ne sont pas connectés à d'autres nœuds périphériques. Les nœuds d'une même classe peuvent être fortement connectés entre eux.
**Structure $\Lambda_{cp}$**

$$ \Lambda_{cp} = \begin{pmatrix} 2 & 5 & 8 \\ 5 & 4 & 6 \\ 8 & 6 & 10 \end{pmatrix} $$

Heatmap d'une matrice simulée

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

**Hiérarchique 2** [FB19] Dans la littérature, nous avons trouvé une seconde structure hiérarchique que nous appelons hiérarchique 2 afin d'éviter la confusion. Dans cette structure, les distances intra classes sont plus faibles que les distances interclasses. De plus, les distances sont organisées de manière hiérarchique "croissante".

C'est-à-dire :

Les distances intra classes vérifient :

$$ \Lambda_{\text{hier}}[3,3] > \Lambda_{\text{hier}}[2,2] > \Lambda_{\text{hier}}[1,1] $$

Les distances interclasses vérifient :

$$ \Lambda_{\text{hier}}[1,2] \le \Lambda_{\text{hier}}[2,3] \le \Lambda_{\text{hier}}[1,3] $$

**Structure $\Lambda_{h2}$**

$$ \Lambda_{h2} = (Th) $$

Heatmap d'une matrice simulée

<figure><img src="image_26.png" /></figure>
## 6.4 Protocole pour l'évaluation de l'approche TT

**Objectif** : Évaluer l'approche TT pour le calcul des marginales unaires et binaires. Pour cela, nous comparons avec les trois méthodes présentées dans la section 6.2.

Dans ce but, nous adopterons la procédure suivante pour chacune parmi les six structures de distances présentées dans la section 6.3 :

Pour une structure de distance donnée :

1. Fixer le nombre de classes $Q = 3$ et $\alpha = \left(\frac{1}{Q} \quad \frac{1}{Q} \quad \frac{1}{Q}\right)$

2. Choisir le nombre d'individus $n$.

3. Affecter à chaque nœud $i$ une classe ($Z_i$ et $q$ tel que $Z_{i,q} = 1$) selon $\alpha$.

4. Simuler $D$ connaissant les classes de chaque nœud et la structure de distance choisie.

5. Calculer les facteurs unaires et binaires.

6. Calculer les marginales unaires et binaires pour chaque méthode.

7. Comparer les résultats obtenus par paire de méthode.

Dans les deux sous-sections suivantes, nous présentons la méthode avec laquelle nous simulons les matrices de distances puis les paramètres utilisés pour nos expérimentations.

### 6.4.1 Simulation des matrices de distances

Pour un nombre d'individus $n$, un nombre de classes $Q$, une structure de distances entre classes $\Lambda$ et un vecteur de proportion d'individus dans les classes $\alpha$, nous simulons la matrice de distances $D$ de cette manière.

**Étape 1** : Simulation des classes des individus :

* Pour chaque individu $i = 1, \dots, n$, nous simulons une valeur $Z_i$ qui sera sa classe par l'intermédiaire d'une loi multinomiale de paramètre $\alpha$.

**Étape 2** : Simulation des matrices de distances :

* Pour chaque $i \in [1, n]$ et $j \in [i, n]$, nous génèrons $D(i, j)$ selon une loi de Poisson de paramètre $\lambda_{Z_i, Z_j}$.

### 6.4.2 Paramètres des méthodes

Voici le paramétrage choisi pour la mise en œuvre des méthodes de calcul des marginales.

* **TT**

* Nous faisons des décompositions en valeurs singulières complètes des matrices des facteurs.

* Nous faisons les calculs des produits des TT matrices avec une précision du rounding $\epsilon$ égale à $10^{-2}$ pour $n = 9$ et $10^{-3}$ pour $n = 18$.

* **Champ moyen** : Nous avons fixé le critère d'arrêt $\epsilon_{pf}$ du schéma du point fixe à $10^{-2}$ pour toutes les expériences.

* **Échantillonneur de Gibbs** : Concernant cette méthode, nous adaptons le nombre d'échantillonneurs de Gibbs ($L$) et d'itérations ($N_{iter}$) selon la structure de distance. En effet, nous avons remarqué que pour certaines structures de distances, nous avons besoin d'augmenter le nombre d'itérations $N_{iter}$ pour simuler une réalisation précise de la loi jointe ainsi que le nombre d'itérations $L$ pour améliorer la précision de l'estimation des marginales.

Dans les tableaux ci-dessous, nous présentons les paramètres choisis.
• $n = 9$

<table><thead><tr><th>Structure de distance</th><th>associative</th><th>hiérarchique</th><th>ordonnée</th><th>dissociative</th><th>core-periphérie</th><th>hiérarchique 2</th></tr></thead><tbody><tr><td>Nombre d'échantillons</td><td>10000</td><td>1000</td><td>1000</td><td>5000</td><td>1000</td><td>5000</td></tr><tr><td>Nombre d'itérations</td><td>1000</td><td>500</td><td>500</td><td>500</td><td>500</td><td>1000</td></tr></tbody></table>

• $n = 18$

<table><thead><tr><th>Structure de distance</th><th>associative</th><th>hiérarchique</th><th>ordonnée</th><th>dissociative</th><th>core-periphérie</th><th>hiérarchique 2</th></tr></thead><tbody><tr><td>Nombre d'échantillons</td><td>7000</td><td>7000</td><td>1000</td><td>12000</td><td>1000</td><td>14000</td></tr><tr><td>Nombre d'itérations</td><td>1000</td><td>1000</td><td>500</td><td>500</td><td>500</td><td>500</td></tr></tbody></table>

### 6.4.3 Visualisation des résultats

Pour la visualisation des différences entre les marginales unaires et binaires pour les quatre méthodes, nous avons choisi d'utiliser les pairplots fournis par la librairie ggplot2 du langage R.

Un pairplot est une matrice 4 par 4 de graphiques :

* **Diagonale** La densité des marginales unaires et binaires pour chaque méthode.

* **Partie triangulaire supérieure** : coefficient de corrélation entre les marginales (unaires ou binaire) pour chaque paire de méthode.

* **Partie triangulaire inférieure** : Le tracé des marginales (unaires ou binaire) d'une méthode en fonction de celle d'une autre méthode. Si les deux méthodes calculent des marginales identiques, ce graphique correspond au tracé de la droite $x = y$.

### 6.4.4 Présentation des résultats

Dans les deux sections suivantes, nous présentons les résultats de comparaisons entre les méthodes via des pairplots entre chaque paire de méthodes à la fois pour les marginales unaires et binaires.

Les résultats seront présentés de cette manière, pour chaque structure des distances, nous présentons les classes simulées, une heatmap de la matrice des distances organisée selon les classes, puis les pairplots des comparaisons des méthodes par paire à la fois pour les marginales unaires à gauche puis les binaires à droite.
# 6.5 Comparaison pour un nombre d'individus égale à 9

## 6.5.1 Associative

* Classes simulées (2 2 2 2 1 2 2 0 2)

* Heatmap de la matrice des distances

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

* Pairplots

<figure><img src="image_28.png" /><figcaption>FIGURE 6.2 – Marginales unaires (structure associative (n=9))</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>FIGURE 6.3 – Marginales binaires (structure associative (n=9))</figcaption></figure>

6.5.2 Hiérarchique

• Classes simulées (0 1 0 0 2 2 1 1 1)

• Heatmap de la matrice des distances

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

• Pair plots
<figure><img src="image_31.png" /><figcaption>FIGURE 6.4 – Marginales unaires (structure hiérarchique (n=9))</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>FIGURE 6.5 – Marginales binaires (structure hiérarchique (n=9))</figcaption></figure>

### 6.5.3 Ordonné

* Classes simulées (2 2 1 2 1 2 1 0 2)

* Heatmap de la matrice des distances
<figure><img src="image_33.png" /></figure>

• Pair plots

<figure><img src="image_34.png" /><figcaption>FIGURE 6.6 – Marginales unaires (structure ordonné (n=9))</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>FIGURE 6.7 – Marginales binaires (structure ordonné (n=9))</figcaption></figure>
### 6.5.4 Dissociative

* Classes simulées (1 1 0 2 2 1 2 0 1)

* Heatmap de la matrice des distances

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

* Pair plots

<figure><img src="image_37.png" /><figcaption>FIGURE 6.8 – Marginales unaires (structure dissociative (n=9))</figcaption></figure>
<figure><img src="image_38.png" /><figcaption>FIGURE 6.9 – Marginales binaires (structure dissociative (n=9))</figcaption></figure>

### 6.5.5 Cores periphery

* Classes simulées (2 0 2 2 0 1 2 0 0)

* Heatmap de la matrice des distances

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

* Pair plots
<figure><img src="image_40.png" /><figcaption>FIGURE 6.10 – Marginales unaires (structure Cores periphery (n=9))</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>FIGURE 6.11 – Marginales binaires (structure Cores periphery (n=9))</figcaption></figure>

### 6.5.6 Hiérarchique 2

* Classes simulées (0 2 2 0 0 0 0 0 0)

* Heatmap de la matrice des distances
<figure><img src="image_42.png" /></figure>

• Pair plots

<figure><img src="image_43.png" /><figcaption>FIGURE 6.12 – Marginales unaires (structure hiérarchique 2 (n=9))</figcaption></figure>

<figure><img src="image_44.png" /><figcaption>FIGURE 6.13 – Marginales binaires (structure hiérarchique 2 (n=9))</figcaption></figure>
### 6.5.7 Remarques sur le temps de calcul

Les expériences ont été réalisées sans parallélisation, sur un ordinateur présentant les caractéristiques suivantes : Processeur : Intel(R) Core(TM) i9-9880H CPU @ 2.30GHz ainsi qu'une mémoire totale de 62 go ( pour **15 go utilisées**). Voici un tableau récapitulatif des temps des calculs :

<table><caption>TABLE 6.1 – Tableau présentant les temps de calcul des marginales unaires et binaires, pour $n = 9$.</caption><thead><tr><th>Méthode</th><th>Point fixe</th><th>Gibbs</th><th>TT</th></tr></thead><tbody><tr><th>Temps de calcul (s)</th><td>autour de 0.02</td><td>entre 189 et 3608</td><td>entre 2.12 et 6</td></tr></tbody></table>
# 6.6 Comparaison pour un nombre d'individus égal à 18

## 6.6.1 Associative

* Classes simulées (2 1 2 1 0 1 0 0 0 1 0 1 0 1 1 0 0 2)

* Heatmap de la matrice des distances

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

* Pair plots

<figure><img src="image_46.png" /><figcaption>FIGURE 6.14 – Marginales unaires (structure associative (n=18))</figcaption></figure>
<figure><img src="image_47.png" /><figcaption>FIGURE 6.15 – Marginales binaires (structure associative (n=18))</figcaption></figure>

6.6.2 Hiérarchique

• Classes simulées (2 0 2 1 0 1 2 0 1 2 0 0 2 1 2 1 2 1)

• Heatmap de la matrice des distances

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

• Pair plots
<figure><img src="image_49.png" /><figcaption>FIGURE 6.16 – Marginales unaires (structure hiérarchique (n=18))</figcaption></figure>

<figure><img src="image_50.png" /><figcaption>FIGURE 6.17 – Marginales binaires (structure hiérarchique (n=18))</figcaption></figure>

### 6.6.3 Ordonné

* Classes simulées (0 0 1 2 2 1 0 0 0 1 0 0 1 0 1 0 1 2)

* Heatmap de la matrice des distances
<figure><img src="image_51.png" /></figure>

• Pair plots

<figure><img src="image_52.png" /><figcaption>FIGURE 6.18 – Marginales unaires (structure ordonné (n=18))</figcaption></figure>

<figure><img src="image_53.png" /><figcaption>FIGURE 6.19 – Marginales binaires (structure ordonné (n=18))</figcaption></figure>
### 6.6.4 Dissociative

* Classes simulées (1 0 0 1 0 1 2 2 0 1 1 2 2 2 1 1 2 1)

* Heatmap de la matrice des distances

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

* Pair plots

<figure><img src="image_55.png" /><figcaption>FIGURE 6.20 – Marginales unaires (structure dissociative (n=18))</figcaption></figure>
<figure><img src="image_56.png" /><figcaption>FIGURE 6.21 – Marginales binaires (structure dissociative (n=18))</figcaption></figure>

6.6.5 Cores periphery

• Classes simulées (2 2 1 2 1 0 0 2 0 1 2 1 0 1 0 1 1 0)

• Heatmap de la matrice des distances

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

• Pair plots
<figure><img src="image_58.png" /><figcaption>FIGURE 6.22 – Marginales unaires (structure Cores periphery (n=18))</figcaption></figure>

<figure><img src="image_59.png" /><figcaption>FIGURE 6.23 – Marginales binaires (structure Cores periphery (n=18))</figcaption></figure>

### 6.6.6 Hiérarchique 2

* Classes simulées (2 1 2 2 1 1 2 0 0 0 1 2 0 0 2 0 0,2)

* Heatmap de la matrice des distances
<figure><img src="image_60.png" /></figure>

• Pair plots

<figure><img src="image_61.png" /><figcaption>FIGURE 6.24 – Marginales unaires (structure hiérarchique 2 (n=18))</figcaption></figure>

<figure><img src="image_62.png" /><figcaption>FIGURE 6.25 – Marginales binaires (structure hiérarchique 2 (n=18))</figcaption></figure>
### 6.6.7 Remarque sur le temps de calcul

Comme pour un nombre d'individus $n = 9$, nous présentons les temps de calcul que chaque méthode a nécessité pour le calcul des marginales pour des SBM avec $n = 18$.

<table><caption>TABLE 6.2 – Tableau présentant les temps de calcul des marginales unaires et binaires, pour $n = 18$.</caption><thead><tr><th>Temps de calcul (s)</th><th>Méthode</th><th>Point fixe</th><th>Gibbs</th><th>TT</th></tr></thead><tbody><tr><td></td><td></td><td>Entre 0.4 et 0.6</td><td>Entre 2834 et 8139</td><td>Entre 50 et 400</td></tr></tbody></table>

Nous observons dans les tableaux 6.1 et 6.2 que le temps est multiplié par un facteur entre 20 et 30 pour la méthode du Point Fixe, il est multiplié par un facteur de 50 pour l'approche TT. Globalement, le temps de calcul augmente énormément entre les deux familles d'expériences.

## 6.7 Discussion

Nous commencerons cette section par un résumé des résultats que nous avons présentés durant ce chapitre, ensuite nous nous intéresserons à la cohérence entre les différentes méthodes de calcul des marginales. Enfin, nous nous pencherons sur la non convergence du schéma point du point fixe utilisé dans la méthode du champs moyen.

### 6.7.1 Résumé des résultats

Nous avons présenté deux expériences où nous avons comparé différentes méthodes de calculs de marginales unaires et binaires d'un SBM sur six structures de distances.

Pour chaque nombre d'individus et chaque structure de distance, nous avons présenté les résultats par l'intermédiaire de pairplots.

* Comparaisons pour un nombre d'individus $n = 9$ : nous prenons comme référence la méthode exacte par énumération.

Dans nos résultats, nous montrons une dépendance entre les précisions des méthodes et la structure de distances à partir de laquelle nous avons simulé la matrice de distance. Mais généralement, voici ce que nous pouvons dire pour les six structures de distances :

* Comparaison entre Gibbs, TT et la méthode exacte par énumération : Les coefficients de corrélations entre l'approche de Gibbs (respectivement l'approche TT) et notre référence varient entre 0.944 et 1 (respectivement entre 0.986 et 1) pour les marginales unaires et 0.931 et 1 (respectivement entre 0.979 et 1) pour les binaires. De plus, les scatterplots confirment cette proximité entre ces deux méthodes et notre référence.

Nous avons en plus établi que l'approche de Gibbs nécessite dans certains cas un nombre d'échantillonneurs $L$ et d'itérations $N_{iter}$ très élevés qui conduisent à des temps de calculs relativement lents pouvant atteindre une heure tandis que l'approche TT ne nécessitent que quelques secondes.

* Comparaison entre le champ moyen et la méthode exacte par énumération : Au vu des scatterplots et des coefficients de corrélation, nous pouvons dire que si nous excluons la structure ordonné où les corrélations sont de -0.499 pour les marginales unaires et -0.218 pour les binaires, la méthode du champ moyen se comporte de manière similaire aux deux autres vis-à-vis de notre référence : les corrélations varient
entre 0.984 et 1 pour les marginales unaires et entre 0.977 et 1 pour les binaires.
L'avantage de cette méthode est qu'elle nécessite des temps de calculs relativement bas par rapport aux deux autres. Cependant, il y a eu plusieurs cas où le schéma du point fixe n'a pas convergé.

* **Comparaison pour un nombre d'individus** $n = 18$ : nous prenons désormais comme référence la méthode de Gibbs. Nous obtenons les résultats suivants :

* **Comparaison entre Gibbs et TT** : Les coefficients de corrélation entre les marginales calculées par Gibbs et l'approche TT varient entre 0.985 et 1 pour les unaires et 0.985 et 1 pour les binaires. De plus, les scatterplot vont vers le sens d'une fonction $X = Y$. Cela indique une très grande proximité entre les deux méthodes.
En ce qui concerne les temps de calcul, celui de l'approche TT est beaucoup plus raisonnable que celui de Gibbs.

* **Comparaison entre Gibbs et le champ moyen** : Les résultats indiquent que pour les structures de distances ordonné et dissociative il y a une absence d'accord entre les deux méthodes qui se manifeste dans les faibles corrélations entre les deux méthodes : pour la structure dissociative (resp. ordonné) les corrélations sont -0.500 (resp. 0.693) pour les marginales unaires et -0.125 (resp. 0.688) pour les binaires.
Pour les quatre autres structures de distances, les corrélations entre les marginales calculées par les deux méthodes varient 0.985 et 1 pour les unaires et entre 0.983 et 1 pour les binaires. Nous pouvons donc en conclure que les méthodes sont en accord pour ces quatre structures.

## 6.7.2 Cohérence entre les différentes méthodes de calculs des marginales

Nous avons constaté que les mêmes patterns pour les deux deux tailles de problème.

### Proximité entre les méthodes

* **Méthode exacte par énumération par rapport aux autres méthodes** : Mise à part le champ moyen dans les structures ordonné et dissociative, les quatre méthodes donnent globalement des résultats très proches.
Comme nous ne pouvons plus utiliser cette méthode quand le nombre d'individus est grand, cela nous permet de valider la méthode de Gibbs comme notre référence.

* **Approche TT et Gibbs** : Pour les six structures de distances étudiées, nous avons le même constat, que ce soit pour les marginales unaires ou binaires.
Nous nous penchons maintenant vers les corrélations entre les deux méthodes, nous remarquons des corrélations au-dessus de 0.999 dans les structures associative, Cores periphery et hiérarchique 2. La corrélation minimale est de 0.886 (atteinte entre les marginales binaires calculées dans le cadre de la structure hiérarchique).
Ces résultats nous permettent de dire que l'approche TT donne de très bons résultats comparés à l'approche de Gibbs qui est réputée très précise en littérature.

* **Approche Point fixe et Gibbs** : L'approche du champ moyen qui est l'approche la plus rapide dans la littérature semble diverger par rapport à la méthode de Gibbs sur deux structures de distances : ordonné et dissociative où nous tombons sur des corrélations très faibles que ce soit pour les marginales binaires ou pour les unaires.
Cette divergence pour les marginales binaires pourrait s'expliquer par le fait que nous faisons l'hypothèse d'indépendance des classes sachant le graphe ce qui nous conduit à faire le calcul des marginales binaires comme produit des unaires.
Sur les quatre autres distances, nous observons une très forte corrélation à la fois avec Gibbs.
### 6.7.3 Temps de calculs/ Précision :

Si nous comparons les temps de calcul, nous réalisons que :

* L'approche de Gibbs nécessite des temps de calculs plus lents étant donné qu'elle exige un nombre important d'itérations pour certaines structures de distances : 12000 pour dissociative et 14000 pour hiérarchique 2.

* L'approche Champ moyen, bien qu'elle soit très rapide, fournit des résultats sont de qualité variable.

* L'approche TT donne des résultats très précis et nécessite des temps de calculs très faibles comparés à Gibbs. Nous pouvons dire que cette méthode donne le meilleur compromis temps de calcul/ précision.

### 6.7.4 Influence de la structure de distance :

Nos résultats montrent que dans les deux tailles de problème, pour certaines structures de distances toutes les méthodes de calcul approché font des calculs précis des marginales. Il s'agit des structures Cores periphery et hiérarchique.

Pour d'autres structures comme dissociative et ordonnée, les calculs des marginales par le champ moyen est peu précis.

Nous avons aussi constaté que pour la structure dissociative, l'échantillonneur de Gibbs nécessite l'utilisation d'un grand nombre d'échantillons afin d'avoir une bonne précision du calcul des marginales.

Cela suggère que les 6 structures s'ordonnent par précision décroissance avec l'ordre suivant : Cores periphery, hiérarchique, hiérarchique 2, associative, ordonné puis dissociative.

**NB** : Les résultats que nous avons présentés sont issus de matrices de distances et des facteurs unaires pour lesquels la méthode du point fixe converge.

## 6.8 Conclusion

### 6.8.1 Qualité de l'approche TT

Les résultats montrent que l'approche TT est très précise lorsque l'on prend pour référence la méthode exacte par énumération ou l'échantillonneur de Gibbs. Nous avons borné l'erreur des marginales comparé à un calcul exacte dans les formules (B.0.11) et (B.0.11) de l'annexe B qui porte sur la précision des calculs au format TT.

De plus, l'approche TT permet de calculer les marginales dans un temps relativement raisonnable comparé à l'approche de Gibbs.

Nous avons également remarqué que cette approche est stable vis-à-vis de la variation du nombre d'individus et de la structure de distance. À partir des résultats obtenus, nous pouvons dire que cette approche est celle qui respecte le plus le compromis précision/temps de calcul.

### 6.8.2 Perspectives

Dans le cadre de ces expériences, nous avons utilisé l'approche TT avec un calcul des marginales comme présentée dans le chapitre 5. Cette approche nécessite beaucoup de stockage et un temps de calcul relativement long. Ceci fait que quand le nombre d'individus augmente, nous ne pourrons plus utiliser cette approche, nous devrons ajouter certaines optimisations que nous allons présenter dans le chapitre 7.
# Chapitre 7

## Vers le passage à l'échelle de l'approche TT

<table><caption>Table des matières</caption><tbody><tr><td>7.1</td><td>Introduction</td><td>110</td></tr><tr><td>7.2</td><td>Procédure d'homothétie</td><td>112</td></tr><tr><td>7.2.1</td><td>Représentation en double précision</td><td>113</td></tr><tr><td>7.2.2</td><td>Motivations</td><td>114</td></tr><tr><td>7.2.3</td><td>Procédure mise en oeuvre</td><td>114</td></tr><tr><td>7.2.4</td><td>Conclusion</td><td>116</td></tr><tr><td>7.3</td><td>Changement de paramètres dans la compression par l'algorithme du rounding</td><td>116</td></tr><tr><td>7.3.1</td><td>Choix de la précision du rounding</td><td>116</td></tr><tr><td>7.3.2</td><td>Explication du problème</td><td>117</td></tr><tr><td>7.3.3</td><td>Présentation du rounding</td><td>120</td></tr><tr><td>7.3.4</td><td>Résultats</td><td>120</td></tr><tr><td>7.3.4.1</td><td>Calcul de la constante de normalisation</td><td>120</td></tr><tr><td>7.3.4.2</td><td>Réitération des calculs</td><td>123</td></tr><tr><td>7.3.4.3</td><td>Comparaison des marginales</td><td>124</td></tr><tr><td>7.3.5</td><td>Discussion</td><td>125</td></tr><tr><td>7.3.6</td><td>Conclusion</td><td>126</td></tr><tr><td>7.4</td><td>Réduction des dimensions de TT-matrices de rang 1 par fusion des cores</td><td>126</td></tr><tr><td>7.4.1</td><td>Motivation</td><td>126</td></tr><tr><td>7.4.2</td><td>Explication de la fusion des cores</td><td>127</td></tr><tr><td>7.4.2.1</td><td>Réduction du nombre des cores d'une TT-matrice de rang 1</td><td>127</td></tr><tr><td>7.4.3</td><td>Fusion des cores de la TT-matrice A<sub>i</sub>[z<sub>i</sub>]</td><td>128</td></tr><tr><td>7.4.3.1</td><td>Fusion avec un nombre uniforme de cores</td><td>129</td></tr><tr><td>7.4.3.2</td><td>Fusion avec volumes uniformes des cores</td><td>131</td></tr><tr><td>7.4.4</td><td>Résultats</td><td>133</td></tr><tr><td>7.4.4.1</td><td>Influence de la méthode de la fusion sur les volumes des cores</td><td>133</td></tr><tr><td>7.4.4.2</td><td>Diminution des rangs des cores lors du calcul de W</td><td>135</td></tr><tr><td>7.4.4.3</td><td>Comparaison des constantes de normalisations obtenues par les trois méthodes</td><td>136</td></tr><tr><td>7.4.5</td><td>Discussion</td><td>138</td></tr><tr><td>7.4.6</td><td>Conclusion</td><td>138</td></tr><tr><td>7.5</td><td>Conclusion du chapitre</td><td>138</td></tr></tbody></table>
## 7.1 Introduction

Dans le chapitre 2, nous avons introduit les SBM qui sont des modèles graphiques très complexes. En effet, comme le graphe des liens entre les sommets est complet, il permet de représenter toutes les interactions entre les individus du graphe, donc deux sommets quelconques sont toujours reliés. Ainsi, si nous avons $n \in \mathbb{N}$ individus, le nombre de liens est $m = \frac{n(n-1)}{2}$. Cela induit une complexité quadratique par rapport au nombre d'individus des calculs d'inférence sur ce type de modèle.

Des difficultés calculatoires apparaissent quand nous appliquons la procédure de Novikov [NROV14] que nous avons expliqué dans le chapitre 5, étant donné que les TT-matrices $A_i[Z_i]$ et $B_i$ sont d'ordre $m$ et donc difficilement manipulables quand $m$ est grand. En effet, à chaque somme où produit de ces TT-matrices, nous devrons manipuler $m$ tenseurs (TT-cores) de petites tailles de manière séquentielle.

Cela nous contraint à des optimisations sur cette procédure afin d'éviter de travailler avec des valeurs numériques trop petites, de limiter l'accroissement des rangs des cores ainsi que de réduire le nombre de cores des TT-matrices. Dans la suite de cette section, nous détaillons les trois limitations que nous avons rencontrées lors de l'implémentation de la procédure de Novikov sur des SBM et comment nous y avons répondu.

**Éviter les problèmes des petites valeurs.** Les valeurs des marginales non normalisées et de la constante de normalisation sont très petites. Nous avons constaté que leur logarithme est proportionnel au nombre de liens $m$ de notre modèle, avec le coefficient de proportionnalité négatif.

Dans l'exemple qui suit, nous avons simulé 10 SBM avec des proportions d'individus dans les classes equi-répartis ($\alpha = [\frac{1}{3}, \frac{1}{3}, \frac{1}{3}]$) et des nombres d'individus variant de 8 à 28 par pas de 2. Pour chaque nombre d'individus, nous simulons une matrice de distance et pour chaque modèle, nous calculons la constante de normalisation $W$.

Nous présentons dans la Figure 7.1 l'évolution de $-\log(W)$ en fonction du nombre de liens $m = \frac{n(n-1)}{2}$ pour deux structures de distances associative et dissociative qui sont des structures classiques pour les SBM et que nous avons présenté dans la section 6.3.

**Structure Associative :**

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

**Structure Hiérarchique :**

$$ \Lambda = \begin{pmatrix} 12 & 3 & 1 \\ 3 & 6 & 3 \\ 1 & 3 & 3 \end{pmatrix} $$

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

<figure><img src="image_64.png" /><figcaption>FIGURE 7.1 – Evolution de $-\log(W)$ par rapport au nombre de liens selon une structure de distance associative (à droite) et hiérarchique (à gauche).</figcaption></figure>
Nous calculons la pente des courbes ($m, -\log(W)$) pour chacun des graphes : la pente est de 1.05 pour la structure associative et 0.90 pour la structure hiérarchique. Nous remarquons que l'ordre de grandeur des valeurs de $W$ diffère selon la structure de distances choisie, mais il reste proportionnel à $10^{-m} = 10^{-\frac{n(n-1)}{2}}$.

À partir de cette remarque, nous pouvons dire que pour des nombres d'individus élevés, les valeurs des constantes de normalisation sont très petites. Cela pose deux problèmes majeurs liés aux limites relatives de la représentation des petites valeurs en arithmétique finie. En effet, quand les valeurs calculées sont très faibles, leur représentation en double précision est difficile et cela conduit à obtenir des valeurs de $-0.0$ ou $0.0$.

En conclusion, à cause des problèmes de représentation des petites valeurs en arithmétique finie, nous nous retrouvons avec des valeurs de la constante de normalisation et des marginales nulles. De plus, lors de la phase $LQ$ de l'algorithme du rounding, nous pouvons obtenir une racine négative ($\sqrt{-0.0}$) ou des logarithmes de valeurs négatives ($\text{Log}(-0.0)$) qui engendre soit des NaN (not-a-number) soit des problèmes de numpy `RuntimeWarning`.

Nous présentons la procédure qui nous permet d'éviter ces problèmes dans la section 7.2.

**Explosion des rangs des cores.** Dans la méthode proposée par Novikov [NROV14], le calcul de la constante de normalisation, qui est très similaire à celui des marginales (détaillé dans la section 5.2.6) est décrit dans l'algorithme 1.

**Algorithm 1 Calcul de W**

**Require:** $n$ nombre d'individus, $[\mathbf{B}_1, \dots, \mathbf{B}_n]$ : liste de TT-matrices intervenant dans le calcul, $\epsilon$: précision du rounding.

**Ensure:** $\mathbf{W}$ : constante de normalisation.

1: $\mathbf{W} = \mathbf{B}_n$

2: **for** $i = n - 1, \dots, 1$ **do**

3: $\mathbf{W} = \mathbf{B}_i \times \mathbf{W}$

4: $\mathbf{W} = \text{rounding}(\mathbf{W}, \epsilon)$

5: **end for**

6: Transformation de $\mathbf{W}$ en scalaire $W$

7: **return** $W$

Dans l'algorithme 1, nous commençons par attribuer le TT-vecteur $\mathbf{B}_n$ à $\mathbf{W}$. Ensuite, chaque itération se décompose en deux phases, une phase de produits de TT-matrices-TT-vecteur (ligne 3 de l'algorithme 1) où nous calculons le produit $\mathbf{B}_n \dots \mathbf{B}_{i+1} \times \mathbf{B}_i$ et une phase de compression (ligne 4 de l'algorithme 1) où nous faisons la compression du TT-vecteur résultant par l'algorithme du rounding.

Nous faisons ces deux phases étant donné que chaque produit entre deux TT-matrices augmente le rang de la TT-matrice résultante et que si nous faisons les calculs sans recompression du résultat (ligne 4) à chaque itération, les rangs des TT-cores qui composent les TT-matrices augmenteront jusqu'à atteindre la valeur théorique de $Q^n$ alors que la constante de normalisation est un scalaire.

Afin d'éviter ces rangs élevés, nous sommes contraints de recompresser le TT-vecteur résultant après chaque itération comme nous présentons dans la ligne 4 de l'algorithme 1 afin de limiter la croissance des rangs des TT-cores.

Dans la procédure de Novikov et al., la compression est faite avec l'algorithme du rounding avec une précision $\epsilon$ fixée, cela nous permet en plus de contrôler la précision du résultat final. Cependant, faire la compression de cette manière n'empêche pas les rangs de croître fortement quand $n$ augmente.
Pour illustrer nos propos, nous calculons la constante de normalisation par l'intermédiaire de l'algorithme 1 pour un modèle SBM avec comme paramètres : $n = 22$, $Q = 3$, $\Lambda = \begin{pmatrix} 2 & 10 & 10 \\ 10 & 3 & 10 \\ 10 & 10 & 4 \end{pmatrix}$ et $\alpha = (\frac{1}{3}, \frac{1}{3}, \frac{1}{3})$.

<figure><img src="image_65.png" /><figcaption>FIGURE 7.2 – Évolution de rang dans le calcul de W en fonction de l'itération.</figcaption></figure>

La Figure 7.2 présente l'évolution du rang maximal des cores de W en fonction des itérations de l'algorithme 1 avec le paramètre de compression $\epsilon$ fixé à $5 \times 10^{-2}$. Nous observons que les rangs augmentent jusqu'à atteindre 140 avant de diminuer pour atteindre 35 lors de la dernière itération. Sans rounding ce calcul n'est pas possible car le TT-rang maximal est : $3^{22} = 31381059609$. Nous présentons la procédure qui nous permet de réduire les rangs des cores des TT-matrices dans la section 7.3.

**Réduction du nombre de cores.** Pour un modèle SBM avec $n$ individus, nous avons $m$ liens. Chacune des TT-matrices (ou TT-vecteurs) $A_i[Z_i]$ et $B_i$ est d'ordre $m$. Cela a pour conséquences : Chaque opération (somme ou produit) sur les TT-matrices $A_i[Z_i]$ et $B_i$ nécessite d'itérer séquentiellement sur les $m$ cores (détaillé dans la section 4.4.2).

De plus, nous avons expliqué dans la section 4.4.4 que pour une TT-matrice composée de $m$ cores, l'algorithme du rounding consiste à faire une orthogonalisation (décomposition QR) puis une décomposition en valeur singulière (SVD) pour chaque core, le tout en parcourant séquentiellement toutes les dimensions.

Ces deux raisons font que nous avons besoin de trouver une manière de diminuer l'ordre du tenseur afin de réduire la taille de la boucle sur les dimensions ainsi que le temps du calcul des produits des TT-matrices (ligne 3 de l'algorithme 1) et de la compression par rounding (ligne 4 de l'algorithme 1) à chaque itération.

Nous présentons la procédure qui nous permet de diminuer l'ordre du tenseur dans la section 7.4.

## 7.2 Procédure d'homothétie

Nous avons présenté auparavant l'évolution des valeurs obtenues pour la constante de normalisation et la figure 7.1 indique qu'elles sont proportionnelles à $10^{-m}$.
Cela ne gêne pas d'un point de vue mathématique, par contre, pour l'implémentation des calculs, les nombres sont stockés soit en simple, soit en double précision. Dans la section qui suit, nous expliquons comment se fait la représentation informatique des valeurs en double précision. Nous présentons ensuite une solution basée sur une homothétie pour s'affranchir de ce problème.

### 7.2.1 Représentation en double précision

En informatique, les nombres sont représentés selon le standard IEEE 754 [46108], celle-ci stipule qu'un réel $v$ est représenté en double précision (resp. simple précision) par sa mantisse, son exposant et son signe en 64 bits (resp. 32 bits) selon la figure a (resp. b) du schéma 7.3 (figure 5 de l'article [NLTK19] ) :

<figure><img src="image_66.png" /><figcaption>FIGURE 7.3 – Standard IEEE 754 pour la double précision.</figcaption></figure>

Ce format de stockage en 64 bits des nombres se décompose en trois parties : 1 bit pour le stockage du signe, 52 bits pour le stockage de la mantisse et 11 bits pour le stockage de l'exposant.

Le réel $v$ est représenté en nombre normalisé de cette manière :

$$v = s \times (1 + ma) \times 2^{(e-d)} \qquad (7.2.1)$$

Avec,

* $s \in \{-1, +1\}$ représente le signe de $v$

* $e$ exposant.

* $d$ décalage.

* $ma$ représente la mantisse $0 \le ma \le 1$.

À partir de la formule (7.2.1), le plus petit nombre positif normalisé qui peut être représenté est : $2^{-(2^{111}-2)} = 2^{-(2^{10}-2)} \approx 2.22 \times 10^{-308}$.

Il existe une autre représentation des réels appelée représentation en nombres dénormalisés, une manière de représenter des nombres ayant une valeur très proche de zéro. Dans cette représentation, les nombres vérifient les conditions suivantes : $ma \neq 0$ et $d = 0$.

L'exposant devient baisé par $2^{111-1}-1$ alors la formule (7.2.1) pour stocker les nombres devient :

$$v = s \times (1 + ma) \times 2^{(-d)} \qquad (7.2.2)$$

À partir de la formule (7.2.2), la plus petite valeur qu'on peut représenter est :

$$v_{min} = 2^{-(2^{se-1}-2+sm)} = 2^{-1074} \approx 5 \times 10^{-324} \qquad (7.2.3)$$
**Remarque :** Au niveau de l'implémentation nous utilisation la représentation en nombre dénor-
malisé [SST03] vue qu'elle nous offre plus de flexibilité et donc le réel le plus petit que nous
pouvons représenter est $5 \times 10^{-324}$.

7.2.2 Motivations

Comme nous travaillons avec des petites valeurs, cela présente deux difficultés majeures pour
l'implémentation des calculs présentés dans la section 5.

* Premièrement, comme les valeurs des marginales et des constantes de normalisation sont très petites, alors il est difficile de les représenter en arithmétique finie et donc une valeur de 0.0 ou -0.0 leur est attribuée.

* Deuxièmement, ses valeurs très basses induisent des erreurs au niveau de la phase d'orthogonalisation lors de la compression via l'algorithme du rounding. Lors de la phase d'orthogonalisation, nous sommes confrontés à des cores dont les normes sont tellement petites qu'ils sont considérés comme nul et d'autres avec une norme égale à 1. Pour illustrer ces problèmes, nous partons d'une TT-matrice A qui s'exprime :

$$
\mathbf{A}[i_1, \dots, i_d; j_1, \dots, j_d] = \mathbf{G}_1(i_1, j_1) \cdot \mathbf{G}_2(i_2, j_2) \cdots \mathbf{G}_d(i_d, j_d)
$$

avec, $k \in [1, d]$, $G_k$ le k-ième core de $\mathbf{A}$.

Si les normes de ces cores sont très faibles, c'est-à-dire proches de $v_{min}$, alors, pour la
phase d'orthogonalisation du rounding, les décompositions $LQ$ des cores $G_k$ donnent des
matrices $Q$ de norme 1 et des matrices $L$ de valeurs très faibles.

Á cause de cela, nous nous retrouvons à faire des opérations comme les sommes, les soustractions, les division et les racines sur des valeurs très faibles. Cela conduit à des problèmes d'instabilités numériques, des divisions par zéros et des NaN.

Afin de remédier à ces problèmes, nous allons faire une homothétie sur les facteurs.

7.2.3 Procédure mise en oeuvre

La procédure d'homothétie consiste à multiplier par un scalaire chaque facteur $\psi'_{i,j}$ défini par l'équation (5.2.5).

Soit $a \in \mathbb{R}$, nous définissons les facteurs normalisés par $\psi_{i,j}^p$ comme suit

$$
\forall i \in [1, n-1], j \in [i+1, n], \quad \psi_{i,j}^{p}(z_i, z_j) = a \psi_{i,j}'(z_i, z_j) \quad (7.2.4)
$$

**Proposition 8 (Homothétie des facteurs).** *Si nous multiplions tous les facteurs par une constante $a$ et que nous mettons ces facteurs à l'entrée de la procédure de Novikov, nous obtenons la nouvelle constante de normalisation $\bar{W}$ agrandie par rapport à $W$ de cette manière :*

$$
\bar{W} = a^m W
$$

et les nouvelles marginales non normalisées

$$
\bar{p}_{i,j}(z_i, z_j) = a^m p_{i,j}(z_i, z_j)
$$

Démonstration. On définit la loi non normalisée $\mathbb{P}_{\theta}(z, D)$ comme le produit de facteurs

$$
\mathbb{P}(z, D) = \left[ \prod_{i=1}^{n} \prod_{j>i} \psi'_{i,j}(z_i, z_j) \right]
$$
et la constante de normalisation s'écrit

$$W = \sum_z \mathbb{P}(z, D)$$

Nous pouvons exprimer $\mathbb{P}(z, D)$ de cette manière :

$$\mathbb{P}(z, D) = \left[ \prod_{i=1}^{n} \prod_{j>i} \frac{\psi_{i,j}^{p}(z_i, z_j)}{a} \right]$$

Alors,

$$\begin{aligned} \mathbb{P}(z, D) &= \left[ \prod_{i=1}^{n} \prod_{j>i} \psi_{i,j}^{p}(z_i, z_j) \right] \prod_{i=1}^{n} \prod_{j>i} \left[ \frac{1}{a} \right] \\ &= \left[ \prod_{i=1}^{n} \prod_{j>i} \psi_{i,j}^{p}(z_i, z_j) \right] a^{\sum_{i=1}^{n} \sum_{i>j} (-1)} \\ \mathbb{P}(z, D) &= \left[ \prod_{i=1}^{n} \prod_{j>i} \psi_{i,j}^{p}(z_i, z_j) \right] a^{-m} \end{aligned}$$

Si nous posons :

$$\bar{\mathbb{P}}(z, D) = \left[ \prod_{i=1}^{n} \prod_{j>i} \psi_{i,j}^{p}(z_i, z_j) \right]$$

Alors,

$$\mathbb{P}(z, D) = a^{-m} \bar{\mathbb{P}}(z, D)$$

Ensuite,

$$W = \sum_z \mathbb{P}(z, D) = a^{-m} \sum_z \bar{\mathbb{P}}(z, D) = a^{-m} \bar{W}$$

Nous aurons la même chose pour les marginales binaires et unaires non normalisées :

$$p_{i,j}(z_i, z_j) = \sum_{z \setminus z_i, z_j} \mathbb{P}(z, D) = a^{-m} \sum_{z \setminus z_i, z_j} \bar{\mathbb{P}}(z, D) = a^{-m} \bar{p}_{i,j}(z_i, z_j)$$

□

Si nous normalisons les marginales par $W$, nous perdons la constante d'homothétie. En effet :

$$\frac{\bar{p}_{i,j}(z_i, z_j)}{\bar{W}} = \frac{a^m p_{i,j}(z_i, z_j)}{a^m W} = \frac{p_{i,j}(z_i, z_j)}{W} \qquad (7.2.5)$$

Pour résumer cette procédure, nous augmentons numériquement les facteurs afin d'éviter les différents problèmes de petites valeurs puis nous retrouvons les vraies valeurs des marginales après normalisation.

Si nous nous basons sur les figures 7.1 présentées durant l'introduction, la pente des courbes nous suggère qu'il faut utiliser une constante d'homothétie $a = 10$.

Cette constante n'est pas la meilleure, mais elle nous permet d'éviter des valeurs très basses pour la constante de normalisation (inférieures à $10^{-324}$).
### 7.2.4 Conclusion

La valeur de la constante de normalisation peut facilement atteindre des ordres de grandeurs inférieurs à $10^{-325}$ au cas où le nombre d'individus dépasse 26, cela conduit à des problèmes de représentation des petits chiffres, ce qui cause des problèmes de stabilités numériques ainsi que des problèmes de NaN (Not a Number) lors de la phase d'orthogonalisation dans le rounding. L'approche par homothétie nous a permis d'atteindre le cap de $n = 26$ pour les calculs des marginales et $n = 33$ pour la constante de normalisation.

Cependant, nous ne pouvons pas utiliser le graphe présenter dans la section 5.3 à cause de l'explosion des rangs des cores des TT-matrices.

Dans ce qui va suivre, nous présentons une procédure qui nous permet de dépasser cette limite.

## 7.3 Changement de paramètres dans la compression par l'algorithme du rounding

Dans la proposition 2, nous avons montré que le produit de TT-matrices (resp. une TT-matrice par un TT-vecteur) augmente le rang de la TT-matrice (resp. le TT-vecteur) résultante. Pour un modèle SBM avec $n$ individus, en utilisant l'algorithme 1, le calcul de la constante de normalisation ou d'une marginale nécessite de faire $n-1$ produits TT-matrices TT-vecteur.

Nous avons indiqué durant l'introduction que l'algorithme classique du rounding ne suffit pas à limiter la croissance des rangs des TT-vecteurs obtenus après chaque itération de l'algorithme 1. Cela a pour effet d'augmenter l'espace nécessaire de stockage du tenseur ainsi que le temps de calcul. Nous montrons dans la section 7.3.2 que la méthode classique du rounding présentée durant la section 4.4.4 est insuffisante quand le nombre d'individus augmente et ne nous permet pas d'avoir un temps de calcul raisonnable.

Puis dans la section 7.3.3, nous regardons comment nous pouvons changer les paramètres de l'algorithme du rounding afin d'accélérer les calculs et nous donnons des résultats numérique dans la section 7.3.4.

### 7.3.1 Choix de la précision du rounding

Afin de définir la précision du rounding que nous utilisons pour le reste du chapitre, nous avons pris cinq valeurs $[1 \times 10^{-1}, 1 \times 10^{-2}, 5 \times 10^{-2}, 1 \times 10^{-3}, 1 \times 10^{-4}]$.

Pour chacune de ces valeurs, nous définissons une structure de distance $\Lambda$ pour chacune des six structures définies dans la section 6.3. Ensuite, nous avons simulé des SBM avec des nombres d'individus $n$ variant entre 8 et 20 par pas de deux. Pour chacun de ces modèles, nous avons défini le nombre de classes $Q=3$ et les proportions d'individus dans les classes $\alpha = (\frac{1}{3}, \frac{1}{3}, \frac{1}{3})$. Pour chaque nombre d'individus, nous simulons la matrice de distances $D$ et nous appliquons la procédure de Novikov, nous calculons les TT-matrices $B_i$. Enfin, nous utilisons l'algorithme 1 avec une précision parmi les cinq valeurs et nous considérons comme référence $\epsilon = 1 \times 10^{-4}$.

Pour la structure associative, le tableau 7.1 présente les constantes de normalisation et le tableau 7.2 présente les temps de calculs.

Nous observons que toutes les valeurs de la constante de normalisation obtenues sont proches de celle obtenue par le $\epsilon$ de référence à part celle obtenue pour un modèle SBM avec $n = 12$ et pour une précision $\epsilon = 1 \times 10^{-1}$.
<table><caption>TABLE 7.1 – Constante de normalisation pour différentes valeurs de $\epsilon$ dans le cas d'une structure associative</caption><thead><tr><th>n</th><th>1 &times; 10<sup>-1</sup></th><th>5 &times; 10<sup>-2</sup></th><th>1 &times; 10<sup>-2</sup></th><th>1 &times; 10<sup>-3</sup></th><th>1 &times; 10<sup>-4</sup></th></tr></thead><tbody><tr><td>8</td><td>1.30 &times; 10<sup>-30</sup></td><td>1.30 &times; 10<sup>-30</sup></td><td>1.30 &times; 10<sup>-30</sup></td><td>1.30 &times; 10<sup>-30</sup></td><td>1.30 &times; 10<sup>-30</sup></td></tr><tr><td>10</td><td>1.18 &times; 10<sup>-52</sup></td><td>1.2 &times; 10<sup>-52</sup></td><td>1.2 &times; 10<sup>-52</sup></td><td>1.2 &times; 10<sup>-52</sup></td><td>1.2 &times; 10<sup>-52</sup></td></tr><tr><td>12</td><td><b>9.54 &times; 10<sup>-74</sup></b></td><td>1.09 &times; 10<sup>-73</sup></td><td>1.09 &times; 10<sup>-73</sup></td><td>1.09 &times; 10<sup>-73</sup></td><td>1.09 &times; 10<sup>-73</sup></td></tr><tr><td>14</td><td>1.37 &times; 10<sup>-102</sup></td><td>1.37 &times; 10<sup>-102</sup></td><td>1.37 &times; 10<sup>-102</sup></td><td>1.37 &times; 10<sup>-102</sup></td><td>1.37 &times; 10<sup>-102</sup></td></tr><tr><td>16</td><td>9.73 &times; 10<sup>-138</sup></td><td>9.72 &times; 10<sup>-138</sup></td><td>9.80 &times; 10<sup>-138</sup></td><td>9.8 &times; 10<sup>-138</sup></td><td>9.8 &times; 10<sup>-138</sup></td></tr><tr><td>18</td><td>4.32 &times; 10<sup>-163</sup></td><td>4.32 &times; 10<sup>-163</sup></td><td>4.35 &times; 10<sup>-163</sup></td><td>4.35 &times; 10<sup>-163</sup></td><td>4.35 &times; 10<sup>-163</sup></td></tr><tr><td>20</td><td>1.94 &times; 10<sup>-203</sup></td><td>1.96 &times; 10<sup>-203</sup></td><td>1.96 &times; 10<sup>-203</sup></td><td>1.96 &times; 10<sup>-203</sup></td><td>1.96 &times; 10<sup>-203</sup></td></tr><tr><td>22</td><td>4.92 &times; 10<sup>-255</sup></td><td>4.92 &times; 10<sup>-255</sup></td><td>4.93 &times; 10<sup>-255</sup></td><td>4.93 &times; 10<sup>-255</sup></td><td>4.93 &times; 10<sup>-255</sup></td></tr></tbody></table>

<table><caption>TABLE 7.2 – Temps de calcul de la constante de normalisation pour différentes valeurs de $\epsilon$ dans le cas d'une structure associative</caption><thead><tr><th></th><th>1 &times; 10<sup>-1</sup></th><th>5 &times; 10<sup>-2</sup></th><th>1 &times; 10<sup>-2</sup></th><th>1 &times; 10<sup>-3</sup></th><th>1 &times; 10<sup>-4</sup></th></tr></thead><tbody><tr><td>8</td><td>0.010712</td><td>0.015826</td><td>0.023264</td><td>0.054599</td><td>0.134134</td></tr><tr><td>10</td><td>0.113050</td><td>0.173764</td><td>0.310835</td><td>1.203731</td><td>3.219305</td></tr><tr><td>12</td><td>0.202726</td><td>0.339975</td><td>2.218524</td><td>5.057278</td><td>12.831984</td></tr><tr><td>14</td><td>0.249002</td><td>0.329584</td><td>0.759825</td><td>2.940815</td><td>8.567263</td></tr><tr><td>16</td><td>0.546148</td><td>0.790832</td><td>1.437784</td><td>3.304456</td><td>8.719653</td></tr><tr><td>18</td><td>0.241273</td><td>0.290420</td><td>1.043211</td><td>5.462873</td><td>19.105516</td></tr><tr><td>20</td><td>0.220903</td><td>0.271335</td><td>0.430281</td><td>1.848418</td><td>5.202058</td></tr><tr><td>22</td><td>1.042942</td><td>1.579870</td><td>14.846940</td><td>15.378039</td><td>63.950157</td></tr></tbody></table>

Nous montrons dans les résultats que le rounding avec $\epsilon = 5 \times 10^{-2}$ est largement suffisant pour avoir une bonne précision de la constante de normalisation. De plus, il nous permet d'avoir un très bon temps de calcul comparé aux autres paramètres de rounding.

### 7.3.2 Explication du problème

Dans cet exemple, nous utilisons un modèle SBM avec comme paramètres : $n = 7$, $Q = 3$
$$ \Lambda = \begin{pmatrix} 2 & 10 & 10 \\ 10 & 3 & 10 \\ 10 & 10 & 4 \end{pmatrix} \text{ et } \alpha = (\frac{1}{3}, \frac{1}{3}, \frac{1}{3}). \text{ Nous simulons les classes des individus et la matrice de distances } D \text{ comme détaillé dans la section 6.4.1. Ensuite, nous calculons la constante de normalisation } W \text{ sans faire la compression par rounding puis avec compression comme présenté dans l'algorithme 1. Pour chaque itération des calculs, nous traçons sur les figures 7.4 et 7.5 les rangs des cores.} $$

Pour les tracés de rangs, à l'itération 1, nous traçons les rangs des 21 cores de $B_n$, pour l'itération $i-1$ nous traçons les rangs des 21 cores de $B_n...B_{n-i}$.

Dans la figure 7.4, nous traçons les rangs des cores des TT-matrices calculées à chaque itération de l'algorithme 1 sans la recompression (ligne 4 de l'algorithme). Nous observons que les rangs de la plupart des cores évoluent de manière exponentielle : ils valent $Q^i$ pour l'itération $i$ et atteignent une valeur de $Q^n$ pour la dernière itération. Dans cet exemple, il s'agit de $3^7 = 2187$.
<figure><img src="image_67.png" /><figcaption>FIGURE 7.4 – Évolution des rangs des cores dans un calcul sans rounding pour $n = 7$.</figcaption></figure>

Nous réitérons l'expérience pour les mêmes paramètres et matrice de distance en recompressant le tenseur W à chaque itération (ligne 4 de l'algorithme 1) avec une précision de $\epsilon = 5 \times 10^{-2}$. La figure 7.5 indique l'évolution du rang de chacun des cores en fonction de l'itération et le tableau 7.3 présente les rangs des cores.

<table><caption>TABLE 7.3 – Tableau présentant les rangs des cores des TT-matrices obtenues à chaque itération</caption><thead><tr><th>core</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th><th>13</th><th>14</th><th>15</th><th>16</th><th>17</th><th>18</th><th>19</th><th>20</th><th>21</th></tr></thead><tbody><tr><th>It</th><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>1</th><td>1</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>3</td><td>1</td></tr><tr><th>2</th><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>3</td><td>9</td><td>9</td><td>9</td><td>9</td><td>9</td><td>9</td><td>9</td><td>9</td><td>9</td><td>9</td><td>9</td><td>9</td><td>9</td><td>3</td><td>1</td></tr><tr><th>3</th><td>1</td><td>1</td><td>1</td><td>1</td><td>3</td><td>8</td><td>19</td><td>19</td><td>19</td><td>19</td><td>19</td><td>19</td><td>19</td><td>19</td><td>19</td><td>19</td><td>9</td><td>3</td><td>1</td><td>1</td><td>1</td></tr><tr><th>4</th><td>1</td><td>1</td><td>1</td><td>3</td><td>9</td><td>19</td><td>30</td><td>30</td><td>30</td><td>31</td><td>32</td><td><b>33</b></td><td>17</td><td>9</td><td>3</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><th>5</th><td>1</td><td>1</td><td>3</td><td>3</td><td>7</td><td>10</td><td>12</td><td>12</td><td>9</td><td>5</td><td>3</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><th>6</th><td>1</td><td>3</td><td>6</td><td>6</td><td>5</td><td>3</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><th>7</th><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr></tbody></table>

Dans le tableau 7.3, chaque ligne correspond à l'itération de l'algorithme 1 et chaque colonne au numéro du core. Nous observons que lors de la première itération, les cores sont de rang 1 ou 3. Puis lors des premières itérations, ils augmentent avant de décroitre à partir de la cinquième itération jusqu'à valoir 1 lors de la dernière.

<figure><img src="image_68.png" /><figcaption>FIGURE 7.5 – Évolution des rangs des cores dans un calcul avec rounding à $\epsilon = 0.05$ pour $n = 7$.</figcaption></figure>

Dance le graphique 7.5, nous constatons qu'après chaque itération de l'algorithme 1 le rang
augmente jusqu'à la quatrième itération jusqu'à atteindre 33, soit à peine plus que $Q^3 = 27$.
Ensuite, les rangs des cores baissent jusqu'à atteindre 1 lors de la dernière itération. De plus,
nous remarquons qu'après chaque itération, le nombre de cores de rangs 1 augmente et qu'à la
dernière itération tous les cores sont de rang 1.

Nous comparons désormais les précisions et les temps de calculs deux expériences présentées
dans les graphes 7.4 et 7.5. Nous observons :

* **Erreur relative.** Nous comparons les résultats obtenus pour chaque calcul et nous obte-
nons une erreur relative de $1.31 \times 10^{-7}$ alors que nous souhaitons $5 \times 10^{-2}$.

* **Temps de calcul.** Le temps de calcul pour la première expérience est de 1.09 secondes, tandis qu'il est de 0.05 secondes pour la deuxième expérience.

Cette expérience justifie l'importance de faire des compressions par rounding à chaque itération de les calculs. De plus, quand le nombre d'individus est égal à 8, le calcul de la constante de normalisation devient difficile étant donné que nous devons manipuler des TT-matrices avec des cores de rangs $3^8 = 6561$.

Nous refaisons la même expérience présentée dans le graphe 7.5, mais nous augmentons le
nombre d'individus *n* à 22. Nous aurons un nombre de liens *m* = 231. Nous traçons ensuite
l'évolution du rang des cores en fonction des itérations dans la figure 7.6.

<figure><img src="image_69.png" /><figcaption>FIGURE 7.6 – Évolution des rangs des cores dans un calcul avec compression avec l’algorithme du rounding à $\epsilon = 5 \times 10^{-2}$ pour un modèle SBM avec $n = 22$.</figcaption></figure>
Dans la figure 7.6, nous avons tracé une surface 3D, une représentation tridimensionnelle des rangs des cores (l'axe des *z*) en fonction de son ordre (l'axe des *y*) et des itérations de l'algorithme 1 (l'axe des *x*). Nous faisons varier les couleurs du bleu (pour les valeurs basses des rangs) au rouge (pour les valeurs hautes des rangs).

Nous constatons que durant les premières itérations, les rangs des cores de la TT-matrice résultante croient jusqu'à atteindre 507 lors de la septième itération avant de décroître au fur et à mesure des itérations jusqu'à atteindre 1 pour la plupart des cores lors de la dernière itération. Á partir de cet exemple, nous pouvons dire que si le nombre d'individus augmente, le calcul de la constante de normalisation devient quasiment impossible vu que le stockage nécessaire est trop grand.

Nous avons illustré dans les figures 7.4 et 7.5 l'intérêt de faire une compression par le rounding. Cependant, dès que nous augmentons le nombre d'individus, les rangs explosent, même avec rounding avec comme paramètre $\epsilon$. Et donc, le calcul de la constante de normalisation à cette précision devient très lent voir impossible dès que le nombre d'individus est supérieur à 33. Dans les sous-sections suivantes, nous présentons comment changer les paramètres du rounding afin d'accélérer les calculs sans pour autant perdre en termes de précision.

### 7.3.3 Présentation du rounding

Dans le cadre de l'approche TT, nous n'utilisons que des TT-vecteurs obtenus par produits de TT-matrice TT-vecteurs comme nous l'avons présenté dans la sous-section 5.2.6.

Á chaque itération, nous appliquons le rounding que nous avons expliqué dans la sous-section 4.4.4 et que nous rappelons dans cette sous-section.

Pour une TT-matrice donnée, nous pouvons réduire les rangs de leurs cores avec deux méthodes de recompression par l'algorithme du rounding présentées dans la section 4.4.4.

* **TT-rounding avec précision prescrite $\epsilon$** où nous pouvons contrôler la précision de l'approximation de la TT-matrice, mais pas le rang de ses cores.

* **TT-rounding avec un rang prescrit** où nous pouvons contrôler les rangs des TT-cores de l'approximation de la TT-matrice, mais pas sa précision.

Dans les différentes librairies python que nous avons utilisé dans le cadre de la manipulation des opérations sur les TT-matrices et TT-vecteur, il existe une méthode alternative aux deux premières. Cette méthode prend comme paramètres une précision $\epsilon$ et un rang prescrits et fait une compression des TT-matrices ou des TT-vecteurs en contrôlant le rang tout en gardant une meilleure précision que le TT-rounding avec un rang prescrit.

Nous appliquons cette méthode pour la compression des TT-matrices et nous la comparons avec les deux autres.

## 7.3.4 Résultats

### 7.3.4.1 Calcul de la constante de normalisation

Pour l'ensemble des expériences que nous présentons, nous avons fixé les paramètres suivants : Le nombre de classes $Q = 3$ et les proportions d'individus dans les classes $\alpha = (\frac{1}{3}, \frac{1}{3}, \frac{1}{3})$. Nous faisons varier la structure de la matrice des distances $\Lambda$ parmi les six présentées dans la section 6.4.

Pour chaque structure, nous faisons varier le nombre d'individus *n* entre 8 et 32 par pas de deux. Pour chaque nombre d'individus, nous simulons la matrice de distances *D* et nous appliquons
la procédure de Novikov avec notre procédure d'homothétie pour calculer les TT-matrices $B_i$. Enfin, nous fixons $\epsilon = 5 \times 10^{-2}$ et nous calculons la constante de normalisation par six méthodes. Premièrement avec $\epsilon$ fixé, puis $r_{max} = 9$ et $\epsilon$ fixé, $r_{max} = 3$ et $\epsilon$ fixé, $r_{max} = 9$, $r_{max} = 27$ et $\epsilon$ fixé et enfin $r_{max} = 81$ et $\epsilon$ fixé. Les résultats présentés sont les suivants :

* Premièrement, nous présentons dans les tableaux 7.5 et 7.7 les résultats du calcul de la constante de normalisation.
Nous présentons aussi le plus grand TT-rang max observé durant les différentes itérations de l'algorithme 1 lorsque le rounding prend comme paramètre uniquement la précision $\epsilon$.

* Deuxièmement, nous présentons les temps de calcul pour les différents calculs dans les tableaux 7.6 et 7.8.

**NB** : Nous ne présentons que les structures associative et cores periphery dans ce chapitre, les quatre autres structures sont présentés dans l'annexe A.

Pour une meilleure lisibilité des tableaux, nous mettrons en gras les résultats qui diffèrent de la valeur de référence qui est celle obtenue lorsqu'on fait un calcul avec rounding avec $\epsilon = 5 \times 10^{-2}$ fixé. Pour les discussions et les commentaires des résultats, nous notons les paramètres du rounding de cette manière :

<table><caption>TABLE 7.4 – Tableau des abréviations pour différentes méthodes de rounding.</caption><thead><tr><th>ε</th><th>r<sub>max</sub> = 3 et ε</th><th>r<sub>max</sub> = 9 et ε</th><th>r<sub>max</sub> = 9</th><th>r<sub>max</sub> = 27 et ε</th><th>r<sub>max</sub> = 81 et ε</th></tr></thead><tbody><tr><td>R<sub>1</sub></td><td>R<sub>2</sub></td><td>R<sub>3</sub></td><td>R<sub>4</sub></td><td>R<sub>5</sub></td><td>R<sub>6</sub></td></tr></tbody></table>

## Associative

<table><caption>TABLE 7.5 – Constante de normalisation calculée ainsi que le $r_{max}$ observé pour la méthode $R_1$ dans le cas associative.</caption><thead><tr><th></th><th>R<sub>1</sub> (r<sub>max</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><th>R<sub>6</sub></th></tr></thead><tbody><tr><td>n=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.26 × 10<sup>-28</sup></td><td>6.21 × 10<sup>-28</sup></td><td>6.21 × 10<sup>-28</sup></td></tr><tr><td>n=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.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>n=12</td><td>1.0 × 10<sup>-75</sup> (31)</td><td><b>3.32 × 10<sup>-77</sup></b></td><td><b>6.73 × 10<sup>-76</sup></b></td><td><b>6.78 × 10<sup>-76</sup></b></td><td>1.0 × 10<sup>-75</sup></td><td>1.0 × 10<sup>-75</sup></td></tr><tr><td>n=14</td><td>5.38 × 10<sup>-101</sup> (8)</td><td><b>1.11 × 10<sup>-101</sup></b></td><td>5.38 × 10<sup>-101</sup></td><td>5.39 × 10<sup>-101</sup></td><td>5.38 × 10<sup>-101</sup></td><td>5.38 × 10<sup>-101</sup></td></tr><tr><td>n=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.34 × 10<sup>-135</sup></td><td>3.30 × 10<sup>-135</sup></td><td>3.30 × 10<sup>-135</sup></td></tr><tr><td>n=18</td><td>2.05 × 10<sup>-161</sup> (15)</td><td><b>1.14 × 10<sup>-161</sup></b></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><td>2.05 × 10<sup>-161</sup></td></tr><tr><td>n=20</td><td>1.18 × 10<sup>-199</sup> (18)</td><td><b>2.65 × 10<sup>-205</sup></b></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><td>1.18 × 10<sup>-199</sup></td></tr><tr><td>n=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><td>1.59 × 10<sup>-247</sup></td></tr><tr><td>n=24</td><td>9.57 × 10<sup>-294</sup> (16)</td><td><b>1.12 × 10<sup>-297</sup></b></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><td>9.57 × 10<sup>-294</sup></td></tr><tr><td>n=26</td><td>1.80 × 10<sup>-353</sup> (16)</td><td><b>8.13 × 10<sup>-375</sup></b></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><td>1.80 × 10<sup>-353</sup></td></tr><tr><td>n=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><td>1.34 × 10<sup>-393</sup></td></tr><tr><td>n=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><td>4.03 × 10<sup>-459</sup></td></tr><tr><td>n=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><td>3.43 × 10<sup>-511</sup></td></tr></tbody></table>

En termes de résultats obtenus et de temps de calculs, la structure associative est la structure la plus simple parmi les six structures utilisées pour simuler les matrices de distances. Nous constatons que seul le rounding $R_2$ donne des résultats éloignés de notre référence ($R_1$). En termes de temps de calcul, le rounding $R_3$ est beaucoup plus rapide si nous excluons $R_2$. Par contre, nous remarquons qu'il fait une erreur dans le calcul des de la constante de normalisation pour $n = 12$.
<table><caption>TABLE 7.6 – Temps de calcul en secondes pour W dans le cas associative.</caption><thead><tr><th></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><th>R<sub>6</sub></th></tr></thead><tbody><tr><td>n=8</td><td>0.05</td><td>0.00</td><td>0.01</td><td>0.01</td><td>0.02</td><td>0.02</td></tr><tr><td>n=10</td><td>0.09</td><td>0.01</td><td>0.04</td><td>0.04</td><td>0.04</td><td>0.04</td></tr><tr><td>n=12</td><td>0.28</td><td>0.01</td><td>0.05</td><td>0.06</td><td>0.14</td><td>0.14</td></tr><tr><td>n=14</td><td>0.33</td><td>0.01</td><td>0.07</td><td>0.10</td><td>0.16</td><td>0.16</td></tr><tr><td>n=16</td><td>1.02</td><td>0.02</td><td>0.06</td><td>0.15</td><td>0.48</td><td>0.50</td></tr><tr><td>n=18</td><td>0.52</td><td>0.02</td><td>0.10</td><td>0.22</td><td>0.26</td><td>0.26</td></tr><tr><td>n=20</td><td>0.33</td><td>0.04</td><td>0.14</td><td>0.32</td><td>0.17</td><td>0.17</td></tr><tr><td>n=22</td><td>2.08</td><td>0.04</td><td>0.17</td><td>0.42</td><td>0.81</td><td>0.93</td></tr><tr><td>n=24</td><td>0.41</td><td>0.06</td><td>0.15</td><td>0.55</td><td>0.19</td><td>0.19</td></tr><tr><td>n=26</td><td>1.79</td><td>0.07</td><td>0.38</td><td>0.72</td><td>0.85</td><td>0.83</td></tr><tr><td>n=28</td><td>0.91</td><td>0.07</td><td>0.32</td><td>0.94</td><td>0.42</td><td>0.42</td></tr><tr><td>n=30</td><td>5.52</td><td>0.08</td><td>0.37</td><td>1.16</td><td>2.35</td><td>2.36</td></tr><tr><td>n=32</td><td>2.65</td><td>0.10</td><td>0.32</td><td>1.38</td><td>1.17</td><td>1.17</td></tr></tbody></table>

Au vu de ces résultats, $R_5$ conduit au meilleur rounding pour obtenir des résultats aussi précis que notre référence bien qu'il soit plus coûteux que les autres.

Nous pouvons aussi remarquer que les temps de calcul entre $R_5$ et $R_6$ sont similaires, cela est lié aux valeurs assez faibles du $r_{max}$ comprises entre 9 et 34, cela fait que faire ces rounding reviennent à appliquer $R_1$.

## Cores periphery

<table><caption>TABLE 7.7 – Constante de normalisation calculée ainsi que le $r_{max}$ observé pour la méthode $R_1$ dans le cas cores periphery.</caption><thead><tr><th></th><th>R<sub>1</sub> (r<sub>max</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><th>R<sub>6</sub></th></tr></thead><tbody><tr><td>n=8</td><td>2.63 × 10<sup>-33</sup> (58)</td><td><b>7.16 × 10<sup>-36</sup></b></td><td><b>6.90 × 10<sup>-34</sup></b></td><td><b>7.33 × 10<sup>-34</sup></b></td><td><b>1.55 × 10<sup>-33</sup></b></td><td>2.62 × 10<sup>-33</sup></td></tr><tr><td>n=10</td><td>1.55 × 10<sup>-46</sup> (99)</td><td><b>3.70 × 10<sup>-47</sup></b></td><td><b>9.87 × 10<sup>-47</sup></b></td><td><b>1.01 × 10<sup>-46</sup></b></td><td>1.33 × 10<sup>-46</sup></td><td>1.54 × 10<sup>-46</sup></td></tr><tr><td>n=12</td><td>2.12 × 10<sup>-70</sup> (34)</td><td>1.67 × 10<sup>-70</sup></td><td>2.01 × 10<sup>-70</sup></td><td>2.01 × 10<sup>-70</sup></td><td>2.12 × 10<sup>-70</sup></td><td>2.12 × 10<sup>-70</sup></td></tr><tr><td>n=14</td><td>1.60 × 10<sup>-92</sup> (30)</td><td><b>5.68 × 10<sup>-102</sup></b></td><td>1.57 × 10<sup>-92</sup></td><td>1.57 × 10<sup>-92</sup></td><td>1.60 × 10<sup>-92</sup></td><td>1.60 × 10<sup>-92</sup></td></tr><tr><td>n=16</td><td>5.83 × 10<sup>-129</sup> (43)</td><td><b>1.36 × 10<sup>-137</sup></b></td><td><b>1.38 × 10<sup>-137</sup></b></td><td><b>2.26 × 10<sup>-137</sup></b></td><td>5.83 × 10<sup>-129</sup></td><td>5.83 × 10<sup>-129</sup></td></tr><tr><td>n=18</td><td>2.47 × 10<sup>-154</sup> (38)</td><td><b>5.63 × 10<sup>-171</sup></b></td><td>2.47 × 10<sup>-154</sup></td><td>2.47 × 10<sup>-154</sup></td><td>2.47 × 10<sup>-154</sup></td><td>2.47 × 10<sup>-154</sup></td></tr><tr><td>n=20</td><td>6.88 × 10<sup>-188</sup> (89)</td><td><b>1.27 × 10<sup>-190</sup></b></td><td>5.66 × 10<sup>-188</sup></td><td>5.55 × 10<sup>-188</sup></td><td>6.53 × 10<sup>-188</sup></td><td>6.62 × 10<sup>-188</sup></td></tr><tr><td>n=22</td><td>4.06 × 10<sup>-238</sup> (83)</td><td><b>1.03 × 10<sup>-253</sup></b></td><td><b>1.71 × 10<sup>-252</sup></b></td><td><b>8.75 × 10<sup>-250</sup></b></td><td>3.71 × 10<sup>-238</sup></td><td>4.06 × 10<sup>-238</sup></td></tr><tr><td>n=24</td><td>2.41 × 10<sup>-280</sup> (135)</td><td><b>4.80 × 10<sup>-308</sup></b></td><td>2.41 × 10<sup>-280</sup></td><td>2.46 × 10<sup>-280</sup></td><td>2.41 × 10<sup>-280</sup></td><td>2.41 × 10<sup>-280</sup></td></tr><tr><td>n=26</td><td>5.17 × 10<sup>-324</sup> (136)</td><td><b>4.45 × 10<sup>-334</sup></b></td><td>5.11 × 10<sup>-324</sup></td><td>5.11 × 10<sup>-324</sup></td><td>5.11 × 10<sup>-324</sup></td><td>5.17 × 10<sup>-324</sup></td></tr><tr><td>n=28</td><td>6.96 × 10<sup>-385</sup> (58)</td><td><b>3.08 × 10<sup>-412</sup></b></td><td><b>4.74 × 10<sup>-392</sup></b></td><td><b>4.74 × 10<sup>-392</sup></b></td><td>6.97 × 10<sup>-385</sup></td><td>6.97 × 10<sup>-385</sup></td></tr><tr><td>n=30</td><td>4.90 × 10<sup>-431</sup> (47)</td><td>4.90 × 10<sup>-431</sup></td><td>4.90 × 10<sup>-431</sup></td><td>4.90 × 10<sup>-431</sup></td><td>4.90 × 10<sup>-431</sup></td><td>4.90 × 10<sup>-431</sup></td></tr><tr><td>n=32</td><td>3.45 × 10<sup>-474</sup> (18)</td><td><b>2.63 × 10<sup>-519</sup></b></td><td>3.45 × 10<sup>-474</sup></td><td>3.46 × 10<sup>-474</sup></td><td>3.45 × 10<sup>-474</sup></td><td>3.45 × 10<sup>-474</sup></td></tr></tbody></table>

En termes de résultats obtenus, la structure cores periphery est la structure plus difficile parmi les six structures utilisées pour simuler les matrices de distances. Nous observons que dans cette structure de distances, les rangs maximaux trouvés augmentent jusqu'à atteindre 136 pour $n = 26$, cela induit des temps de calcul qui peuvent attendre à leur tour 379 secondes pour le rounding $R_1$ soit plus de 6 minutes pour ce calcul.

Nous constatons aussi que cette augmentation du rang n'a pas nécessairement d'influence sur les résultats. En effet, le calcul de la constante de normalisation donnent des résultats précis quand nous appliquons une compression par rounding avec les paramètres $R_3$ et $R_4$ pour $n = 24$ (resp. $n = 26$) où
<table><caption>TABLE 7.8 – Temps de calcul en secondes pour W dans le cas cores periphery.</caption><thead><tr><th></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><th>R<sub>6</sub></th></tr></thead><tbody><tr><td>n=8</td><td>0.81</td><td>0.00</td><td>0.01</td><td>0.02</td><td>0.07</td><td>0.23</td></tr><tr><td>n=10</td><td>6.28</td><td>0.01</td><td>0.03</td><td>0.04</td><td>0.33</td><td>2.72</td></tr><tr><td>n=12</td><td>0.88</td><td>0.01</td><td>0.05</td><td>0.07</td><td>0.30</td><td>0.42</td></tr><tr><td>n=14</td><td>1.15</td><td>0.02</td><td>0.09</td><td>0.11</td><td>0.54</td><td>0.59</td></tr><tr><td>n=16</td><td>2</td><td>0.02</td><td>0.13</td><td>0.17</td><td>0.70</td><td>0.84</td></tr><tr><td>n=18</td><td>3.57</td><td>0.03</td><td>0.12</td><td>0.27</td><td>1.33</td><td>2.09</td></tr><tr><td>n=20</td><td>79.80</td><td>0.04</td><td>0.25</td><td>0.36</td><td>3.46</td><td>31.90</td></tr><tr><td>n=22</td><td>79.20</td><td>0.05</td><td>0.31</td><td>0.53</td><td>4.36</td><td>35.80</td></tr><tr><td>n=24</td><td>302</td><td>0.06</td><td>0.41</td><td>0.67</td><td>5.65</td><td>79</td></tr><tr><td>n=26</td><td>350</td><td>0.07</td><td>0.54</td><td>0.79</td><td>7.08</td><td>74.70</td></tr><tr><td>n=28</td><td>379</td><td>0.08</td><td>0.48</td><td>1.02</td><td>5.52</td><td>93.10</td></tr><tr><td>n=30</td><td>11.40</td><td>0.09</td><td>0.34</td><td>1.19</td><td>1.97</td><td>4.28</td></tr><tr><td>n=32</td><td>1.73</td><td>0.12</td><td>0.43</td><td>1.50</td><td>0.77</td><td>0.77</td></tr></tbody></table>

les rangs sont de 135 (resp. 136) tandis que, alors que ces roudings donnent des résultats très éloignés du calcul de référence pour $n = 16$ où $r_{max} = 43$. En ce qui concerne les autres méthodes de calcul, nous remarquons que $R_5$ donne des résultats très précis à part pour $n = 8$. De plus, avec cette méthode, nous faisons chuter le temps de calcul de 379 secondes à 5.52 secondes pour $n = 28$ contre 93.10 secondes (1 minute 30 seconde) pour la méthode $R_6$.

### 7.3.4.2 Réitération des calculs

Dans le cadre des expériences précédentes, nous avons fixé $\alpha$ et $\Lambda$ et $Q$. De ce fait, la constante de normalisation $W$ ne dépend plus que de la matrice distance $D$ qui est générée aléatoirement et peut conduire entre 2 matrices différentes à de fortes variations sur $W$ et les $r_{max}$ observés lors des itérations de l'algorithme 1. Nous avons remarqué que pour certaines expériences, il y a de grandes différences entre la constante de normalisation calculée avec un $\epsilon$ et celles calculées avec $\epsilon$ et $r_{max}$. Afin de visualiser cela, nous générons 50 matrices de distance par structure de distance pour deux nombres d'individus $n = 20$ et $n = 30$ et nous comparons les résultats obtenus pour différents choix du rounding avec le calcul de la méthode $R_1$ (rounding avec précision) qui devient la référence. Dans les deux tableaux 7.10 et 7.11, nous présenterons en ligne les structures de distances et en colonnes les paramètres du rounding. Chaque colonne donne le nombre de fois où le calcul avec la méthode choisie donne un résultat différent de celui obtenu avec $R_1$ avec une erreur relative supérieur à 0.1.

Á noter que quand l'erreur relative est supérieure à 0.1, elle est en fait de l'ordre de 1. Donc les résultats sont très loin de la référence.

**NB** : Pour une meilleure lisibilité des tableaux, nous mettrons en gras le nombre de fois où une méthode donne un résultat différent de notre référence. Pour les discussions et les commentaires des résultats, nous notons les paramètres du rounding de cette manière.

<table><caption>TABLE 7.9 – Tableau des abréviations pour différentes méthodes de rounding 2.</caption><thead><tr><th></th><th>r<sub>max</sub> = 9 et &epsilon;</th><th>r<sub>max</sub> = 9</th><th>&epsilon; et r<sub>max</sub> = 27</th><th>r<sub>max</sub> = 50 et &epsilon;</th><th>r<sub>max</sub> = 81 et &epsilon;</th></tr></thead><tbody><tr><td></td><td>R<sub>3</sub></td><td>R<sub>4</sub></td><td>R<sub>5</sub></td><td>R<sub>5bis</sub></td><td>R<sub>6</sub></td></tr></tbody></table>
<table><caption>TABLE 7.10 – Résultats pour les calculs de la constante de normalisation pour de 50 matrices de distances par structure de SBM avec un nombre d'individus $n = 20$.</caption><thead><tr><th></th><th colspan="5">Paramètre du rounding</th></tr><tr><th></th><th>R<sub>4</sub></th><th>R<sub>3</sub></th><th>R<sub>5</sub></th><th>R<sub>5bis</sub></th><th>R<sub>6</sub></th></tr></thead><tbody><tr><th>associative</th><td><b>1</b></td><td><b>1</b></td><td>0</td><td>0</td><td>0</td></tr><tr><th>Hiérarchique</th><td><b>2</b></td><td><b>2</b></td><td>0</td><td>0</td><td>0</td></tr><tr><th>ordered</th><td><b>2</b></td><td><b>2</b></td><td>0</td><td>0</td><td>0</td></tr><tr><th>dissociative</th><td><b>9</b></td><td><b>7</b></td><td>0</td><td>0</td><td>0</td></tr><tr><th>cores periphery</th><td><b>3</b></td><td><b>2</b></td><td><b>1</b></td><td>0</td><td>0</td></tr><tr><th>Hiérarchique 2</th><td><b>7</b></td><td><b>3</b></td><td>0</td><td>0</td><td>0</td></tr></tbody></table>

<table><caption>TABLE 7.11 – Résultats pour les calculs de la constante de normalisation pour de 50 matrices de distances par structure de SBM avec un nombre d'individus $n = 30$.</caption><thead><tr><th></th><th colspan="5">Paramètre du rounding</th></tr><tr><th></th><th>R<sub>4</sub></th><th>R<sub>3</sub></th><th>R<sub>5</sub></th><th>R<sub>5bis</sub></th><th>R<sub>6</sub></th></tr></thead><tbody><tr><th>associative</th><td><b>3</b></td><td><b>2</b></td><td>0</td><td>0</td><td>0</td></tr><tr><th>Hiérarchique</th><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><th>Ordered</th><td><b>2</b></td><td><b>1</b></td><td>0</td><td>0</td><td>0</td></tr><tr><th>dissociative</th><td><b>5</b></td><td><b>3</b></td><td>0</td><td>0</td><td>0</td></tr><tr><th>cores periphery</th><td><b>8</b></td><td><b>5</b></td><td>0</td><td>0</td><td>0</td></tr><tr><th>Hiérarchique 2</th><td><b>2</b></td><td><b>1</b></td><td>0</td><td>0</td><td>0</td></tr></tbody></table>

Dans un premier temps, nous observons que le nombre de calculs égaux au calcul avec $R_1$ pour les diverses structures de distances est supérieur à 95% à part pour les cas dissociative pour $n = 20$ et cores périphérie pour $n = 30$ où la précision des calculs avec un rounding ayant comme paramètres $r_{max} = 9$ et $r_{max} = 9$ et $\epsilon$ où le nombre de calculs exacts est compris entre 80% et 90%.

Si nous comparons les tableaux 7.10 et 7.11, nous constatons que les résultats pour $n = 30$ sont légèrement meilleurs que ceux pour $n = 20$. Nous pouvons dire que plus le nombre d'individus augmente, plus l'approximation de la constante de normalisation a tendance à être proche de notre référence.

Dans un second temps, nous constatons que l'approximation avec rounding avec comme paramètres $\epsilon$ et $r_{max} = 27$ suffit largement pour obtenir un résultat précis presque tout le temps. De plus, il est plus rapide que ceux considérons comme paramètres $\epsilon$ et $r_{max} = 50$ et $\epsilon$ et $r_{max} = 81$.

En vue des résultats présentés dans cette section, nous pouvons présumer que le meilleur choix des paramètres afin d'améliorer le rapport précision du résultat et temps de calcul est le rounding avec comme paramètres $r_{max} = Q^3 = 27$ et $\epsilon$.

### 7.3.4.3 Comparaison des marginales

Nous avons montré durant les sections 6.5 et 6.6 que le calcul des marginales unaires et binaires avec l'approche TT est très précis lorsque nous faisons les calculs avec une précision $\epsilon = 5 \times 10^{-2}$. Dans cette section, nous comparons les calculs des marginales avec des compressions avec un rounding de paramètre $\epsilon$ ($R_1$) et une autre avec un rounding de paramètres $r_{max} = Q^3$ et $\epsilon$ ($R_3$).

Nous commençons par choisir deux structures associative et dissociative puis $n = 18$ puis nous simulons un modèle SBM de la même manière qu'auparavant. Nous calculons ensuite les marginales et la constante de normalisation avec un rounding $R_1$ puis avec un rounding $R_3$. Enfin, dans les figures 7.7 et 7.8 nous présentons les résultats pour le cas associative avec un pairplot (expliqué dans la section 6.4.3). Dans les figures 7.9 et 7.10 nous présentons les résultats pour le cas dissociative.
<figure><img src="image_70.png" /><figcaption>FIGURE 7.7 – Comparaison des marginales unaires calculées avec rounding $R_1$ et $R_3$ pour le cas associative.</figcaption></figure>

<figure><img src="image_71.png" /><figcaption>FIGURE 7.8 – Comparaison des marginales binaires calculées avec rounding $R_1$ et $R_3$ pour le cas associative.</figcaption></figure>

<figure><img src="image_72.png" /><figcaption>FIGURE 7.9 – Comparaison des marginales unaires avec rounding $R_1$ et $R_3$ pour le cas dissociative</figcaption></figure>

<figure><img src="image_73.png" /><figcaption>FIGURE 7.10 – Comparaison des marginales binaires avec rounding $R_1$ et $R_3$ pour le cas dissociative</figcaption></figure>

Dans ces deux exemples, les calculs de la constante de normalisation donnent les mêmes résultats et les marginales sont très proches (la corrélation est de 1). Nous présentons dans le tableau 7.12 le temps de calcul des marginales unaires et binaires pour les deux modèles.

<table><caption>TABLE 7.12 – Temps de calcul des marginales en secondes pour un calcul avec rounding $R_3$ comparé à un calcul avec rounding $R_1$</caption><thead><tr><td>Paramètre du rounding</td><td>R<sub>3</sub></td><td>R<sub>1</sub></td></tr></thead><tbody><tr><td>associative</td><td>35</td><td>82</td></tr><tr><td>dissociative</td><td>91</td><td>201</td></tr></tbody></table>

### 7.3.5 Discussion

**Limite du rounding avec précision.** Pour un modèle SBM, le calcul avec des compressions de TT-matrices par un rounding classique (avec une précision fixée $\epsilon$) ne permet ni de faire un calcul de la constante de normalisation si le nombre d'individus est élevé ($n > 33$) ni celui des marginales
$n > 26$.

Cela peut s'expliquer par le fait qu'utiliser le rounding classique permet de contrôler la précision $\epsilon$ mais ne permet pas de contrôler les rangs des cores durant les calculs intermédiaires, ce qui induit a des stockages de cores de rangs élevés, ce qui soit augmente les temps de calcul soit rend les calculs impossibles.

Pour palier à cette difficulté, nous avons changé les paramètres en entrées du rounding ce qui modifie la ligne 4 de l'algorithme 1 ainsi :

$$W = \text{rounding}(W, \epsilon, r_{max})$$

**Influence des paramètres du rounding** À partir des différents tableaux 7.7, 7.8, 7.5, 7.6, 7.10, 7.11 et ceux présentés dans l'annexe A, nous pouvons dire que pour les six structures de distances choisies, à part cores periphery 7.8 et dissociative A.2, le meilleur choix des paramètres du rounding en termes de temps de calcul est celui prenant $r_{max} = 9$ et $\epsilon = 5 \times 10^{-2}$ étant donné qu'il donne des résultats assez précis et le deuxième meilleur temps de calcul. Cependant, nous remarquons qu'il peut donner parfois des résultats peu précis pour certains calculs de la constante de normalisation comme nous l'avons présenté dans les tableaux 7.10 et 7.11.

Ainsi, afin d'avoir le meilleur rounding en termes de rapport précision temps de calcul pour les six structures, il est plus adéquat d'utiliser un $r_{max} = 27$ et $\epsilon = 5 \times 10^{-2}$ qui permet par exemple de faire le calcul de la constante de normalisation avec une excellente précision avec un temps de calcul assez faible (5.52 secondes au lieu de 379 pour un calcul avec une précision $\epsilon$).

**Choix optimal des paramètres du rounding** A travers les différents résultats que nous avons présentés, nous avons remarqué que plus le $r_{max}$ augmente plus la précision des calculs est meilleure jusqu'à atteindre un point de saturation quand $r_{max} = Q^3$.

Le meilleur choix en termes de rapport temps de calcul/précision est de prendre un $\epsilon = 0.05$ et $r_{max} = Q^3$.

### 7.3.6 Conclusion

Durant cette section, nous avons présenté comment évoluent les rangs dans les différentes itérations du calcul des marginales et de la constante de normalisation. Nous avons ensuite expliqué comment limiter cette évolution par l'intermédiaire d'une compression par l'algorithme classique du rounding puis comment changer les paramètres de cet algorithme pour limiter le stockage et accélérer les calculs. Enfin, nous avons présenté des résultats expérimentaux de comparaison de plusieurs calculs avec différents paramètres en entrées sur différentes structures de distances.

Nous avons ajouté cette procédure dans la méthode de calcul présentée dans la section 5.3 et nous l'avons associée à l'homothétie que nous avons présenté dans la section 7.2 permis de faire les calculs de la constante de normalisation et des marginales pour un modèle avec un nombre d'individus égale à 44.

Cependant, nous avons présenté dans les tableaux 7.10 et 7.11 que dans certains cas, les calculs de la constante de normalisation donne des résultats peu précis.

## 7.4 Réduction des dimensions de TT-matrices de rang 1 par fusion des cores

### 7.4.1 Motivation

Nous avons expliqué précédemment que pour un modèle SBM avec $n$ individus, nous manipulons des TT-matrices avec $m = \frac{n(n-1)}{2}$ cores. Cela a plusieurs répercussions sur les calculs.
Tout d'abord, au niveau de la complexité, plus le nombre de cores augmente plus la longueur de la boucle sur les dimensions et le temps des calculs augmentent, que ce soit lors de la construction des $\mathbf{B}_i$ (formule (5.2.33)), pour les produits des TT-matrices TT-vecteurs $\mathbf{B}_i$ et pour le rounding nécessaires au calcul de la constante des marginales et de la constate de normalisation.

De plus, dans la section 5.3, nous montrons que le calcul des marginales binaires nécessite un nombre quadratique ($O(Q^2n^2)$, équation 5.3.16) de produits TT-matrice-TT-matrice.

Chacun de ces produits de TT-matrices impose $O(n^2)$ factorisations $LQ$ et décompositions SVD, donc, le calcul de la constante de normalisation (resp. le calcul les marginales binaires) nécessite $O(n^3)$ (resp. $O(Q^2n^4)$) factorisations $LQ$ et décompositions SVD.

D'autre part, au niveau de l'implémentation, la majorité des bibliothèques qui implémentent ces factorisations sont optimisées pour des matrices de grandes tailles. Or, les cores des TT-matrices $\mathbf{A}_i[z_i]$ sont de très petites tailles comme nous avons expliqués dans la formule (5.2.24).

De plus, nous avons implémenté les calculs en utilisant la librairie ttpy (Tenseur Train (TT) toolbox en Python) développée en 2015 par Oseledets [OTZ11]. La bibliothèque impose une limite de 1024 pour la dimension des tenseurs au format TT. Le nombre d'individus du modèle SBM est donc limité à $n = 45$ pour pouvoir construire les TT-matrices $\mathbf{A}_i[z_i]$.

Afin de résoudre ce problème, nous proposons plusieurs méthodes de fusion des cores pour réduire les dimensions des TT-matrices $\mathbf{A}_i[z_i]$.

## 7.4.2 Explication de la fusion des cores

### 7.4.2.1 Réduction du nombre des cores d'une TT-matrice de rang 1

Pour expliquer la procédure de fusion des cores, nous commençons par l'illustrer sur un exemple où nous considérons une TT-matrice $\mathbf{A}$ d'ordre 6 ayant pour dimensions $(n_0, n_1, \dots, n_6) \times (m_0, m_1, \dots, m_6)$ et dont les cores sont de rangs $(r_0, r_1, \dots, r_6)$. Les éléments de $\mathbf{A}$ peuvent s'exprimer :

$$ \mathbf{A}[i_1, \dots, i_6; j_1, \dots, j_6] = \mathbf{G}_1(i_1, j_1) \cdot \mathbf{G}_2(i_2, j_2) \dots \mathbf{G}_6(i_6, j_6) \quad (7.4.1) $$

avec $\mathbf{G}_k \in \mathbb{R}^{r_{k-1} \times n_k \times m_k \times r_k}$ le $k^e$ core de $\mathbf{A}$.

Si $\mathbf{A}$ est de TT-rang 1, les rangs de ces TT-cores sont toujours égaux à 1, En développant les produits matriciels, $\mathbf{A}$ s'écrit :

$$ \mathbf{A}[i_1, \dots, i_6; j_1, \dots, j_6] = \sum_{\alpha_1=1}^{r_1} \dots \sum_{\alpha_5=1}^{r_5} \mathbf{G}_1(1, i_1, j_1, \alpha_1) \mathbf{G}_2(\alpha_1, i_2, j_2, \alpha_2) \dots \mathbf{G}_6(\alpha_5, i_6, j_6, 1) $$

qui se simplifie en un produit de scalaire

$$ \mathbf{A}[i_1, \dots, i_6; j_1, \dots, j_6] = \mathbf{G}_1(1, i_1, j_1, 1) \mathbf{G}_2(1, i_2, j_2, 1) \dots \mathbf{G}_6(1, i_6, j_6, 1) \quad (7.4.2) $$

Supposons que nous souhaitons fusionner les cores de $\mathbf{A}$ par groupe de deux, c'est-à-dire trouver une TT-matrice $\mathbf{B}$ d'ordre 3 qui permet d'exprimer les éléments de $\mathbf{A}$ :

$$ \mathbf{B}[\overline{i_1 i_2}, \overline{i_3 i_4}, \overline{i_5 i_6}; \overline{j_1 j_2}, \overline{j_3 j_4}, \overline{j_5 j_6}] = \mathbf{A}[i_1, \dots, i_6; j_1, \dots, j_6] \quad (7.4.3) $$

Nous utilisons la notion de multiplexage présentée dans la section 4.4.1 et la fonction $\rho$ définie dans la formule (4.4.3) pour exprimer les indices de $\mathbf{B}$ de cette manière :

$$ k_1 = \rho((i_1, i_2)) = \overline{i_1 i_2}, \quad k_2 = \rho((i_3, i_4)) = \overline{i_3 i_4}, \quad k_3 = \rho((i_5, i_6)) = \overline{i_5 i_6}, \quad l_1 = \rho((j_1, j_2)) = \overline{j_1 j_2}, \quad l_2 = \rho((j_3, j_4)) = \overline{j_3 j_4}, \quad l_3 = \rho((j_5, j_6)) = \overline{j_5 j_6} $$

En développant $\mathbf{B}$ nous obtenons :

$$ \mathbf{B}[k_1, k_2, k_3; l_1, l_2, l_3] = \sum_{\alpha'_1=1}^{r'_1} \sum_{\alpha'_2=1}^{r'_2} \mathbf{C}_1(1, k_1, l_1, \alpha'_1) \mathbf{C}_2(\alpha'_1, k_2, l_2, \alpha'_2) \mathbf{C}_3(\alpha'_2, k_3, l_3, 1) \quad (7.4.4) $$
où les TT-cores **C**, de **B** s'expriment par produit de Kronecker des TT-cores initiaux **G** :

$$
\begin{align}
\mathbf{C}_1(k_1, l_1) &= \mathbf{G}_1(i_1, j_1) \otimes_k \mathbf{G}_2(i_2, j_2), \nonumber \\
\mathbf{C}_3(k_3, l_3) &= \mathbf{G}_5(i_5, j_5) \otimes_k \mathbf{G}_6(i_6, j_6) \tag{7.4.5}
\end{align}
$$

Puisque les cores de la TT-matrice **B** sont calculés par produits de Kronecker des cores de la TT-matrice **A**. Alors, ils vérifient les propriétés suivantes :

**Proposition 9.** Les TT-matrices **B** sont de TT-rang 1 :

$$
\mathrm{rang}(\mathbf{C}_1) = \mathrm{rang}(\mathbf{C}_2) = \mathrm{rang}(\mathbf{C}_3) = 1
$$

Démonstration. Comme les TT-matrices **A** sont de TT-rang 1 alors :

$$
\mathbf{G}_1(i_1, j_1) \cdot \mathbf{G}_2(i_2, j_2) \cdots \mathbf{G}_6(i_6, j_6) \in \mathbb{R}
$$

Les cores de **B** sont calculés par produits de Kronecker des cores de **A**. Nous utilisons la propriété de Kronecker qui stipule que pour deux matrices A et B nous avons : $rang(A \otimes_k B) = rang(A) \times rang(B)$. Donc,

$$
\left\{
\begin{array}{l}
\mathrm{rang}(\mathbf{C}_1(k_1, l_1)) = \mathrm{rang}(\mathbf{G}_1(i_1, j_1)) \times \mathrm{rang}(\mathbf{G}_2(i_2, j_2)) = 1 \\
\mathrm{rang}(\mathbf{C}_2(k_2, l_2)) = \mathrm{rang}(\mathbf{G}_3(i_3, j_3)) \times \mathrm{rang}(\mathbf{G}_4(i_4, j_4)) = 1 \\
\mathrm{rang}(\mathbf{C}_1(k_3, l_3)) = \mathrm{rang}(\mathbf{G}_5(i_5, j_5)) \times \mathrm{rang}(\mathbf{G}_6(i_6, j_6)) = 1
\end{array}
\right.
\qquad (7.4.6)
$$

Donc les TT-matrices **B** sont de TT-rang 1.

Cette procédure peut facilement se généraliser pour des TT-matrices de TT-rang 1 d'ordre quelconque.

**Proposition 10 (Réduction des dimensions TT-matrices).** *Soit A une TT-matrice d'ordre m et de TT-rang 1. Nous pouvons réduire l'ordre de cette TT-matrice vers un ordre d < m .*

$$
\mathbf{B}[k_1, \dots, k_d; l_1, \dots, l_d] = \mathbf{A}[i_1, \dots, i_m; j_1, \dots, j_m] \quad (7.4.7)
$$

Les cores de la TT-matrices **B** est fait par produit de Kronecker des cores de **A** comme présenté dans
la formule (7.4.5).

**Proposition 11** (Réduction des dimensions de TT-matrices de rang 1). *Soit A une TT-matrice de rang 1, alors, si nous réduisons ses dimensions, la TT-matrice résultante est de rang 1.*

La démonstration de cette propriété vient de la généralisation de la formule (7.4.6).

Dans le cadre des modèles graphiques, la procédure de Novikov [NROV14] génère des TT-matrices
**A**<sub>*i*</sub>[*z*<sub>*i*</sub>] sont des TT-matrices de TT-rang 1. Donc la procédure de fusion peut être appliquée sans aucune
réserve.

Comme les TT-matrices **B**<sub>*i*</sub> sont construites par somme des TT-matrices **A**<sub>*i*</sub>[*z*<sub>*i*</sub>]. Alors, les rangs des
cores sont *r*<sub>1</sub> = ... = *r*<sub>*d*−1</sub> = *Q* et *r*<sub>0</sub> = *r*<sub>*d*</sub> = 1. Nous pouvons donc appliquer la procédure de fusion sur
les TT-matrices **B**<sub>*i*</sub> de manière à vérifier la condition (2.). Dans la suite de cette section, nous allons
présenter cette procédure sur des TT-matrices de TT-rang 1.

7.4.3 Fusion des cores de la TT-matrice Aᵢ[zᵢ]

L'idée derrière cette méthode de fusion consiste à transformer la TT-matrice $\mathbf{A}_i[z_i]$ d'ordre $m$ en une
TT-matrice d'ordre $d$ avec $d < m$.

Pour cela, nous découpons l'intervalle $[1, m]$ en $d$ sous paquets $(P_k)_{k \in [1, d]}$. Puis, nous construisons
les nouveaux cores comme produits de Kronecker des cores d'un paquet. Pour la TT-matrice qui
s'exprime :

$$
\mathbf{A}_i[z_i][i_1, \dots, i_m; j_1, \dots, j_m] = \mathbf{G}_i^1[z_i](i_1, j_1) \mathbf{G}_i^2[z_i](i_2, j_2) \dots \mathbf{G}_i^m[z_i](i_m, j_m) \quad (7.4.8)
$$
avec $\mathbf{G}_i^\ell[z_i] \in \mathbb{R}^{1 \times n_\ell \times m_\ell \times 1}$ le $\ell^e$ core de $\mathbf{A}_i[z_i]$, la transformation permet de la réécrire :

$$ \mathbf{A}_i[z_i][l_1, \dots, l_d; k_1, \dots, k_d] = \mathbf{C}_i^1[z_i](l_1, k_1) \mathbf{C}_i^2[z_i](l_2, k_2) \dots \mathbf{G}_i^{d+1}[z_i](l_d, k_d) \quad (7.4.9) $$

avec,

$$ \forall f \in [1, d+1], \quad \mathbf{C}_i^f[z_i] = \bigotimes_{\ell \in P_f} \mathbf{G}_i^\ell[z_i] \quad (7.4.10) $$

avec : $\mathbf{C}_i^f[z_i] \in \mathbb{R}^{1 \times (\Pi_{\ell \in P_f} n_\ell) \times (\Pi_{\ell \in P_f} m_\ell) \times 1}$

**Définition 2 (Volume d'un core).** Le volume d'un core est défini comme le logarithme du produit des dimensions des cores.

Si nous considérons un core $\mathbf{G}_i^\ell[z_i] \in \mathbb{R}^{r_{\ell-1} \times n_\ell \times m_\ell \times r_\ell}$, son volume s'exprime :

$$ V_i^\ell = \log(r_{\ell-1} \times m_\ell \times m_\ell \times r_\ell) = \log(r_{\ell-1}) + \log(n_\ell) + \log(m_\ell) + \log(r_\ell) \quad (7.4.11) $$

Alors $\forall \ell \in [1, m]$ le volume de chaque cores de $\mathbf{A}_i[z_i]$ s'exprime,

$$ V_i^\ell = \log(n_\ell) + \log(m_\ell) \quad (7.4.12) $$

### 7.4.3.1 Fusion avec un nombre uniforme de cores

La fusion avec un nombre uniforme de cores consiste en un découpage uniforme de l'intervalle $[1, m]$ en $d$ intervalles égaux et un dernier intervalle avec le reste.

Soit $p$ la taille d'un paquet, alors le nombre de paquets est lié à $m$ par la relation suivante :

$$ m = dp + s \quad \text{avec} \quad s \in [0, p-1] \quad (7.4.13) $$

Dans ce cas, la nouvelle TT-matrice contiendra un nombre de cores égale à $d = \lceil \frac{m}{p} \rceil + 1$ si $s \neq 0$ et $d = \lceil \frac{m}{p} \rceil$ si $s = 0$.

Dans la figure 7.11, nous indiquons la distribution du nombre de cores, $P_f$ que contiendra chaque nouveau core $\mathbf{C}_i^k[z_i]$ dans le cadre $d = 6$.

<figure><img src="image_74.png" /><figcaption>FIGURE 7.11 – Répartition uniforme en 6 paquets des $m$ cores</figcaption></figure>

**Procédure de fusion :** Considérons $\mathbf{A}_i[z_i]$ une TT-matrice d'ordre $m$ et de TT-rang 1. Soit $p \in [1, m]$.

La procédure de la fusion avec un nombre uniforme de paquet consiste à :

1. Définir l'ensemble des paquets $P_f$

$$ \forall f \in [1, d], \quad P_f = [(f-1)p+1, \min(fp, m)] \quad (7.4.14) $$

où $d$ est le nouvel ordre de la TT-matrice $\mathbf{A}_i[z_i]$.
2. Chaque nouveau core est formé par le produit de Kronecker de $p$ cores, le dernier core est formé par produit de Kronecker des $s$ derniers cores.

Construire les cores $C_i^f[z_i]$ de la TT-matrice $A_i[z_i]$ d'ordre $d$ qui s'exprime :

$$A_i[z_i][l_1, \dots, l_d; k_1, \dots, k_d] = C_i^1[z_i](l_1, k_1) C_i^2[z_i](l_2, k_2) \dots G_i^{d}[z_i](l_d, k_d) G_i^{d+1}[z_i](l_{d+1}, k_{d+1}) \quad (7.4.15)$$

où $C_i^f[z_i]$ s'exprime selon la formule (7.4.10)

**NB** : Le choix de $p$ doit être fait de sorte que la taille des cores ne soit pas trop grande. Pour ce choix, nous préconisons un entier $p \in [2, 7]$ si $Q = 3$. Si nous considérons $p = 8$, la taille des cores peut atteindre dans le pire des cas une taille de $6561 \times 6561$ ce qui est difficilement manipulable car cela demandera beaucoup de mémoire.

**Avantages et inconvénients de cette méthode** : L'avantage majeur de cette méthode est qu'elle nous permet de passer d'une TT-matrice d'ordre $m$ à une TT-matrice d'ordre $d$ plus petit avec une procédure très simple.

L'inconvénient majeur est que les cores qui composent les TT-matrices ne sont pas tous de mêmes tailles. Donc si nous appliquons cette procédure, nous aurons une mauvaise distribution des volumes des cores. En effet, certains cores auront un grand volume tandis que d'autres auront un volume nul. Nous illustrons les difficultés qui apparaissent dans cette méthode par un exemple où nous considérons un modèle SBM<fnref n="1" /> avec $n = 60$ (soit un nombre de liens $m = 1770$). Nous traçons dans le graphique 7.12 le nombre de cores identités de tailles $1 \times Q \times Q \times 1$ en bleu et le nombre de cores de tailles $1 \times 1 \times 1 \times 1$ en orange.

<figure><img src="image_75.png" /><figcaption>FIGURE 7.12 – Évolution du nombre de cores par rapport à l'évolution de $i$ l'indice des TT-matrices $A_i[z_i]$ pour un modèle SBM avec $n = 60$ ($m = 1770$)</figcaption></figure>

Nous constatons dans ce graphique que les cores de grandes tailles sont concentrés dans les TT-matrices $A_i[z_i]$ où $i$ est proche de $\frac{n}{2}$. Ce nombre atteint $870 \sim \frac{m}{2}$ pour la TT-matrice $A_{\frac{n}{2}}[z_{\frac{n}{2}}]$.

Nous remarquons aussi le grand déséquilibre pour les premières et dernières TT-matrices $A_i[z_i]$ où nous avons très peu de cores identités mais beaucoup de cores $1 \times 1 \times 1 \times 1$. Le maximum vaux 1711 et est atteint pour $A_1[z_1]$ et $A_n[z_n]$.

Cela engendre des problèmes lorsque nous faisons une fusion avec un nombre uniforme de cores car nous serons amené à fusionner soit des cores de très petites tailles soit des cores de très grandes tailles.

<footnote n="1">1. A noté que les tailles des cores des TT-matrices $A_i[z_i]$ ne dépendent pas de la structure de distance choisie</footnote>
### 7.4.3.2 Fusion avec volumes uniformes des cores

Nous avons expliqué durant la procédure de fusion dans la sous-section 7.4.3.1 que nous devons choisir judicieusement le nombre de cores à fusionner afin d'éviter les limites de dépassement de mémoire. Dans cette sous-section, nous allons désormais nous intéresser à la taille des cores des TT-matrices $A_i[z_i]$.

Nous avons précisé les tailles des cores dans la formule (5.2.24). À partir de cette formule, nous pouvons pressentir que pour un modèle SBM avec $n$ individus, la TT-matrice qui contiendra les plus grands cores est la TT-matrice $A_{\frac{n}{2}}$.

Nous nous sommes intéressés au nombre exact des cores de dimensions $1 \times Q \times Q \times 1$. Pour cela, nous proposons la propriété 12.

**Proposition 12.** Soit $n \in \mathbb{N}$, la TT-matrice $A_i[z_i]$, $i = 1, \dots, n$ est constituée de :

* $(P_1)$ $i-1$ cores de tailles $1 \times Q \times 1 \times 1$

* $(P_2)$ $n-i$ cores de tailles $1 \times 1 \times Q \times 1$

* $(P_3)$ $(i-1)(n-i)$ cores de tailles $1 \times Q \times Q \times 1$

* $(P_4)$ $\frac{(n-1)(n-2)}{2} - (i-1)(n-i)$ cores de tailles $1 \times 1 \times 1 \times 1$

Cette proposition se démontre de la même manière que la proposition 5.

*Démonstration.* Soit $n \in \mathbb{N}$, $m = \frac{n(n-1)}{2}$ et $i \in [1, n]$. La matrice $A_i[z_i]$ s'exprime :

$$A_i[z_i] = G_i^{1,2}[z_i] \bullet_{\alpha_1} G_i^{1,3}[z_i] \bullet_{\alpha_2} \dots \bullet_{\alpha_{d-2}} G_i[z_i]^{n-2,n} \bullet_{\alpha_{d-1}} G_i[z_i]^{n-1,n} \quad (7.4.16)$$

$$ \text{Avec, } \left| \begin{array}{llll} i < k & \Rightarrow & G_i^{k,l}[z_i] = 1 & \in \mathbb{R}^{1 \times 1 \times 1 \times 1} \\ i = k & \Rightarrow & G_i^{k,l}[z_i] = G_i^{i,l}[z_i] & \in \mathbb{R}^{1 \times 1 \times Q \times 1} \\ k < i < \ell & \Rightarrow & G_i^{k,l}[z_i] = \mathbb{I}_{Q,Q} & \in \mathbb{R}^{1 \times Q \times Q \times 1} \\ i = l & \Rightarrow & G_i^{k,l}[z_i] = G_i^{k,l}[z_i] & \in \mathbb{R}^{1 \times Q \times 1 \times 1} \\ i > \ell & \Rightarrow & G_i^{k,l}[z_i] = 1 & \in \mathbb{R}^{1 \times 1 \times 1 \times 1} \end{array} \right. \quad (7.4.17) $$

Le nœud $i$ apparaît dans les liens suivants :

* Par la droite : $1 \to i, 2 \to i, \dots, i-1 \to i$

* Par la gauche : $i \to i+1, i \to i+2, \dots, i \to n$

À partir de ces constats, nous pouvons déterminer les emplacements des différents cores.

* Quand $i$ apparaît à droite d'une connexion, cela donne des cores de taille $1 \times Q \times 1 \times 1$ : Nous observons que nous avons $i-1$ connexions par la droite et donc $i-1$ cores de taille $1 \times Q \times 1 \times 1$. Ce qui donne $(P_1)$.

* Quand $i$ apparaît à gauche d'une connexion, cela donne des cores de taille $1 \times 1 \times Q \times 1$ : Comme nous avons $n-i$ connexions à gauche donc $n-i$ cores de taille $1 \times 1 \times Q \times 1$. D'où $(P_2)$.

* À partir des deux premières propriétés, nous avons $m - (n-1) = \frac{(n-1)(n-2)}{2}$ cores scalaires ou des identités. Nous utilisons la formule (7.4.17) pour remarquer que les cores identités sont :

$$ \underbrace{G_i^{1,i+1}[z_i], G_i^{1,i+2}[z_i], \dots, G_i^{1,n}[z_i]}_{n-i}, \underbrace{G_i^{2,i+1}[z_i], G_i^{2,i+2}[z_i], \dots, G_i^{2,n}[z_i]}_{n-i}, \dots \underbrace{G_i^{i-1,i+1}[z_i], \dots, G_i^{i-1,n}[z_i]}_{n-i} $$

Donc, il y aurait tout simplement $(i-1)(n-i)$ cores identités. Ce qui conduit à $(P_3)$.

* À partir de $(P_1)$, $(P_2)$ et $(P_3)$, nous avons $\frac{(n-1)(n-2)}{2} - i(n-i)$ cores de tailles $1 \times 1 \times 1 \times 1$.

□
**Proposition 13.** La TT-matrice $\mathbf{A}_{\frac{n}{2}}[z_{\frac{n}{2}}]$ est celle qui contient le plus de de cores de tailles $1 \times Q \times Q \times 1$.

**Démonstration.** Les matrices $\mathbf{A}_i[z_i]$ contiennent $(i-1)(n-i)$ cores de tailles $1 \times Q \times Q \times 1$. Soit :

$$f : \llbracket 1, n \rrbracket^n \xrightarrow{f} \mathbb{N} \\ i \xrightarrow{\phantom{f}} (i-1)(n-i) \tag{7.4.18}$$

$f$ est une fonction dérivable et sa dérivée est : $\forall i \in \llbracket 1, n \rrbracket f'(i) = n - 2i$.

Alors,

$$f'(i) = 0 \Leftrightarrow i = \frac{n}{2}$$

Le maximum de $f$ est atteint en $\frac{n}{2}$ donc la TT-matrice $\mathbf{A}_{\frac{n}{2}}[z_{\frac{n}{2}}]$ est celle qui contient le plus de de cores de tailles $1 \times Q \times Q \times 1$ □

Si nous nous basons sur le graphique 7.12 et la propriété 12, nous pouvons dire qu'une meilleure fusion des cores se ferait en fusionnant les cores de sorte à avoir des cores de mêmes volumes.

Contrairement à la fusion des cores avec un nombre uniforme de cores, la fusion des cores avec des volumes de cores uniformes nécessite un pré-traitement qui permettra de définir le nombre de cores ($P_k$) dans chaque paquet et donc nous n'avons pas nécessairement l'égalité des nombres cores dans chaque paquet qui vérifient $P_k \in [1, m]$ tel que $\sum_k P_k = m$. De plus, nous aurons une contrainte qui consiste à commencer par trouver le meilleur regroupement des cores selon $\mathbf{A}_{\frac{n}{2}}[z_{\frac{n}{2}}]$ puis le propager dans les autres cores.

Dans cette méthode de fusion, chaque nouveau core est composé d'un nombre de cores $P_k$. Dans la figure 7.13, nous présentons la distribution du nombre de cores $P_k$ que contiendra chaque nouveau core $\mathbf{C}_i^k[z_i]$ dans le cadre ou $d = 5$.

<figure><img src="image_76.png" /><figcaption>FIGURE 7.13 – Répartition des cores dans les paquets de volumes uniformes de cores</figcaption></figure>

**Procédure de fusion.** Considérons $\mathbf{A}_i[z_i]$ une TT-matrice d'ordre $m$ et de TT-rang égal à 1. Cette procédure se découpe en quatre phases : premièrement, nous faisons une analyse symbolique sur les facteurs afin de trouver la meilleure manière d'équilibrer les cores pour l'ensemble des $\mathbf{A}_i[z_i]$. Deuxièmement, nous faisons une phase de rééquilibrage, ensuite une phase de propagation du rééquilibrage et enfin une phase de fusion des cores par produits de kronecker.

Nous détaillons les différentes étapes de cette procédure :

1. **Phase d'initialisation.** Nous faisons la procédure de Novikov de manière symbolique afin de déterminer les volumes des cores et des TT-matrices $\mathbf{A}_i[z_i]$
(a) **Traitement symbolique des facteurs.** Construction symbolique des cores des TT-matrices $\mathbf{A}_i[z_i]$.

Nous déterminons dans un premier temps les emplacements (c'est-à-dire les indices $\ell$ et $i$) des cores de taille $1 \times r \times 1 \times 1$ et ceux de taille $1 \times 1 \times r \times 1$ obtenus après la phase de décomposition des facteurs par une SVD.

Nous déterminons ensuite les cores de taille $1 \times r \times r \times 1$ et $1 \times 1 \times 1 \times 1$ obtenus lors de la phase d'ajout des dimensions non essentielles.

(b) **Calcul des volumes des différents cores.** Nous calculons le volume $V_i^\ell$ de chaque core $\mathbf{G}_i^\ell[z_i]$ qui compose les différentes TT-matrices $\mathbf{A}_i[z_i]$.

(c) **Calcul des volumes des $\mathbf{A}_i[z_i]$**. Nous calculons la somme $V_i^\ell$ pour obtenir les volumes des TT-matrices $\mathbf{A}_i[z_i]$.

$$V^{\mathbf{A}_i} = \sum_{\ell=1}^{m} V_i^{\ell}$$

2. **Phase de rééquilibrage.** Nous avons remarqué que les matrices identités se concentrent dans les TT-matrices du milieu et principalement la TT-matrice $\mathbf{A}_{\frac{n}{2}}$. Nous allons utiliser $V^{\mathbf{A}_{\frac{n}{2}}}$ pour faire un premier équilibrage des cores. Cet équilibrage se fait de la manière suivante :

(a) Choisir une contrainte $V_{\max}$, le volume maximal des cores.

(b) Choisir $d \in [1, m]$ tel que $d+1$ est notre nombre de cores.

(c) Calcul de la taille moyenne optimale que devraient avoir les nouveaux corps :

$$size = \frac{V^{\mathbf{A}_{\frac{n}{2}}}}{d}$$

(d) Regrouper les cores en $d$ groupes de sorte à ne pas dépasser la taille moyenne optimale.

(e) Mettre les cores restant dans un seul core puis calculer $V^{d+1}$ leur volume total. Si $V^{d+1} > V_{\max}$, alors morceler ce nouveau core en autant de paquets de cores qu'il faut afin de respecter cette contrainte. Nous notons ce nombre de paquets $b$.

À la fin de cette phase, nous obtenons une liste avec les indices $\ell$ des cores à fusionner afin d'avoir un bon équilibrage pour la TT-matrice $\mathbf{A}_{\frac{n}{2}}[z_{\frac{n}{2}}]$.

3. **Phase de propagation du rééquilibrage sur les autres TT-matrices :** Cette phase consiste à faire une propagation du rééquilibrage présenté auparavant en respectant la contrainte des cores à équilibrer.

Pour faire cette propagation, nous parcourons les indices de 1 à $n$ en vérifiant que si nous fusionnons les cores des autres $\mathbf{A}_i[z_i]$ nous respectons la contrainte.

4. **Construction des TT-matrices $\mathbf{A}_i[z_i]$ :** Calculer les $p+b$ nouveaux cores par produit de Kronecker.

$$\forall 1 \le k \le d+b-1, \quad \mathbf{C}_i^k[z_i] = \bigotimes_{\ell=(k-1)\times P_k+1}^{k \times P_{k+1}} \mathbf{G}_i^\ell[z_i] \quad (7.4.19)$$

À la fin de cette étape, nous obtenons les TT-matrices $\mathbf{A}_i[z_i]$ :

$$\mathbf{A}_i[z_i][i_1, \dots, i_{d+1}; j_1, \dots, j_{d+1}] = \mathbf{C}_i^1[z_i](i_1, j_1) \mathbf{C}_i^2[z_i](i_2, j_2) \dots \mathbf{G}_i^{d+b}[z_i](i_{d+b}, j_{d+b}) \quad (7.4.20)$$

### 7.4.4 Résultats

#### 7.4.4.1 Influence de la méthode de la fusion sur les volumes des cores

Nous nous intéresserons aux logarithmes des volumes des cores des TT-matrices. Pour cela, nous considérons un modèle SBM avec $n = 15$ et une structure de distance associative et nous traçons la heatmap des volumes des cores pour différentes fusions.
<figure><img src="image_77.png" /><figcaption>FIGURE 7.14 – Volumes des cores des $A_i[z_i]$ sans fusion ($n = 15$ et $m = 105$)</figcaption></figure>

<figure><img src="image_78.png" /><figcaption>FIGURE 7.15 – Volumes des cores des $B_i$ sans fusion ($n = 15$ et $m = 105$)</figcaption></figure>

<figure><img src="image_79.png" /><figcaption>FIGURE 7.16 – Volumes des cores des $A_i[z_i]$ pour la fusion par un nombre uniforme de cores ($n = 15$ et $d = 35$)</figcaption></figure>

<figure><img src="image_80.png" /><figcaption>FIGURE 7.17 – Volumes des cores des $B_i$ pour la fusion par un nombre uniforme de cores ($n = 15$ et $d = 35$)</figcaption></figure>

<figure><img src="image_81.png" /><figcaption>FIGURE 7.18 – Volumes des cores des $A_i[z_i]$ avec fusion par un volume uniforme de cores ($n = 15$ et $d = 58$)</figcaption></figure>

<figure><img src="image_82.png" /><figcaption>FIGURE 7.19 – Volumes des cores des $B_i$ avec fusion par un volume uniforme de cores ($n = 15$ et $d = 58$)</figcaption></figure>
Les trois figures 7.14, 7.16 et 7.18 représentent les heatmaps des volumes des cores des TT-matrices $\mathbf{A}_i[z_i]$. Les trois figures 7.15, 7.17 et 7.19 représentent les heatmaps des volumes des cores des TT-matrices $\mathbf{B}_i$.

Pour chacune de ces figures, nous traçons en ligne les indices $i$ des TT-matrices et en colonnes les indices de leurs cores.

- Dans les graphes 7.14 et 7.15 nous traçons les volumes des TT-matrices sans fusion le nombre de cores est 105.

- Dans les graphes 7.16 et 7.17 nous traçons les volumes des TT-matrices avec fusion avec un nombre uniforme de cores par paquet de trois, le nombre de cores est 35.

- Dans les graphes 7.18 et 7.19 nous traçons les volumes des TT-matrices avec fusion avec un volume uniforme de cores par paquet de trois, le nombre de cores est 58.

Nous constatons que les mêmes motifs qui apparaissent dans heatmaps des cores $\mathbf{A}_i[z_i]$ se retrouvent dans les heatmaps des cores $\mathbf{B}_i$. Nous allons juste faire une comparaison des figures 7.14, 7.16 et 7.18. Tout d'abord, nous nous intéressons aux similarités entre les trois heatmaps. Nous remarquons qu'une grande partie des motifs sont partagés par les trois figures (les pixels noirs). Cela vient du fait que les TT-matrices $\mathbf{A}_i[z_i]$ pour $i = 1, \dots, 5$ et $i = 12, \dots, 15$ contiennent beaucoup de cores de volumes nuls. Cela vient du fait que seuls $(i-1)(n-i)$ cores sont de tailles $1 \times Q \times Q \times 1$. Dans le cas où $i = 1$ ou $15$, nous n'aurons aucun core de taille $1 \times Q \times Q \times 1$. Dans le cas où $i = 2$ ou $14$, nous aurons 14 cores de cette taille.

Nous observons aussi que les heatmaps 7.14 et 7.16 sont très similaires. Cela est dû au fait que comme nous fusionnons les cores de manière consécutive par paquet de 3, les mêmes patterns de volumes des cores se retrouvent sans fusion et avec fusion avec nombre uniforme de cores.

Nous observons aussi que pour la fusion de cores en paquets uniformes, il y a beaucoup de cores avec des volumes proches de 6.5, tandis que pour la fusion par volumes uniformes des cores, très peu de cores ont un grand volume : seuls les six derniers cores des TT-matrices $\mathbf{A}_i[z_i]$ ont des cores avec de très grandes tailles (le volume maximal est de 10.4).

À partir de ces figures, nous pouvons dire que la fusion avec volumes uniformes des cores nous donne une meilleure distribution des volumes des cores.

### 7.4.4.2 Diminution des rangs des cores lors du calcul de W

Nous avons remarqué que lors du calcul de la constante de normalisation sans fusion et avec fusion, il y a une différence pour les TT-rangs observés lors des calculs intermédiaires, et cela, quelque soit la structure de distances choisie. Afin d'illustrer ce phénomène, nous avons fait trois expériences.

Nous définissons un nombre de classes $Q = 3$ et les proportions d'individus dans les classes $\alpha = (\frac{1}{3}, \frac{1}{3}, \frac{1}{3})$ et une structure de distances $\Lambda$ associative (pour 7.20) et cores periphery (pour les figures 7.21 et 7.22).

Nous choisissons trois nombres d'individus : $n = 15$, $n = 33$ et $n = 50$. Puis, pour les deux graphes 7.20 et 7.21, nous appliquons la procédure de Novikov à la fois sans fusion des cores (courbe orange), avec une fusion avec nombres un uniforme de cores par paquets de trois (courbe bleue) puis, avec un volume uniforme des cores (courbe verte). Ensuite, nous calculons la constante de normalisation pour chaque découpage avec une compression via l'algorithme du rounding avec comme paramètres $\epsilon = 5 \times 10^{-2}$.

Enfin, dans le graphique 7.22, nous appliquons la procédure de Novikov avec fusion des cores d'abord, avec nombre uniforme de cores (par paquet de trois) puis avec un volume uniforme de cores étant donné que la procédure sans fusion ne peut être appliquée. De plus, nous compressons les résultats intermédiaires via un rounding ayant pour paramètres $\epsilon = 5 \times 10^{-2}$ et $r_{max} = 27$.

Dans les graphiques 7.20, 7.21 et 7.22, nous traçons l'évolution des TT-rangs maximaux des cores obtenus après chaque itération du calcul de la constante de normalisation.
<figure><img src="image_83.png" /></figure>

FIGURE 7.20 – Évolution des rangs des cores en fonction des itérations de l’algorithme 1 pour les 3 procédures ($n = 15, \epsilon = 5 \times 10^{-2}$ et $r_{max} = 27$)

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

FIGURE 7.21 – Évolution des rangs des cores en fonction des itérations de l’algorithme 1 pour les 3 procédures ($n = 33, \epsilon = 5 \times 10^{-2}$ et $r_{max} = 27$)

Dans les deux graphiques 7.20 et 7.21, nous constatons que la fusion des cores permet de diminuer le rang lors du calcul de la constante de normalisation. En effet, dans la figure 7.20 le rang maximal des cores observé atteint 48 lors de l'itération 4 de l'algorithme 1 pour le calcul sans fusion, tandis qu'il n'est que de 40 (resp. 39) pour le calcul avec fusion avec nombre de cores uniforme (resp. pour le calcul avec fusion avec volumes uniformes).

Concernant le temps de calcul pour la constante de normalisation, il est relativement le même pour les trois méthodes. Nous exprimons les temps de calculs dans le tableau 7.4.4.2.

Dans la figure 7.22, les rangs maximaux des cores des TT-matrices par itération de calcul par les deux méthodes sont les mêmes quelle que soit la procédure de fusion adoptée.

### 7.4.4.3 Comparaison des constantes de normalisations obtenues par les trois méthodes

**Protocole.** Nous définissons un nombre de classes $Q = 3$ et les proportions d'individus dans les classes $\alpha = (\frac{1}{3}, \frac{1}{3}, \frac{1}{3})$ et une structure de distances $\Lambda$ associative.

Puis, nous faisons varier le nombre d'individus $n$ entre 12 et 42 par pas de 3 et nous appliquons la
<table><thead><tr><th>Méthode</th><th>Sans fusion</th><th>Fusion avec nombres uniformes</th><th>Fusion avec volumes uniformes</th></tr></thead><tbody><tr><td>n = 15</td><td>0.83 s</td><td>0.79 s</td><td>0.81 s</td></tr><tr><td>n = 33</td><td>3.35 s</td><td>3.42 s</td><td>3.40 s</td></tr></tbody></table>

<figure><img src="image_85.png" /><figcaption>FIGURE 7.22 – Évolution des rangs des cores en fonction des itérations de l’algorithme 1 en changeant les paramètres du rounding ($n = 50, \epsilon = 5 \times 10^{-2}$ et $r_{max} = 27$)</figcaption></figure>

procédure de Novikov sur le modèle à la fois sans fusion des cores, avec une fusion avec un nombre uniforme de cores (par paquets de trois) et avec un volume uniforme de cores.

Enfin, nous calculons de la constante de normalisation avec une compression via un rounding avec comme paramètres $\epsilon = 5 \times 10^{-2}$ et $r_{max} = 27$.

**Résultats** : Nous présentons les résultats dans le tableau 7.13 avec en ligne le nombre d’individus et en colonnes la constante de normalisation. Nous observons que les erreurs relatives entre les différentes procédures de fusion sont très faibles voir nulles, le maximum de ces erreurs est atteint en $n = 27$ où il vaut 0.008.

Les erreurs relatives entre les calculs sans fusion et les deux méthodes de fusion et le calcul sans fusion est très faible lui aussi. Le maximum est de 0.08 et est atteint pour $n = 27$.

<table><caption>TABLE 7.13 – Tableau présentant les résultats de la constante de normalisation avec et sans fusion</caption><thead><tr><th>n</th><th>Sans fusion</th><th>Fusion avec nombres uniformes (paquet de 5 cores)</th><th>Fusion avec volumes uniformes</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.4.5 Discussion

Durant cette section, nous nous sommes intéressés à la fusion des cores, nous avons expliqué comment réduire les dimensions des TT-matrices grâce à une procédure de fusion des cores par produits de Kronecker. Nous avons expliqué deux procédures de fusion, la première avec un nombre uniforme de cores et une deuxième avec un volume uniforme de cores.

**Résumé des résultats** Durant la sous-section 7.4.4.1, nous avons présenté l'influence de la méthode de fusion des cores adoptées sur les nouveaux cores résultants. Nous avons montré que la fusion avec des volumes uniformes des cores conduits à une meilleure distribution des nombres d'éléments dans les cores dans les TT-matrices $A_i[z_i]$ et $B_i$.

Durant la sous-section 7.4.4.2, nous nous sommes intéressés au calcul de la constante de normalisation $W$ par l'intermédiaire de la fusion des cores. Dans les résultats, nous avons remarqué que ce calcul avec fusion permet de diminuer les rangs des TT-matrices à chaque itération de l'algorithme 1.

Enfin, durant la section 7.4.4.3, nous nous sommes intéressé aux résultats de la constante de normalisation obtenus avec fusion et sans fusions des cores. Nous avons expliqué que les différences entre les différents résultats obtenus sont très faibles.

**Temps de calcul.** Nous avons remarqué que le temps de calcul est plus faible lorsque nous faisons la procédure sans fusion avec rounding à précision prescrite. Ce temps de calcul devient quasiment identique pour les trois méthodes lorsque le calcul se fait avec une précision et un rang max prescrits. Cependant, comme précisé lors de la sous-section 7.3.6 nous ne pouvons pas faire les calculs, que ce soit des marginales ou de la constante de normalisation pour un nombre d'individus supérieur à 44, grâce à cette procédure, nous pouvons maintenant le faire.

**Passage à l'échelle** Avec la procédure de fusion des cores, nous arrivons à faire les calculs des marginales et de la constante de normalisation pour des nombres d'individus relativement plus élevés. En effet, nous avons étendu la procédure de Novikov pour un nombre d'individus jusqu'à $n = 66$. Soit un modèle graphique complet avec un nombre de liens $m = 2145$.

### 7.4.6 Conclusion

Dans cette section, nous avons proposé deux méthodes de réduction des dimensions des TT-matrices $A_i[z_i]$ utilisées dans l'approche de Novikov [NROV14]. La première consiste à faire une fusion avec un nombre uniforme de cores et la seconde avec un volume uniforme de cores.

Nous avons ensuite présenté les résultats pour chacune de ces méthodes de fusion, à partir de ces résultats, nous pouvons dire que ses deux méthodes de fusion sont très intéressantes étant donné qu'elles fournissent des résultats très précis, un rang plus faible et un temps de calcul quasiment identique que la procédure sans fusion.

## 7.5 Conclusion du chapitre

Durant ce chapitre, nous nous sommes interressés à trois limitations que nous avons rencontrées lors de la déclinaison de la procédure de Novikov sur des SBM que nous avons détaillé durant l'introduction. Nous résumons les procédures adoptées pour résoudre ces limitations :

**Éviter les problèmes des petites valeurs.** Nous avons présenté durant la section 7.2 comment faire une homothétie pour agrandir les facteurs binaires $\psi'_{i,j}$ à l'entrée de la procédure de Novikov pour éviter des valeurs nulles pour la constante de normalisation et des marginales. Nous avons ensuite expliqué quelle constante choisir pour faire ces calculs.

Grâce à cette procédure, nous avons réussi à calculer la constante de normalisation pour des SBM avec
un nombre d'individus $n = 33$ et de calculer des marginales pour des SBM avec un nombre d'individus $n = 26$.

**Explosion des rangs des cores.** Nous avons expliqué durant l'introduction que chaque produit de TT-matrices augmente le rang de la TT-matrice résultante et que sans rounding sur la constante de normalisation $W$ lors de chaque itération augmentera jusqu'à atteindre $Q^n$ alors que $W$ est un scalaire. Nous avons expliqué durant la section 7.3, que le rounding avec une précision $\epsilon$ ne stabilise pas ce rang et qu'il faudrait modifier les paramètres de cet algorithme en le faisant avec $r_{max}$ et $\epsilon$. Grâce à cette procédure, nous avons réussi à accélérer les calculs.

En associant cette procédure à l'homothétie, nous avons réussi à faire les calculs de la constante de normalisation et des marginales pour un modèle avec nombre d'individus égale à 44.

**Réduction du nombre de cores.** Nous avons présenté dans la section 7.4 deux procédures de réduction de l'ordre des TT-matrices $A_i[z_i]$ et $B_i$. Une première est basée sur une fusion avec un nombre uniforme de cores et une deuxième, consiste à faire une fusion avec un volume uniforme des cores. Nous avons ensuite expliqué les avantages et les inconvénients de chaque approche.

Nous avons conclu cette section par une présentation des résultats que nous avons décomposé en trois parties, influence de la méthode sur la distribution des volumes des cores, la diminution des rangs des produits de TT-matrices TT-vecteurs lors du calcul de $W$ et enfin, nous avons comparé les résultats de la constante de normalisation obtenus avec et sans fusion.

À partir des résultats présentés, cette procédure, combinée aux deux autres, nous permettent d'étendre la procédure de Novikov pour des SBM avec un nombre d'individus égale à $n = 66$ soit un modèle graphique avec un nombre de liens égal à $m = 2145$.
# Chapitre 8

# Conclusion et perspectives

## Table des matières

<table><tr><td>8.1</td><td>Rappel des contributions</td><td>141</td></tr><tr><td>8.2</td><td>Perspectives</td><td>141</td></tr><tr><td>8.2.1</td><td>Perspectives numériques</td><td>141</td></tr><tr><td>8.2.2</td><td>Réduction du temps de calcul</td><td>142</td></tr><tr><td>8.2.3</td><td>Proposition d'un algorithme TT-EM pour l'estimation d'un SBM</td><td>143</td></tr><tr><td>8.2.4</td><td>Extension de la procédure de Novikov à d'autres modèles graphiques</td><td>143</td></tr></table>
## 8.1 Rappel des contributions

Nous nous sommes intéressés au calcul des marginales binaires des SBM pondérés par un tableau de distances par une approche de type Tensor Train et nous avons présenté les contributions suivantes.

Dans le chapitre 3 et en annexe B.2, nous avons réalisé un clustering de matrices de distances de séquences de marqueurs d'intérêt taxonomique (barcodes) d'arbres guyanais. Nous avons montré les similitudes et complémentarités entre les classifications obtenues par SBM et par les méthodes plus classiques de CAH. Nous avons aussi montré que la classification SBM permet de détecter, par le biais de la distribution estimée des distances intra-classe, des situations où les classifications moléculaires peuvent être en désaccord avec la classification botanique.

Nous avons présenté dans le chapitre 5 une écriture exacte de la loi jointe d'un SBM conditionnellement à des distances observées comme un tenseur au format TT, sans approximation. Cela mène naturellement à la séparation des variables dans la loi jointe, qui permet le calcul des marginales. Nous avons expliqué les matrices intervenant dans le calcul des marginales puis décrit le processus de tensorisation (quantisation) de ces matrices qui permet de les écrire sous la forme de TT-matrices ainsi que la méthode d'implémentation des calculs des marginales.
Nous avons ensuite présenté un algorithme de type programmation dynamique pour le calcul des marginales binaires par une mutualisation des calculs ainsi que sa complexité.

Dans le chapitre 6, nous avons comparé l'approche TT pour le calcul des marginales unaires et binaires d'un SBM avec trois autres méthodes présentes dans la littérature (énumération, approximation champ moyen, échantillonneur de Gibbs). Nous avons montré dans nos résultats que l'approche TT est stable vis-à-vis de la variation du nombre d'individus et de la structure de distance. De plus, elle permet de calculer les marginales dans un temps relativement raisonnable et avec une bonne précision comparée à l'approche de Gibbs.
Cette approche soulève des difficultés numériques redoutables dès que le nombre d'individus est significatif (contrairement à l'approche par le champ moyen par exemple).

Dans le chapitre 7, nous avons étudié les verrous numériques qui surviennent lors du passage à l'échelle de l'approche TT puis nous avons présenté une procédure opérationnelle qui intègre quelques solutions à ces verrous. Premièrement, les valeurs des marginales non normalisées et de la constante de normalisation sont très petites, nous avons proposé une procédure d'homothétie pour lever ce verrou. Un deuxième problème vient de l'explosion des rangs des cores, nous avons proposé un changement de paramètres dans la compression par l'algorithme du rounding.
Un troisième verrou que nous avons rencontré est que la dimension des TT-matrices intermédiaires $A_i[z_i]$ et $B_i$ est quadratique par rapport au nombre d'individus. Nous avons proposé une procédure de réduction des dimensions de TT-matrices $A_i[z_i]$ par fusion de ses cores.
Grâce à ces trois procédures, nous avons étendu la procédure de Novikov pour des SBM avec des nombres d'individus pouvant atteindre à $n = 66$.

Ce domaine numérique reste cependant en grande partie inexploré et mal compris.

## 8.2 Perspectives

### 8.2.1 Perspectives numériques

Lors de l'application de la procédure de Novikov [NROV14] pour le calcul des marginales des SBM, nous nous sommes confrontés à des limitations numériques quand le nombre d'individus augmente. Afin d'y remédier, nous avons proposé trois procédures qui nous ont permis de traiter des problèmes de plus grande taille. Nous avons ainsi résolu le problème des petites valeurs. Cependant, en dépit
de ces procédures, nous sommes toujours confrontés à des limitations liées à l'évolution des rangs des cores et leur nombre durant les calculs des marginales ou de la constante de normalisation des modèles.

Nous avons essayé de lever le premier verrou par le biais d'un changement de paramètres du rounding. Nous avons réussi à faire des calculs plus rapidement, mais nous ne pouvons pas contrôler la précision (les calculs du rounding à rang prescrit sont plus rapides qu'avec une précision prescrite, mais au détriment justement de la précision). Une des perspectives est de repenser un autre algorithme de compression de TT-matrices ou de TT-vecteurs. Al Daas & al. ont proposé une librairie [DBB20] pour faires des opérations parallélisées au format TT. Dans un autre article [DBC+21], Al Daas & al. ont proposé un nouvel algorithme pour faire le rounding. Une deuxième possibilité pour éviter l'explosion du rang des cores consiste à remplacer la SVD exacte pour la séparation des variables par une approximation de rang $r < Q$ (où $Q$ est le nombre de classe). Si $r = 1$ on se retrouve avec une approximation de type champ moyen.

Le deuxième verrou auquel nous avons été confrontés est l'augmentation du nombre de cores qui composent les TT-matrices $A_i[z_i]$ et $B_i$ de manière quadratique par rapport aux nombres d'individus. Nous avons essayé de résoudre ce problème par une procédure de fusion des cores des TT-matrices $A_i[z_i]$ par produit de kronecker. Grâce à cette procédure, nous avons réussi à traiter des problèmes de plus grande taille, cependant nous n'avons pas réussi à dépasser le cap des SBM avec un nombre d'individus $n = 66$.

Une perspective possible pour dépasser cette limite consiste à trouver une méthode de fusion des cores plus optimale que celles abordées dans le chapitre 7.

### 8.2.2 Réduction du temps de calcul

Une première perspective de ce travail serait de réduire les temps de calcul nécessaires pour le calcul des marginales grâce à la parallélisation, avec deux possibilités. En effet, la plupart des calculs sont indépendants.

Premièrement, il est possible de distribuer les calculs des feuilles de l'arbre présenté dans la section 5.3 sur plusieurs nœuds d'un cluster. Pour expliquer une méthode possible pour faire cette parallélisation, nous reprenons la formule de récurrence selon le parcours de l'arbre des feuilles à la racine :

$$ \left\{ \begin{array}{ll} \mathbf{G}_{i,i} = \mathbf{B}_i & \forall i \in [1, n-1] \\ \mathbf{G}_{i,j+1} = \mathbf{G}_{i,j} \times \mathbf{B}_{j+1} & \forall j \in [i+1, n] \end{array} \right. \qquad (8.2.1) $$

Nous reprenons la figure 5.2 présentée dans la section 5.3 (où le nombre d'individus est $n = 5$) en coloriant de la même couleur les ensembles de calculs indépendants.

Dans la figure 8.1, la construction de cet arbre nécessite 10 calculs. Ces derniers peuvent être rassemblés en quatre ensembles indépendants : le calcul $G_{12}$, $G_{13}$, $G_{14}$ et $G_{15}$ (noir), celui de $G_{23}$, $G_{24}$, $G_{2F5}$ (orange), celui de $G_{34}$, $G_{35}$ (bleu) et celui de $G_{45}$ (vert). Nous remarquons que chaque ensemble ne dépend que de $B_1, \dots, B_5$. Nous pouvons donc distribuer la totalité des calculs de cet arbre sur 4 nœuds d'un cluster, chaque nœud correspond à une couleur : 4 calculs pour le premier, 3 pour le second, 2 pour le troisième et un pour le dernier. Cette méthode peut se généraliser pour la distribution des calculs des feuilles d'un SBM avec un nombre d'individus $n \in \mathbb{N}$ quelconque sur $n-2$ nœuds d'un cluster de calcul. Pour le faire nous devons procéder de cette manière : pour chaque $i \in [1, n-1]$, faire la boucle pour le calcul de $\forall j \in [i+1, n]$, $G_{i,j}$ sur un nouveau nœud du cluster.

Une deuxième phase de parallélisation peut être faite au niveau du calcul des marginales. En effet, pour un modèle SBM nous devenons calculer $\frac{n(n-1)Q^2}{2}$. Chaque calcul d'une marginale est indépendant des autres. Il est donc facile de distribuer les calculs des marginales sur plusieurs nœuds d'un cluster de calcul.
<figure><img src="image_86.png" /><figcaption>FIGURE 8.1 – Décomposition des calculs des parties centrales de l’expression des marginales binaires</figcaption></figure>

### 8.2.3 Proposition d'un algorithme TT-EM pour l'estimation d'un SBM

Nous avons présenté durant la section 2.3 l’algorithme EM [DLR77] où chaque itération se décompose en deux étapes : une étape de calcul d’espérance (E-step) qui nécessite le calcul de toutes les marginales binaires conditionnelles et une étape de maximisation (M-step).

Nous pouvons proposer une alternative à cet algorithme via une approche TT : l’algorithme TT-EM prendra comme entrées : *D* matrice des distances, *Q* nombre de classes, *n* nombre d’individus, *ϵ* précision du rounding, *τ* critère d’arrêt de l’algorithme. Il retournera les classes des individus, $\alpha^{opt}$ et $\Lambda^{opt}$, les paramètres du modèle SBM.

L’avantage de cet algorithme est que l’approche TT permet de calculer les marginales avec une bonne précision comme nous avons expliqué dans le chapitre 6. Cependant, il reste limité comparé à des algorithmes concurrents comme le VEM [DPR08] utilisés dans diverses applications comme [TBC17, APF22]. En effet, même s’il repose sur une moins bonne approximation des marginales, le VEM donne des résultats de très bonnes qualités en terme d’estimation des paramètres d’un SBM, avec un temps de calcul raisonnable.

### 8.2.4 Extension de la procédure de Novikov à d'autres modèles graphiques

Nous avons travaillé sur un modèle SBM, dont tous les facteurs sont binaires et pour le calcul des marginales binaires. Chacun de ces trois choix peut être étendu comme perspective ouverte par ce travail. Nous les détaillons ci-dessous.

Appliquer l’approche pour le calcul des marginales binaires d’un modèle graphique dont les facteurs sont binaires au plus, comme le modèle d’Ising par exemple [BRU67, Nis09]. Notre approche peut se transposer "nec varietur" à de tels modèles. L’approche sera plus simple, car le modèle SBM est le cas le pire quant au nombre de facteurs (en $n(n-1)/2$ alors que pour le modèle d’Ising il est en $2n$).

L’approche de Novikov pour le calcul de la constante de normalisation ou des marginales unaires s’applique à tout modèle graphique. Si les facteurs sont d-naires (avec $d > 2$), on perd l’exactitude de
la décomposition TT : on remplace l'exactitude de la SVD pour la séparation des variables des facteurs binaires présentée dans la section 5.2.3 par une approximation par TT-SVD sur les facteurs d-naires si $d > 2$. On obtient alors une meilleure approximation par un tenseur en format TT de chaque facteur. Au-delà de ce changement, la suite des calculs se décline sans changement, et notre approche pour le calcul des marginales binaires peut se décliner également. On calcule les TT-matrices intermédiaires $\mathbf{A}_i[z_i]$ et $\mathbf{B}_i$ par l'approche de Novikov puis on utilise l'algorithme présenté dans la section 5.3 pour le calcul des marginales des binaires. Celles-ci s'expriment :

$$ \forall i,j \in [1,n]^2, p_{i,j}(z_i, z_j) = \underbrace{\mathbf{B}_1 \times \dots \times \mathbf{B}_{i-1}}_{\mathbf{L}_i} A_i[z_i] \underbrace{\mathbf{B}_{i+1} \times \dots \times \mathbf{B}_{j-1}}_{\mathbf{C}_{i,j}} A_j[z_j] \underbrace{\mathbf{B}_{j+1} \times \dots \times \mathbf{B}_n}_{\mathbf{R}_i} \quad (8.2.2) $$

Cette approche peut également se décliner pour le calcul des marginales k-naires (où $k > 2$) d'un modèle graphique avec des facteurs d-naires. Le calcul des marginales ternaires par exemple s'exprime :
$$ \forall i,j,k \in \llbracket 1,n \rrbracket^3, \quad p_{i,j,k}(z_i, z_j, z_k) $$

$$ = \underbrace{\mathbf{B}_1 \times \dots \times \mathbf{B}_{i-1}}_{\mathbf{L}_i} A_i[z_i] \underbrace{\mathbf{B}_{i+1} \times \dots \times \mathbf{B}_{j-1}}_{\mathbf{C}_{i,j}} A_j[z_j] \underbrace{\mathbf{B}_{j+1} \times \dots \times \mathbf{B}_{k-1}}_{\mathbf{C}_{j,k}} A_k[z_k] \underbrace{\mathbf{B}_{k+1} \times \dots \times \mathbf{B}_n}_{\mathbf{R}_k} \quad (8.2.3) $$

La formule (8.2.3) peut être généralisée pour des marginales k-naires en rajoutant plusieurs parties centrales de cette manière : Pour $\forall i_1, i_2, \dots, i_k \in \llbracket 1, n \rrbracket^k$

$$ p_{i_1, i_2, \dots, i_k}(z_{i_1}, z_{i_2}, \dots, z_{i_k}) = \mathbf{L}_{i_1} \mathbf{A}_{i_1}[z_{i_1}] \mathbf{C}_{i_1, i_2} \mathbf{A}_{i_2}[z_{i_2}] \dots \mathbf{A}_{i_{k-1}}[z_{i_{k-1}}] \mathbf{C}_{i_{k-1}, i_k} \mathbf{A}_{i_k}[z_{i_k}] \mathbf{R}_{i_k} \quad (8.2.4) $$

Les matrices $\mathbf{C}_{i_{k-1}, i_k}$ peuvent être calculées via la mutualisation des calculs de la même manière que pour des modèles graphiques avec des facteurs binaires.

À partir de cela, la procédure présentée dans le chapitre 5 peut être transposée pour le calcul des marginales de n'importe quel modèle graphique.
# Appendices
# Annexe A

## Complément des résultats sur les méthodes de passage à l'échelle de l'approche TT

Dans cette annexe, nous complétons les résultats présentés dans la sous-section résultats de la section 3 du chapitre. Nous présentons les résultats du calcul de la constante de normalisation pour les quatre structures Hiérarchique, ordonnée, dissociative et hiérarchique 2 que nous avons omis. Nous rappelons le tableau des paramètres du rounding utilisé pour les discussions et les commentaires des résultats.

<table><caption>TABLE A.1 – Tableau des abréviations pour différentes méthodes de rounding.</caption><thead><tr><th>ε</th><th>r<sub>max</sub> = 3 et ε</th><th>r<sub>max</sub> = 9 et ε</th><th>r<sub>max</sub> = 9</th><th>r<sub>max</sub> = 27 et ε</th><th>r<sub>max</sub> = 81 et ε</th></tr></thead><tbody><tr><td>R<sub>1</sub></td><td>R<sub>2</sub></td><td>R<sub>3</sub></td><td>R<sub>4</sub></td><td>R<sub>5</sub></td><td>R<sub>6</sub></td></tr></tbody></table>

### Hiérarchique

<table><caption>TABLE A.2 – Temps de calcul en secondes pour la constante de normalisation dans le cas Hiérarchique</caption><thead><tr><th></th><th>R<sub>1</sub> (r<sub>max</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><th>R<sub>6</sub></th></tr></thead><tbody><tr><td>n=8</td><td>1.35 × 10<sup>-25</sup>(8)</td><td>1.35 × 10<sup>-25</sup></td><td>1.33 × 10<sup>-25</sup></td><td>1.35 × 10<sup>-25</sup></td><td>1.35 × 10<sup>-25</sup></td><td>1.35 × 10<sup>-25</sup></td></tr><tr><td>n=10</td><td>3.97 × 10<sup>-43</sup>(21)</td><td><b>4.43 × 10<sup>-85</sup></b></td><td>3.96 × 10<sup>-43</sup></td><td>3.98 × 10<sup>-43</sup></td><td>3.97 × 10<sup>-43</sup></td><td>3.97 × 10<sup>-43</sup></td></tr><tr><td>n=12</td><td>4.39 × 10<sup>-73</sup>(13)</td><td>4.39 × 10<sup>-73</sup></td><td>4.39 × 10<sup>-73</sup></td><td>4.41 × 10<sup>-73</sup></td><td>4.39 × 10<sup>-73</sup></td><td>4.39 × 10<sup>-73</sup></td></tr><tr><td>n=14</td><td>1.63 × 10<sup>-81</sup>(14)</td><td><b>2.20 × 10<sup>-84</sup></b></td><td>1.63 × 10<sup>-81</sup></td><td>1.63 × 10<sup>-81</sup></td><td>1.63 × 10<sup>-81</sup></td><td>1.63 × 10<sup>-81</sup></td></tr><tr><td>n=16</td><td>6.38 × 10<sup>-111</sup>(37)</td><td><b>2.44 × 10<sup>-140</sup></b></td><td>6.38 × 10<sup>-111</sup></td><td>6.43 × 10<sup>-111</sup></td><td>6.38 × 10<sup>-111</sup></td><td>6.38 × 10<sup>-111</sup></td></tr><tr><td>n=18</td><td>4.18 × 10<sup>-137</sup>(57)</td><td><b>6.78 × 10<sup>-188</sup></b></td><td>4.37 × 10<sup>-137</sup></td><td>4.49 × 10<sup>-137</sup></td><td>4.18 × 10<sup>-137</sup></td><td>4.18 × 10<sup>-137</sup></td></tr><tr><td>n=20</td><td>2.61 × 10<sup>-164</sup>(59)</td><td><b>7.13 × 10<sup>-168</sup></b></td><td>2.25 × 10<sup>-164</sup></td><td><b>9.69 × 10<sup>-168</sup></b></td><td>2.61 × 10<sup>-164</sup></td><td>2.61 × 10<sup>-164</sup></td></tr><tr><td>n=22</td><td>4.34 × 10<sup>-190</sup>(5)</td><td>4.35 × 10<sup>-190</sup></td><td>4.35 × 10<sup>-190</sup></td><td>4.34 × 10<sup>-190</sup></td><td>4.34 × 10<sup>-190</sup></td><td>4.34 × 10<sup>-190</sup></td></tr><tr><td>n=24</td><td>5.01 × 10<sup>-229</sup>(16)</td><td><b>1.63 × 10<sup>-263</sup></b></td><td>5.04 × 10<sup>-229</sup></td><td>5.02 × 10<sup>-229</sup></td><td>5.01 × 10<sup>-229</sup></td><td>5.01 × 10<sup>-229</sup></td></tr><tr><td>n=26</td><td>3.15 × 10<sup>-285</sup>(9)</td><td><b>1.55 × 10<sup>-408</sup></b></td><td>3.15 × 10<sup>-285</sup></td><td>3.15 × 10<sup>-285</sup></td><td>3.15 × 10<sup>-285</sup></td><td>3.15 × 10<sup>-285</sup></td></tr><tr><td>n=28</td><td>1.09 × 10<sup>-333</sup>(100)</td><td>1.09 × 10<sup>-333</sup></td><td>1.09 × 10<sup>-333</sup></td><td>1.09 × 10<sup>-333</sup></td><td>1.09 × 10<sup>-333</sup></td><td>1.09 × 10<sup>-333</sup></td></tr><tr><td>n=30</td><td>2.64 × 10<sup>-384</sup>(18)</td><td>2.64 × 10<sup>-384</sup></td><td>2.64 × 10<sup>-384</sup></td><td>2.64 × 10<sup>-384</sup></td><td>2.64 × 10<sup>-384</sup></td><td>2.64 × 10<sup>-384</sup></td></tr><tr><td>n=32</td><td>2.64 × 10<sup>-439</sup>(11)</td><td>2.64 × 10<sup>-439</sup></td><td>2.64 × 10<sup>-439</sup></td><td>2.64 × 10<sup>-439</sup></td><td>2.64 × 10<sup>-439</sup></td><td>2.64 × 10<sup>-439</sup></td></tr></tbody></table>
<table><caption>TABLE A.3 – Temps de calcul en secondes pour la constante de normalisation dans le cas Hiérarchique</caption><thead><tr><th></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><th>R<sub>6</sub></th></tr></thead><tbody><tr><td>n=8</td><td>0.22</td><td>0.01</td><td>0.01</td><td>0.02</td><td>0.11</td><td>0.10</td></tr><tr><td>n=10</td><td>0.15</td><td>0.01</td><td>0.02</td><td>0.03</td><td>0.06</td><td>0.06</td></tr><tr><td>n=12</td><td>2.55</td><td>0.01</td><td>0.04</td><td>0.06</td><td>0.40</td><td>1.00</td></tr><tr><td>n=14</td><td>0.02</td><td>0.01</td><td>0.01</td><td>0.10</td><td>0.01</td><td>0.01</td></tr><tr><td>n=16</td><td>1.46</td><td>0.02</td><td>0.06</td><td>0.17</td><td>0.28</td><td>0.49</td></tr><tr><td>n=18</td><td>0.57</td><td>0.02</td><td>0.08</td><td>0.23</td><td>0.24</td><td>0.25</td></tr><tr><td>n=20</td><td>55.60</td><td>0.03</td><td>0.22</td><td>0.34</td><td>2.71</td><td>22.40</td></tr><tr><td>n=22</td><td>3.49</td><td>0.04</td><td>0.14</td><td>0.45</td><td>0.46</td><td>0.97</td></tr><tr><td>n=24</td><td>0.43</td><td>0.04</td><td>0.07</td><td>0.56</td><td>0.10</td><td>0.10</td></tr><tr><td>n=26</td><td>2.68</td><td>0.05</td><td>0.12</td><td>0.70</td><td>0.63</td><td>0.81</td></tr><tr><td>n=28</td><td>0.75</td><td>0.06</td><td>0.16</td><td>0.94</td><td>0.20</td><td>0.20</td></tr><tr><td>n=30</td><td>9.67</td><td>0.08</td><td>0.25</td><td>1.17</td><td>2.23</td><td>4.07</td></tr><tr><td>n=32</td><td>3.73</td><td>0.10</td><td>0.35</td><td>1.47</td><td>1.36</td><td>1.36</td></tr></tbody></table>

Le calcul de la constante de normalisation fait avec le rounding $R_2$ donne des valeurs très éloignées de la référence que nous avons choisi $R_1$.

Les calculs avec $R_3$, $R_5$ et $R_6$ donne des résultats très proches de $R_1$. De plus, les temps de calculs sont très bas lorsque nous utilisons la méthode $R_3$, nous pouvons donc dire que ce sont les meilleurs paramètres en termes de compromis temps de calcul précision pour cette structure. Dans ce qui va suivre, nous éliminons la méthode $R_2$ de la discussion étant donné qu'elle fournit des résultats peu précis.

## Ordonnée

<table><caption>TABLE A.4 – Constante de normalisation calculées ainsi que les $r_{max}$ trouvée dans le cas ordonnés</caption><thead><tr><th></th><th>R<sub>1</sub> (r<sub>max</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><th>R<sub>6</sub></th></tr></thead><tbody><tr><td>8</td><td>1.06 × 10<sup>-26</sup> (7)</td><td>1.04 × 10<sup>-26</sup></td><td>1.04 × 10<sup>-26</sup></td><td>1.06 × 10<sup>-26</sup></td><td>1.06 × 10<sup>-26</sup></td><td>1.06 × 10<sup>-26</sup></td></tr><tr><td>10</td><td>1.27 × 10<sup>-49</sup> (15)</td><td><b>7.78 × 10<sup>-50</sup></b></td><td>1.27 × 10<sup>-49</sup></td><td>1.26 × 10<sup>-49</sup></td><td>1.27 × 10<sup>-49</sup></td><td>1.27 × 10<sup>-49</sup></td></tr><tr><td>12</td><td>1.89 × 10<sup>-66</sup> (12)</td><td>1.89 × 10<sup>-66</sup></td><td>1.89 × 10<sup>-66</sup></td><td>1.89 × 10<sup>-66</sup></td><td>1.89 × 10<sup>-66</sup></td><td>1.89 × 10<sup>-66</sup></td></tr><tr><td>14</td><td>6.47 × 10<sup>-93</sup> (12)</td><td>6.47 × 10<sup>-93</sup></td><td>6.47 × 10<sup>-93</sup></td><td>6.47 × 10<sup>-93</sup></td><td>6.47 × 10<sup>-93</sup></td><td>6.47 × 10<sup>-93</sup></td></tr><tr><td>16</td><td>1.39 × 10<sup>-122</sup> (17)</td><td><b>5.87 × 10<sup>-123</sup></b></td><td>1.39 × 10<sup>-122</sup></td><td>1.39 × 10<sup>-122</sup></td><td>1.39 × 10<sup>-122</sup></td><td>1.39 × 10<sup>-122</sup></td></tr><tr><td>18</td><td>2.44 × 10<sup>-159</sup> (22)</td><td>2.44 × 10<sup>-159</sup></td><td>2.44 × 10<sup>-159</sup></td><td>2.45 × 10<sup>-159</sup></td><td>2.44 × 10<sup>-159</sup></td><td>2.44 × 10<sup>-159</sup></td></tr><tr><td>20</td><td>1.30 × 10<sup>-191</sup> (22)</td><td><b>2.00 × 10<sup>-200</sup></b></td><td>1.30 × 10<sup>-191</sup></td><td>1.30 × 10<sup>-191</sup></td><td>1.30 × 10<sup>-191</sup></td><td>1.30 × 10<sup>-191</sup></td></tr><tr><td>22</td><td>2.51 × 10<sup>-240</sup> (29)</td><td><b>1.19 × 10<sup>-290</sup></b></td><td>2.51 × 10<sup>-240</sup></td><td>2.51 × 10<sup>-240</sup></td><td>2.51 × 10<sup>-240</sup></td><td>2.51 × 10<sup>-240</sup></td></tr><tr><td>24</td><td>5.82 × 10<sup>-280</sup> (33)</td><td><b>1.51 × 10<sup>-305</sup></b></td><td>5.82 × 10<sup>-280</sup></td><td>5.82 × 10<sup>-280</sup></td><td>5.82 × 10<sup>-280</sup></td><td>5.82 × 10<sup>-280</sup></td></tr><tr><td>26</td><td>1.03 × 10<sup>-323</sup> (15)</td><td><b>5.31 × 10<sup>-340</sup></b></td><td>1.03 × 10<sup>-323</sup></td><td>1.03 × 10<sup>-323</sup></td><td>1.03 × 10<sup>-323</sup></td><td>1.03 × 10<sup>-323</sup></td></tr><tr><td>28</td><td>4.12 × 10<sup>-383</sup> (36)</td><td><b>2.06 × 10<sup>-383</sup></b></td><td>4.12 × 10<sup>-383</sup></td><td>4.12 × 10<sup>-383</sup></td><td>4.12 × 10<sup>-383</sup></td><td>4.12 × 10<sup>-383</sup></td></tr><tr><td>30</td><td>2.54 × 10<sup>-442</sup> (48)</td><td><b>7.47 × 10<sup>-485</sup></b></td><td>2.53 × 10<sup>-442</sup></td><td>2.54 × 10<sup>-442</sup></td><td>2.54 × 10<sup>-442</sup></td><td>2.54 × 10<sup>-442</sup></td></tr><tr><td>32</td><td>5.6 × 10<sup>-496</sup> (14)</td><td>5.6 × 10<sup>-496</sup></td><td>5.6 × 10<sup>-496</sup></td><td>5.6 × 10<sup>-496</sup></td><td>5.6 × 10<sup>-496</sup></td><td>5.6 × 10<sup>-496</sup></td></tr></tbody></table>
<table><thead><tr><th></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><th>R<sub>6</sub></th></tr></thead><tbody><tr><td>8</td><td>0.24</td><td>0</td><td>0.01</td><td>0.02</td><td>0.16</td><td>0.12</td></tr><tr><td>10</td><td>0.03</td><td>0.01</td><td>0.01</td><td>0.04</td><td>0.01</td><td>0.01</td></tr><tr><td>12</td><td>0.69</td><td>0.01</td><td>0.06</td><td>0.08</td><td>0.33</td><td>0.33</td></tr><tr><td>14</td><td>0.29</td><td>0.01</td><td>0.06</td><td>0.11</td><td>0.15</td><td>0.14</td></tr><tr><td>16</td><td>0.05</td><td>0.02</td><td>0.02</td><td>0.17</td><td>0.02</td><td>0.03</td></tr><tr><td>18</td><td>0.38</td><td>0.02</td><td>0.09</td><td>0.25</td><td>0.19</td><td>0.19</td></tr><tr><td>20</td><td>1.95</td><td>0.04</td><td>0.15</td><td>0.35</td><td>0.80</td><td>1.06</td></tr><tr><td>22</td><td>1.05</td><td>0.05</td><td>0.15</td><td>0.48</td><td>0.53</td><td>0.54</td></tr><tr><td>24</td><td>0.75</td><td>0.06</td><td>0.16</td><td>0.65</td><td>0.33</td><td>0.33</td></tr><tr><td>26</td><td>10.10</td><td>0.06</td><td>0.20</td><td>0.85</td><td>2.18</td><td>4.92</td></tr><tr><td>28</td><td>1.37</td><td>0.08</td><td>0.31</td><td>1.05</td><td>0.56</td><td>0.56</td></tr><tr><td>30</td><td>15.60</td><td>0.11</td><td>0.64</td><td>1.40</td><td>7.11</td><td>8.09</td></tr><tr><td>32</td><td>17.10</td><td>0.12</td><td>0.36</td><td>1.52</td><td>2.92</td><td>7.84</td></tr></tbody></table>

TABLE A.5 – Temps de calcul en secondes pour la constante de normalisation dans le cas ordonnée

Les quatre méthodes de rounding $R_3$, $R_4$, $R_5$ et $R_6$ donnent le même résultat, qui est le même que celui de la méthode de référence.

De plus, $R_3$ nous offre un meilleur temps de calcul que les autres méthodes.

## Dissociative

<table><thead><tr><th></th><th>R<sub>1</sub> (r<sub>max</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><th>R<sub>6</sub></th></tr></thead><tbody><tr><td>n= 8</td><td>3.57 × 10<sup>-32</sup> (12)</td><td><strong>1.68 × 10<sup>-32</sup></strong></td><td>3.57 × 10<sup>-32</sup></td><td>3.57 × 10<sup>-32</sup></td><td>3.57 × 10<sup>-32</sup></td><td>3.57 × 10<sup>-32</sup></td></tr><tr><td>n= 10</td><td>2.08 × 10<sup>-44</sup> (12)</td><td><strong>5.96 × 10<sup>-45</sup></strong></td><td>2.08 × 10<sup>-44</sup></td><td>2.08 × 10<sup>-44</sup></td><td>2.08 × 10<sup>-44</sup></td><td>2.08 × 10<sup>-44</sup></td></tr><tr><td>n= 12</td><td>1.77 × 10<sup>-63</sup> (12)</td><td>1.26 × 10<sup>-63</sup></td><td>1.77 × 10<sup>-63</sup></td><td>1.77 × 10<sup>-63</sup></td><td>1.77 × 10<sup>-63</sup></td><td>1.77 × 10<sup>-63</sup></td></tr><tr><td>n= 14</td><td>6.79 × 10<sup>-89</sup> (13)</td><td><strong>8.13 × 10<sup>-90</sup></strong></td><td>6.79 × 10<sup>-89</sup></td><td>6.79 × 10<sup>-89</sup></td><td>6.79 × 10<sup>-89</sup></td><td>6.79 × 10<sup>-89</sup></td></tr><tr><td>n= 16</td><td>2.03 × 10<sup>-117</sup> (14)</td><td><strong>8.42 × 10<sup>-120</sup></strong></td><td>2.03 × 10<sup>-117</sup></td><td>2.03 × 10<sup>-117</sup></td><td>2.03 × 10<sup>-117</sup></td><td>2.03 × 10<sup>-117</sup></td></tr><tr><td>n= 18</td><td>1.08 × 10<sup>-151</sup> (21)</td><td><strong>8.63 × 10<sup>-152</sup></strong></td><td><strong>7.79 × 10<sup>-152</sup></strong></td><td><strong>7.79 × 10<sup>-152</sup></strong></td><td>1.08 × 10<sup>-151</sup></td><td>1.08 × 10<sup>-151</sup></td></tr><tr><td>n= 20</td><td>3.88 × 10<sup>-179</sup> (44)</td><td><strong>8.61 × 10<sup>-189</sup></strong></td><td>2.18 × 10<sup>-179</sup></td><td>2.21 × 10<sup>-179</sup></td><td>3.88 × 10<sup>-179</sup></td><td>3.88 × 10<sup>-179</sup></td></tr><tr><td>n= 22</td><td>1.17 × 10<sup>-231</sup> (42)</td><td>1.03 × 10<sup>-231</sup></td><td>1.17 × 10<sup>-231</sup></td><td>1.17 × 10<sup>-231</sup></td><td>1.17 × 10<sup>-231</sup></td><td>1.17 × 10<sup>-231</sup></td></tr><tr><td>n= 24</td><td>6.07 × 10<sup>-276</sup> (30)</td><td><strong>3.51 × 10<sup>-280</sup></strong></td><td><strong>1.83 × 10<sup>-278</sup></strong></td><td><strong>1.83 × 10<sup>-278</sup></strong></td><td>6.07 × 10<sup>-276</sup></td><td>6.07 × 10<sup>-276</sup></td></tr><tr><td>n= 26</td><td>2.03 × 10<sup>-301</sup> (19)</td><td><strong>4.97 × 10<sup>-304</sup></strong></td><td>2.03 × 10<sup>-301</sup></td><td>2.03 × 10<sup>-301</sup></td><td>2.03 × 10<sup>-301</sup></td><td>2.03 × 10<sup>-301</sup></td></tr><tr><td>n= 28</td><td>3.25 × 10<sup>-361</sup> (50)</td><td>1.06 × 10<sup>-361</sup></td><td>2.07 × 10<sup>-361</sup></td><td>2.07 × 10<sup>-361</sup></td><td>3.23 × 10<sup>-361</sup></td><td>3.25 × 10<sup>-361</sup></td></tr><tr><td>n= 30</td><td>7.19 × 10<sup>-404</sup> (23)</td><td><strong>1.79 × 10<sup>-425</sup></strong></td><td>7.13 × 10<sup>-404</sup></td><td>7.13 × 10<sup>-404</sup></td><td>7.19 × 10<sup>-404</sup></td><td>7.19 × 10<sup>-404</sup></td></tr><tr><td>n= 32</td><td>2.80 × 10<sup>-476</sup> (36)</td><td><strong>4.99 × 10<sup>-495</sup></strong></td><td><strong>4.99 × 10<sup>-495</sup></strong></td><td><strong>4.99 × 10<sup>-495</sup></strong></td><td>2.80 × 10<sup>-476</sup></td><td>2.80 × 10<sup>-476</sup></td></tr></tbody></table>

TABLE A.6 – Constante de normalisation calculées ainsi que les $r_{max}$ trouvée dans le cas disso- ciative
<table><thead><tr><th></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><th>R<sub>6</sub></th></tr></thead><tbody><tr><td>8</td><td>0.02</td><td>0</td><td>0.01</td><td>0.02</td><td>0.01</td><td>0.01</td></tr><tr><td>10</td><td>0.09</td><td>0.01</td><td>0.03</td><td>0.04</td><td>0.06</td><td>0.05</td></tr><tr><td>12</td><td>0.11</td><td>0.01</td><td>0.05</td><td>0.07</td><td>0.05</td><td>0.05</td></tr><tr><td>14</td><td>0.22</td><td>0.02</td><td>0.08</td><td>0.12</td><td>0.12</td><td>0.12</td></tr><tr><td>16</td><td>0.34</td><td>0.02</td><td>0.14</td><td>0.18</td><td>0.17</td><td>0.18</td></tr><tr><td>18</td><td>1.05</td><td>0.03</td><td>0.16</td><td>0.27</td><td>0.52</td><td>0.52</td></tr><tr><td>20</td><td>4.03</td><td>0.04</td><td>0.20</td><td>0.38</td><td>1.15</td><td>2.05</td></tr><tr><td>22</td><td>5.28</td><td>0.05</td><td>0.29</td><td>0.52</td><td>1.37</td><td>2.78</td></tr><tr><td>24</td><td>3.88</td><td>0.07</td><td>0.54</td><td>0.69</td><td>1.50</td><td>1.55</td></tr><tr><td>26</td><td>1.70</td><td>0.07</td><td>0.27</td><td>0.84</td><td>0.70</td><td>0.69</td></tr><tr><td>28</td><td>27.10</td><td>0.07</td><td>0.43</td><td>0.96</td><td>4.42</td><td>9.79</td></tr><tr><td>30</td><td>3.45</td><td>0.10</td><td>0.30</td><td>1.21</td><td>1.63</td><td>1.64</td></tr><tr><td>32</td><td>15.90</td><td>0.14</td><td>0.71</td><td>1.57</td><td>5.16</td><td>7.07</td></tr></tbody></table>

TABLE A.7 – Temps de calcul en secondes pour la constante de normalisation dans le cas dissociative

Les deux méthodes $R_5$ et $R_6$ sont les méthodes les plus précises. Cependant, $R_5$ est un peu plus rapide pour le calcul.

## Hiérarchique 2

<table><thead><tr><th></th><th>R<sub>1</sub> (r<sub>max</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><th>R<sub>6</sub></th></tr></thead><tbody><tr><td>8</td><td>8.09 × 10<sup>-34</sup> (19)</td><td><b>3.25 × 10<sup>-46</sup></b></td><td>8.09 × 10<sup>-34</sup></td><td>8.09 × 10<sup>-34</sup></td><td>8.09 × 10<sup>-34</sup></td><td>8.09 × 10<sup>-34</sup></td></tr><tr><td>10</td><td>2.20 × 10<sup>-51</sup> (18)</td><td><b>1.60 × 10<sup>-51</sup></b></td><td>2.20 × 10<sup>-51</sup></td><td>2.20 × 10<sup>-51</sup></td><td>2.20 × 10<sup>-51</sup></td><td>2.20 × 10<sup>-51</sup></td></tr><tr><td>12</td><td>8.51 × 10<sup>-65</sup> (29)</td><td><b>5.18 × 10<sup>-113</sup></b></td><td>8.51 × 10<sup>-65</sup></td><td>8.51 × 10<sup>-65</sup></td><td>8.51 × 10<sup>-65</sup></td><td>8.51 × 10<sup>-65</sup></td></tr><tr><td>14</td><td>2.13 × 10<sup>-111</sup> (8)</td><td><b>1.44 × 10<sup>-120</sup></b></td><td>2.13 × 10<sup>-111</sup></td><td>1.44 × 10<sup>-102</sup></td><td>2.13 × 10<sup>-111</sup></td><td>2.13 × 10<sup>-111</sup></td></tr><tr><td>16</td><td>4.91 × 10<sup>-135</sup> (4)</td><td>4.91 × 10<sup>-135</sup></td><td>4.91 × 10<sup>-135</sup></td><td>4.91 × 10<sup>-135</sup></td><td>4.91 × 10<sup>-135</sup></td><td>4.91 × 10<sup>-135</sup></td></tr><tr><td>18</td><td>2.27 × 10<sup>-175</sup> (14)</td><td>2.27 × 10<sup>-175</sup></td><td>2.27 × 10<sup>-175</sup></td><td>2.27 × 10<sup>-175</sup></td><td>2.27 × 10<sup>-175</sup></td><td>2.27 × 10<sup>-175</sup></td></tr><tr><td>20</td><td>8.31 × 10<sup>-228</sup> (13)</td><td>8.31 × 10<sup>-228</sup></td><td>8.31 × 10<sup>-228</sup></td><td>8.31 × 10<sup>-228</sup></td><td>8.31 × 10<sup>-228</sup></td><td>8.31 × 10<sup>-228</sup></td></tr><tr><td>22</td><td>1.36 × 10<sup>-265</sup> (17)</td><td><b>1.61 × 10<sup>-545</sup></b></td><td>1.36 × 10<sup>-265</sup></td><td>1.36 × 10<sup>-265</sup></td><td>1.36 × 10<sup>-265</sup></td><td>1.36 × 10<sup>-265</sup></td></tr><tr><td>24</td><td>9.83 × 10<sup>-328</sup> (46)</td><td><b>1.60 × 10<sup>-353</sup></b></td><td>9.83 × 10<sup>-328</sup></td><td>9.83 × 10<sup>-328</sup></td><td>9.83 × 10<sup>-328</sup></td><td>9.83 × 10<sup>-328</sup></td></tr><tr><td>26</td><td>1.79 × 10<sup>-369</sup> (37)</td><td>1.79 × 10<sup>-369</sup></td><td>1.79 × 10<sup>-369</sup></td><td>1.79 × 10<sup>-369</sup></td><td>1.79 × 10<sup>-369</sup></td><td>1.79 × 10<sup>-369</sup></td></tr><tr><td>28</td><td>1.68 × 10<sup>-418</sup> (37)</td><td><b>6.96 × 10<sup>-549</sup></b></td><td>1.68 × 10<sup>-418</sup></td><td>1.68 × 10<sup>-418</sup></td><td>1.68 × 10<sup>-418</sup></td><td>1.68 × 10<sup>-418</sup></td></tr><tr><td>30</td><td>2.05 × 10<sup>-498</sup> (17)</td><td><b>1.11 × 10<sup>-532</sup></b></td><td>2.05 × 10<sup>-498</sup></td><td>2.05 × 10<sup>-498</sup></td><td>2.05 × 10<sup>-498</sup></td><td>2.05 × 10<sup>-498</sup></td></tr><tr><td>32</td><td>1.52 × 10<sup>-586</sup> (21)</td><td>1.38 × 10<sup>-586</sup></td><td>1.52 × 10<sup>-586</sup></td><td>1.52 × 10<sup>-586</sup></td><td>1.52 × 10<sup>-586</sup></td><td>1.52 × 10<sup>-586</sup></td></tr></tbody></table>

TABLE A.8 – Constante de normalisation calculées ainsi que les $r_{max}$ trouvé dans le cas hiérarchique 2
<table><thead><tr><th></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><th>R<sub>6</sub></th></tr></thead><tbody><tr><td>8</td><td>0.07</td><td>0</td><td>0.01</td><td>0.01</td><td>0.04</td><td>0.04</td></tr><tr><td>10</td><td>0.08</td><td>0.01</td><td>0.02</td><td>0.04</td><td>0.03</td><td>0.03</td></tr><tr><td>12</td><td>0.37</td><td>0.01</td><td>0.03</td><td>0.06</td><td>0.20</td><td>0.20</td></tr><tr><td>14</td><td>0.05</td><td>0.01</td><td>0.02</td><td>0.10</td><td>0.03</td><td>0.03</td></tr><tr><td>16</td><td>0.03</td><td>0.01</td><td>0.02</td><td>0.14</td><td>0.02</td><td>0.02</td></tr><tr><td>18</td><td>0.19</td><td>0.02</td><td>0.07</td><td>0.21</td><td>0.10</td><td>0.10</td></tr><tr><td>20</td><td>0.49</td><td>0.03</td><td>0.14</td><td>0.30</td><td>0.24</td><td>0.24</td></tr><tr><td>22</td><td>0.65</td><td>0.04</td><td>0.15</td><td>0.41</td><td>0.32</td><td>0.32</td></tr><tr><td>24</td><td>6.73</td><td>0.05</td><td>0.12</td><td>0.53</td><td>1.20</td><td>2.20</td></tr><tr><td>26</td><td>4.95</td><td>0.06</td><td>0.30</td><td>0.73</td><td>1.66</td><td>2.10</td></tr><tr><td>28</td><td>4.73</td><td>0.07</td><td>0.17</td><td>0.93</td><td>1.36</td><td>2.16</td></tr><tr><td>30</td><td>0.69</td><td>0.08</td><td>0.17</td><td>1.29</td><td>0.34</td><td>0.37</td></tr><tr><td>32</td><td>3.21</td><td>0.12</td><td>0.40</td><td>1.49</td><td>1.59</td><td>1.62</td></tr></tbody></table>

TABLE A.9 – Temps de calcul en secondes pour la constante de normalisation dans le cas hiérarchique 2

L'ensemble des résultats suggèrent que les six méthodes de rounding s'ordonnent par précision dé-
croissance avec l'ordre suivant : $R_6 \backslash R_5$, $R_3$, $R_4$ et enfin $R_2$
# Annexe B

## Précision des calculs des marginales et la constante de normalisation

Nous présentons dans cette annexe un théorème pour borner une TT-matrice **B** obtenue par la compression d'une TT-matrice **A** via l'algorithme du rounding avec précision prescrite $\epsilon$. Nous allons utiliser ce théorème pour présenter un corollaire qui porte sur l'erreur de la compression des produits de TT-matrices par l'algorithme du rounding avec précision prescrite.

Nous utiliserons ce corollaire afin de présenter la précision des marginales et de la constante de normalisation calculées par l'approche TT.

### B.1 Compression d'une TT-matrice par le rounding

**Théorème 14.** Soit **A** une TT-matrice. Soit **B** une approximation de **A** obtenue par l'algorithme rounding avec précision prescrite $\epsilon$, c'est-à-dire :

$$ \mathbf{B} = \text{rounding}(\mathbf{A}, \epsilon) $$

alors :

$$ (1 - \epsilon) \|\mathbf{A}\| \le \|\mathbf{B}\| \le (1 + \epsilon) \|\mathbf{A}\| \quad (B.1.1) $$

*Démonstration.* Soit **A** une TT-matrice, alors d'après Oseledets [Ose11], si on a $\mathbf{B} = \text{rounding}(\mathbf{A}, \epsilon)$ alors l'erreur s'exprime :

$$ \|\mathbf{A} - \mathbf{B}\| \le \epsilon \|\mathbf{A}\| \quad (B.1.2) $$

Comme $\|\cdot\|$ est une norme, alors elle vérifie l'inégalité triangulaire. Donc,

$$ \|\mathbf{A}\| - \|\mathbf{B}\| \le \|\mathbf{A} - \mathbf{B}\| $$

D'après B.1.2, nous obtenons

$$ \|\mathbf{A}\| - \|\mathbf{B}\| \le \epsilon \|\mathbf{A}\| \quad (B.1.3) $$

Donc,

$$ (1 - \epsilon) \|\mathbf{A}\| \le \|\mathbf{B}\| \quad (B.1.4) $$

Nous procédons de la même manière pour avoir la borne supérieure de **B**.

$$ \|\mathbf{A} - \mathbf{B}\| \le \epsilon \|\mathbf{A}\| $$

Donc,

$$ \|\mathbf{B} - \mathbf{A}\| \le \epsilon \|\mathbf{A}\| $$

Nous développons,

$$ \|\mathbf{A}\| + \|\mathbf{B} - \mathbf{A}\| \le \|\mathbf{A}\| + \epsilon \|\mathbf{A}\| $$
D'après la propriété de l'inégalité triangulaire : $||\mathbf{A} + \mathbf{B} - \mathbf{A}|| \le ||\mathbf{A}|| + ||\mathbf{B} - \mathbf{A}||$. Donc,

$$||\mathbf{A} + \mathbf{B} - \mathbf{A}|| \le (1 + \epsilon)||\mathbf{A}||$$

Donc,

$$||\mathbf{B}|| \le (1 + \epsilon)||\mathbf{A}|| \qquad (B.1.5)$$

Enfin, nous pouvons borner la norme de l'erreur à partir des formules (B.1.5) et (B.1.4) :

$$(1 - \epsilon)||\mathbf{A}|| \le ||\mathbf{B}|| \le (1 + \epsilon)||\mathbf{A}|| \qquad (B.1.6)$$

**Corollaire 15** (Produit de TT-matrices). Soient $\mathbf{A}, \mathbf{B}, \mathbf{C}$ trois TT-matrices et, $\epsilon \in \mathbb{R}$.

Soient $\mathbf{D} = \mathbf{A} \times \mathbf{B} \times \mathbf{C}$ et $\hat{\mathbf{D}} = \text{rounding}(\text{rounding}(\mathbf{A} \times \mathbf{B}, \epsilon) \times \mathbf{C}, \epsilon)$

Alors,

$$(1 - \epsilon)^2 \mathbf{D} \le \hat{\mathbf{D}} \le (1 + \epsilon)^2 \mathbf{D} \qquad (B.1.7)$$

La démonstration de ce corollaire vient du théorème 14. Ce corollaire peut se généraliser facilement pour un produit de $n \in \mathbb{N}$ TT-matrices.

## B.2 Erreur relative dans le calcul des marginales et la constante de normalisation

Nous avons expliqué dans le chapitre 5 comment appliquer la procédure de Novikov [NROV14] pour des SBM et comment faire le calcul des marginales par produits de TT-matrices. Nous rappelons la formule exacte pour le calcul des marginales unaires d'un modèle avec $n$ individus.

$$\mathbf{p}_i(z_i) = \frac{\mathbf{B}_1 \times \dots \times \mathbf{B}_{i-1} \mathbf{A}_i[z_i] \mathbf{B}_{j+1} \times \dots \times \mathbf{B}_n}{W}$$

Nous avons expliqué dans la fin de ce chapitre ainsi que dans le chapitre 7 l'intérêt de faire des compressions par rounding à chaque itération. À partir de cela, le calcul d'une marginale unaire se fait :

$$\hat{\mathbf{p}}_i(z_i) = \frac{\text{rounding}(\dots\text{rounding}(\mathbf{B}_1 \times \mathbf{B}_2, \epsilon) \times \dots \times \mathbf{B}_{i-1}, \epsilon) \times \mathbf{A}_i[z_i], \epsilon) \times \mathbf{B}_{j+1}, \epsilon) \times \dots \times \mathbf{B}_n, \epsilon)}{W}$$

Si nous généralisons la formule (B.1.7), nous obtienons une de l'erreur :

$$(1 - \epsilon)^{n-2} W \times ||\mathbf{p}_i(z_i)|| \le W ||\hat{\mathbf{p}}_i(z_i)|| \le (1 + \epsilon)^{n-2} W \times ||\mathbf{p}_i(z_i)|| \qquad (B.2.1)$$

Ainsi,

$$(1 - \epsilon)^{n-2} ||\mathbf{p}_i(z_i)|| \le ||\hat{\mathbf{p}}_i(z_i)|| \le (1 + \epsilon)^{n-2} ||\mathbf{p}_i(z_i)|| \qquad (B.2.2)$$

Alors, la formule (B.2.2) peut s'exprimer :

$$(1 - \epsilon)^{n-2} ||p_i(z_i)|| \le ||\hat{p}_i(z_i)|| \le (1 + \epsilon)^{n-2} ||p_i(z_i)|| \qquad (B.2.3)$$

Ainsi, l'erreur relative des marginales unaires s'exprime :

$$(1 - \epsilon)^{n-2} \le \frac{||\hat{p}_i(z_i)||}{||p_i(z_i)||} \le (1 + \epsilon)^{n-2}$$

De la même manière, nous pouvons borner l'erreur des marginales binaires de cette manière :

$$(1 - \epsilon)^{n-2} \le \frac{||\hat{p}_{i,j}(z_i, z_j)||}{||p_{i,j}(z_i, z_j)||} \le (1 + \epsilon)^{n-2} \qquad (B.2.4)$$
Article

Nous présentons dans cette annexe l'article "Does clustering of DNA barcodes agree with botanical classification directly at high taxonomic levels ? Trees in French Guiana as a case study" publié dans *Molecular Ecology Resources*.
# Evaluating the adequacy between morphological-based and molecular-based inventories at high taxonomic level

Mohamed Anwar Abouabdallah¹,³, Nathalie Peyrard²,*, and Alain Franc ¹,³

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

²Université de Toulouse, INRAE, UR MIAT, F-31320, Castanet-Tolosan, France

³Pleiade, EPC INRIA-INRAE-CNRS, Université de Bordeaux 33405

Talence France

*corresponding author: nathalie.peyrard@inrae.fr

**Running title:** Barcoding angiosperms directly at order level

## Abstract

Characterising biodiversity is one of the big challenges for the coming decades.

Most diversity has not been morphologically described and metabarcoding is taking over from morphological-based taxonomy to further develop inventories. The approaches have been cross-validated at the level of species and OTUs. However, many known species are not listed in reference databases. One path to speed up inventories using metabarcoding could be to directly identify individuals at coarser taxonomic levels. We therefore studied whether morphology-based and molecular-based approaches are in adequacy at genus, family and order levels. The comparison
was performed on a dataset of approximately 1,500 trees from a French Guianese
experimental plot, well-established both botanically and molecularly. We used Ag-
glomerative Hierarchical Clustering (with Ward, Complete and Single Linkage) and
Stochastic Block Models (SBM), with two dissimilarity measures, one based on Smith-
Waterman local alignment scores, and one based on kmers. The adequacy between
morphological-based and molecular-based classifications can be good to very good
at taxonomic levels beyond species, even though it decreases when taxonomic levels
increase, or when using the tetramer-based distance. Adequacy is correlated with
the entropy of morphological-based taxonomy and with the ratio of the mean within-
and mean between-groups dissimilarities. Ward method globally leads to the best
adequacy whereas Single Linkage can show poor behaviours. SBM provides a useful
signal as to whether or not the dissimilarities are structured according to the botan-
ical groups. We assert that it is time for scaling up the Ward and SBM methods to
analyse massive datasets established by metabarcoding.

**Keywords:** taxonomy; barcoding; Stochastic Block Model; clustering; Ward method; French Guianan Trees

# 1 Introduction

Numerical taxonomy and hierarchical clustering have coevolved since the 1960s' (Cole, 1969; Sneath and Sokal, 1973). Both approaches rely on the assumption that the diversity of life for taxonomy, or patterns in distances between some items in clustering, are organized as a nested hierarchy, modelled as a tree. This approach has survived the revolution of molecular-based taxonomy (Hillis et al., 1996) and molecular phylogenies (Felsenstein, 2004; Yang, 2006), with a current revival due to barcoding (Floyd et al., 2002; Hebert et al., 2003), and, at the same time, metabarcoding (López-García et al., 2001; Sogin et al., 2006; Hajibabaei et al., 2011; Taberlet et al., 2012; Kermarrec et al., 2013). As far as morphological-based taxonomy is concerned, most of the diversity in many clades of
organisms is still unknown (Leray and Knowlton, 2015, point out that between 33% and 91% of all marine biodiversity has never been named), and much efforts has been devoted to speeding up the process of producing large inventories with metabarcoding by bypassing identified blocking points (Bik et al., 2012).

Therefore, the notion of OTU (which stands for Operational Taxonomic Unit) has been coined (Floyd et al., 2002; Blaxter et al., 2005). Such units are produced by clustering sets of barcodes by aggregation at a level assumed to be similar to the level of species in morphological-based classifications. The authors in Blaxter et al. (2005) emphasize that they are "agnostic" as to whether OTU are species or not. Identifying OTUs in an environmental sample and organising molecular diversity as the frequency of OTUs make it possible to produce molecular-based inventories at previously unparalleled speed.

A classical approach is therefore to build OTUs and to map them on reference databases that contain reference barcodes. A standard tool for mapping is BLAST (Altschul et al., 1990), but other and more sophisticated solutions exist (e.g., the use of Bayesian Phylogenetics, Munch et al., 2008). When taxonomic expertise and references exist at the species level, the adequacy between molecular and morphological-based classification can be excellent (Ji et al., 2013). It may happen that such a comparison is not feasible when morphological-based taxonomy is unknown or when only partial references exist. Leray and Knowlton (2015) report in their study that less than 12% of their OTUs matched with GenBank or BOLD. The same observation was made in White et al. (2010) regarding intestinal microbial flora. Hence, most of inventories with supervised learning are made at a grain often much coarser than the genus/species level when the reference databases on which a match is looked for are annotated at high taxonomic levels like that of phylum.

Trying to complete databases at the species level is known to be highly time-consuming. We privilege another approach here and suggest directly identifying coarse levels in taxonomy present in the inventory by clustering<fnref n="1" /> the barcodes in a limited number of groups.

<footnote n="1">¹In this article, the term *clustering* makes reference to any numerical method for the unsupervised</footnote>
Supervised learning can then be used in a second step to annotate the groups found as taxa at deep taxonomic levels. This raises the question of adequacy between morphological-based and molecular-based taxonomy when clusters of sequences are built at a level coarser than species, like, e.g., class or order. Comparing morphological-based classifications and OTU produced by barcodes clustering has been thoroughly studied (see, e.g., White et al., 2010). Several methods have been recently designed and widely used for delineating species on the basis of barcodes (Pons et al., 2006; Fontaneto et al., 2008; Puillandre et al., 2012; Talavera et al., 2013; Zhang et al., 2013). However, to our knowledge, such a question has seldom been addressed at coarser levels like orders.

Our objective here is to study whether the clustering of barcodes in molecular-based taxonomy makes it possible to directly recover all the taxonomic levels beyond the species present in a sample, and, if so, with which tool, accuracy and robustness.

We performed this study on a well-known dataset where deep knowledge of the taxonomy is available for several orders of angiosperms. We selected flowering plants because the botanical classification is well known, both morphologically (it is organised as a nested system of different taxonomic levels as a classification system) and molecularly with the APG initiative, even if it is under continuous revision (The Angiosperm Phylogeny Group et al., 2016). The dataset itself (Caron et al., 2019) is composed of some 1500 trees from French Guiana which have been botanically identified and sequenced (chloroplastic marker trnH-psbA). An advantage is that we do not have to deal with the computational burden of treatment of metabarcoding data, and we can therefore concentrate on the analysis of adequacy. The question of the scaling of those methods which will show itself be effective on this well-known and manageable dataset, will be the object of further studies, and some leads are given in the Discussion section of this article.

It can be expected that there is not a clear answer to the degree of adequacy between the grouping of the individuals, while the term classification designates the method's output, i.e. the partition of individuals into classes.
two visions of classification (morphology-based or molecular-based), especially over several
taxonomic levels. There may be favourable situations where the adequacy is strong, and
others where the two classifications are surprisingly quasi-independent of each other. To
identify potential factors that may explain variations in adequacy in our study: (i) we
varied the taxonomic level at which the classification is performed (order, family, genus,
species), (ii) we used two definitions of dissimilarity between sequences; and, finally, (iii)
we considered four numerical methods for the clustering of the molecular data.

More specifically, we first worked with 30 distinct sub-samples of the whole dataset,
each sub-sample being either all the individuals of an order or of a family, and for each of
them, we compared the botanical classification of the individuals at the next taxonomic
level with the molecular-based classifications. In a second step, we studied whether the
mean behaviour observed from these replicates is recovered when the set of individuals to
be classified carries more diversity, by comparing the botanical classification of the whole
dataset (into orders, families, genera and species) with the molecular-based classifications.

Regarding the computation of dissimilarity between sequences of unequal length, the
Smith-Waterman score is a classical choice (Smith and Waterman, 1981; Gusfield, 1997).
Even if this algorithm relies on dynamic programming, thus making it very efficient (and
accurate), its complexity is in $O(n^2l^2)$ if n is the number of barcodes or reads, and l their
length. This becomes prohibitive for large datasets. A classical way to circumvent this
difficulty is to use kmer-based distances (Sun et al., 2009; Mahé et al., 2014), a priori with
a decrease in the quality of the estimation of the dissimilarity, but much faster to compute.
A comparison between Smith-Waterman scores and kmer-based distances can be found
in Sun et al. (2009). The question here is to know whether the loss in quality remains
acceptable and does not lead to a decrease in adequacy between the botanical and the
molecular-based classifications.

If the morphology-based taxonomic classification is a priori unique, this is not true
for a molecular-based classification. A diversity of softwares for implementing hierarchical
clustering has been proposed for a decade or so in metabarcoding with the objective of efficient scaling with respect to the growing size of environmental datasets. This includes Uclust (Edgar, 2010), ESPRIT (Sun et al., 2009) and SWARM (Mahé et al., 2014) which make it possible to cluster millions of barcodes on a laptop. All of the Hierarchical Clustering algorithms mentioned above rely on heuristics (like computing kmer-based distances, considering short distances only i.e. working with sparse distance matrices) to make computation feasible within a reasonable time with reasonable memory. In this study, we focus on understanding the adequacy (or not) between molecular-based numerical classification and botanical classification, without computational constraints. We therefore consider Aggregative Hierarchical Clustering (AHC), whose above-mentioned algorithms can be seen as heuristic versions for scaling up, with three different aggregation methods (Single Linkage, Complete Linkage, Ward (Murtagh, 1983; Müllner, 2013)). Statistical models like Bayesian classifiers with mixture models have also been considered to cluster sequences (Hao et al., 2011). To extend the scope of statistical modeling in molecular-based taxonomy, we explore here the potential interest of a model-based clustering method, the Stochastic Block Model (SBM, Holland et al., 1983; Daudin et al., 2008; Lee and Wilkinson, 2019) as an alternative to AHC. SBM are already widely applied with success in domains like the social sciences (Barbillon et al., 2017), the analysis of ecological interaction networks (Miele and Matias, 2017) and neurology (Faskowitz et al., 2018). They rely on a more flexible definition of a cluster than AHC (searching for general groups and not just communities), and we hypothesized that SBM and AHC could be complementary in their capacity to distinguish meaningful groups of individuals in an inventory.

In the following section, we provide a brief description of the dataset. We also describe the method. Results on the quality of the adequacy between molecular-based and morphological-based classifications obtained on replicates are presented in Section 3.2, the results obtained on the whole dataset are presented in Section 3.3.
141

# 2 Materials and methods

142

## 2.1 Dataset and computation of dissimilarities

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

This study was conducted on a dataset built from a collection of some 1,500 trees located in the "Piste de Saint-Elie" experimental plot in French Guiana, mainly composed of low-land tropical rainforest (Sabatier et al., 1997). The data used here are part of a dataset gathered for the study published in Caron et al. (2019), which focused on adequacy or not between botanical-based and molecular-based classification at the species level over a wide range of diversity along the angiosperm tree. The principal result in Caron et al. (2019) is that molecular-based clustering is highly consistent with species delineation in a majority of cases, and that introgression or incomplete lineage sorting are the most likely explanations in the case of non-adequacy. We focus here on a similar question, at the level of genera, families and orders. The main elements for the material are recalled here, and the reader can refer to Caron et al. (2019) for details. Among this dataset, 1458 individual trees were selected for this study. For each tree, we used the botanical name as given by field botanists working with the Cayenne IRD Herbarium, and a sequence of chloroplastic marker trnH-psbA. These trees encompass 20 orders, 56 families, 182 genera and 428 species. Two 1458 × 1458 matrices of pairwise distances or dissimilarities between sequences were built, a first one using the Smith-Waterman algorithm for local sequence alignment (Smith and Waterman, 1981), and a second one for the distance between kmers distributions. The local alignment score is the most precise quantification of genetic dissimilarities between sequences, but it is time consuming. Several methods for building OTUs therefore rely on alternatives to local alignment scores. A classical way to circumvent this computational burden is to build kmer histograms for each sequence, and then compute the distance between the histograms. A kmer is a contiguous sub-sequence of length *k* in a given sequence. We selected short kmers here with *k* = 4 (there are 4⁴ = 64 different tetramers which is a good compromise between longer ones with more resolution, but too
sparse histograms, or smaller ones with coarse resolution and less empty categories). We designed an efficient algorithm that counts the frequency of each tetramer in each sequence, and a short program that computes a distance between any pair of frequency distributions as the $ℓ^1$ norm, i.e., the sum of absolute values of difference per tetramer.

The dataset used in the rest of this paper is composed of three files (see Frigerio et al. (2021)):

* a csv file of botanical names for each sequence for order, family, genus and species;

* a csv file of pairwise dissimilarities computed with the Smith-Waterman algorithm;

* a csv file of pairwise distances based on the comparison of tetramer histograms (same format as Smith-Waterman dissimilarities)

## 2.2 Visualisation of the whole dataset using MDS

A preliminary step is to propose a global picture of the dataset based on the dissimilarity matrices, without a classification objective. Multidimensional Scaling (MDS) is a method that, once a dissimilarity matrix is given, builds a point cloud in a Euclidean space of prescribed dimension where each point corresponds to a sequence, and such that the Euclidean distance between any two points is as close as possible to the dissimilarity given in the matrix (Cox and Cox, 2001; Izenman, 2008). In our case, we selected the so-called Classical Scaling, as proposed initially in Torgerson (1952). It is expected that the projection of the point cloud on the first axis encompasses much of the information about the structure of the point cloud. MDS was run with the dissimilarity matrices built with the Smith-Waterman algorithm and as distances between kmer histograms. To perform MDS, we have filtered the whole data set by keeping sequences in orders with 15 specimen or more only.