# Model Based Co-Clustering: High Dimension and Estimation Challenges

Christophe Biernacki, Julien Jacques, Christine Keribin

► To cite this version:

Christophe Biernacki, Julien Jacques, Christine Keribin. Model Based Co-Clustering: High Dimension and Estimation Challenges. RMR 2024: Modèles statistiques pour des données dépendantes et applications, Jun 2024, Rouen, France. hal-04867840

HAL Id: hal-04867840

https://inria.hal.science/hal-04867840v1

Submitted on 6 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# MODEL BASED CO-CLUSTERING: HIGH DIMENSION & ESTIMATION CHALLENGES

Christophe BIERNACKI<fnref n="1" />, Julien JACQUES<fnref n="2" />, Christine KERIBIN<fnref n="3" />

RMR 2024 : Modèles statistiques pour des données dépendantes et applications - 19-21 juin 2024, Rouen, France

<footnote n="1">¹Inria, CNRS, Laboratoire de Mathématiques Painlevé, U. Lille, F59650 V. d'Ascq, France</footnote>

<footnote n="2">²U. Lyon 2, ERIC, UR 3083, 5 Avenue Pierre Mendès-France, F69676 Bron Cedex, France</footnote>

<footnote n="3">³U. Paris-Saclay, CNRS, Inria, Laboratoire de mathÂ©matiques dÂOrsay, F91405 Orsay, France</footnote>
<figure><img src="image_2.png" /></figure>

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

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

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

# Outline

1. Motivation

2. Latent Block Model

3. Estimation challenges

4. Co-clustering for HD

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

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

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

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

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

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

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

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

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

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

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

<figure><img src="image_17.png" /></figure>
<figure><img src="image_18.png" /></figure>

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

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

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

# Clustering

## Unsupervised ML framework

* data set $\mathbf{x} = (\mathbf{x}_1, \dots, \mathbf{x}_n) \in \mathcal{X}^n$ with $dim(\mathcal{X}) = d$

* to partition into $K$ clusters $G_1, \dots, G_K$

* classification matrix $\mathbf{z} = (z_1, \dots, z_n) : z_{ik} = \mathbb{1}_{i \in G_k}$

## Model based clustering (MBC) reformulates cluster analysis in a well-posed estimation problem

* $\mathbf{x}_i$ are iid observations of a mixture pdf

$$f(\cdot; \boldsymbol{\theta}) = \sum_{k=1}^{K} \pi_k \varphi(\cdot; \boldsymbol{\alpha}_k); \quad \sum_{i=1}^{K} \pi_K = 1, \pi_k \ge 0; \quad \mathbb{P}(z_{ik} = k) = \pi_k$$

* estimation with EM algorithm, partition is deduced from the conditional probability $\rho(z_i|\mathbf{x}_i; \hat{\boldsymbol{\theta}})$, typically from the MAP principle, as soon as we have an estimate $\hat{\boldsymbol{\theta}}$ of $\boldsymbol{\theta}$

<figure><img src="image_22.png" /></figure>
High dimension (HD) settings: $d \sim n, d \gg n$

More and more variables (d) and not so much observations (n)

► marketing: $d \sim 10^2$ for $n = 10^4 - 10^6$

► microarray gene expression: $d \sim 10^2 - 10^4$ for $n = 20 - 50$ stress conditions

► fMRI images: $d \sim 10^4 - 10^5$ voxels for $n = 50 - 100$ patients

► text mining, curve, internet logs...

Moreover, variables can be of mixed types

► different modes: qualitative, quantitative, functional

► different semantic

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

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

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

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

Clustering in HD: curse (and blessing?)

► **curse**: density estimation quality decreases with *d* for a same *n*

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

Kullback($f, \hat{f}$) for a two-component *d*-variate Gaussian mixture

$$x_j | z_j \sim \mathcal{N}_d(\mu_{z_j}, \Sigma); \pi_1 = \pi_2 = 0.5$$

Basic case: $\mu_1 = 0_d; \mu_2 = 1_d; \Sigma = I_d$

► **blessing**: but components are (in that case) more and more separated

► which impact on the **misclassification** error?

- when the variables are **correlated**, redundant, non informative?

- when $d \gg n$?

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Basic case: HD is a blessing for clustering

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

