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

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

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

Clustering : une vision unifiée pour une utilisation éclairée

—

Partie 4

Traitement de la grande dimension & co-clustering
(leçon plus avancée)

C. Biernacki

Ateliers statistiques de la SFdS, 15 et 16 juin 2023, IHP, Paris
<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>

## Outline

1. HD clustering

2. Modeling

3. Estimating

4. Selecting

5. Experiments

6. To go further
<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>

# Motivation

High dimensional (HD) data sets are now frequent:

* Marketing: $d \sim 10^2$

* microarray gene expression: $d \sim 10^2-10^4$

* SNP data: $d \sim 10^6$

* Curves: depends on discretization but can be very high

* Text mining

* ...

Clustering has to be applied for HD datasets for the same reasons as the lower dimensional datasets:

* Data summary

* Data exploratory

* Preprocessing for more flexibility of a forthcoming prediction step

But clustering is even more important since visualization in the HD setting can be hazardous...
<figure><img src="image_19.png" /></figure>

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

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

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

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

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

Today': exponential growing of dimension<fnref n="1" />

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

<footnote n="1">¹S. Alelyani, J. Tang and H. Liu (2013). Feature Selection for Clustering: A Review. *Data Clustering: Algorithms and Applications*, 29</footnote>
<figure><img src="image_26.png" /></figure>

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

# HD data: definition (1/2)

## An attempt in the non-parametric case

Dataset $\mathbf{x} = (\mathbf{x}_1, \dots, \mathbf{x}_n)$, $\mathbf{x}_j$ described by $d$ variables, where $n = o(e^d)$

## Justifications:

* To approximate within error $\epsilon$ a (Lipschitz) function of $d$ variables, about $(1/\epsilon)^d$ evaluations on a grid are required [Bellman, 61]

* Approximate a Gaussian distribution with fixed Gaussian kernels and with approximate error of about 10% [Silverman, 86]

$$ \log_{10} n(d) \approx 0.6(d - 0.25) $$

For instance, $n(10) \approx 7.10^5$
<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>

HD data: definition (2/2)

An attempt in the parametric case

Dataset **x** = (x₁, ..., xₙ), **x**ᵢ described by **d** variables and a model **m** with **v**
parameters, where **n** = o(g(**v**)), with **g** a given function

Justification:

■ We consider the heteroscedastic Gaussian mixture with of true parameter θ* with K* components. We note θ̂ the Gaussian MLE with K* components. We have g linear from the following result [Michel, 08]: it exists constants κ, A and C such that

$$
\mathrm{Ex}[\mathrm{Hellinger}^2(\mathrm{p}_{\theta^*}, \mathrm{p}_{\hat{\theta}_K})] \le C \left[ \kappa \frac{\nu}{n} \left\{ 2A \ln d + 1 - \ln \left( 1 \wedge \left[ \frac{\nu}{n} A \ln d \right] \right) \right\} + \frac{1}{n} \right].
$$

But $\nu$ can be high since $\nu \sim d^2/2$, combined with potentially large constants.
<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>

HD density estimation: curse

A two-component d-variate Gaussian mixture:

$$ \pi_1 = \pi_2 = \frac{1}{2}, \quad \mathbf{X}_1|z_{11} = 1 \sim N_d(0, I), \quad \mathbf{X}_1|z_{12} = 1 \sim N_d(1, I) $$

Components are more and more separated when d grows: $||\mu_2 - \mu_1||_1 = \sqrt{d}...$

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

... but density estimation quality decreases with d
<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>

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

HD clustering: blessing (1/2)

A two-component d-variate Gaussian mixture:

$$ \pi_1 = \pi_2 = \frac{1}{2}, \quad \mathbf{X}_1|z_{11} = 1 \sim N_d(\mathbf{0}, \mathbf{I}), \quad \mathbf{X}_1|z_{12} = 1 \sim N_d(\mathbf{1}, \mathbf{I}) $$

Each variable provides equal and own separation information

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

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

... and empirical error rate decreases also with d!
<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>

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

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

HD clustering: blessing (2/2)

FDA

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

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

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

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

HD clustering: curse (1/2)

Many variables provide no separation information

Same parameter setting except:

