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

Christine Keribin, Christophe Biernacki, Julien Jacques

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

## To cite this version:

Christine Keribin, Christophe Biernacki, Julien Jacques. Model Based Co-Clustering: High Dimension and Estimation Challenges. 2024. hal-04862826

HAL Id: hal-04862826

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

Submitted on 3 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

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

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

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

# MODEL BASED CO-CLUSTERING: HIGH DIMENSION AND ESTIMATION CHALLENGES

Christine KERIBIN<fnref n="1" />

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

## Séminaire parisien de statistique, IHP, 11 mars 2024

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

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

<footnote n="1">¹Université Paris-Saclay, CNRS, Inria, Laboratoire de mathématiques d'Orsay, 91405, Orsay, France. - Equipe
CELESTE</footnote>

<footnote n="2">²Université de Lille
INRIA - Lille - MODAL</footnote>

<footnote n="3">³Université de Lyon, Lyon 2, ERIC UR 3083</footnote>
<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>

## Outline

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

### 1 Motivation

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

#### 2 Latent Block Model

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

#### 3 Estimation challenges

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

#### 4 CC for HD

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

## Unsupervised ML framework

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

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

* classification matrix $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

* $x_i$ are iid observations of a mixture pdf

$$f(.; \theta) = \sum_{k=1}^{K} \pi_k \varphi(.; \alpha_k) ; \sum_{i=1}^{K} \pi_K = 1, \pi_K \ge 0; \mathbb{P}(z_{ik} = k) = \pi_k$$

* estimation with EM algorithm, partition is deduced from the conditional probability $p(z_i|x_i; \hat{\theta})$

<figure><img src="image_22.png" /></figure>
<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>

HD settings: *d* ~ *n*, *d* >> *n*

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

* marketing: *d* ~ 10² for *n* = 10⁴ – 10⁶

* microarray gene expression: *d* ~ 10² – 10⁴ for *n* = 20 – 50 stress conditions

* fMRI images: *d* ~ 10⁴ – 10⁵ voxels for *n* = 50 – 100 patients

* textmining, curve, internet logs ...

Variables with mixed types

* different modes: qualitative, quantitative, functional

* different semantic

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

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

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

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

# Clustering in HD: curse (and blessing?)

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

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

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

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

Basic case: $\mu_1 = 0d$, $\mu_2 = 1d$, $\Sigma = ld$

► **blessing:** but components could be more and more separated.

► Which impact on the misclassification error ?

- what if the variables are *correlated, redundant, non informative* ?

- what if $d >> n$ ?

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

Basic case: HD is a blessing for clustering

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

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 nb_ob= 30

<figure><img src="image_44.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_45.png" /></figure>
<figure><img src="image_46.png" /></figure>

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

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

<figure><img src="image_49.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_50.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$$

less informative: $\Sigma = I_d$; $\mu_1 = 0_d$

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

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

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

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

<figure><img src="image_55.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 nb_obss = 30

<figure><img src="image_56.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+ correlated variables: $\Sigma = \sum_{d}(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, too much variance**

* **error decreases with the biased model which does not dammaged the separation**

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

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

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

► Dimension reduction in the canonical space (variable selection)
for Gaussian mixture, 3 variable-clusters (informative, non informative, linearly dependent) [Maugis et al (2009)], for latent class analysis [Fop et al (2018)]