<figure><img src="image_50.png" /></figure>
<figure><img src="image_51.png" /></figure>

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

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

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

**Basic case: HD is a blessing for clustering**

Each new variable is informative, sufficiently large $n \sim d$

► theoretical classification error decreases when $d$ grows: $err_{theo} = \Phi(-\sqrt{d}/2)$

basic, c= 0 eps= 1 nbobs= 30

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

A two-component $d$-variate Gaussian mixture

$$x_i | z_i \sim N_d(\mu z_i, \Sigma); \pi_1 = \pi_2 = 0.5$$

Basic case: $\Sigma = I_d; \mu_1 = 0_d; \mu_2 = 1_d$

► and empirical error rate also decreases with $d$

► no model bias, classification ($\hat{z}$) variance decreases

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

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

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

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

But HD can be less favorable for clustering

(Many) less informative variables

► theoretical error slowly decreases when *d* grows: $err_{theo} = \Phi(-||\mu_1 - \mu_2||/2)$

less informative, c=0 eps= 1 nbobs= 100

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

A two-component *d*-variate Gaussian mixture

$$x_j|z_j \sim N_d(\mu_{z_j}, \Sigma); \pi_1 = \pi_2 = 0.5$$

Less informative: $\Sigma = I_d$; $\mu_1 = 0_d$.

$$\mu_2 = \left(1, \frac{1}{2^2}, \dots, \frac{1}{d^2}\right)$$

► but empirical error rate does not decrease any more...

► no model bias, classification ($\tilde{z}$) variance increases

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

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

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

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

# Bias/variance trade-off for clustering

## Correlated variables

* theoretical error decreases $err_{theo} = \Phi(-||\mu_1 - \mu_2||_{\Sigma} - 1/2)$

basic, c= 0.2 eps= 1 nbobs= 30

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

A two-component *d*-variate Gaussian mixture

$$x_j | z_j \sim N_d(\mu_{z_j}, \Sigma); \pi_1 = \pi_2 = 0.5$$

Basic case+ correlated variables (correlation c):

$$\Sigma = \Sigma(c); \mu_1 = 0_d; \mu_2 = 1_d$$

* high error when clustered with the **real** (corr. G.) mixture model **higher** than with the **biased** (diag. G) mixture model, since corr. G. has too much variance

* error decreases with the **biased** model diag. G, which does not dammaged the separation

<figure><img src="image_68.png" /></figure>
<figure><img src="image_69.png" /></figure>

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

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

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

# Some alternatives for reducing variance

* Dimension reduction in non canonical space (PCA-like typically)

* Dimension reduction in the canonical space (variable selection)

- Gaussian mixture with 3 variable-clusters (informative, non informative, linearly dependent) [Maugis et al (2009)]

- latent class analysis [Fop et al (2018)]

* Model parsimony in the initial HD space (constraints on model parameters)