$$X_1|z_{12} = 1 \sim N_d((1 \ 0 \ \dots \ 0)', I)$$

Groups are not separated more when d grows: $||\mu_2 - \mu_1||_1 = 1...$

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

...thus theoretical error is constant (= $\Phi(-\frac{1}{2})$) and empirical error increases with d
<figure><img src="image_66.png" /></figure>

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

HD clustering: curse (2/2)

Many variables provide redundant separation information

Same parameter setting except:

$$ \mathbf{X}_1^j = \mathbf{X}_1^1 + N_1(0, 1) \quad (j = 2, \dots, d) $$

Groups are not separated more when d grows: $||\mu_2 - \mu_1||_\Sigma = 1...$

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

...thus err<sub>theo</sub> is constant (= $\Phi(-\frac{1}{2})$) and empirical error increases (less) with d
<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>

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

# The trade-off bias/variance

## The fundamental statistical principle

Always minimize an error err between truth (z) and estimate (ŷ)

* Gap between true (z) and model-based ($\mathcal{Z}_p$) partitions: $z^* = \arg\min_{\hat{z} \in \mathcal{Z}_p} \Delta(z, \hat{z})$

* Estimation ŷ of z* in $\mathcal{Z}_p$: any relevant method (bias, consistency, efficiency...)

* Fundamental decomposition of the observed error err(z, ŷ):

$$ 
\begin{align*}
\text{err}(z, \hat{z}) &= \{\text{err}(z, z^*) - \text{err}(z, z)\} + \{\text{err}(z, \hat{z}) - \text{err}(z, z^*)\} \\
&= \{\text{bias}\} + \{\text{variance}\} \\
&= \{\text{error of approximation}\} + \{\text{error of estimation}\}
\end{align*}
$$
<figure><img src="image_79.png" /></figure>

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

Bias/variance in HD: reduce variance, accept bias

A two-component *d*-variate Gaussian mixture with intra-dependency:

$$ \pi_1 = \pi_2 = \frac{1}{2}, \quad \mathbf{X}_1|z_{11} = 1 \sim N_d(0, \Sigma), \quad \mathbf{X}_1|z_{12} = 1 \sim N_d(1, \Sigma) $$

* Each variable provides equal and own separation information

* Theoretical error decreases when *d* grows: $err_{theo} = \Phi(-||\mu_2 - \mu_1||_{\Sigma}^{-1/2})$

* Empirical error rate with the (true) intra-correlated model worse with *d*

* Empirical error rate with the (false) intra-independent model better with *d*!

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

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

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

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

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

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

## Some alternatives for reducing variance

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

* Dimension reduction in the canonical space (variable selection)

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

### But which kind of parsimony?

* Remember that clustering is a way for dealing with large *n*

* Why not reusing this idea for large *d*?

### Co-clustering

It performs parsimony of row clustering through variable clustering
<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>

# Gaussian mixture of factor analysers (1/2)

## Definition

[Ghahramani and Hinton, 97], [McLachlan et al., 03]

$$ \Sigma_k = B_k B'_k + \omega_k \Lambda_k $$

where

* $B_k$ is a loadings $d \times q$ non-square real matrix ($1 \le q \le q_{max}$, $q_{max} < d$)

* $\omega_k$ is a positive real number

* $\Lambda_k$ is a $d \times d$ diagonal positive definite matrix such that $|\Lambda_k| = 1$

Interpretation $X_1 \in \mathbb{R}^d$ is generated by a latent variable $Y_1 \in \mathbb{R}^q$

$$ X_{1|Y_1, Z_{1k}=1} = B_k \underbrace{Y_1}_{\text{factor}} + \mu_k + \varepsilon_k $$

where $Y_1 \perp \varepsilon_k$ and

$Y_1 \sim N_q(0, I)$ and $\varepsilon_k \sim N_d(0, \omega_k \Lambda_k)$
<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>

Gaussian mixture of factor analysers (2/2)

Complexity

$$ \nu = (K - 1) + Kd + Kq(d - (q - 1)/2) + Kd $$

Model selection

- Models in competition: $m = (q, K) + 12$ parsimonious versions [McNicholas and Murphy, 10]

- Classical criteria can be used

Package

PGMM: http://cran.r-project.org/web/packages/pgmm/index.html
<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>

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

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

HD Gaussian models (1/2)

Definition

[Bouveyron et al., 07]

$$ \Sigma_k = D_k \Delta_k D'_k $$

where

- $D_k$ is the orthogonal matrix of the eigenvectors of $\Sigma_k$

- $\Delta_k$ is a diagonal matrix containing the eigenvalues of $\Sigma_k$

$$ \Delta_k = \left( \begin{array}{cc} a_{k1} & 0 \\ \vdots & \vdots \\ 0 & a_k \delta_k \end{array} \right) \left. \begin{array}{c} 0 \\ b_k \\ \vdots \\ 0 \end{array} \right\} \delta_k \quad (d - \delta_k) $$

with $a_{kj} \ge b_k$, for $j = 1, \dots, \delta_k$ and $\delta_k < d$

Complexity

$$ \nu = (K - 1) + Kd + \sum_{k=1}^{K} \delta_k \left[ d - \left( \delta_k + 1 \right) / 2 \right] + \sum_{k=1}^{K} \delta_k + 2K $$
<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>

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

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

HD Gaussian models (2/2)

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

Model selection

■ Models in competition: $\mathbf{m} = ((\delta_1, \dots, \delta_K), K) + 8$ parsimonious versions
[Bouveyron et al., 07]

■ Since it could be **combinatorics**, detect **empirically** a break in the eigenvalue scree of the empirical covariance matrix of each group

■ Classical criteria can be used for the remaining models

Package

R MIXMOD: http://cran.r-project.org/web/packages/Rmixmod/index.html
<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>

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

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

Gaussian “variable selection” (1/2)

Definition

[Raftery and Dean, 06], [Maugis et al., 09a], [Maugis et al., 09b]

$$ p(\mathbf{x}_1; \theta) = \underbrace{\left\{ \sum_{k=1}^{K} \pi_k p(\mathbf{x}_1^S; \mu_k, \Sigma_k) \right\}}_{\text{clustering variables}} \times \underbrace{\left\{ p(\mathbf{x}_1^U; \mathbf{a} + \mathbf{x}_1^R \mathbf{b}, \mathbf{C}) \right\}}_{\text{redundant variables}} \times \underbrace{\left\{ p(\mathbf{x}_1^W; \mathbf{u}, \mathbf{V}) \right\}}_{\text{independent variables}} $$

where

■ all parts are Gaussians

■ S: set of variables useful for clustering

■ U: set of redundant clustering variables, expressed with $R \subseteq S$

■ W: set of variables independent of clustering

Trick