► Model parsimony in the initial HD space (constraints on model parameters)
GMM [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 perfect modeling and what can be estimated in practice

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

# 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

* 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_67.png" /></figure>
<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>

[Govaert (2008)]

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

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

► PCA for the i-th data individual $\mathbf{x}_i$ :

$$ \mathbf{x}_i = \sum_{j=1}^{d} x_i^j \mathbf{e}^j = \sum_{j=1}^{d} a_i^j \mathbf{u}^j \approx \sum_{j=1}^{J} a_i^j \mathbf{u}^j. $$

► specific CC case reduced to K = 1

$$ \mathbf{x}_i = \sum_{j=1}^{d} x_i^j \mathbf{e}^j = \sum_{j=1}^{d} (\mu_{\bar{w}_j} + \varepsilon_i^j) \mathbf{e}^j = \sum_{\ell=1}^{L} \mu_\ell \mathbf{v}^\ell + \mathbf{r}_i \approx \sum_{\ell=1}^{L} \mu_\ell \mathbf{v}^\ell. $$

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

► From sequential method PCA then MBC to simultaneous clustering with the dimension reduction process

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

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

<footnote n="1">¹[C. Biernacki, C. Keribin, J. Jacques, A survey on Model-Based Co-Clustering, JOC 2023]</footnote>
<figure><img src="image_80.png" /></figure>

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

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

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

# Outline

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

1 Motivation

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

2 Latent Block Model

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

3 Estimation challenges

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

4 CC for HD

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

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

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

* **Co-clustering:** 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$

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

* **In our case:** dimension reduction for HD clustering

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

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

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

$H_3$: the $n \times d$ cells $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_{jk} = 1, w_{j\ell} = 1 \sim \varphi(x_{ij}; \alpha_{k\ell})$$

Observed likelihood

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

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

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

Parametric model

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

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

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

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

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

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

# Many appealing extensions

* variable type diversity: ordinal data [Bieracki 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_110.png" /><figcaption>Figure 4.1: Graphical representation of the Zero-Inflated dLBM model.</figcaption></figure>

<figure><img src="image_111.png" /></figure>
<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>

Particular case: graph clustering

► **Stochastic Block Model (SBM):** $x$ is the adjacency matrix of a graph
  $\Leftrightarrow n = d, K = L, \mathbf{z} = \mathbf{w}$

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

<figure><img src="image_117.png" /></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>

# Estimation

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

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

* likelihood for 2 × 2 blocks and 20 × 20 matrix: ≈ 2²⁰ × 2²⁰ ≈ 10¹² terms: intractable

* E-step: intractable

* ↪→ Approximation with Variational EM [Govaert and Nadif '08]

* ↪→ Estimation with SEM [K. et al '10]

* ↪→ Stochastic gradient descent [Baey et al. '23]

* ↪→ Bayesian versions: Gibbs and V-Bayes algorithms [K., Brault, Celeux, Govaert, '15]

* Consistency and asymptotic normality of ML- and V- estimators (when $log(n)/d \to +\infty$ and $log(d)/n \to +\infty$) [Brault, K., Mariadassou, '20]

* model selection (ICL)

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

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

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

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

# EM for LBM

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

## Expectation Maximization (iteration c)

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

$$ Q(\theta|\theta^{(c)}) = \sum_{i,k} \underbrace{P(z_{ik} = 1|\theta^{(c)}, \mathbf{x}) \log \pi_k}_{s_{ik}^{(c)}} + \sum_{j,\ell} \underbrace{P(w_{j\ell} = 1|\theta^{(c)}, \mathbf{x}) \log \rho_\ell}_{t_{j\ell}^{(c)}} + \sum_{i,j,k,\ell} \underbrace{P(z_{ik} w_{j\ell} = 1|\theta^{(c)}, \mathbf{x}) \log \varphi(x_{ij}; \alpha_k \ell)}_{e_{i,j,k,\ell}^{(c)}} $$

$\Leftrightarrow s_{ik}^{(c)}, t_{j\ell}^{(c)}, e_{i,j,k,\ell}^{(c)}$ are intractable

* **M-step:** maximizing $Q(\theta|\theta^{(c)})$ in $\theta$, classical

$\Leftrightarrow \theta^{(c+1)} = \arg\max_{\theta} Q(\theta|\theta^{(c)})$

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

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

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

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

# VEM for LBM

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

## VEM (iteration c)

### ► approximated E-step:

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

$$ \Leftrightarrow s_{ik}^{(c)} \approx q_{zik}^{(c)}, t_{j\ell}^{(c)} \approx q_{wj\ell}^{(c)}, e_{ijk\ell}^{(c)} \approx q_{zik}^{(c)} q_{wj\ell}^{(c)} $$

$$ Q(\theta|\theta^{(c)}) \approx \sum_{i,k} q_{zik}^{(c)} \log \pi_k + \sum_{j,\ell} q_{wj\ell}^{(c)} \log \rho_\ell + \sum_{i,j,k,\ell} q_{zik}^{(c)} q_{wj\ell}^{(c)} \log \varphi(x_{jj}; \alpha_k) $$

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

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

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

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

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

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

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

SEM-Gibbs (iteration c)

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

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

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

$$ \leftarrow (\mathbf{z}^{(c)}, \mathbf{w}^{(c)}) $$

► step M: maximize $\log p(\mathbf{x}, \mathbf{z}^{(c)}, \mathbf{w}^{(c)}; \theta)$ in $\theta$

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

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

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

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

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

# Stochastic gradient descent for latent variables <fnref n="2" />

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

## SGD-Gibbs for LBM (iteration c)

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

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

$$J_{ij}^{(c), \text{obs}} = \nabla_{\theta} \log p(x_{ij} | z_{i}^{(c)}, w_{j}^{(c)}, \theta^{(c)}); \Delta_{ij}^{(c), \text{obs}} = (1 - \gamma^{(c)}) \Delta_{ij}^{(c), \text{obs}} + \gamma^{(c)} J_{ij}^{(c), \text{obs}}$$

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

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

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

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

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

<footnote n="2">²[Baey, Delattre, Kuhn, Leger, Lemler '23</footnote>
<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>

# Outline

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

1 Motivation

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

2 Latent Block Model

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

3 Estimation challenges

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

4 CC for HD

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

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

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

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

# Estimation challenges

in model based clustering context:

* "classical" local maxima

* empty cluster solutions

* degenerate solutions

* worse for LBM ?

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

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

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

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

► linked to the number of latent variables

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

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

<figure><img src="image_162.png" /></figure>
<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>

Empty clusters in LBM

Really frequent (and reported) !

► a component such that $\pi_k^{(c)} \simeq 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

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

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

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

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

# Degenerate and spurious solutions in MBC

► Degenerate solutions lay on the border parameter space

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

► spurious local maximizers due to one or several relatively small (but not null) generalized variance

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

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

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

# Degenerate and spurious solutions in LBM

no reported such problems ?

* frequency 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}. $$

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

$$ \frac{\#(Z \times W)^1_{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

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

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

<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 Gaussian LBM: illustration<fnref n="3" />

Figure: *Left:* Gaussian $2 \times 2$ LBM; e.c.d.f. of ELBO values obtained from $B = 100$ initializations, standard precision (*center*) and high precision (*right*)

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

* some potential slow convergence effect

* different parameters but the same partition
  $\leftrightarrow$ clustering is easier than estimation, no need to estimate with too much precision

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

<footnote n="3">³R package blockcluster</footnote>
<figure><img src="image_186.png" /></figure>

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

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

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

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

<figure><img src="image_190.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_191.png" /></figure>
<figure><img src="image_192.png" /></figure>

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

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

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

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

<figure><img src="image_196.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>> 10<sup>4</sup></td><td>> 10<sup>4</sup></td><td>> 100</td></tr></table>

<figure><img src="image_197.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_198.png" /></figure>

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

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

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

<figure><img src="image_202.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_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>

## Outline

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

1. Motivation

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

2. Latent Block Model

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

3. Estimation challenges

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

4. CC for HD

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

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

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

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

# Properties of LBM in view of HD clustering

* parcimonious model

* the true partition is recovered when there is a consistent estimator

[Mariadassou and Matias, '15]

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

↔ naturally regularized candidate for HD clustering

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

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

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

<table><thead><tr><th></th><th><i>I</i><sub>d</sub></th><th><i>Σ</i><sub>d</sub>(<i>c</i>)</th></tr></thead><tbody><tr><td>$\mu_1 = \mathbf{0}_d, \mu_2 = \mathbf{1}_d$</td><td>(<i>M</i><sub>1</sub>)</td><td>(<i>M</i><sub>2</sub>)</td></tr><tr><td>$\mu_1 = \mathbf{0}_d, \mu_2 = (1, 2^{-2}, \dots, d^{-2})$</td><td>(<i>M</i><sub>3</sub>)</td><td>(<i>M</i><sub>6</sub>)</td></tr><tr><td>($M_1$) of size $n \times d/2$ duplicated twice</td><td>(<i>M</i><sub>4</sub>)</td><td></td></tr><tr><td>$\mu_1 \sim \mathcal{N}_d(\mathbf{0}, \mathbf{I}_d), \mu_2 \sim \mathcal{N}_d((\underbrace{1, \dots, 1}_{\sqrt{d}}, 0, \dots, 0), \mathbf{I}_d)$</td><td>(<i>M</i><sub>5</sub>)</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_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>

Numerical illustration

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

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

Resists to HD

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

<figure><img src="image_229.png" /></figure>
<figure><img src="image_230.png" /></figure>

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

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

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

Numerical illustration

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

non info. variables

$$ \mu_1 \sim \mathcal{N}_d(0, ld) $$

$$ \mu_2 \sim \mathcal{N}_d((1, \dots, 1, 0, \dots, 0), ld) \qquad \sqrt{d} $$

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

redundancy: duplicated variables

<figure><img src="image_236.png" /></figure>
<figure><img src="image_237.png" /></figure>

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

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

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

Numerical illustration: too many noninformative variables

$$ \mu_1 = \mathbf{0}_d, \mu_2 = (1, 2^{-2}, \dots, d^{-2}) $$

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

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

* 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_243.png" /></figure>
<figure><img src="image_244.png" /></figure>

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

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

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

# Take home message

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

* co-clustering 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_248.png" /></figure>
<figure><img src="image_249.png" /></figure>

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

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

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