- parsimonious Gaussian mixture [Celeux and Govaert, '95]

- mixture of factor analyzers [[Ghahramani and Hinton, 97], [McLachlan et al], 03]

- HD Gaussian models [[Bouveyron et al, 07]]

↔ trade-off between “meaningful” modeling and what can be estimated in practice

<figure><img src="image_73.png" /></figure>
<figure><img src="image_74.png" /></figure>

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

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

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

# An alternative for reducing the variance

## Main message:

* clustering is a way for dealing with large *n*
* why not also reusing this idea for large *d*?

## Co-clustering (CC)

* simultaneously clusters rows and columns
* performs parsimony of row clustering through variable clustering
* parsimonious, model-based, independent of data type, available

<figure><img src="image_78.png" /></figure>
<figure><img src="image_79.png" /></figure>

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

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

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

From clustering to co-clustering

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

(1)

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

(2)

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

(3)

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

(4)

$$ (n = 10) \times (d = 7) $$

$$ (K = 3) \times (d = 7) $$

$$ (K = 3) \times (L = 3) $$

[Govaert (2008)]

<figure><img src="image_87.png" /></figure>
<figure><img src="image_88.png" /></figure>

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

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

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

Interpreting LBM as an MBC dimension reduction method<fnref n="1" />

(See the Latent Block Model (LBM) definition in the next section...)

► PCA for the *i*-th data individual $\mathbf{x}_i$: acts as a dimension reduction method

$$ \mathbf{x}_i = \sum_{j=1}^{d} x_j^l \mathbf{e}^l = \sum_{j=1}^{d} a_j^l \mathbf{u}^l \approx \sum_{j=1}^{J} a_j^l \mathbf{u}^l $$

► specific CC case reduced to $K = 1$: again, acts as a dim. reduction method

$$ \mathbf{x}_i = \sum_{j=1}^{d} x_j^l \mathbf{e}^l = \sum_{j=1}^{d} (\mu_{\tilde{w}_j} + \varepsilon_j^l) \mathbf{e}^l = \sum_{\ell=1}^{L} \mu_{\ell} \mathbf{v}^{\ell} + \mathbf{r}_i \approx \sum_{\ell=1}^{L} \mu_{\ell} \mathbf{v}^{\ell} $$

with $\tilde{w}_j = \ell \Leftrightarrow w_j = 1$, $\varepsilon_j^l \sim \mathcal{N}(0, \sigma_{\tilde{w}_j}^2)$, $\mathbf{v}^l = \sum_{\{j: \tilde{w}_j = \ell\}} \mathbf{e}^l$ and $\mathbf{r}_i = \sum_{j=1}^{d} \varepsilon_j^l \mathbf{e}^l$

↔ return now to the clustering task:

- sequential clustering method: (1) PCA and then (2) MBC

- simultaneous clustering method: only CC, since dim. reduction is included

<footnote n="1">¹ [C. Biernacki, C. Keribin, J. Jacques, A survey on Model-Based Co-Clustering, JOC 2023] □ ▶ ◀ ▶ ◁ ▶ ◁ ▶ ◁ ▶ ▶ ▶ ▶ ▶ ▶</footnote>
## Outline

1. Motivation

2. Latent Block Model

3. Estimation challenges

4. Co-clustering for HD

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

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

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

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

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

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

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

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

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

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

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

<figure><img src="image_103.png" /></figure>
# Unsupervised block clustering framework

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

* **co-clustering (CC):** to find a block clustering structure simultaneously on rows and columns

* **$H_1$: blocks** form a Cartesian product of a row-partition in $K$ row-groups by a column-partition in $L$ column-groups

$z = (z_{ik})$ where $z_{ik} = 1$ if row $i$ belongs to row-group $k$

$w = (w_{j\ell})$ where $w_{j\ell} = 1$ if column $j$ belongs to column-group $\ell$

* **traditional applications:** huge data sets arising in recommendation systems (binary), genomic data analysis (Gaussian), text mining (Poisson),...

* **our proposed application:** dimension reduction for HD clustering

<figure><img src="image_105.png" /></figure>
Latent Block Model (LBM): a mixture model<fnref n="2" />

$H_2$: row and column labels are independently assigned :

$$
z_j^{i,l,d} \sim \mathcal{M}(1, \pi) \perp w_j^{i,l,d} \sim \mathcal{M}(1, \rho)
$$

$H_3$: the $n \times d$ blocks $x_{ij}$ are conditionally independent given $z$ and $w$ and follow the
same conditional distribution $\varphi$ which parameter $\alpha$ only depends on the block:

$$
x_{ij} | z_{ik} = 1, w_{j\ell} = 1 \sim \varphi(x_{ij}; \alpha_{k\ell})
$$

Observed likelihood

$$
\rho(x; \theta) = \sum_{(z,w) \in \mathcal{Z} \times \mathcal{W}} \underbrace{\prod_{i,k} \pi_k^{z_{ik}} \prod_{j,l} \rho_l^{w_{jl}}}_{\rho(z,w;\theta)=\rho(z,\theta)\rho(w;\theta)} \underbrace{\prod_{i,j,k,l} \varphi(x_{ij}; \alpha_{k\ell})^{z_{ik}w_{jl}}}_{\rho(x|z,w;\theta)}
$$

with $\theta = (\pi, \rho, \alpha) \in \Theta$

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

<footnote n="2">² | (Govaert and Nadif (2008))</footnote>
<figure><img src="image_107.png" /></figure>

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

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

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

Parametric model

► very parsimonious model, good candidate to deal with HD settings...

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

<table><thead><tr><th>Model</th><th>Number of parameters</th></tr></thead><tbody><tr><td>Binary</td><td>dim(π) + dim(ρ) + KL</td></tr><tr><td>Categorical</td><td>dim(π) + dim(ρ) + KL(m - 1)</td></tr><tr><td>Contingency</td><td>dim(π) + dim(ρ) + KL</td></tr><tr><td>Continuous</td><td>dim(π) + dim(ρ) + 2KL</td></tr></tbody></table>

$$ \text{nb. param}_{\text{HD}} = \text{nb. param}_{\text{classic}} \times \frac{L}{d} $$

<figure><img src="image_112.png" /></figure>
<figure><img src="image_113.png" /></figure>

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

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

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

Many appealing extensions

* variable type diversity: ordinal data [Biernacki and Jacques, '18], functional data [Bouveyron et al., '18], mixed-type data [Selosse et al., '20], textual interaction data [Bergé et al., '19]

* dynamic, e.g. [Marchello, 23] (zero inflated dLBM)

<figure><img src="image_117.png" /><figcaption>Figure 4.1: Graphical representation of the Zero-Inflated dLBM model.</figcaption></figure>

<figure><img src="image_118.png" /></figure>
<figure><img src="image_119.png" /></figure>

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

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

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

Particular case: graph clustering

<figure><img src="image_123.png" /><figcaption>► Stochastic Block Model (SBM): **x** is the adjacency matrix of a graph
$\leftrightarrow n = d, K = L, \mathbf{z} = \mathbf{w}$</figcaption></figure>

<figure><img src="image_124.png" /></figure>
<figure><img src="image_125.png" /></figure>

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

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

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

# Estimation

$$ \rho(x; \theta) = \sum_{(z, w) \in \mathcal{Z} \times \mathcal{W}} \prod_{i,k} \pi_k^{z_{ik}} \prod_{j,\ell} \rho_\ell^{w_{j\ell}} \prod_{i,j,k,\ell} \varphi(x_{ij}; \alpha_{k\ell})^{z_{ik}} w_{j\ell} $$

* generic identifiability [Keribin et al., '15]

* likelihood intractable: $2 \times 2$ blocks and $20 \times 20$ matrix: $\approx 2^{20} \times 2^{20} \approx 10^{12}$ terms

* E-step intractable

* $\leftrightarrow$ approximation with Variational EM (VEM) [Govaert and Nadif '08]

* $\leftrightarrow$ estimation with Stochastic EM (SEM) [Keribin et al., '10]

* $\leftrightarrow$ stochastic gradient descent [Baey et al., '23]

* $\leftrightarrow$ Bayesian versions: Gibbs and V-Bayes algorithms [Keribin, Brault, Celeux, Govaert, '15]

* consistency and asymptotic normality of ML- and V- estimators (when $\log(n)/d \to 0$ and $\log(d)/n \to 0$) [Brault, Keribin, Mariadassou, 20]

* model selection (ICL)

<figure><img src="image_129.png" /></figure>
<figure><img src="image_130.png" /></figure>

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

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

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

# EM for LBM

$$ \log p(\mathbf{x}; \boldsymbol{\theta}) = \underbrace{\mathbf{E}[\log p(\mathbf{x}, \mathbf{z}, \mathbf{w}; \boldsymbol{\theta})|\mathbf{x}; \boldsymbol{\theta}^{(r)}]}_{Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(r)})} - \underbrace{\mathbf{E}[\log p(\mathbf{z}, \mathbf{w}|\mathbf{x}; \boldsymbol{\theta})|\mathbf{x}; \boldsymbol{\theta}^{(r)}]}_{\boldsymbol{\theta}^{(r)}} $$

## Expectation Maximization (iteration $r$)

* **E-step:** computing $Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(r)})$

$$ \begin{aligned} Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(r)}) = & \sum_{i,k} \underbrace{\mathbf{P}(z_{ik} = 1|\mathbf{x}; \boldsymbol{\theta}^{(r)}) \log \pi_k}_{s_{ik}^{(r)}} + \sum_{j,\ell} \underbrace{\mathbf{P}(w_{j\ell}|\mathbf{x}; \boldsymbol{\theta}^{(r)}) \log \rho_\ell}_{t_{j\ell}^{(r)}} \\ & + \sum_{i,j,k,\ell} \underbrace{\mathbf{P}(z_{ik} w_{j\ell} = 1|\mathbf{x}; \boldsymbol{\theta}^{(r)}) \log \varphi(x_{ij}; \alpha_{k\ell})}_{e_{i,j,k,\ell}^{(r)} \text{ intractable}} \end{aligned} $$