Variable selection is recasted as a particular variable role
<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>

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

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

Gaussian “variable selection” (2/2)

Model selection

■ Models in competition: **m** = (S, R, U, W, K) → combinatorics

■ Use a backward stepwise algorithm guided by a model selection criterion: *d* ≈ 30

■ Use alternatively a lasso-like procedure for ranking quickly different sets of clustering related and clustering independent variables [Sedki et al., 14]

$$
\operatorname{crit}_{\lambda, \rho} = \ell(\boldsymbol{\theta}; \bar{\boldsymbol{x}}) - \lambda \sum_{k=1}^{K} \sum_{j=1}^{d} |\mu_{kj}| - \rho \sum_{k=1}^{K} \sum_{(j,j') \neq j'j}^{d} |(\Sigma_k^{-1})_{jj'}|
$$

where θ full Gaussian parameters, x̄ is x centered and (λ, ρ) are on a grid
A variable j is considered independent of clustering if ſ̂_{kj}(λ, ρ) = 0 for all k

■ Classical criteria are available

Package

■ SELVARCLUST:
http://www.math.univ-toulouse.fr/~maugis/SelvarClustHomepage.html

■ SELVARMIX:
http://cran.r-project.org/web/packages/SelvarMix/index.html
<figure><img src="image_125.png" /></figure>

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

From clustering to co-clustering

<figure><img src="image_127.png" /><figcaption>(1)</figcaption></figure>

<figure><img src="image_128.png" /><figcaption>(2)</figcaption></figure>

<figure><img src="image_129.png" /><figcaption>(3)</figcaption></figure>

<figure><img src="image_130.png" /><figcaption>(4)</figcaption></figure>

[Govaert, 2011]
<figure><img src="image_131.png" /></figure>

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

Bi-clustering

■ A generalization of co-clustering

■ Look for submatrices of x which are homogeneous

■ We do not consider bi-clustering here

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

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

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

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

Link between co-clustering and PCA/kmeans<fnref n="2" />

**Interpreting LBM as an MBC dimension reduction method** PCA is certainly the most emblematic reduction dimension method for numerical data sets **x** (we consider in the following that **x** is centered in columns). It proceeds in two steps. First, it expresses each data unit **x**ᵢ of the initial data set **x** in a new vector basis (u¹, ..., uᵈ), ordered in decreasing value of the preserved variance, with respective coordinates (a₁ᵢ, ..., aₚᵢ) and where each **u**ᵢ is defined by a linear combination of the canonical vector basis (e¹, ..., eᵈ), namely **u** = ∑<sub>j=1</sub><sup>d</sup> bⱼ**e**<sup>j</sup>. Second, it selects a reduced number of these new coordinates only (say J), resulting in a new data set of smaller dimension (J < d). This PCA sequential procedure for the i-th data individual **x**ᵢ can be expressed as follows:

$$ \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_j^j \mathbf{u}^j. \quad (7) $$

Consider now the (very) specific co-clustering case reduced to K = 1, thus meaning that only a variable clustering in L column clusters is performed. Let ε<sup>j</sup> ~ N(0, σ<sup>2</sup><sub>wj</sub>) in an i.i.d manner, where ȳ<sub>j</sub> is the cluster index such that ȳ<sub>j</sub> = ℓ ⇔ w<sub>jℓ</sub> = 1. We can write, with **v**<sup>ℓ</sup> = ∑<sub>{j:ȳ<sub>j</sub>=ℓ}</sub> **e**<sup>j</sup> and **r**<sub>i</sub> = ∑<sub>j=1</sub><sup>d</sup> ε<sup>j</sup> **e**<sup>j</sup>:

$$ \mathbf{x}_i = \sum_{j=1}^{d} x_i^j \mathbf{e}^j = \sum_{j=1}^{d} (\mu_{\tilde{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. \quad (8) $$

It reinforces the interest of co-clustering for the HD case...

<footnote n="2">²C. Biernacki, J. Jacques, C. Keribin (2023). A Survey on Model-Based Co-Clustering: High Dimension and</footnote>
<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>

## Outline

1. HD clustering

2. Modeling

3. Estimating

4. Selecting

5. Experiments

6. To go further
<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>

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

Notations

■ $z_i$: the cluster of the row $i$

■ $w_j$: the cluster of the column $j$

■ $(z_i, w_j)$: the block of the element $x_{ij}$ (row $i$, column $j$)

■ $z = (z_1, \dots, z_n)$: partition of individuals in K clusters of rows

■ $w = (w_1, \dots, w_d)$: partition of variables in L clusters of columns

■ $(z, w)$: bi-partition of the whole data set $x$

■ Both space partitions are respectively denoted by $\mathcal{Z}$ and $\mathcal{W}$

Restriction

All variables are of the same kind (see discussion at the end)
<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>

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

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

# A geometric approach

* Example in the continuous case: $\mathbf{x} \in \mathbb{R}^{n \times d}$

* It could be possible to define a within-block inertia criterion

$$W(\mathbf{z}, \mathbf{w}) = \frac{\sum_{k=1}^{K} \sum_{l=1}^{L} \sum_{i=1}^{n} \sum_{j=1}^{d} z_{ik} w_{jl} \|x_{ij} - \mu_{kl}\|^2}{\sum_{i,j,k,l} \mu_{kl}}$$

with $\mu_{kl}$ the center of the block $(k, l)$

$$\mu_{kl} = \frac{1}{n_{kl}} \sum_{i,j} z_{ik} w_{jl} x_{ij}$$

where $n_{kl} = \sum_{ij} z_{ik} w_{jl}$ is the sample size of the block $(k, l)$

But we know now that it hides some model-based assumptions...
<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>

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

The latent block model (LBM)

■ Generalization of some existing non-probabilistic methods

■ Extend the latent class principle of local (or conditional) independence

■ Thus $x_{ij}$ is assumed to be independent once $z_i$ and $w_j$ are fixed ($\alpha = (\alpha_{kl})$):

$$p(\mathbf{x}|\mathbf{z}, \mathbf{w}; \boldsymbol{\alpha}) = \prod_{i,j} p(x_{ij}; \boldsymbol{\alpha}_{z_i w_j})$$

■ $\pi = (\pi_k)$: vectors of proba. $\pi_k$ that a row belongs to the kth row cluster

■ $\rho = (\rho_k)$: vectors of proba. $\rho_k$ that a row belongs to the lth column cluster

■ Independence between all $z_i$ and $w_j$

■ Extension of the traditional mixture model-based clustering ($\alpha = (\alpha_{kl})$):

$$p(\mathbf{x}; \boldsymbol{\theta}) = \sum_{(\mathbf{z}, \mathbf{w}) \in \mathcal{Z} \times \mathcal{W}} \prod_{i,j} \pi_{z_i} \rho_{w_j} p(x_{ij}; \boldsymbol{\alpha}_{z_i w_j})$$
<figure><img src="image_169.png" /></figure>

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

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

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

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

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

Distribution for different kinds of data

[Govaert and Nadif, 2014] The pdf $p(\cdot; \alpha_{z_i w_j})$ depends on the kind of data $x_{ij}$:

■ Binary data: $x_{ij} \in \{0, 1\}, p(\cdot; \alpha_{kl}) = B(\alpha_{kl})$

■ Categorical data with $m$ levels:

$x_{ij} = \{x_{ijh}\} \in \{0, 1\}^m$ with $\sum_{h=1}^m x_{ijh} = 1$ and $p(\cdot; \alpha_{kl}) = M(\alpha_{kl})$ with $\alpha_{kl} = \{\alpha_{kjh}\}$

■ Count data: $x_i^j \in \mathbb{N}, p(\cdot; \alpha_{kl}) = P(\mu_k \nu_l \gamma_{kl})^3$

■ Continuous data: $x_i^j \in \mathbb{R}, p(\cdot; \alpha_{kl}) = N(\mu_{kl}, \sigma_{kl}^2)$

<footnote n="3">³The Poisson parameter is here split into $\mu_k$ and $\nu_l$ the effects of the row $k$ and the column $l$ respectively and $\gamma_{kl}$ the effect of the block $kl$. Unfortunately, this parameterization is not identifiable. It is therefore not possible to estimate simultaneously $\mu_k$, $\nu_l$ and $\gamma_{kl}$ without imposing further constraints. Constraints</footnote>

<footnote>$\sum_k \pi_k \gamma_{kl} = \sum_l p_l \gamma_{kl} = 1$ and $\sum_k \gamma_{kl} = 1$, $\sum_l \mu_k = 1$, $\sum_l \nu_l = 1$ are a possibility.</footnote>
<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>

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

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

Extreme parsimony ability

<table><thead><tr><td>Model</td><td>Number of parameters</td></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>

Very parsimonious so well suitable for the (ultra) HD setting

$$nb. \ param_{HD} = nb. \ param_{classic} \times \frac{L}{d}$$

Other advantage: stay in the canonical space thus meaningful for the end-user
<figure><img src="image_181.png" /></figure>

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

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

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

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

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

Binary illustration: easy interpretation

[Govaert, 2011]

<table><thead><tr><th></th><th>abcdefghij</th></tr></thead><tbody><tr><td>Y1</td><td>1010001101</td></tr><tr><td>Y2</td><td>0101110011</td></tr><tr><td>Y3</td><td>1000001100</td></tr><tr><td>Y4</td><td>1010001100</td></tr><tr><td>Y5</td><td>0111001100</td></tr><tr><td>Y6</td><td>0101110101</td></tr><tr><td>Y7</td><td>0111110111</td></tr><tr><td>Y8</td><td>1100111101</td></tr><tr><td>Y9</td><td>0100110000</td></tr><tr><td>Y10</td><td>1010101101</td></tr><tr><td>Y11</td><td>1010001100</td></tr><tr><td>Y12</td><td>1010000110</td></tr><tr><td>Y13</td><td>1010001101</td></tr><tr><td>Y14</td><td>0010011100</td></tr><tr><td>Y15</td><td>0010010100</td></tr><tr><td>Y16</td><td>1111001100</td></tr><tr><td>Y17</td><td>0101110011</td></tr><tr><td>Y18</td><td>1010011101</td></tr><tr><td>Y19</td><td>1000010000</td></tr><tr><td>Y20</td><td>1100101100</td></tr></tbody></table>

Raw data

<table><thead><tr><th></th><th>acgh</th><th>bdefij</th></tr></thead><tbody><tr><td>Y2</td><td>0000</td><td>1 1 1 1 1</td></tr><tr><td>Y6</td><td>0001</td><td>1 1 1 1 0 1</td></tr><tr><td>Y7</td><td>0101</td><td>1 1 1 1 1 1</td></tr><tr><td>Y8</td><td>1010</td><td>1 0 1 1 1 0</td></tr><tr><td>Y9</td><td>0000</td><td>1 0 1 1 0 1</td></tr><tr><td>Y17</td><td>0000</td><td>1 1 1 1 1 1</td></tr><tr><td>Y1</td><td>1 1 1 1</td><td>0 0 0 0 0 1</td></tr><tr><td>Y3</td><td>1 0 1 1</td><td>0 0 0 0 0 0</td></tr><tr><td>Y4</td><td>1 1 1 1</td><td>0 0 0 0 0 0</td></tr><tr><td>Y5</td><td>1 0 1 1 1</td><td>1 1 0 0 0 0 0</td></tr><tr><td>Y10</td><td>1 1 1 1 1</td><td>0 0 1 0 0 1</td></tr><tr><td>Y11</td><td>1 1 1 1</td><td>0 0 0 0 0 0</td></tr><tr><td>Y12</td><td>1 1 0 1</td><td>0 0 0 0 0 0</td></tr><tr><td>Y13</td><td>1 1 1 1</td><td>0 0 0 0 0 1</td></tr><tr><td>Y14</td><td>1 1 1 1</td><td>0 0 0 0 0 0</td></tr><tr><td>Y15</td><td>1 1 1 1</td><td>0 0 0 1 0 0</td></tr><tr><td>Y16</td><td>0 1 1 1</td><td>0 0 0 1 0 0</td></tr><tr><td>Y17</td><td>0 1 0 1</td><td>0 0 1 0 0 0</td></tr><tr><td>Y18</td><td>1 1 1 1</td><td>1 0 0 0 0 0</td></tr><tr><td>Y19</td><td>1 1 1 1</td><td>0 0 0 1 0 0</td></tr><tr><td>Y20</td><td>0 1 1 1</td><td>0 0 0 0 0 0</td></tr></tbody></table>

Permuted data
 rows/columns

mode

<table><tr><td>0</td><td>1</td></tr><tr><td>1</td><td>0</td></tr></table>

Summary

<table><tr><td>0.86</td><td>0.79</td></tr><tr><td>0.83</td><td>0.86</td></tr></table>

Homogeneity

probamode

iid Bin(0.83)
<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>

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

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

Binary illustration: user-friendly visualization

[Govaert, 2011]

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

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

$$n = 500, d = 10, K = 6, L = 4$$
<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>

Other kind of data: ordinal

[Jacques and Biernacki, 2018]

<figure><img src="image_201.png" /><figcaption>Figure 11: Top 100 Amazone Fine Food Review data (left) and co-clustering result (right).</figcaption></figure>
<figure><img src="image_202.png" /></figure>

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

Other kind of data: functional

[Jacques, 2016]

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

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

Other kind of data: image

Original Data

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

Co-Clustered Data

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

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

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

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

Particular case: graph clustering

Stochastic Block Model (SBM): adjacency matrix with $n = d$ and $K = L$

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

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

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

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

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

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

## Outline

1. HD clustering

2. Modeling

3. Estimating

4. Selecting

5. Experiments

6. To go further
<figure><img src="image_231.png" /></figure>

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

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

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

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

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

MLE estimation: log-likelihood(s)

■ Remember Lesson 3: first estimate $\theta$, then deduce estimate of $(z, w)$

■ Observed log-likelihood: $\ell(\theta; x) = \ln p(x; \theta)$

■ MLE:

$$ \hat{\theta} = \arg\max_{\theta} \ell(\theta; x) $$

■ Complete log-likelihood:

$$ \ell_c(\theta; x, z, w) = \ln p(x, z, w; \theta) \\ = \sum_{i,k} z_{ik} \log \pi_k + \sum_{k,l} w_{jl} \log \rho_l + \sum_{i,j,k,l} z_{ik} w_{jl} \log p(x_i^j; \alpha_{kl}) $$

Be careful with asymptotics...

If $\ln(d)/n \to 0$, $\ln(n)/d \to 0$ when $n \to \infty$ and $d \to \infty$, then the MLE is consistent
[Brault et al., 2017]
<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>

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

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

MLE estimation: EM algorithm

■ E-step of EM (iteration $q$):

$$ 
\begin{aligned}
Q(\theta, \theta^{(q)}) &= E[\ell_c(\theta; \mathbf{x}, \mathbf{z}, \mathbf{w}) | \mathbf{x}; \theta^{(q)}] \\
&= \sum_{i,k} \underbrace{p(z_i = k | \mathbf{x}; \theta^{(q)}) \ln \pi_k}_{t_{ik}^{(q)}} + \sum_{j,l} \underbrace{p(w_j = l | \mathbf{x}; \theta^{(q)}) \ln \rho_l}_{s_{jl}^{(q)}} \\
&\quad + \sum_{i,j,k,l} \underbrace{p(z_i = k, w_j = l | \mathbf{x}; \theta^{(q)}) \ln p(x_{ij}; \alpha_{kl})}_{e_{ijkl}^{(q)}}
\end{aligned}
$$

■ M-step of EM (iteration $q$): classical. For instance, for the Bernoulli case, it gives

$$ 
\pi_k^{(q+1)} = \frac{\sum_j t_{ik}^{(q)}}{n}, \quad \rho_l^{(q+1)} = \frac{\sum_j s_{jl}^{(q)}}{d}, \quad \alpha_{kl}^{(q+1)} = \frac{\sum_{i,j} e_{ijkl}^{(q)} x_{ij}}{\sum_{i,j} e_{ijkl}^{(q)}}
$$
<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>

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

MLE: intractable E step

$e_{ijkl}^{(q)}$ is usually intractable...

* Consequence of dependency between $x_{ij}$s (link between rows and columns)

* Involve $K^n L^d$ calculus (number of possible blocks)

* Example: if $n = d = 20$ and $K = L = 2$ then $10^{12}$ blocks

* Example (cont'd): 33 years with a computer calculating 100,000 blocks/second

Alternatives to EM

* Variational EM (numerical approx.): conditional independence assumption

$$p(z, w|x; \theta) \approx p(z|x; \theta)p(w|x; \theta)$$

* SEM-Gibbs (stochastic approx.): replace E-step by a S-step approx. by Gibbs

$$z|x, w; \theta \text{ and } w|x, z; \theta$$
<figure><img src="image_249.png" /></figure>

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

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

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

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

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

MLE: variational EM (1/2)

* Use a general variational result from [Hathaway, 1985]

* Maximizing $\ell(\theta; x)$ on $\theta$ is equivalent to maximize $\tilde{\ell}_c(\theta; x, e)$ on $(\theta, e)$

$$ \tilde{\ell}_c(\theta; x, e) = \sum_{i,k} t_{ik} \ln \pi_k + \sum_{j,l} s_{jl} \ln \rho_l + \sum_{i,j,k,l} e_{ijkl} \ln p(x_{ij}; \alpha_{kl}) $$

where $e = (e_{ijkl})$, $e_{ijkl} \in \{0, 1\}$, $\sum_{k,l} e_{ijkl} = 1$, $t_{ik} = \sum_{j,l} e_{ijkl}$, $s_{jl} = \sum_{i,k} e_{ijkl}$

* Of course maximizing $\ell(\theta; x)$ or $\tilde{\ell}_c(\theta; x, e)$ are both intractable

* Idea: restriction on $e$ to obtain tractability $e_{ijkl} = t_{ik}s_{jl}$

* New variables are thus now $t = (t_{ik})$ and $s = (s_{jl})$

* As a consequence, it is a maximization of a lower bound of the max. likelihood

$$ \max_{\theta} \ell(\theta; x) \geq \max_{\theta,t,s} \tilde{\ell}_c(\theta; x, e) $$
<figure><img src="image_255.png" /></figure>

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

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

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

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

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

MLE: variational EM (2/2)

Approximated E-step

$$Q(\theta, \theta^{(q)}) \approx \sum_{i,k} t_{ik}^{(q)} \ln \pi_k + \sum_{j,l} s_{jl}^{(q)} \ln \rho_l + \sum_{i,j,k,l} t_{ik}^{(q)} s_{jl}^{(q)} \ln p(x_{ij}; \alpha_{kl})$$

* We called it now VEM

* Also known as mean field approximation

* Consistency of the variational estimate [Brault et al., 2017]
<figure><img src="image_261.png" /></figure>

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

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

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

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

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

MLE: local maxima

■ More local maxima than in classical mixture models

■ It is a consequence of many more latent variables (blocks)

■ Thus: either many VEM runs, or use the SEM-Gibbs algorithm
<figure><img src="image_267.png" /></figure>

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

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

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

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

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

MLE: SEM-Gibbs

* We have already seen the SEM algorithm in Lesson 3 (thus we do not detail more)

* It limits dependency to starting point, so it limits local maxima

* The S-step: a draw ($z^{(q)}, w^{(q)}) \sim p(z, w|x; \theta^{(q)})$ instead an expectation

* But it is still intractable, thus use a Gibbs algorithm to approximate this draw

Approximated S-step

Two easy draws

$$z^{(q)} \sim p(z|w^{(q-1)}, x; \theta^{(q)})$$

and

$$w^{(q)} \sim p(w|z^{(q)}, x; \theta^{(q)})$$

* Rigorously speaking, many draws within the S-step should be performed

* Indeed, Gibbs has to reach a stochastic convergence

* In practice it works well while saving computation time
<figure><img src="image_273.png" /></figure>

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

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

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

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

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

MLE: degeneracy

■ More degenerate situations than in classical mixture models

■ It is again a consequence of many more latent variables (blocks)

■ The Bayesian regularization (instead MLE) can be an answer
<figure><img src="image_279.png" /></figure>

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

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

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

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

<figure><img src="image_284.png" /><figcaption>Illustration of a degenerate situation</figcaption></figure>

<figure><img src="image_285.png" /></figure>
<figure><img src="image_286.png" /></figure>

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

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

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

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

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

Bayesian estimation: pitch

■ Everything passes by the **posterior distribution of θ**

$$p(\theta|x) \propto \underbrace{p(x|\theta)}_{\text{log-likelihood}} \underbrace{p(\theta)}_{\text{prior}}$$

■ Then, take (for instance) the **MAP** as a **θ** estimate (use a VEM like algo...)

$$\hat{\theta} = \arg \max_{\theta} p(\theta|x)$$
<figure><img src="image_292.png" /></figure>

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

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

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

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

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

Bayesian estimation: limiting degeneracy

■ Interest for avoiding degeneracy is the prior: it acts as a **penalization term**

■ Typical choices are **Dirichlet** for $\pi$ and $\rho$ (with independence between $\pi, \rho, \alpha$)

$$p(\theta) = \underbrace{p(\pi)}_{D_K(a, \dots, a)} \times \underbrace{p(\rho)}_{D_L(a, \dots, a)} \times \underbrace{p(\alpha)}_{\text{model dependent}}$$

■ The **Dirichlet** distribution is conjugate, thus easy calculus

■ Control degeneracy frequency with the a value:

■ a = 1: uniform prior, so $\hat{\theta}$ is strictly the MLE (no regularisation)

■ a = 1/2: Jeffreys prior, classical (no informative prior) but may favor degeneracy

■ a = 4: a rule of thumb working well for limiting degeneracy frequency
<figure><img src="image_298.png" /></figure>

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

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

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

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

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

Bayesian estimation: prior overview

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

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

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

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

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

<figure><img src="image_309.png" /></figure>
<figure><img src="image_310.png" /></figure>

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

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

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

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

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

Block estimation: estimate

■ Once we have a parameter estimate $\hat{\theta}$, we need to have an block estimate ($\hat{z}, \hat{w}$)

■ But MAP not directly available because of the following maximization difficulty

$$ (\hat{z}, \hat{w}) = \arg\max_{(z,w)} \underbrace{p(z, w|x; \hat{\theta})}_{\text{intractable}} $$

■ Instead the following (easily, as classical mixtures) estimates are usually retained

$$ \hat{z} = \arg\max_{z} p(z|x; \hat{\theta}) \quad \text{and} \quad \hat{w} = \arg\max_{w} p(w|x; \hat{\theta}) $$
<figure><img src="image_316.png" /></figure>

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

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

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

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

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

Block estimation: evaluation

■ Empirical error rate between blocks:

$$ \text{err}_{\text{blocks}}(\underbrace{(z, w)}_{\text{"True" blocks}}, \underbrace{(\hat{z}, \hat{w})}_{\text{Estimated blocks}}) = \text{err}(z, \hat{z}) + \text{err}(w, \hat{w}) - \text{err}(z, \hat{z}) \times \text{err}(w, \hat{w}) $$

■ Rand index between blocks: it exists also a recent definition...
<figure><img src="image_322.png" /></figure>

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

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

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

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

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

Block estimation: consistency

[Mariadassou and Matias, 12]

$$ \underbrace{\hat{\theta}^{n,d \to \infty} \to \theta^*}_{\text{we have seen that...}} \Rightarrow \underbrace{p(\hat{z} = z^*, \hat{w} = w^* | x; \hat{\theta})^{n,d \to \infty}}_{\text{exact bi-partition retrieval!}} \mathbf{1} $$

Thus we retrieve the HD clustering blessing...
<figure><img src="image_328.png" /></figure>

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

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

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

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

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

Block estimation: non asymptotic properties (1/2)

Binary case: marginals seems so simple mixtures! [Brault, 14]

Matrice initiale

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

Lignes

Colonnes

Histogrammes
des semaines

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

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

Separations

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

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

Matrice reorga-
nise

<figure><img src="image_339.png" /></figure>
<figure><img src="image_340.png" /></figure>

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

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

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

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

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

Block estimation: non asymptotic properties (2/2)

[Brault, 14]

■ Probability of $x_{ij}$ with no regard to the column membership is Bernoulli

$$p(x_{ij} = 1 | z_{ik} = 1) = \tau_k = \sum_{l=1}^{L} \alpha_{kl} \rho_l$$

■ Thus marginal distribution of $x_{ij}$ is a mixture (indep. of $x_{ij}$ cond. $z_{ik} = 1$)

$$\left( \sum_j x_{ij} \right) |z_{ik} = 1 \sim B(d, \tau_k)$$

■ Control of error on this partition mixture estimate $\hat{z}^{mix}$ of binomial distributions

$$p(\hat{z}^{mix} \neq z^*) \leq 2n \exp \left\{ -\frac{1}{8} d \left[ \min_{k \neq k'} |\tau_k - \tau_{k'}| \right] \right\} + K \left(1 - \min_k \pi_k\right)^n$$

■ We retrieve also consistency for very high dimension with constraint

$$\ln(n) = o(d)$$
<figure><img src="image_346.png" /></figure>

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

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

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

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

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

## Illustration: document clustering (1/2)

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

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

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

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

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

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

## Illustration: document clustering (2/2)

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

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

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

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

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

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

## Outline

1. HD clustering

2. Modeling

3. Estimating

4. Selecting

5. Experiments

6. To go further
<figure><img src="image_366.png" /></figure>

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

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

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

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

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

Models in competition

$\mathbf{m} = (K, L)$ typically, but not restricted to
<figure><img src="image_372.png" /></figure>

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

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

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

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

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

BIC criterion: two difficulties

■ Difficult 1: which BIC definition because of the double asymptotic on n and d?

■ Difficult 2: the observed log-likelihood value is intractable

$$ \ell(\theta; \mathbf{x}) = \sum_{(z, \mathbf{w}) \in \mathcal{Z} \times \mathcal{W}} p(\mathbf{x}, z, \mathbf{w}; \theta) $$

Could be estimated by harmonic mean but time consuming and high variance
ICL criterion: overcome both difficulties

■ ICL uses complete likelihood thus no intractability

$$ \text{ICL} = \ln p(\mathbf{x}, \hat{\mathbf{z}}, \hat{\mathbf{w}}) = \ln p(\mathbf{x}|\hat{\mathbf{z}}, \hat{\mathbf{w}}) + \ln p(\hat{\mathbf{z}}) + \ln p(\hat{\mathbf{w}}) $$

■ Multinomial case (r levels): [Keribin et al., 2014]

■ Derive an exact (non-asymptotic) ICL version

■ Deduce an asymptotic approximation of ICL

$$ \text{ICLbic} = \ell_c(\hat{\theta}; \mathbf{x}, \hat{\mathbf{z}}, \hat{\mathbf{w}}) - \frac{K-1}{2} \ln(n) - \frac{L-1}{2} \ln(d) - \frac{KL(r-1)}{2} \ln(n) $$

■ We can make a conjecture for the general case

$$ \text{ICLbic} = \ell_c(\hat{\theta}; \mathbf{x}, \hat{\mathbf{z}}, \hat{\mathbf{w}}) - \frac{K-1}{2} \ln(n) - \frac{L-1}{2} \ln(d) - \frac{KL\nu_{\alpha kl}}{2} \ln(n) $$
<figure><img src="image_378.png" /></figure>

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

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

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

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

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

ICL criterion: consistency

■ We can obtain a BIC expression from ICLbic

$$
\begin{align*}
\text{BIC} &= \text{ICLbic} - \ln p(\hat{z}, \hat{w}|\mathbf{x}; \hat{\theta}) \\
&= \underbrace{\ell(\hat{\theta}; \mathbf{x})}_{\text{difficult}} - \frac{K-1}{2} \ln(n) - \frac{L-1}{2} \ln(d) - \frac{KL(m-1)}{2} \ln(nd)
\end{align*}
$$

■ [Brault et al., 2017] establish that asymptotically on *n* and *d*

$$
“\ell(\hat{\theta}; \mathbf{x}) = \ell_c(\hat{\theta}; \mathbf{x}, \hat{z}, \hat{w})”
$$

■ Thus, since BIC is consistent, ICL is also consistent

Again the HD clustering blessing is here!
<figure><img src="image_384.png" /></figure>

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

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

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

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

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

Strategy to smart browsing of (K, L)

<figure><img src="image_390.png" /><figcaption>[Robert, 2017] Algorithm Bi-KM1</figcaption></figure>
<figure><img src="image_391.png" /></figure>

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

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

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

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

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

Illustration: discuss the dimension (1/2)

■ SPAM E-mail Database<fnref n="4" />

■ $n = 4601$ e-mails composed by 1813 “spams” and 2788 “good e-mails”

■ $d = 48 + 6 = 54$ continuous descriptors<fnref n="5" />

■ 48 percentages that a given word appears in an e-mail ("make", "you"...)

■ 6 percentages that a given char appears in an e-mail ("a", "b", "c"...)

■ Transformation of continuous descriptors into binary descriptors

$$x_{ij} = \begin{cases} 1 & \text{if word/char } j \text{ appears in e-mail } i \\ 0 & \text{otherwise} \end{cases}$$

<footnote n="4">⁴https://archive.ics.uci.edu/ml/machine-learning-databases/spambase/</footnote>

<footnote n="5">⁵There are 3 other continuous descriptors we do not use</footnote>
<figure><img src="image_397.png" /></figure>

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

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

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

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

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

Illustration: discuss the dimension (2/2)

■ Perform co-clustering with $K = 2$ and $L = 5$: ICL=-92,682, err=0.1984

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

■ Perform clustering<fnref n="6" /> with $K = 2$: ICL=-89,433, err=0.1837

Thus use preferably co-clustering in the HD setting

<footnote n="6">⁶Equivalent to co-clustering with $L = 54$</footnote>
<figure><img src="image_404.png" /></figure>

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

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

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

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

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

## Outline

1. HD clustering

2. Modeling

3. Estimating

4. Selecting

5. Experiments

6. To go further
<figure><img src="image_410.png" /></figure>

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

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

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

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

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

# The Blockcluster package

## BlockCluster

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

Visit notebooks repo

BlockCluster can estimate the parameters of co-clustering models for binary, contingency and continuous data. Simply put, when considering a set of data as rows and columns, BlockCluster will make simultaneous permutations of rows and columns in order to organise the data into homogeneous blocks.

Here is a related publication for this software. Here is the related package on the CRAN.

## Two ways to be used

* On the CRAN:
  https://cran.r-project.org/web/packages/blockcluster/index.html

* As a note book (nothing to install):
  https://team.inria.fr/modal/software/notebooks/
<figure><img src="image_417.png" /></figure>

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

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

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

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

# A synthetic example

## Blockcluster basic example

### Data

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

### Clustering with Blockcluster

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

Define the strategy.

Specify the number of group clusters to create in rows and columns.

<figure><img src="image_424.png" /></figure>
HD clustering

Modeling

Estimating

Selecting

Experiments
●

To go further

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

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

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

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

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

# A synthetic example

## Output's Analysis

### Model criterion

This chart represents the criterion value for each model that was built. The lower the value (close to 0) the better the model.

In [7]:

res@ICLvalue

-45557.0741834614

### Cluster Plot

Draw the original data matrix, and the matrix obtained after performing co-clustering

In [8]:

plot(res)

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

<figure><img src="image_431.png" /></figure>
<figure><img src="image_432.png" /></figure>

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

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

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

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

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

## Outline

1. HD clustering

2. Modeling

3. Estimating

4. Selecting

5. Experiments

6. To go further
<figure><img src="image_438.png" /></figure>

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

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

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

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

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

# Co-clustering of mixed data

* Same partitions in lines, disjoint partitions in columns

* Example: data set TED talks, with talks × (terms,scores)

<figure><img src="image_444.png" /><figcaption>Poisson</figcaption></figure>

<figure><img src="image_445.png" /><figcaption>Gaussian</figcaption></figure>