* **M-step:** $\boldsymbol{\theta}^{(c+1)} = \arg\max_{\boldsymbol{\theta}} Q(\boldsymbol{\theta}|\boldsymbol{\theta}^{(r)})$, classical

<figure><img src="image_134.png" /></figure>
<figure><img src="image_135.png" /></figure>

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

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

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

VEM for LBM

$$ \log p(\mathbf{x}; \boldsymbol{\theta}) = \underbrace{\mathbf{E}_{q_z q_w} \left[ \log \frac{p(\mathbf{x}, \mathbf{z}, \mathbf{w}; \boldsymbol{\theta})}{q_z q_w} | \mathbf{x} \right]}_{\mathcal{F}(\boldsymbol{\theta}; q_z, q_w)} + KL(q_z q_w, \log p(\mathbf{z}, \mathbf{w} | \mathbf{x}; \boldsymbol{\theta})] $$

VEM (iteration $r$)

► approximated E-step:

$$ q_z^{(r)}, q_w^{(r)} = \arg \max_{q_z, q_w} \mathcal{F}(\boldsymbol{\theta}; q_z, q_w) $$

$$ \leftarrow s_{ik}^{(r)} \approx q_{zik}^{(r)}, t_{\ell}^{(r)} \approx q_{w\ell}^{(r)}, e_{ijk\ell}^{(r)} \approx q_{zik}^{(r)} q_{w\ell}^{(r)} $$

$$ Q(\boldsymbol{\theta} | \boldsymbol{\theta}^{(r)}) \approx \sum_{i,k} q_{zik}^{(r)} \log \pi_k + \sum_{j,\ell} q_{w\ell}^{(r)} \log \rho_{\ell} + \sum_{i,j,k,\ell} q_{zik}^{(r)} q_{w\ell}^{(r)} \log \varphi(x_{ij}; \alpha_{k\ell}) $$

► M-step: $\theta^{(c+1)} = \arg \max_{\theta} \mathcal{F}(\boldsymbol{\theta}; q_z^{(r)}, q_w^{(r)})$ in $\boldsymbol{\theta}$, classical

$$ \max_{\boldsymbol{\theta}} \log p(\mathbf{x}; \boldsymbol{\theta}) \geq \max_{\boldsymbol{\theta}, q_z, q_w} \mathcal{F}(\boldsymbol{\theta}; q_z, q_w) $$

<figure><img src="image_139.png" /></figure>
SEM for LBM

$$ \log p(\mathbf{x}; \boldsymbol{\theta}) = \underbrace{\mathbb{E}[\log p(\mathbf{x}, \mathbf{z}, \mathbf{w}; \boldsymbol{\theta}) | \mathbf{x}; \boldsymbol{\theta}^{(r)}]}_{Q(\boldsymbol{\theta} | \boldsymbol{\theta}^{(r)})} - \underbrace{\mathbb{E}[\log p(\mathbf{z}, \mathbf{w}|\mathbf{x}; \boldsymbol{\theta}) | \mathbf{x}; \boldsymbol{\theta}^{(r)}]}_{\boldsymbol{\theta}^{(r)}} $$

## SEM-Gibbs (iteration $r$)

* **SE-step:** impute missing values $(\mathbf{z}, \mathbf{w}) \sim p(\mathbf{z}, \mathbf{w}|\mathbf{x}; \boldsymbol{\theta}^{(r)})$ | Celeux and Diebolt (1986) with Gibbs-sampling

* $\mathbf{z}^{(t+1)} \sim p(\mathbf{z}|\mathbf{w}^{(t)}, \mathbf{x}; \boldsymbol{\theta}^{(r)})$

* $\mathbf{w}^{(t+1)} \sim p(\mathbf{w}|\mathbf{z}^{(t+1)}, \mathbf{x}; \boldsymbol{\theta}^{(r)})$

$\leftarrow$ $(\mathbf{z}^{(r)}, \mathbf{w}^{(r)})$, after several Gibbs iterations

* **M-step:** maximize $\log p(\mathbf{x}, \mathbf{z}^{(r)}, \mathbf{w}^{(r)}; \boldsymbol{\theta})$ in $\boldsymbol{\theta}$, classical

$\leftarrow$ stochastic, result fluctuates around the optimal value

<figure><img src="image_140.png" /></figure>
# Stochastic gradient descent for latent variables<fnref n="3" />

* Fisher Identity ([Cappé et al. '05]: $\nabla_{\theta} \log p(x;\theta) = E(\nabla_{\theta} \log p(x,\mathbf{z}|x;\theta))$)

## SGD-Gibbs for LBM (iteration $r$)

* **S-step:** Gibbs-sampling: $(\mathbf{z}^{(r)}, \mathbf{w}^{(r)}) \sim p(\cdot, |\mathbf{x}; \theta^{(r)})$

* gradient for observations $i = 1, \dots, n$; $j = 1, \dots, d$

$$J_{ij}^{(r), \text{obs}} = \nabla_{\theta} \log p(x_{ij} | z_j^{(r)}, w_j^{(r)}, \theta^{(r)}); \Delta_{ij}^{(r), \text{obs}} = (1-\gamma^{(r)}) \Delta_{ij}^{(r), \text{obs}} + \gamma^{(r)} J_{ij}^{(r), \text{obs}}$$

* gradient for latent variables $i = 1, \dots, n$; $j = 1, \dots, d$

$$J_{i}^{(r),z} = \nabla_{\theta} \log p(z_{i}^{(r)}, \theta^{(r)}); \Delta_{i}^{(r),z} = (1-\gamma^{(r)}) \Delta_{i}^{(c-1)}, z_{+\gamma^{(r)}} J_{i}^{(r),z}$$

$$J_{j}^{(r),w} = \nabla_{\theta} \log p(w_{j}^{(r)}, \theta^{(r)}); \Delta_{j}^{(r),w} = (1-\gamma^{(r)}) \Delta_{j}^{(r-1)}, z_{+\gamma^{(r)}} J_{j}^{(r),w}$$

* $\theta^{(r+1)} = \theta^{(r)} + \gamma^{(k)} P_k v_k$

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

<footnote n="3">³[Baey, Delatre, Kuhn, Leger, Lemler '23]</footnote>
# Outline

1 Motivation

2 Latent Block Model

3 Estimation challenges

4 Co-clustering for HD

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

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

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

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

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

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

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

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

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

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

<figure><img src="image_152.png" /></figure>
# Estimation challenges

## Challenges in the (classical) MBC context:

* local maxima solutions

* empty cluster solutions

* degenerate solutions

**Question to be addressed here**

What about LBM?

<figure><img src="image_153.png" /></figure>
<figure><img src="image_154.png" /></figure>

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

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

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

Local maxima

► not just linked to the number of parameters!

► linked to the number of latent variables (here a toy two-param. Gaussian mixture)

40% of the partition is unknown

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

85% of the partition is unknown

70% of the partition is unknown

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

100% of the partition is unknown

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

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

► additional order of magnitude: from $K^n$ to $K^n L^d \leftarrow$ pb for LBM?
but no way to draw a likelihood map...

<figure><img src="image_162.png" /></figure>
Empty clusters in LBM

Really frequent (and reported)!

► a component such that $\pi_K^{(r)} \approx 0$

► number of empty clusters in MBC: $\#Z^0$

► number of empty clusters in LBM: $\#(Z \times W)^0$

$$ \frac{\#(Z \times W)^0}{\#Z^0} = \frac{K^n L^d - S(n, K)S(d, L)K!L!}{K^n - S(n, K)K!} \rightarrow \infty \text{ when } n \rightarrow \infty \text{ or } d \rightarrow \infty $$

$S(n, K)$ is the Stirling number of second kind

► exponential speed even with very small sample sizes

$$ \frac{\#(Z \times W)^0}{\#Z^0} = \begin{cases} 62 & \text{when } n = d = 5 \text{ and } K = L = 2 \\ 710\,768 & \text{when } n = d = 9 \text{ and } K = L = 4 \end{cases} $$

↔ initialization often achieved by performing independent clustering of the rows and columns but the best initialization strategy stays an open question

<figure><img src="image_163.png" /></figure>
<figure><img src="image_164.png" /></figure>

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

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

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

Degenerate and spurious solutions in MBC

<figure><img src="image_168.png" /><figcaption>► spurious local maximizers due to one or several relatively small (but not null) generalized variance</figcaption></figure>

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

<figure><img src="image_170.png" /></figure>
# Degenerate and spurious solutions in LBM

## No reported such problems?

* freq. of obtaining a cluster with a single element in the whole partition latent space

$$ \frac{\#Z_1^1}{\#Z} = \left(\frac{K-1}{K}\right)^{n-1} $$

* freq. of obtaining a block with a single element in the whole partition latent space

$$ \frac{\#(Z \times W)^{1,1}}{\#(Z \times W)} = \frac{\#Z_1^1}{\#Z} \left(\frac{L-1}{L}\right)^{d-1} $$

→ degeneracy situations are expected to be significantly less present in co-clustering than in clustering

Ex with $L = 4$ and $d = 50$: $(1 - 1/L)^{d-1} = 7.5510^{-7}$

* spurious case: same combinatorial arguments (but more tedious computations)

<figure><img src="image_171.png" /></figure>
Local maxima in Gaussian LBM: illustration<fnref n="4" />

Figure: Gaussian 2 × 2 LBM (**left**) ; e.c.d.f. of ELBO values obtained from *B* = 100 initializations, at both standard precision (**center**) and high precision (**right**)

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

► some potential slow convergence effect

► different parameters but the same partition

→ clustering is easier than estimation, no need to estimate with too much precision

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

<footnote n="4">⁴R package blockcluster</footnote>
<figure><img src="image_174.png" /></figure>

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

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

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

Local maxima in Gaussian LBM: illustration (cont'd)

<figure><img src="image_178.png" /><figcaption>Projection of the solutions on the PCA planes</figcaption></figure>

► *blockcluster* initialization performs several initial small VEM steps, that can explain these results

<figure><img src="image_179.png" /></figure>
<figure><img src="image_180.png" /></figure>

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

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

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

Local maxima in binary LBM: illustration<fnref n="5" />

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

<table><tr><td>ELBO</td><td>-1615.48</td><td>-1615.43</td><td>-1596.67</td></tr><tr><td>count</td><td>9</td><td>38</td><td>254</td></tr><tr><td>nb steps</td><td>&gt; 10<sup>4</sup></td><td>&gt; 10<sup>4</sup></td><td>&lt; 100</td></tr></table>

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

* high ELBO values (red): 36 relabelled configurations

* low ELBO values (black): 30 co-clustering configurations, all with an empty row or cluster

↔ very lazy convergence

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

<footnote n="5">⁵R package bikm1</footnote>
<figure><img src="image_187.png" /></figure>

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

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

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

Initialization strategies: crucial!

► two independent clustering (row / column)

► using several initializations with other algorithms (CEM, SEM, EM-VBayes) on few iterations

► resampling in order to avoid getting empty blocks during the first iterations

► *blockmodels* [Leger, '16] smart and robust initialization, combining *spectral clustering* adapted for LBM and a *forward/backward* reinitialization strategy:

- forward exploration of the space of models (K, L) by splitting already existing clusters [Robert '17]

- backward exploration by merging groups

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

<figure><img src="image_192.png" /></figure>
# Outline

1 Motivation

2 Latent Block Model

3 Estimation challenges

4 Co-clustering for HD

<figure><img src="image_193.png" /></figure>
## Properties of LBM in view of HD clustering

* parsimonious model

* the true partition is recovered when there is a consistent estimator [Mariadassou and Matias, '15]

$$ \hat{\theta}^{n,d \to \infty} \theta^* \Rightarrow p(\hat{\mathbf{z}} = \mathbf{z}^*, \hat{\mathbf{w}} = \mathbf{w}^* | \mathbf{x}; \hat{\theta})^{n,d \to \infty} 1, $$

* non asymptotic properties (binary case) [Braault, '14]

$$ \mathbb{P}(\hat{\mathbf{z}} \neq \mathbf{z}^*) \leq 2n \exp \left\{ -\frac{1}{8} d \left[ \underbrace{\min_{k \neq k'} |\tau_k - \tau_{k'}|}_{\text{overlap}} \right] \right\} + K(1 - \min_k \pi_k)^n $$

→ LBM is a natural regularized candidate for HD clustering

<figure><img src="image_194.png" /></figure>
Document clustering: exact recovery

► Mixture of 1033 medical summaries and 1398 aeronautics summaries

► Lines: 2431 documents

► Columns: present words (except stop), thus 9275 unique words

► Data matrix: cross counting document×words

► Poisson model

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

<figure><img src="image_196.png" /></figure>
<figure><img src="image_197.png" /></figure>

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

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

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

Document clustering: exact recovery (cont'd)

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

Results with 2×2 blocks

<table><thead><tr><th></th><th>Medline</th><th>Cranfield</th></tr></thead><tbody><tr><th>Medline</th><td>1033</td><td>0</td></tr><tr><th>Cranfield</th><td>0</td><td>1398</td></tr></tbody></table>

Experiment illustrates previous theory: HD clustering is blessing

<figure><img src="image_202.png" /></figure>
Numerical illustrations: experimental design

Generating models, with two balanced ($\pi_1 = \pi_2$) row clusters in $\mathbb{R}^d$

<table><thead><tr><th></th><th>$\mathbf{I}_d$</th><th>$\Sigma(c)$</th></tr></thead><tbody><tr><td>$\mu_1 = \mathbf{0}_d, \mu_2 = \mathbf{1}_d$</td><td>$(M_1)$</td><td>$(M_2)$</td></tr><tr><td>$\mu_1 = \mathbf{0}_d, \mu_2 = (1, 2^{-2}, \dots, d^{-2})$</td><td>$(M_3)$</td><td>$(M_6)$</td></tr><tr><td>$(M_1)$ of size $n \times d/2$ duplicated twice</td><td>$(M_4)$</td><td></td></tr><tr><td>$\mu_1 \sim \mathcal{N}_d(\mathbf{0}, \mathbf{I}_d), \mu_2 \sim \mathcal{N}_d((1, \dots, \underbrace{1, \dots, 0, \dots, 0}_{\sqrt{d}}), \mathbf{I}_d)$</td><td>$(M_5)$</td><td></td></tr></tbody></table>

↔ only $(M_1)$ is a nominal LBM (2x1)

Row clustering is performed using four different methods

* clustering with a mixture of two spherical *d*-dim. Gaussian distributions

* clustering with a mixture of two full-covariance *d*-dim. Gaussian distr.

* co-clustering with a Gaussian ($2 \times 1$) LBM

* co-clustering with a Gaussian ($2 \times 2$) LBM

Classification error averaged over 30 samples of size $n = 30$

<figure><img src="image_203.png" /></figure>
<figure><img src="image_204.png" /></figure>

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

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

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

Numerical illustrations: results

<figure><img src="image_208.png" /><figcaption>Model (M1)</figcaption></figure>

LBM resists to HD

<figure><img src="image_209.png" /><figcaption>Model (M2)</figcaption></figure>

LBM performs better than simple
mixture, although the underlying
model overrides the co-clustering
conditional independence
assumptions

<figure><img src="image_210.png" /></figure>
<figure><img src="image_211.png" /></figure>

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

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

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

## Numerical illustrations: results (cont'd)

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

non informative variables

redundancy: duplicated variables

Same conclusion for LBM, which resists to such HD scenarios

<figure><img src="image_216.png" /></figure>
<figure><img src="image_217.png" /></figure>

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

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

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

Numerical illustrations: results (cont'd)

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

too many noninformative variables

* too much bias, classification error with CC no longer lines up with Bayes risk

* but CC methods still perform better than simple mixture

<figure><img src="image_222.png" /></figure>
<figure><img src="image_223.png" /></figure>

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

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

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

# Take home message

* in HD clustering, accept model bias to reduce variance on the classification

* CC acts as a **regularized tool** to perform clustering
  → very parsimonious model, with model bias but often offers better classification error

* LBM also offers a level of **flexibility** despite its parsimony
  → using LBM for HD clustering should be more emphasized
  → estimation is challenging, but there are solutions
  → less cases of spurious / degenerate solutions than in simple mixture models

## For further research...

* interpret the groups of variables in co-clustering

* extend the empirical results by theoretical guidelines

* robust estimation

* model selection: to be fixed by theoretical results

<figure><img src="image_227.png" /></figure>
# A Survey on Model-Based Co-Clustering: High Dimension and Estimation Challenges

C. Biernacki¹ · J. Jacques² · C. Keribin³

Accepted: 9 May 2023

© The Author(s) under exclusive licence to The Classification Society 2023

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