# Multiway latent block model for pca tensor decomposition

Rodolphe Priam

► To cite this version:

Rodolphe Priam. Multiway latent block model for pca tensor decomposition. 2025. hal-04803423v5

HAL Id: hal-04803423

https://hal.science/hal-04803423v5

Preprint submitted on 23 Mar 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
# Multiway latent block model for pca tensor decomposition*

R.Priam†

March 23, 2025

## Abstract

A multiway latent block model (WLBM) is proposed for modeling the reduction and the clustering of m-modes tensors via material latent variables. BEM-like algorithms allow the parameters inference of the models for tensors with multiple modes. This results into summarizing high-dimensional tensors for pca reduction and visualization in a generic multidimensional setting. The parameterizations are -related or approximated- Canonical Polyadic and Tucker decompositions with dramatically less dimensions than their usual full versions for the latent components.

## 1 Introduction

When data matrices are large, a clustering can give a quicker access to the data contents than a projection method for reducing the dimensionality of the features. Combining clustering and reduction for the projection of the clusters rather than the projection of the rows or the columns is therefore a promising purpose for data analysis. In visualization, one way to reach this requirement is to show the clusters -plus the rows or columns- on a map after clustering the data by an ad hoc algorithm and reducing the dimensions of the matrix, both separately. For reduction or decomposition/factorization, a related usual way is to construct reduced vector or views of the rows and columns in order to summarize with few dimensions. From a data matrix to a data tensor, one or several modes are added such that the method needs to manage a different algebraic space at three or more modes, otherwise some generative models may be more suitable here. The family of co-clustering methods makes possible to reveal the hidden association between the rows and the columns of a data table and the modes in a tensor, after that k-means [1] has been improved from intensive researches in clustering. A simultaneous partitioning treats symmetrically the data table or tensor on the contrary to a clustering of just only one mode. Algorithms were introduced earlier in the literature in [2] and [3, 4]. There is also the information-theoretic co-clustering method [5], its generalization [6] for contingency matrices. Other approaches [7] are for instance the general method for prediction [8] or the non negative matricial decompositions [9, 10, 11] for non negative matrice. In [12, 13, 14], the Latent Block Model (LBM) and its inference by

a variational generalized EM algorithm [15] are introduced. Reparameterization of these models via latent variables [16, 17, 18] bring new properties to these models for: rows and columns visualizations, dimensions reductions or even tensor decomposition as proposed herein. Clustering decreases the size of the latent vectors for the reduced views from the modes, this approach is expected to be more scalable than the usual decompositions thanks to the homogeneous blocks summarizing the tensor and dividing the computation burden. Including latent variables for reduction in LBM brings new properties to this model for analysis while complementing recent works to deal with big or huge tensors. Current libraries available as open sources in Python (or R) language are for instance found in [19, 20, 21, 22, 23] for clustering or decomposition. Note that the vocabury in literature on tensor is diverse with diverse terms like clustering, co-clustering, tri(-)clustering, biclustering, partitioning for groups but also related terms like reordering with coordinates in reduced views of the different modes. Next, distributions are introduced for tensors in order be more general because k-means is known for its non optimality due to its spherical equiprobable clusters only relevant for continuous data despite its good results on some datasets. The plan of the paper is as follows. Section 1 is for introducing the problematic. Section 2 is for presenting the foundations for the proposed approach. Section 3 is presents the new multiway latent block model with also the latent variables for the decomposition. Section 4 illustrates the model with the Gaussian cases by presenting three variants, one for co-clustering, one for sparse co-clustering and one for the decomposition or principal component analysis. Section 5 is for the experiments with tensor data for reduction and visualization plus a discussion and comparison with alternative decompositions. Section 6 is the conclusion with the perspectives.

## 2 Review of foundation models

This section revise several models for reduction or clustering, they are principal component analysis, decomposition/factorization, and tri-clustering. This introduced the basis for the following proposed approach which proceeds via hidden latent variables within a probabilistic model.

### Principal component analysis (PCA) for matrix (m=2)

Factorization [24] is a living activity for reducing the dimen-

<footnote>*Fifth version of the draft: some notations are improved or fixed and some sections are changed for better presentation while new experiments are included.</footnote>

<footnote>†rpriam@gmail.com.</footnote>
sionality of matrices. For helping at denoising and clustering, principal component analysis (PCA) [25] aims at the construction of synthetic variables from the original variables or columns of the data matrix. These new variables maximize an explained variance via a projection from the full data cloud to a subspace spanning the columns space. The limited number of new variables or components in the new vectors that replace the (possibly high-dimensional) original one leads to the property of reducing the dimensionality. A new matrix is available with less columns (or rows) than the original matrix. For all these reasons, this method is useful in numerous research areas from an applied statistical analysis to a pre-reduction for another method such that clustering. More formally for $m = 2$ modes, let $\mathbf{x}$ the numerical matrix with $n$ rows $d$-dimensional $\mathbf{x}_i$ and thus $n \times d$ cells $x_{ij}$. PCA has several definitions actually in the literature according to diverse purposes. There can be bayesian priors for the reduced vectors, orthogonality and/or positivity/non negativity for the loading matrix, or even a factorization without any other hypothesis. For a matrix with cells $(x_{ij})_{ij}$, the following definitions are widely considered:

* For centered data, PCA is a matricial method [26] which is originally defined as solving for:

$$ \hat{\mathbf{R}} = \underset{\mathbf{R}}{\operatorname{argmin}} \sum_i \|\mathbf{x}_i - \mathbf{R}\mathbf{R}^T \mathbf{x}_i\|^2 $$

This is equivalent to define as minimizing the reconstruction error of the projected data points into the new subspace, a linear transformation for reducing the space of the columns. Eventually, orthogonality supposes $\mathbf{R}^T \mathbf{R} = \mathbf{I}$ which often improves the results.

* A factorization [24] is written as follows,

$$ (\xi^{\mathbf{x}^{(1)}}, \hat{\xi}^{\mathbf{x}^{(2)}}) = \underset{\xi^{\mathbf{x}^{(1)}}, \xi^{\mathbf{x}^{(2)}}}{\operatorname{argmin}} \sum_{i,j} (x_{ij} - \xi_i^{\mathbf{x}^{(1)}T} \xi_j^{\mathbf{x}^{(2)}})^2 . $$

A solution is found by regression for solving for $\xi^{\mathbf{x}^{(1)}}$ and keeping $\xi^{\mathbf{x}^{(2)}}$ constant, and similarly for $\xi^{\mathbf{x}^{(2)}}$ when $\xi^{\mathbf{x}^{(1)}}$ is constant. This may be the definition the more convenients for modeling a PCA(-like) method within non Gaussian models, with the distribution function for $x_{ij}$ as depending on these latent vectors or instead cluster related versions in the case of a partitioning as herein.

* In Probabilistic PCA (PPCA [27], it may be written also for non centered data in a cluster:

$$ (\hat{\mu}, \hat{W}) = \underset{\mu, W}{\operatorname{argmin}} \sum_i \|\mathbf{x}_i - \mu - \mathbf{W}\mathbf{y}_i\|^2 . $$

This is the approach used in the probabilistic setting of PCA where this notation is more usual in this family of models. Note that it may not be required the hypothesis of orthogonality anymore for some authors. Without this constraint, The matrix $\mathbf{W}$ may not correspond to the solution of PCA $\mathbf{R}$ for an orthogonal projection resulting into principal latent variables $\mathbf{y}_i = \mathbf{R}^T \mathbf{x}_i$, like in an usual autoencoder [28] without nonlinearities and deep learning.

Note that the three forms of PCA listed above may be not yet directly converted into a co-clustering setting while the second one is intensively studied for tensors.

## Decomposition for m-modes tensor (m=3)

For PCA of a tensor data, there exists also several definitions - called multilinear PCA - which may be not always related to the former ones above and which may be out of the scope herein. When PCA is identified to a factorization, as often in generative models, the nearest and direct implementation of this idea is the **Canonical Polyadic decomposition (CP-D)** and **Tucker decomposition (TK-D)**, see [29, 30, 31, 32, 33, 34, 35, 36, 37], which finds reduced latent matrices $\xi^{\mathbf{x}^{(1)}}, \xi^{\mathbf{x}^{(2)}}, \xi^{\mathbf{x}^{(3)}}$ for each mode of the tensor associated to scaling constants aggregated in $\Lambda$. Here for three modes:

$$ (\hat{\Lambda}, \hat{\xi}^{\mathbf{x}^{(1)}}, \hat{\xi}^{\mathbf{x}^{(2)}}, \hat{\xi}^{\mathbf{x}^{(3)}}) = \left\{ \begin{array}{ll} \underset{\Lambda, \xi^{\mathbf{x}^{(1)}}, \xi^{\mathbf{x}^{(2)}}, \xi^{\mathbf{x}^{(3)}}}{\operatorname{argmin}} \sum_{ijk} \left( x_{ijk} - \sum_s \lambda_s^\mathbf{x} \xi_{is}^{\mathbf{x}^{(1)}} \xi_{js}^{\mathbf{x}^{(2)}} \xi_{ks}^{\mathbf{x}^{(3)}} \right)^2, & \text{CP-D} \\ \\ \underset{\Lambda, \xi^{\mathbf{x}^{(1)}}, \xi^{\mathbf{x}^{(2)}}, \xi^{\mathbf{x}^{(3)}}}{\operatorname{argmin}} \sum_{ijk} \left( x_{ijk} - \sum_{s_1 s_2 s_3} \lambda_{s_1 s_2 s_3}^\mathbf{x} \xi_{is_1}^{\mathbf{x}^{(1)}} \xi_{js_2}^{\mathbf{x}^{(2)}} \xi_{ks_3}^{\mathbf{x}^{(3)}} \right)^2, & \text{TK-D.} \end{array} \right. $$

Note that in order to induce a better fitting, some data may need additional scalars depending on the rows and columns, $\gamma_i$ and $\gamma_j$, to add just after the overall products within the squares above. The algorithm k-means has been already extended to co-clustering before too, hence extending k-means [38] to more modes - thus $m > 2$ - for tensors is natural. By adding the classification variables $\mathbf{z}_\ell = (z_{i_\ell k_\ell})$ for each dimension, and when aggregating the means as $\mu = (\mu_{k_1 k_2 k_3})$, it may be written:

$$ (\hat{\mu}, \hat{z}_1, \hat{z}_2, \hat{z}_3) = \underset{\mu, \mathbf{z}_1, \mathbf{z}_2, \mathbf{z}_3}{\operatorname{argmin}} C_{tri}^{123}(\mu, \mathbf{z}_1, \mathbf{z}_2, \mathbf{z}_3). $$

When this criterion - see Figure 1 for the fulll expression - is solved, this leads to a tri-clustering of the tensor for three modes -also called more generally co-clustering or clustering or eventually multi-clustering - in the spirit of k-means. An extension comes from adding to the criterion, entropies and replacing the euclidean distance by the log of a density function, which allows to define a new criterion with a probabilistic setting. This extends the former classical mixture model with only one classification variable for the rows or the columns separately, with several classification variables for all the modes simultaneously.

Herein, this is the latent block model which is extended for tensors clustering while including latent variables for the decomposition. The two decompositions presented above are modeled in a more economical way next sections after introducing the multiway latent block model as the backbone of the proposed reduction models and the algorithms for their parameters inference for tensors. Several variants are also considered in order to improve the decomposition and the visualization because this sparse setting can be replaced in order to add more constraints for different results to be compared. While tensors are in large dimensions,
instead of $C_{tri}^{123}$ as follows, with now $z_{i\ell k\ell}$ not anymore binary but fuzzy replaced by $c_{i\ell k\ell}^{\ell}$ in $[0; 1[$ and aggregated in $c^{\ell}$, with,

$$C_{tri,f}^{123}(\theta, \mathbf{c}^1, \mathbf{c}^2, \mathbf{c}^3) = \sum_{i_1 i_2 i_3} \sum_{k_1 k_2 k_3} c_{i_1 k_1}^1 c_{i_2 k_2}^2 c_{i_3 k_3}^3 f(\mathbf{x}, \theta_{k_1 k_2 k_3}) - \sum_{\ell} \sum_{i_\ell k_\ell} c_{i_\ell k_\ell}^{\ell} \log c_{i_\ell k_\ell}^{\ell} + \sum_{i_1 i_2 i_3} \sum_{k_1 k_2 k_3} c_{i_\ell k_\ell}^{\ell} p_{\ell k_\ell}.$$

The idea of tri-clustering next after will be to optimize in an alternated way for $c_{i_1 k_1}^1$, $c_{i_2 k_2}^2$, $c_{i_3 k_3}^3$ at a first step and all the parameters in $\theta$ at the second step. The latent block model comes to this solution thanks to a probabilistic setting via a statistical parametric likelihood approximated by the function above for the inference of the parameters. But if under this model, the function $f(.)$ is generally a density probability function or a mass probability function, it is clear that other function may be used here, like robust distances for instance. Next, mostly the Bernoulli, Gaussian and Poisson probability functions for $f(.)$ are discussed in order to deal with tensors with binary, continuous or count values.

# 3 Multiway latent block model (WLBM)

A review of the latent block model and its general definition summarizes briefly the modeling foundation before presenting our proposal in this section.

## 3.1 Former 2-modes model and algorithm ($m=2$)

### Model LBM (or $W_2$LBM)

Within the context of the classical mixture model for a matrix with $n$ rows and $d$ columns, a partition of $I$ into $g$ clusters is represented by the binary classification matrix $\mathbf{z} = (z_{ik})_{n \times g}$ such that $\sum_{k=1}^g z_{ik} = 1$ and $z_{ik} = 1$ indicates the component of the row $i$. Just as $I$ is partitioned into $g$ clusters, columns can be partitioned into $m$ clusters by the binary classification matrix $\mathbf{w} = (w_{j\ell})_{d \times m}$. If the most usual clustering methods deal with partitioning of only a set of rows $(x_i)_{i=1}^n$ with indices $I = \{1, 2, \dots, n\}$ or eventually a set of columns $(x^j)_{j=1}^d$ with indices $J = \{1, 2, \dots, d\}$, co-clustering is interested in the clustering of both simultaneously. The set of all possible assignments $\mathbf{w}$ of $J$ (resp. $\mathbf{z}$ of $I$) is denoted $\mathcal{W}$ (resp. $\mathcal{Z}$). The $n \times d$ random variables generating the observed $x_{ij}$ cells of the data matrix are assumed to be independent in LBM, once $\mathbf{z}$ and $\mathbf{w}$ are fixed. For the parameters optimization, the completed data are taken to be the vector $(\mathbf{x}, \mathbf{z}, \mathbf{w})$ where the latent variables $\mathbf{z}, \mathbf{w}$ are the random labels for the rows and the columns respectively. Let denote the different density and mass functions as follows:

$$P(\mathbf{z}, \mathbf{w}) = \prod_{i,k} (p_k)^{z_{ik}} \prod_{j,\ell} (q_\ell)^{z_{j\ell}} = P(\mathbf{z})P(\mathbf{w})$$

$$P(\mathbf{x}|\mathbf{z}, \mathbf{w}) = \prod_{i,i,k,\ell} \varphi(x_{ij}; \theta_{k\ell})^{z_{ik}w_{j\ell}}.$$

The classification or completed log-likelihood $L_C$ can then be written conditionally to all the labels at each mode:

$$L_C(\mathbf{z}, \mathbf{w}; \mathbf{x}, \theta) = \log P(\mathbf{x}|\mathbf{z}, \mathbf{w}) + \log P(\mathbf{z}) + \log P(\mathbf{w}).$$

Hence, with the data $\mathbf{x}$ seen as a set of cells $(x_{11}, x_{12}, \dots, x_{ij}, \dots, x_{nd})$, these observations are independently generated conditionally to the row and columns labels. The two sets of possible assignments $\mathbf{w}$ and $\mathbf{z}$ aggregate the cells of the matrix $\mathbf{x}$ into a number of contiguous, non-overlapping blocks. The following sum is obtained [12] by independence of $\mathbf{z}$ and $\mathbf{w}$, by summing over all the assignments $\mathcal{Z} \times \mathcal{W}$ of rows and columns clusters:

$$f_{LBM}(\mathbf{x}; \theta) = \sum_{(\mathbf{z}, \mathbf{w})} e^{L_C(\mathbf{z}, \mathbf{w}; \mathbf{x}, \theta)} = \sum_{(\mathbf{z}, \mathbf{w})} \prod_{i,k} p_k^{z_{ik}} \prod_{j,\ell} q_\ell^{w_{j\ell}} \prod_{i,j,k,\ell} \varphi(x_{ij}; \theta_{k\ell})^{z_{ik}w_{j\ell}}.$$

Here $\varphi(\cdot; \theta_{k\ell})$ is a probability function defined on the real line $\mathbb{R}$ or the set of integers $\mathbb{N}$, while $\theta = \{\theta_{k\ell}; 1 \le k \le g, 1 \le \ell \le m\}$ are the unknown parameters. The vectors of the probabilities $p_k$ and $q_\ell$ that a row and a column belong to the $k$th component and to the $\ell$th component are respectively denoted $\mathbf{p} = (p_1, \dots, p_g)$ and $\mathbf{q} = (q_1, \dots, q_m)$. The set of parameters is denoted $\theta$ and is compound of $\mathbf{p}$ and $\mathbf{q}$ plus $\alpha$ which aggregates the $g \times m$ scalar $\alpha_{k\ell}$ via the matrix $(\alpha_{k\ell})$. Hereafter, to simplify the notation, the sums and the products relating to rows, columns or clusters will be subscripted respectively by the letters $i, j, k$, or $\ell$ for matrices and more generally $i_1, i_1, \dots, i_m$ and $k_1, k_1, \dots, k_m$ for tensors, without indicating the limits of variation, which are implicit.

The set of parameters $\theta$ of the model can be estimated by maximizing the log-likelihood:

$$L(\theta; \mathbf{x}) = \log f_{LBM}(\mathbf{x}; \theta).$$

A particular distribution is needed for the cells, for continuous, binary and contingency matrices with originally Gaussian, Bernoulli and Poisson distributions where the parameters are expectations plus sometimes other parameters, say $(\mu_{k\ell}$ or $\alpha_{k\ell}, \sigma_{k\ell})$, $\alpha_{k\ell}$ and $\lambda_{k\ell}^{ij} = \mu_i v_j \alpha_{k\ell}$ respectively. Note that the effects $\mu = (\mu_1, \dots, \mu_n)$ and $v = (v_1, \dots, v_d)$ are assumed equal to the following constant margin totals by rows and by columns, $\mu_i = \sum_j x_{ij}$ for $i \in I$ and $v_j = \sum_i x_{ij}$ for $j \in J$ for contingency matrices.

### Objective function and optimization algorithm

For the proposed model with the introduced constraints, we aim to address the problem of the estimation of the parameters by a maximum likelihood (ML) approach such that:

$$\hat{\theta} = \operatorname*{argmax}_{\theta} L(\theta; \mathbf{x}).$$

Let us focus on the estimation of a value of $\theta$ by the maximum likelihood approach associated to the block mixture model. For this model, the complete data are $(\mathbf{x}, \mathbf{z}, \mathbf{w})$ where the unobservable vectors $\mathbf{z}$ and $\mathbf{w}$ are the labels. The complete log-likelihood of $(\mathbf{x}, \mathbf{z}, \mathbf{w})$ leads to an Expectation-maximization algorithm (EM) [46] in two steps, where the E-step is intractable, thus an alternative approach has been proposed in the literature.

An approach based on a Generalized EM and a variational approximation by the product $c_{ik}^{(t)} d_{j\ell}^{(t)}$ has been proposed previously in the literature, and named Block EM (BEM). Here it is denoted the variational probabilities $c_{ik}$ such that $\sum_k c_{ik} = 1$, and
<table><caption>Table 2: Examples of p.d.f. or p.m.f. with ($\beta_{1}^{m} = \mu_{i1} \cdots \mu_{im}$ and $\mu_{i\ell} = \sum_{l \neq \ell} x_{i1 \cdots i_m}$), from $\phi_E(x_{i1}^m; \theta) = e^{(x_{i1}^m A(\alpha_{k_1}^{i1}) - B(\alpha_{k_1}^{i1}) + C(x_{i1}^m))}$.</caption><thead><tr><th>Model</th><th>pdf of cell (k<sub>1</sub><sup>m</sup>)</th><th>&phi;</th><th>A(&alpha;<sub>k<sub>1</sub><sup>m</sup></sub><sup>2</sup>)</th><th>B(&alpha;<sub>k<sub>1</sub><sup>m</sup></sub><sup>2</sup>)</th><th>&beta;<sub>1</sub><sup>m</sup></th><th>&alpha;<sub>k<sub>1</sub><sup>m</sup></sub></th><th>x<sub>1</sub><sup>m</sup></th></tr></thead><tbody><tr><td>G-W<sub>m</sub>LBM</td><td>N(&alpha;<sub>k<sub>1</sub><sup>m</sup></sub>; &sigma;<sub>k<sub>1</sub><sup>m</sup></sub>)</td><td>exp(&minus;(x<sub>1</sub><sup>m</sup> &minus; &alpha;<sub>k<sub>1</sub><sup>m</sup></sub>)<sup>2</sup>/2&sigma;<sub>k<sub>1</sub><sup>m</sup></sub><sup>2</sup>) / &radic;(2&pi;&sigma;<sub>k<sub>1</sub><sup>m</sup></sub>)</td><td>&alpha;<sub>k<sub>1</sub><sup>m</sup></sub>/&sigma;<sub>k<sub>1</sub><sup>m</sup></sub><sup>2</sup></td><td>&alpha;<sub>k<sub>1</sub><sup>m</sup></sub><sup>2</sup>/2&sigma;<sub>k<sub>1</sub><sup>m</sup></sub><sup>2</sup></td><td>1</td><td>&Re;</td><td>&Re;</td></tr><tr><td>B-W<sub>m</sub>LBM</td><td>B(&alpha;<sub>k<sub>1</sub><sup>m</sup></sub>)</td><td>(&alpha;<sub>k<sub>1</sub><sup>m</sup></sub>)<sup>x<sub>1</sub><sup>m</sup></sup>(1 &minus; &alpha;<sub>k<sub>1</sub><sup>m</sup></sub>)<sup>1 &minus; x<sub>1</sub><sup>m</sup></sup></td><td>log (1/&alpha;<sub>k<sub>1</sub><sup>m</sup></sub>)</td><td>&minus;log(1 &minus; &alpha;<sub>k<sub>1</sub><sup>m</sup></sub>)</td><td>1</td><td>[0; 1]</td><td>{0, 1}</td></tr><tr><td>P-W<sub>m</sub>LBM</td><td>P(&beta;<sub>1</sub><sup>m</sup>; &alpha;<sub>k<sub>1</sub><sup>m</sup></sub>)</td><td>exp(&minus;&beta;<sub>1</sub><sup>m</sup>&alpha;<sub>k<sub>1</sub><sup>m</sup></sub>)(&beta;<sub>1</sub><sup>m</sup>&alpha;<sub>k<sub>1</sub><sup>m</sup></sub>)<sup>x<sub>1</sub><sup>m</sup></sup>/x<sub>1</sub><sup>m</sup>!</td><td>log &beta;<sub>1</sub><sup>m</sup>&alpha;<sub>k<sub>1</sub><sup>m</sup></sub></td><td>&beta;<sub>1</sub><sup>m</sup>&alpha;<sub>k<sub>1</sub><sup>m</sup></sub></td><td>&Pi;<sub>&ell;</sub> &mu;<sub>i<sub>&ell;</sub></sub></td><td>[0; 1]</td><td>&Nplus;</td></tr><tr><td>P<sub>1</sub>-W<sub>m</sub>LBM</td><td>P(&alpha;<sub>k<sub>1</sub><sup>m</sup></sub>)</td><td>exp(&minus;&alpha;<sub>k<sub>1</sub><sup>m</sup></sub>)(&alpha;<sub>k<sub>1</sub><sup>m</sup></sub>)<sup>x<sub>1</sub><sup>m</sup></sup>/x<sub>1</sub><sup>m</sup>!</td><td>log &alpha;<sub>k<sub>1</sub><sup>m</sup></sub></td><td>&alpha;<sub>k<sub>1</sub><sup>m</sup></sub></td><td>1</td><td>&Re;</td><td>&Nplus;</td></tr><tr><td>M-W<sub>m</sub>LBM</td><td>M(&alpha;<sub>k<sub>1</sub><sup>m</sup></sub>)</td><td>&Pi;<sub>h=1</sub><sup>H</sup> &alpha;<sub>k<sub>1</sub><sup>h</sup>,h</sub><sup>x<sub>1</sub><sup>h</sup></sup></td><td>[log &alpha;<sub>k<sub>1</sub><sup>m</sup>,1</sub>/&alpha;<sub>k<sub>1</sub><sup>m</sup>,H</sub>, &middot;&middot;&middot; , log &alpha;<sub>k<sub>1</sub><sup>m</sup>,H&minus;1</sub>/&alpha;<sub>k<sub>1</sub><sup>m</sup>,H</sub>, 0]</td><td>&minus;log(1 &minus; &Pi;<sub>h=1</sub><sup>H&minus;1</sup> &alpha;<sub>k<sub>1</sub><sup>h</sup>,h</sub>)</td><td>1</td><td>[0; 1]</td><td>[0; 1]<sup>&otimes;H</sup></td></tr></tbody></table>

also $d_{j\ell}$ such that $\sum_{\ell} d_{j\ell} = 1$. Their matricial representations are respectively $\mathbf{c} = (c_{ik})$, and $\mathbf{d} = (d_{j\ell})$, such that the variational distribution for the clustering is defined as $Q_{(c,d)}(\mathbf{z}, \mathbf{w}) = \prod_{i,k} (c_{ik})^{z_{ik}} \prod_{j,\ell} (d_{j\ell})^{w_{j\ell}}$. Then, by the Jensen inequality a bound $\mathcal{F}(\mathbf{c}, \mathbf{d}; \theta)$ can be defined for the intractable loglikelihood function. After defining this lower bound of the log-likelihood (see [15]), the algorithm maximizes, instead of $L(\theta; \mathbf{x})$, this alternative surrogate function:

$$ \mathcal{F}(\mathbf{c}, \mathbf{d}; \theta) = \sum_{i,j,k,\ell} c_{ik} d_{j\ell} \log \phi(x_{ij}; \theta_{k\ell}) + \sum_{i,k} c_{ik} (p_k - \log c_{ik}) + \sum_{j,\ell} d_{j\ell} (q_\ell - \log d_{j\ell}). $$

Then, it repeats until convergence the two following steps:

* **E-step** The posterior probabilities $\mathbf{e} = (\mathbf{c}, \mathbf{d})$ are found at the current time step $(t)$ by maximizing $\mathcal{F}$ with respect to $c_{ik}$ and $d_{j\ell}$ with the normalizing constraint of summing to one while $\theta = \theta^{(t)}$, hence leading:

$$ (\mathbf{c}^{(t)}, \mathbf{d}^{(t)}) = \operatorname*{argmax}_{\theta} \mathcal{F}(\mathbf{c}, \mathbf{d}; \theta^{(t)}). $$

* **M-step** A temporary value of the parameters is found for the next current time step $(t+1)$ by maximizing $\mathcal{F}$ with respect to $\theta$ while $\mathbf{c} = \mathbf{c}^{(t)}$ and $\mathbf{d} = \mathbf{d}^{(t)}$ from previous step. Say:

$$ \hat{\theta}^{(t+1)} = \operatorname*{argmax}_{\theta} \tilde{Q}_{\text{LBM}}(\theta, \theta^{(t)}). $$

Where, $\tilde{Q}_{\text{LBM}}(\theta, \theta^{(t)}) = \mathcal{F}(\mathbf{c}^{(t)}, \mathbf{d}^{(t)}; \theta)$. The optimization provides the updated solutions $\alpha_{k\ell}^{(t+1)}$, $p_k^{(t+1)}$, $q_\ell^{(t+1)}$ and eventually other parameters in $\theta$, after derivatives equated to zero.

Hence, the parameters are estimated by an iterative way, the BEM algorithm proceeds by an alternated maximization of $\tilde{Q}$ and converges to a final solution which maximizes (locally) the log-likelihood of the latent block model. A hat is added to each parameter or statistics which is estimated and found at the final stage of the inference algorithm.

## 3.2 Multiway version and parameterization ($m \ge 2$)

Co-clustering methods are dedicated to the detection of the main homogeneous sub-structures in a data cloud when the variables can be clustered. At the end, one has a small matrix $\alpha_{k\ell}$ which can be considered as a very synthetic view of the dataset. The new multiway version is defined in this section for changing a 2-dimensional problem into a $m$-dimensional problem. The (co-)clustering tackles the modes/sides of a tensor $\mathbf{x} = x_{i_1, \cdots, i_m} \in \mathbb{R}^{n_1 \times \cdots \times n_m}$ with $m$ the number of modes of the tensor, from now until the end. Thus, one may denote that:

$$ \begin{aligned} I_m &= \{i_1, \cdots, i_m\} \in \{1, \cdots, n_1\} \times \cdots \times \{1, \cdots, n_m\} \\ \mathcal{K}_m &= \{k_1, \cdots, k_m\} \in \{1, \cdots, g_1\} \times \cdots \times \{1, \cdots, g_m\}, \end{aligned} $$

Thus, it is denoted here $i_\ell \in \{1, \cdots, n_\ell\}$ for the mode $\ell$ of the tensor, and $k_\ell \in \{1, \cdots, g_\ell\}$ for the number of clusters $g_\ell$ of the mode $\ell$ of the tensor. One gets for the new random variables/matrices for the clustering along of each mode, $\mathbf{z}_\ell = (z_{i_\ell k_\ell})$ for $\ell \in \{1, \cdots, m\}$. To summarize, from the 2-modes to the $m$-modes, it is replaced by $i_1$, $j$ by $i_2$, and added the new indices $i_3$, $\cdots$, $i_m$ for the additional modes. Such that also $z$ becomes $z_1$, $w$ becomes $z_2$ to which are added $z_3$, $\cdots$, $z_m$ and so on for the other quantites, with $(p_k)_k$ replaced by $(p_{1k_1})_{k_1}$, $q_\ell$ replaced by $(p_{2k_2})_{k_2}$ to which are also added the new priors $(p_{3k_3})_{k_3}$, $\cdots$, $(p_{mk_m})_{k_m}$. Here, the probabilities $p_{\ell k_\ell}$ that a mode of the tensor belongs to the $k_\ell^{th}$ component are aggregated into,

$$ \mathbf{p}_{\ell} = (p_{\ell k_{\ell}})_{1 \le k_{\ell} \le g_{\ell}} = (p_{\ell 1}, \cdots, p_{\ell g_{\ell}}) \text{ for } \ell = 1, \cdots, m. $$

Also the posterior probabilities are $c_{i_1 k_1}^1$ and $c_{i_2 k_2}^2$ replacing the previous $c_{ik}$ and $d_{j\ell}$ such that it is added for the new modes $c_{i_3 k_3}^3$, $\cdots$, $c_{i_m k_m}^m$, with the matrices,

$$ \mathbf{C}^{\ell} = (c_{i_{\ell} k_{\ell}}^{\ell})_{i_{\ell} \in \{1, \cdots, n_{\ell}\}, k_{\ell} \in \{1, \cdots, g_{\ell}\}} \text{ for } \ell = 1, \cdots, m. $$

Similarly, the criteria are rewritten by summing to the new variables and new modes, in order to get the extension to a tensor of any dimensions. The corresponding expressions are presented in
<table><caption>Table 3: Posterior probabilities for the E-step from BEM with the multiway model WLBM.</caption><thead><tr><th>Model</th><th>log <sup>ℓ(t)</sup><sub>c<sub>i<sub>ℓ</sub>k<sub>ℓ</sub></sub></sub> ∝ &sum;<sub>I<sub>m</sub><sup>ℓ</sup>, K<sub>an</sub><sup>ℓ</sup></sub> c<sub>i<sub>ℓ</sub>k<sub>m</sub></sub><sup>ℓ(t)</sup> log φ(x<sub>i<sub>ℓ</sub>m</sub>; α<sub>k<sub>m</sub></sub><sup>(t)</sup>)</th><th>log <sup>ℓ(t)</sup><sub>p<sub>i<sub>ℓ</sub>k<sub>ℓ</sub></sub></sub> ∝ &sum;<sub>I<sub>m</sub><sup>ℓ</sup>, K<sub>an</sub><sup>ℓ</sup></sub> c<sub>i<sub>ℓ</sub>k<sub>m</sub></sub><sup>ℓ(t)</sup> log φ<sub>N</sub>(x<sub>i<sub>ℓ</sub>m</sub>; θ<sub>k<sub>m</sub></sub><sup>(t)</sup>)</th><th>log <sup>ℓ(t)</sup><sub>I<sub>m</sub><sup>ℓ</sup>, K<sub>an</sub><sup>ℓ</sup></sub> ∝ &sum;<sub>I<sub>m</sub><sup>ℓ</sup>, K<sub>an</sub><sup>ℓ</sup></sub> <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">x<sub>k<sub>m</sub>;2</sub><sup>ℓ(t)</sup> - 2x<sub>k<sub>m</sub>;1</sub><sup>ℓ(t)</sup> α<sub>k<sub>m</sub></sub> + c<sub>k<sub>m</sub></sub><sup>ℓ(t)</sup> (α<sub>k<sub>m</sub></sub><sup>2</sup> + 2σ<sub>k<sub>m</sub></sub><sup>2</sup> log σ<sub>k<sub>m</sub></sub><sup>(t)</sup>)</div><div style="padding: 0 5px;">2σ<sub>k<sub>m</sub></sub><sup>2</sup></div></div></td></tr></thead><tbody><tr><td>G-W<sub>m</sub>LBM</td><td>∞ &sum;<sub>I<sub>m</sub><sup>ℓ</sup>, K<sub>an</sub><sup>ℓ</sup></sub> c<sub>i<sub>ℓ</sub>k<sub>m</sub></sub><sup>ℓ(t)</sup> log φ<sub>N</sub>(x<sub>i<sub>ℓ</sub>m</sub>; θ<sub>k<sub>m</sub></sub><sup>(t)</sup>)</td><td>∞ &sum;<sub>K<sub>an</sub><sup>ℓ</sup></sub> <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">x<sub>k<sub>m</sub>;2</sub><sup>ℓ(t)</sup> - 2x<sub>k<sub>m</sub>;1</sub><sup>ℓ(t)</sup> α<sub>k<sub>m</sub></sub> + c<sub>k<sub>m</sub></sub><sup>ℓ(t)</sup> (α<sub>k<sub>m</sub></sub><sup>2</sup> + 2σ<sub>k<sub>m</sub></sub><sup>2</sup> log σ<sub>k<sub>m</sub></sub><sup>(t)</sup>)</div><div style="padding: 0 5px;">2σ<sub>k<sub>m</sub></sub><sup>2</sup></div></div></td><td>∞ &sum;<sub>I<sub>m</sub><sup>ℓ</sup>, K<sub>an</sub><sup>ℓ</sup></sub> <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">x<sub>k<sub>m</sub>;1</sub><sup>ℓ(t)</sup> log(θ<sub>k<sub>m</sub></sub><sup>(t)</sup>) + (c<sub>k<sub>m</sub></sub><sup>ℓ(t)</sup> - x<sub>k<sub>m</sub>;1</sub><sup>ℓ(t)</sup>) log(1 - θ<sub>k<sub>m</sub></sub><sup>(t)</sup>)</div><div style="padding: 0 5px;">2σ<sub>k<sub>m</sub></sub><sup>2</sup></div></div></td></tr><tr><td>B-W<sub>m</sub>LBM</td><td>∞ &sum;<sub>I<sub>m</sub><sup>ℓ</sup>, K<sub>an</sub><sup>ℓ</sup></sub> c<sub>i<sub>ℓ</sub>k<sub>m</sub></sub><sup>ℓ(t)</sup> log φ<sub>B</sub>(x<sub>i<sub>ℓ</sub>m</sub>; θ<sub>k<sub>m</sub></sub><sup>(t)</sup>)</td><td>∞ &sum;<sub>K<sub>an</sub><sup>ℓ</sup></sub> <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">x<sub>k<sub>m</sub>;1</sub><sup>ℓ(t)</sup> log(θ<sub>k<sub>m</sub></sub><sup>(t)</sup>) + (c<sub>k<sub>m</sub></sub><sup>ℓ(t)</sup> - x<sub>k<sub>m</sub>;1</sub><sup>ℓ(t)</sup>) log(1 - θ<sub>k<sub>m</sub></sub><sup>(t)</sup>)</div><div style="padding: 0 5px;">2σ<sub>k<sub>m</sub></sub><sup>2</sup></div></div></td><td>∞ &sum;<sub>I<sub>m</sub><sup>ℓ</sup>, K<sub>an</sub><sup>ℓ</sup></sub> <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">x<sub>k<sub>m</sub>;1</sub><sup>ℓ(t)</sup> log(θ<sub>k<sub>m</sub></sub><sup>(t)</sup>) + (c<sub>k<sub>m</sub></sub><sup>ℓ(t)</sup> - x<sub>k<sub>m</sub>;1</sub><sup>ℓ(t)</sup>) log(1 - θ<sub>k<sub>m</sub></sub><sup>(t)</sup>)</div><div style="padding: 0 5px;">2σ<sub>k<sub>m</sub></sub><sup>2</sup></div></div></td></tr><tr><td>P-W<sub>m</sub>LBM</td><td>∞ &sum;<sub>I<sub>m</sub><sup>ℓ</sup>, K<sub>an</sub><sup>ℓ</sup></sub> c<sub>i<sub>ℓ</sub>k<sub>m</sub></sub><sup>ℓ(t)</sup> log φ<sub>P</sub>(x<sub>i<sub>ℓ</sub>m</sub>; θ<sub>k<sub>m</sub></sub><sup>(t)</sup>)</td><td>∞ &sum;<sub>K<sub>an</sub><sup>ℓ</sup></sub> <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">- ∑<sub>I<sub>m</sub><sup>ℓ</sup></sub> c<sub>i<sub>ℓ</sub>k<sub>m</sub></sub><sup>ℓ(t)</sup> λ<sub>k<sub>m</sub></sub><sup>im</sup> + ∑<sub>I<sub>m</sub><sup>ℓ</sup></sub> c<sub>i<sub>ℓ</sub>k<sub>m</sub></sub><sup>ℓ(t)</sup> x<sub>i<sub>ℓ</sub>m</sub> log(λ<sub>k<sub>m</sub></sub><sup>im</sup>) - x<sub>i<sub>ℓ</sub>m</sub><sup>ℓ(t)</sup></div><div style="padding: 0 5px;">I<sub>m</sub><sup>ℓ</sup></div></div></td><td>∞ &sum;<sub>I<sub>m</sub><sup>ℓ</sup></sub> <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">x<sub>k<sub>m</sub>;1</sub><sup>ℓ(t)</sup> log(λ<sub>k<sub>m</sub></sub><sup>im</sup>) - x<sub>k<sub>m</sub>;1</sub><sup>ℓ(t)</sup></div><div style="padding: 0 5px;">I<sub>m</sub><sup>ℓ</sup></div></div></td></tr><tr><td>M-W<sub>m</sub>LBM</td><td>∞ &sum;<sub>I<sub>m</sub><sup>ℓ</sup>, K<sub>an</sub><sup>ℓ</sup></sub> c<sub>i<sub>ℓ</sub>k<sub>m</sub></sub><sup>ℓ(t)</sup> log φ<sub>M</sub>(x<sub>i<sub>ℓ</sub>m</sub>; θ<sub>k<sub>m</sub></sub><sup>(t)</sup>)</td><td>∞ &sum;<sub>K<sub>an</sub><sup>ℓ</sup></sub> <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">∑<sub>h=1</sub><sup>H</sup> x<sub>k<sub>m</sub>;1,h</sub><sup>ℓ(t)</sup> log α<sub>k<sub>m</sub>;h</sub></div><div style="padding: 0 5px;">H</div></div></td><td>∞ &sum;<sub>K<sub>an</sub><sup>ℓ</sup></sub> <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">x<sub>k<sub>m</sub>;1,h</sub><sup>ℓ(t)</sup> log α<sub>k<sub>m</sub>;h</sub></div><div style="padding: 0 5px;">H</div></div></td></tr></tbody></table>

the following just after before the optimization which introduces new latent vectors to be estimated.

### Model W<sub>m</sub>LBM for m-modes tensor

The new required expressions are for the likelihood $f_{W_mLBM}$, the update formula of the corresponding variational probabilities $c_{i_ℓk_ℓ}^ℓ$, and the maximization for the M-step optimizing w.r.t the parameters $\Theta_{K_m}$, is defined by following the previous 2-dimensional matrix extended into a m-dimensional tensor, with $x_{i_1}^m = x_{i_1 \dots i_m}$, $\Theta_{k_1}^m = \Theta_{k_1 \dots k_m}$, and also $Z_m$ for the indices from $(\mathbf{z}_1, \dots, \mathbf{z}_m)$,

$$f_{W_mLBM}(\mathbf{x}; \mathbf{\theta}) = \sum_{Z_m} \prod_{\ell=1}^{m} \prod_{i_\ell k_\ell} p_{\ell k_\ell}^{z_{i_\ell k_\ell}} \prod_{I_m, K_m} \varphi(x_{i_1}^m; \theta_{k_1}^m)^{z_{i_1 k_1} \dots z_{i_m k_m}}.$$

This leads also to the variational bound which is maximized w.r.t. m clustering variables instead of just 2 ones before, when adding the index $ℓ$ for the posterior and prior of each mode $ℓ$ in order to a better readability, while $ℓ = 1, \dots, m$ and $i_ℓ = 1, \dots, n_ℓ$:

$$c_{i_\ell k_\ell}^{\ell(t)} \propto p_{\ell k_\ell}^{(t)} \exp \left( \sum_{I_m \setminus i_\ell, K_m \setminus k_\ell} \left[ \prod_{\sigma \neq \ell} c_{i_\ell k_\sigma}^{o(t)} \right] \log \varphi(x_{i_1}^m; \theta_{k_1}^{(t)}) \right).$$

Similarly, just the maximum may be preferred for a classifying version of BEM, called BCEM, for binary variables instead of fuzy ones. For maximizing the new function $\mathcal{F}$ with respect to $\Theta$, the objective function becomes:

$$\tilde{Q}_{W_mLBM}(\theta, \theta^{(t)}) = \sum_{I_m, K_m} c_{i_1}^{1(t)} \cdots c_{i_m k_m}^{m(t)} \log \varphi(x_{i_1}^m; \theta_{k_1}^m) + \sum_{\ell=1}^{m} \sum_{i_\ell, k_\ell} c_{i_\ell k_\ell}^{\ell(t)} p_{\ell k_\ell} - \sum_{\ell=1}^{m} \sum_{i_\ell, k_\ell} c_{i_\ell k_\ell}^{\ell(t)} \log c_{i_\ell k_\ell}^{\ell(t)}$$

There was just changed the indices, as the algebra remains the same: this defines our multiway latent block model and its corresponding parameters inference. The new M-step solves for the optimization:

$$\theta^{(t+1)} = \tilde{Q}_{W_mLBM}(\theta, \theta^{(t)}).$$

The expressions for the variational expectation step with the posterior probabilities of any cluster of index $k\ell$ for mode $ℓ$ in

a m-modes tensor for three case of data are given in Table 2. Here it is denoted $i_1^m = i_1 \dots i_m$, $k_1^m = k_1 \dots k_m$ as indices and $I_m^\ell = I_m \setminus i_\ell$, $K_m^\ell = K_m \setminus k_\ell$ as sets. While $x_{k_1}^{\ell(t)} = \sum_{I_m^\ell} c_{i_1 k_1}^{\ell(t)} x_{i_1}^m$, $x_{k_1;2}^{\ell(t)} = \sum_{I_m^\ell} c_{i_1 k_1}^{\ell(t)} x_{i_1}^2$, $x_{k_1;L_1}^{\ell(t)} = \sum_{I_m^\ell} c_{i_1 k_1}^{\ell(t)} \log x_{i_1}^m!$, $c_{k_\ell}^{\ell(t)} = \sum_{i_\ell} c_{i_\ell k_\ell}^{\ell(t)}$, and $c_{i_1 k_1}^{\ell(t)} = \prod_{o \neq \ell} c_{i_\ell k_\sigma}^{o(t)}$, while $c_{k_1}^{\ell(t)} = \prod_{o \neq \ell} c_{k_\sigma}^{o(t)}$. For the solution for the maximization step, let denote, $x_{k_1}^{(t)} = \sum_{I_m} c_{i_1 k_1}^{(t)} x_{i_1}^m$, $x_{k_1;2}^{(t)} = \sum_{I_m} c_{i_1 k_1}^{(t)} x_{i_1}^2$, $c_{k_1}^{(t)} = \prod_{\ell} c_{k_\ell}^{\ell(t)}$, and $c_{k_\ell}^{\ell(t)} = \sum_{i_\ell} c_{i_\ell k_\ell}^{\ell(t)}$, for the mode $ℓ$ of the tensor and $k_\ell \in \{1, \dots, g_\ell\}$, then:

$$\alpha_{k_1}^{(t+1)} = x_{i_1}^{(t)}/c_{k_1}^{(t)} \quad \text{and} \quad p_{\ell k_\ell}^{(t+1)} = c_{k_\ell}^{\ell(t)}/n_\ell \\ \sigma_{k_1}^{(t+1)} = \sqrt{(x_{k_1}^{(t)} - 2x_{k_1;1}^{(t)}\alpha_{k_1}^{(t+1)} + c_{k_1}^{(t)}\alpha_{k_1}^{(t+1)2})/c_{k_1}^{(t)}} \\ \sigma^{(t+1)} = \sqrt{\sum_{K_m} (x_{k_1;2}^{(t)} - 2x_{k_1;1}^{(t)}\alpha_{k_1}^{(t+1)} + c_{k_1}^{(t)}\alpha_{k_1}^{(t+1)2})/\sum_{K_m} c_{k_1}^{(t)}}$$

Here, the updates are for $\varphi_B$ and $\varphi_G$ mostly. Note that some constraints may be better for $\Sigma$ too, as for instance all variances equal to $\sigma_{k_1}^m = \sigma$ or eventually only for clusters from one given mode, leading to a similar but different expression. When constraints are added, multidimensional latent vectors are estimated instead of directly the centers $\mu_{k_1}^m$ or $\alpha_{k_1}^m$ from the blocks.

### Related and extended models

As a generalized model, the multiway latent block model counts as particular cases several former models from the literature. With diverse distribution functions for the cells in the blocks - see the Table 2 revisited from [18] for tensors - it results different expressions for the posterior probabilities in Table 3, estimation step and maximization step with particular regularizations on the parameters such as diverse model were introduced early in the literature when the LBM was created for co-clustering of matrices.

* The more related researches are a multi-view LBM [47], a multivariate LBM [42], and a 3-modes Poisson LBM [41] illustrating the relevance of BEM for tensors.

* Constrained model with a reparameterization [18] of the parameters with reductions of the dimensions. See the Gaus-
<table><caption>Table 4: Models from the likelihood of WLBM with $m = 3$ and $\varphi = \varphi_N$, with parameters $\Sigma = (\sigma_{k_1}^2) \in \mathbb{R}_+^{g_1 \times g_2 \times g_3}$ and $\alpha = (\alpha_{k_1}^2) \in \mathbb{R}^{g_1 \times g_2 \times g_3}$.</caption><thead><tr><th>Model</th><th>Criterion</th></tr></thead><tbody><tr><td>Unconstrained</td><td>Q&#773;<sub>GW<sub>3</sub>LBM</sub>((α<sub>k<sub>1</sub><sup>3</sup></sub>)<sub>k<sub>3</sub></sub>, (σ<sub>k<sub>1</sub><sup>3</sup></sub>)<sub>k<sub>3</sub></sub>, p)|θ<sup>(t)</sup>) = Q&#773;<sub>W<sub>m</sub>LBM</sub>(θ, θ<sup>(t)</sup>)</td></tr><tr><td>Sparse</td><td>Q&#773;<sub>SGW<sub>3</sub>LBM</sub>(θ|θ<sup>(t)</sup>) = Q&#773;<sub>GW<sub>3</sub>LBM</sub>((τ<sup>T</sup>ξ<sub>k<sub>1</sub>k<sub>2</sub>k<sub>3</sub></sub>), (σ<sub>k<sub>1</sub>k<sub>2</sub>k<sub>3</sub></sub>), p)|θ<sup>(t)</sup>)</td></tr><tr><td>Constrained</td><td>Q&#773;<sub>SGW<sub>3</sub>LBM</sub>(θ|θ<sup>(t)</sup>) = Q&#773;<sub>GW<sub>3</sub>LBM</sub>((τ<sup>T</sup>ξ<sub>k<sub>1</sub>k<sub>2</sub>k<sub>3</sub></sub>), (σ<sub>k<sub>1</sub>k<sub>2</sub>k<sub>3</sub></sub>), p)|θ<sup>(t)</sup>) + λF(ξ&#773;<sub>k<sub>1</sub>k<sub>2</sub>k<sub>3</sub></sub>)</td></tr><tr><td>Decomposition</td><td>Q&#773;<sub>PCAW<sub>3</sub>LBM</sub>(θ|θ<sup>(t)</sup>) = Q&#773;<sub>GW<sub>3</sub>LBM</sub>((Λ ⊕ s ξ<sub>k<sub>1</sub></sub> ⊕ s ξ<sub>k<sub>2</sub></sub> ⊕ s ξ<sub>k<sub>3</sub></sub>), (σ<sub>k<sub>1</sub>k<sub>2</sub>k<sub>3</sub></sub>), p)|θ<sup>(t)</sup>)</td></tr></tbody></table>

sian case for instance in Table 4 as explained at the next section.

* Constrained model with a sparsification of the parameters, with constraints on parameters. The diagonal version [48, 49] is written in the proposed model with parameterizing vectors: one latent $\tau$ and ones observed for the cells in the tensor $\xi_{k_1k_2k_3}$ for more generality. The purpose is to induce constraints in a generic way on the contrary to the "sparse" diagonal case for instance which supposes that the vocabulary is separated for each cluster in the case of textual classes, while an additional perspective may be to find automatically non binary constraints or at least binary ones via a relaxation for instance. See in the Gaussian case for instance in Table 4, leading to an equivalence between the diagonal approach and the proposed parameterization with well chosen<fnref n="1" /> binary vectors $\xi_{k_1k_2k_3}$. Note that an alternative, for this equality of the means in the blocks with a regression, is to keep a full model with $\xi_{k_1k_2k_3}$ equal to one of the basis vectors of $\mathbb{R}^{g_1 \times g_2 \times g_3}$, and to add the constraint externally similarly than an usual restricted regression, such that, one may write when keeping only $r$ clusters, as defined via the

matrix $R$:

$$ \operatorname{argmin}_{\tau} \tilde{Q}_{\mathrm{SGW}_{3}\mathrm{LBM}}(\theta|\theta^{(t)}) \text{ s.c. } \Psi R = O_r, $$

leads to the new estimate $\tau^{(t+1)}$ with the constraints but this approach supposes the constraints known on the contrary to the previous model which is more flexible. In both cases, it becomes possible just before the regression to optimize for the vector $\xi_{k_1k_2k_3}$ during the iterative loop by relaxing the hypothesis of binary values. For instance with $\tilde{\xi}_{k_1k_2k_3}$ a vector of softmax functions with common denominators and unknown parameters, while $F(.)$ is an usual function for the norm penalty such as entropy or $L_1$-norm, see Table 4.

In summary, if there are some correspondences of the previous models with the multiway variant, each particular model has its own properties or implementations without embedding latent variables. In the following, the focus is on the low dimensional reductions of a tensor with the proposed clustering model.

# 4 PCA Multiway Latent Block Model

This section discussed a reparameterization with low dimensional matrices of the center in the blocks.

## 4.1 Decomposition/factorization for m-modes tensor

Following the classical tensor algebra [33], the two main decompositions may be in stake via a cell-wise function $\phi(.)$ such as indentity, exponential or sigmoidal ones, with the reduced latent variables when denoting $\lambda_{s_1}^m = \lambda_{s_1 \cdots s_m}$ and $S_m$ for all the indices from $(s_1, \cdots, s_m)$,

$$ \begin{gathered} \xi^l = (\xi_{k_l}^l)_{k_l \in \mathbb{R}^{g_l \times r_l}} \text{ for } l = 1, \cdots, m \\ \Lambda = \begin{cases} (\lambda_s)_s \in \mathbb{R}^{r_1 \times 1}, & \text{CP-D} \\ (\lambda_{s_1}^m)_{s_1 \in \mathbb{R}^{r_1 \times \cdots \times r_m}}, & \text{TK-D.} \end{cases} \end{gathered} $$

Thus, with eventually $r_m = r$, as the tensor to reduce is now $\alpha$ instead of $\mathbf{x}$, one may write the decomposition with $\Lambda, \xi^{(1)}, \cdots, \xi^{(m)}$ unknown and:

$$ \alpha_{k_1 \cdots k_m} \approx \phi (\Lambda \oplus_s \xi_{k_1}^1 \oplus_s \cdots \oplus_s \xi_{k_m}^m). $$

Where:

$$ \Lambda \oplus_s \xi_{k_1}^1 \oplus_s \cdots \oplus_s \xi_{k_m}^m = \begin{cases} \sum_s \lambda_s \xi_{k_1 s}^1 \xi_{k_2 s}^2 \cdots \xi_{k_m s}^m, & \text{CP-D} \\ \sum_s \lambda_{s_1}^m \xi_{k_1 s_1}^1 \xi_{k_2 s_2}^2 \cdots \xi_{k_m s_m}^m, & \text{TK-D.} \end{cases} $$

$$ \begin{aligned} \tau^{(t+1)} &= \left( \sum_{\mathcal{K}_3} \frac{c_{k_1}^{(t)}}{\sigma_{k_1}^{(t)2}} \xi_{k_1}^T \xi_{k_1}^T \right)^{-1} \left( \sum_{\mathcal{K}_3} \frac{x_{k_1}^{(t)}}{\sigma_{k_1}^{(t)2}} \xi_{k_1}^T \right) \\ &= \left[ \frac{h_{111}^x}{h_{111}} \frac{h_{222}^x}{h_{222}} \cdots \frac{h_{ggg}^x}{h_{ggg}} \frac{\sum_{\mathcal{K}_{\text{off}}} h_{k_1k_2k_3}^x}{\sum_{\mathcal{K}_{\text{off}}} h_{k_1k_2k_3}} \right]^T \\ &= \left[ \alpha_{111}^{(t+1)} \quad \alpha_{222}^{(t+1)} \quad \cdots \quad \alpha_{ggg}^{(t+1)} \quad \alpha_{ooo}^{(t+1)} \right]^T. \end{aligned} $$

$$ \alpha_{ooo}^{(t+1)} = \left( x_{\bullet\bullet} - \sum_{l_3, \mathcal{K}_{\text{dia}}} c_{i_1}^{(t)} x_{l_3}^{i_3} \right) / \sum_{\mathcal{K}_{\text{off}}} c_{i_1}^{(t)}. $$

<footnote n="1">¹For $m = 3$ and $g_1 = g_2 = g_3 = g$, let define the vectors $\xi_{k_1k_2k_3}$ with $\xi_{kkk}$ equal to the basis vector of $\mathbb{R}^g$ completed by a $g+1$ component equal to zero, and with $\xi_{k_1k_2k_3}$ equal to a vector with $g+1$ components where only the last one $g+1$ component equal to one and equal to zero otherwise: $[1, 0, \dots, 0, 0], [0, 1, \dots, 0, 0], \dots, [0, 0, \dots, 1, 0], [0, 0, \dots, 0, 1]$ such that the solution for $\tau$, when denoting $h_{k_1k_2k_3} = c_{k_1}^{(1)} c_{k_2}^{(2)} c_{k_3}^{(3)} / \sigma_{k_1k_2k_3}^{(t+1)}$ and $h_{k_1k_2k_3}^x = \sum_{l,m} c_{i_1k_1}^{(1)} c_{i_2k_2}^{(2)} c_{i_3k_3}^{(3)} / \sigma_{k_1k_2k_3}^{(t+1)}$ becomes the same solution than for the diagonal case:</footnote>

<footnote>Then the different binary vectors $\xi_{k_1k_2k_3}$ selects one of these components. Hence, with $\mathcal{K}_{\text{dia}} = \{k_1, k_2, k_3\}_{k_1=k_2=k_3}$, $\mathcal{K}_{\text{off}} = \{k_1, k_2, k_3\}_{k_1 \neq k_2 \neq k_3}$ and $x_{\bullet\bullet} = \sum_{l} x_{l1}^i$, as $\sum_{\mathcal{K}_3} \alpha_{k_1}^{(t+1)} c_{k_1}^{(t)} = x_{\bullet\bullet}$, the equality constraint induces that the common value is equal to:</footnote>
This expression extends the scalar product between vectors [17] (as previously proposed) to a tensor decomposition with latent matrices $\xi_{k_\ell s_\ell}^\ell$ to be estimated for all $\ell$ and $s_\ell$. The family of new models is related to a CP decomposition and for more generality to a Tucker decomposition from theory in tensor algebra. The difference here is that this is the summarizing tensor from the latent block model which is decompoonded as illustrated for $m = 3$ in Figure 2 where the sum is from the posterior probabilities.

<figure><img src="image_2.png" /><figcaption>Figure 2: (CP/TK) decomposition + WLBM for a 3-modes tensor.</figcaption></figure>

This is useful in many domains such as visualization, reduction, and factorization or even predictive methods because of the reducing matrices along each mode. To study further the parameterization, four variants may be implemented as presented in the Table 5 below, with $\xi_{k_2^m} = \xi_{k_2 \cdots k_m}$.

<table><caption>Table 5: Latent variables for tensor parameters α in WLBM.</caption><thead><tr><th>Name</th><th>α<sub>k_1</sub><sup>3</sup></th><th>constr. ξ<sub>k_1</sub></th><th>constr. ξ<sub>k_2</sub> ... ξ<sub>k_m</sub></th></tr></thead><tbody><tr><td>D-F-m</td><td>Λ ⊕<sub>s</sub> ξ<sub>k_1</sub> ⊕<sub>s</sub> ... ⊕<sub>s</sub> ξ<sub>k_m</sub></td><td>no</td><td>no</td></tr><tr><td>D-F-1</td><td>Λ ⊕<sub>s</sub> ξ<sub>k_1</sub> ⊕<sub>s</sub> ξ<sub>k_2</sub><sup>m</sup></td><td>grid</td><td>no</td></tr><tr><td>PCA-m</td><td>Λ ⊕<sub>s</sub> ξ<sub>k_1</sub> ⊕<sub>s</sub> ... ⊕<sub>s</sub> ξ<sub>k_m</sub></td><td>no</td><td>no</td></tr><tr><td>PCA-1</td><td>Λ ⊕<sub>s</sub> ξ<sub>k_1</sub> ⊕<sub>s</sub> ξ<sub>k_2</sub><sup>m</sup></td><td>grid</td><td>no</td></tr></tbody></table>

<table><caption>Table 6: Example of size for the reduced views of the tensors</caption><thead><tr><th></th><th>Sizes for <b>x</b></th><th>Sizes for α</th><th>Ratio</th></tr></thead><tbody><tr><td>CP-D</td><td>Σ<sub>ℓ</sub> n<sub>ℓ</sub> r</td><td>Σ<sub>ℓ</sub> g<sub>ℓ</sub> r</td><td>= <sup>g<sub>0</sub></sup>/<sub>n<sub>0</sub></sub></td></tr><tr><td>TK-D</td><td>Σ<sub>ℓ</sub> n<sub>ℓ</sub> r<sub>ℓ</sub> + r<sub>1</sub> r<sub>2</sub> ... r<sub>ℓ</sub></td><td>Σ<sub>ℓ</sub> g<sub>ℓ</sub> r<sub>ℓ</sub> + r<sub>1</sub> r<sub>2</sub> ... r<sub>ℓ</sub></td><td>≈ <sup>g<sub>0</sub> + r<sub>0</sub><sup>m-1</sup></sup>/<sub>n<sub>0</sub> + r<sub>0</sub><sup>m-1</sup></sub></td></tr></tbody></table>

The four cases are related to an approximated decompositions CP-D or TK-D for D-F-m and D-F-1 while PCA-m and PCA-1 are even more closely for an approximated PCA decomposition. The parameters are not any more completely free, for this version of principal component analysis, as instead this is a regular grid which allows to have an ordering on the plane [17, 16, 18] as wanted for visualization purposes.

## 4.2 Post-treatments

After the inference of the parameters, diverse approaches are available for data analysis from the obtained final latent variables.

### Reconstruction of cells

From the decomposition of the tensor to the dimensions $r_1 \times \cdots \times r_m$, the reconstruction of the cells remains as:

$$ \begin{aligned} \hat{x}_{i_1}^m &= \sum_{\mathcal{K}_m} \hat{c}_{i_1 k_1}^1 \hat{c}_{i_2 k_2}^2 \cdots \hat{c}_{i_m k_m}^m \beta_{i_1}^{jm} \hat{\alpha}_{k_1 k_2 \cdots k_m} \\ &\approx \sum_{\mathcal{K}_m} \hat{c}_{i_1 k_1}^1 \hat{c}_{i_2 k_2}^2 \cdots \hat{c}_{i_m k_m}^m \beta_{i_1}^{jm} \phi \left( \hat{\Lambda} \oplus_s \hat{\xi}_{k_1}^1 \oplus_s \cdots \oplus_s \hat{\xi}_{k_m}^m \right). \end{aligned} $$

This expression is useful for missing data for filling values in unknown cells in an incomplete data tensor or also for running a cross-validation in model selection for instance.

### Reconstruction of reductions

When $\hat{C}^\ell = (\hat{c}_{i_\ell k_\ell})_{i_\ell k_\ell} \in \mathbb{R}^{n_\ell \times g_\ell}$ is the matrix aggregating the final posterior probabilities, the reduction of the matrix $\alpha$ induces also a reduction for each mode of the data tensor, with:

$$ \hat{\xi}_{\text{WLBM}}^{\mathbf{x}^{(\ell)}} \stackrel{\text{def}}{=} \hat{C}^{\ell} \hat{\xi}^{\ell}. $$

Let denote $n_0$ and $r_0$ the average values of $n_\ell$ and $r_\ell$ in Table 6 where is recalled the number of parameters to be estimated with the re-parameterized WLBM for factorization, the sizes of the latent variables are dramatically reduced in comparison with a full decomposition of the data tensor $\mathbf{x}$ without clustering. Say for each mode $\ell$, $\hat{\xi}^\ell \in \mathbb{R}^{g_\ell \times r_\ell}$ instead of $\hat{\xi}^{\mathbf{x}^{(\ell)}} \in \mathbb{R}^{n_\ell \times r_\ell}$ previously, while the ranks of the involved tensors w.r.t. $r_\ell$ must be discussed. Moreover, for big tensors, the matrices $\hat{c}^\ell$ may be kept sparse by preferring the maximum a posteriori for hard clustering with binary values instead of soft clustering with probability values.

### Direct alternative

One may prefer to focus on the clustering from WLBM without the proposed parameterization via latent vectors. Hence, it may be performed a decomposition $\tilde{\Lambda}, \tilde{\xi}_{k_1}, \cdots, \tilde{\xi}_{k_m}$ only after the inference of the co-clustering model from the final estimates $\tilde{\alpha}_{k_1 \cdots k_m}$. This alternative approach may be written:

$$ \phi^{-1}(\tilde{\alpha}_{k_1, \cdots, k_m}) \approx \tilde{\Lambda} \oplus_s \tilde{\xi}_{k_1}^1 \oplus_s \cdots \oplus_s \tilde{\xi}_{k_m}^m. $$

Note that the purpose may be $\hat{\xi}_{\text{WLBM}}^{\mathbf{x}^{(\ell)}} \approx \hat{\xi}^{\mathbf{x}^{(\ell)}}$ but the dimensions of the tensor $\alpha_{k_1}^m$ are dramatically smaller than the clustered observed tensor $\mathbf{x}$. Hence the obtained reduction $\hat{\xi}_{\text{WLBM}}^{\mathbf{x}^{(\ell)}}$ may be of smaller dimensions than $\hat{\xi}^{\mathbf{x}^{(\ell)}}$ and a quantization required by increasing $g_\ell$. But including the decomposition within the model offers different and eventually better modeling for data analysis. This may happen for varying variances among blocks or overlapping with fuzzy reconstructions.

## 4.3 Algorithm for $W_3LBM+(CP/TK)-D$

For the estimation of the parameters, several algorithms may be tested for the decomposition. First a Newton-Raphson procedure
(NRP) which was successful in [18] may be implemented. In all cases, a BEM for WLBM is implemented in order to increase the log-likelihood as only the parameters have changed, hence the variational framework remains. First, the posterior probabilities are estimated at the E-step. Second, the maximization is performed at the M-step with a (non)linear algorithm for the parameter optimization. The two steps, E and M, are iterated until a stopping criterion is reached with $||$|| a matricial norm such that the Frobenius one, and $\epsilon$ a small positive scalar. Many algorithms are possible for the M-step, but not with all same results. For the algorithm Algo NRP, it is computed the Hessian matrix $\tilde{H}_{\xi}^{(t)}$ and the gradient vector $\nabla \tilde{Q}_{\xi}^{(t)}$ because it is based on a quadratic approximation of the criterion in a vicinity of the current solution at step $t$,

$$ \xi^{(t+1)} = \xi^{(t)} - [\tilde{H}_{\xi}^{(t)}]^{-1} \nabla Q_{\xi}^{(t)}. $$

Considering the bad local minima in some cases, a constraint may be added here, with a sequential quadratic procedure (SQP). The chosen constraint is the weighted sum of the blocks centers which is equal to the sum of the cells of the data tensors, as explained above at the paragraph on sparse latent block model. One may prefer an algorithm Algo SQP when it is defined the constraint $h_{\xi} = h(\xi)$ which enters the optimization criterion with a Lagrangian $\eta$ to be optimized simultaneously via its derivatives. Here, the maximization step is kept to change. The difference with the first algorithm is to insure that the constraint is reached because the nonlinear optimization is prone to local minimal this is a more secure procedure. Experiments must validate this parameterization in order to check the convergence of the algorithms and the final quality of the eventual resulting reductions. Note that in the case of binary posterior probabilities, the decompositions are computed in each block separately. To finish, extensions to non Gaussian model are available by optimizing a bound of the criterion such as found in [50, 51] for binary regression models.

## 4.4 Related models and methods

The latent block model is closely related to a $\chi^2$ criterion also found in correspondence analysis (CA) [52, 53, 26], such that the two methods share several similarities as explained next paragraphs.

* First for $m=2$, co-clustering with CA can be seen as a co-factorization of the rows and columns before clustering of the rows and the column from the components of the decomposition. Following [54] with the multiway notation herein, first one computes the row margins $\mu_1$ and column margins $\mu_2$ to construct the matrix of pseudo-residuals $Z_x$. It is denoted with cells $z_{i_1i_2} = \frac{x_{i_1i_2}/x_{\bullet\bullet} - \mu_{i_1}\mu_{i_2}}{\sqrt{\mu_{i_1}\mu_{i_2}}}$ or matricially $Z_x$, while $D_{\mu_1} = \text{diag}(\mu_{i_11}, \dots, \mu_{i_1n_1})$ and $D_{\mu_2} = \text{diag}(\mu_{i_21}, \dots, \mu_{i_2n_2})$ with the vectors $\mu_1$ and $\mu_2$ for their diagonals. The method computes the r-rank truncated SVD of $Z_x$, such that:

$$ Z_x = D_{\mu_1}^{-1/2} \left( \frac{\mathbf{x}}{x_{\bullet\bullet}} - \mu_1 \mu_2^T \right) D_{\mathbf{c}}^{-1/2} \approx \tilde{\xi}_{\text{CA}}^{\mathbf{x}(1)} D_{\mathbf{r}} \tilde{\xi}_{\text{CA}}^{\mathbf{x}(2)T}. $$

Thus, its standard row coordinates and standard column coordinates are computed:

$$ \xi_{\text{CA}}^{\mathbf{x}(1)} = D_{\mu_1}^{-1/2} \tilde{\xi}_{\text{CA}}^{\mathbf{x}(1)} \quad \text{and} \quad \xi_{\text{CA}}^{\mathbf{x}(2)} = D_{\mu_2}^{-1/2} \tilde{\xi}_{\text{CA}}^{\mathbf{x}(2)}. $$

This induces:

$$ \begin{aligned} \xi_{\text{CA}}^{\mathbf{x}(1)} D_{\mu_1} \xi_{\text{CA}}^{\mathbf{x}(2)T} &\approx D_{\mu_1}^{-1/2} Z_x D_{\mu_2}^{-1/2} \\ &= D_{\mu_1}^{-1} \left( \frac{\mathbf{x}}{x_{\bullet\bullet}} - \mu_1 \mu_2^T \right) D_{\mu_2}^{-1}. \end{aligned} $$

Finally, it follows a matricial expression (for the one in the second point below) by changing the number of modes from two to three.

* Secondly for $m=3$, models for (multiple) correspondence analysis (M)CA are related too, with the right parameterization, see for instance [54] and rewritten for $m=3$ as in [55] with tensor count data:

$$ \begin{aligned} \tilde{x}_{i_1}^{3} &= \frac{x_{i_1}^{3} / x_{\bullet\bullet} - \beta_{i_1}^{3}}{\sqrt{\beta_{i_1}^{3}}} \\ \text{svd} \quad \Rightarrow \quad \frac{x_{i_1}^{3}}{x_{\bullet\bullet}} &\approx \beta_{i_1}^{3} \left( 1 + \sum_{s} \lambda_{s} \xi_{i_1 s}^{\mathbf{x}(1)} \xi_{i_2 s}^{\mathbf{x}(2)} \xi_{i_3 s}^{\mathbf{x}(3)} \right). \end{aligned} $$

This expression is a direct extension to three modes from CA defined for contingency tables $\mathbf{x}$, as recalled in [54]. The equality recalls the probabilistic models herein with the Poisson and Bernoulli distributions. For count data, with $e^u \approx 1+u$, one may write for $x_{i_1}^3$ in the block $z_{i_1i_2}z_{i_2i_3}$, say $\frac{x_{i_1}^3}{x_{\bullet\bullet}} \approx \beta_{i_1}^3 e^{\sum_s \lambda_s \xi_{i_1 s}^1 \xi_{i_2 s}^2 \xi_{i_3 s}^3}$. This is the main purpose herein, to approximate or improve the decomposition of the tensor $\mathbf{x}$ with latent variables from blocks which are estimated jointly or before.

* Thirdly, for $m=3$, the tensor decomposition is also closely related to linear-bilinear models. See for instance [54] as written for $m=3$ with tensor data:

$$ x_{i_1}^3 \approx \mu_{i_1} + \mu_{i_2} + \mu_{i_3} + \sum_{s} \lambda_s \xi_{i_1 s}^{\mathbf{x}(1)} \xi_{i_2 s}^{\mathbf{x}(2)} \xi_{i_3 s}^{\mathbf{x}(3)}. $$

Note that in the expression just above some additional terms for the effects from the sides of the tensors are added while they are ignored herein.

Finally, main other related models are tensor decomposition with exponential family such as often called generalized[56] in the literature. Next the clustering method with CA is compared for (co-)clustering of texts with state of art before experimental projections with multiway latent block model.

## 5 Numerical experiments

In this section, the proposed multiway latent block model and its parameterization are illustrated and validated with real data for a reduction associated to a co-clustering by comparing to current state of art. Correspondence analysis is also involved in order to validate the reduction before clustering when compared with more usual co-clustering methods.
<table><caption>Table 7: Numerical results for (co-)clustering. (*) Only best score kept, while for CA one unique score available from the clustered coordinates.</caption><thead><tr><th></th><th>NMI<br>PLBM</th><th>NMI<br>SPLBM</th><th>NMI<br>CA<sub>km</sub></th><th>NMI<br>CA<sub>gmm</sub></th><th>NMI<sup>(*)</sup><br>CA<sub>lbm</sub></th><th>NMI<br>(PB−)τCC</th><th>ARI<br>PLBM</th><th>ARI<br>SPLBM</th><th>ARI<br>CA<sub>km</sub></th><th>ARI<br>CA<sub>gmm</sub></th><th>ARI<sup>(*)</sup><br>CA<sub>lbm</sub></th></tr></thead><tbody><tr><td>C3SMALL</td><td></td><td></td><td>0.54</td><td>0.74</td><td>(0.79)</td><td></td><td></td><td></td><td>0.41</td><td>0.80</td><td>(0.85)</td></tr><tr><td>TR41</td><td></td><td></td><td>0.57</td><td>0.71</td><td></td><td>0.49 ± 5%</td><td></td><td></td><td>0.40</td><td>0.62</td><td></td></tr><tr><td>CSTR</td><td>0.67 ± 3%</td><td>0.84 ± 5%</td><td>0.72</td><td>0.74</td><td>(0.76)</td><td>0.75 ± 2%</td><td>0.59 ± 3%</td><td>0.73 ± 4%</td><td>0.73</td><td>0.80</td><td>(0.77)</td></tr><tr><td>CLASSIC3</td><td>0.93 ± 1%</td><td>0.93 ± 1%</td><td>0.93</td><td>0.95</td><td></td><td>0.90 ± 2%</td><td>0.94 ± 1%</td><td>0.94 ± 1%</td><td>0.96</td><td>0.97</td><td></td></tr><tr><td>SPORTS</td><td>0.55 ± 3%</td><td>0.56 ± 6%</td><td>0.56</td><td>0.69</td><td></td><td>0.59 ± 4%</td><td>0.42 ± 3%</td><td>0.53 ± 9%</td><td>0.41</td><td>0.58</td><td></td></tr><tr><td>TDT2</td><td>0.75 ± 1%</td><td>0.80 ± 2%</td><td>0.80</td><td>0.76</td><td></td><td></td><td>0.48 ± 3%</td><td>0.79 ± 4%</td><td>0.60</td><td>0.56</td><td></td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Figure 3: First principal map (from coordinates 1 and 2) from CA with from left to right: C3SMALL, TR41, CSTR, CLASSIC3, SPORTS, TDT2.</figcaption></figure>

## 5.1 Experimental settings

The datasets are described in Table 8 with the number of rows, columns, number of clusters and sparsity. One may perform a

<table><caption>Table 8: Several real datasets.</caption><thead><tr><th>name</th><th>n</th><th>d</th><th>g<sub>1</sub></th><th>g<sub>2</sub></th><th>%zero</th><th>x<sub>••</sub></th></tr></thead><tbody><tr><td>C3SMALL</td><td>450</td><td>171</td><td>3</td><td>3</td><td>6.63</td><td>8976</td></tr><tr><td>TR41</td><td>878</td><td>7453</td><td>10</td><td>10</td><td>2.61</td><td>355524</td></tr><tr><td>CSTR</td><td>475</td><td>1000</td><td>4</td><td>4</td><td>3.40</td><td>1754375</td></tr><tr><td>CLASSIC3</td><td>3891</td><td>4303</td><td>3</td><td>3</td><td>1.05</td><td>256348</td></tr><tr><td>SPORTS</td><td>8580</td><td>14870</td><td>7</td><td>7</td><td>0.86</td><td>1722405</td></tr><tr><td>TDT2</td><td>9394</td><td>36771</td><td>30</td><td>30</td><td>0.35</td><td>1725683</td></tr></tbody></table>

clustering from the reduced views of the rows and columns, say $\xi^{x(1)}$ and $\xi^{x(2)}$ from correspondence analysis. Two approaches are available as follows.

**(a) Two sides clustering:** The pseudo algorithm is with nbc for the number of clusters while KM stands for K-means as the clustering algorithm for instance,

$$1. \text{ CA Reduction : } (\xi_{\text{CA}}^{x(1)}, \xi_{\text{CA}}^{x(2)}) = \text{CA}(\mathbf{x}, \text{rank} = r),$$

$$2. \text{ CA Clustering : } \hat{\mathbf{z}}_{\ell}^{\text{CA}} = \text{KM}(\xi_{\text{CA}}^{x(\ell)}, \text{nbc} = g_{\ell}), \ell = 1, 2.$$

Here the approach leads to co-clusters but the two clusterings procedures perform separately: alternative algorithms are simultaneous reduction and clustering, see for instance [57] which is out of the scope.

**(b) Tensor clustering:** For a 3-modes tensor with three sides, the two first ones may correspond to a matrix, such as the third

mode may be informations from the cells of this matrix. For co-clustering the cells $x_{ij}$ from a matrix, a way is to have for the third mode, the information from the row $i$ and column $j$, say both from $x_i = [x_{i1}, x_{i2}, \cdots, x_{id}]^T \in \mathbb{R}^n$ and $x^j = [x_{1j}, x_{2j}, \cdots, x_{nj}]^T \in \mathbb{R}^d$. The corresponding tensor is then available for clustering by comparing the reduced views of the rows and columns simultaneously.

$$1. \text{ Tensorization : } [x_{ij1}, x_{ij2}, \cdots, x_{ijt_r}] = [\Phi(x_i)^T, \Psi(x^j)^T],$$

$$2. \text{ T Clustering : } (\hat{\mathbf{z}}_{\ell})^{\text{WLBM}} = \text{GW}_3\text{LBM}(T, \text{nbc} = ((g_{\ell})_{\ell}).)$$

Here $\hat{x}_i \in \mathbb{R}^{n_r}$ and $\hat{x}^j \in \mathbb{R}^{d_r}$ denote reductions from a projection method such that an autoencoder or pca. When denoting $t_r = n_r + d_r$, the tensor is defined from these low dimensional reductions $\Phi(\cdot)$ and $\Psi(\cdot)$, for the rows and columns of the data matrix when denoting $\hat{x}_i = \Phi(x_i)$ and $\hat{x}^j = \Psi(x^j)$. In the tensor $T$, each tube corresponding to a cell from the contingency table appends the reduction of the corresponding row and column.

The proposed results are computed with Python code when the running time is under a given maximum. Next correspondence analysis is tested as a reducing method for the clustering when defining the tensor for several datasets while standardizing the variance to one for each slice corresponding to a previous row in a contingency table after keeping the $g_1$ first main CA components, assuming that the first side is for the $n_1$ texts and second one for the $n_2$ terms.

## 5.2 Numerical results for clustering

A visualization from first map of CA is in Figure 3 where the separation degree is partially retrieved as other dimensions are informative too. The numerical results (with NMI and ARI for quality criteria) are in Table 7 for several methods: PLBM and SPLBM from [58, 48], $CA_{km}$ for CA+K-means, $CA_{gmm}$ for CA+GMM
and CA<sub>lbm</sub> for CA+GW<sub>3</sub>LBM after tensorization and (PB−)τCC from [59] (page 2171, Table 3). For each method from the literature, the published numbers were extracted such that the obtained Table 7 is not completely filled. As model selection is not in stake for CA<sub>lbm</sub>, the best model with best clustering criteria was kept when were computed the quality criteria. The clustering with CA followed by a clustering algorithm or with a tensorization from CA followed by WLBM are both able to reach a better score for several datasets but selection of the final best obtained model is wanted as a perspective.

## 5.3 Numerical results for a PCA after CA+CP-D

The tensorization is able to improve the projection of CA according to one indicator, called Calinski-Harabasz score (CH) [60] which computes the ratio of the inter versus intra variances, but also the silhouette score (SL) [61] related to averages dissimilarities normalized by their minimum values. For the three datasets for which the results were obtained with the available computer, this leads to compute the indicators for the first plane from CA and the first plane from PCA with the CP-decomposition of the tensorized original data matrix (CA+CP+PCA), see the graphics in Figure 4 for comparison with previous maps. With higher

<figure><img src="image_4.png" /><figcaption>Figure 4: First map CA+CP+PCA for CSTR, C3SMALL, CLASSIC3.</figcaption></figure>

<table><caption>Table 9: Quality criteria for CA map before and after tensorization.</caption><thead><tr><th></th><th>SL<sub>CA</sub></th><th>SL<sub>CA+CP+PCA</sub></th><th>CH<sub>CA</sub></th><th>CH<sub>CA+CP+PCA</sub></th></tr></thead><tbody><tr><th>CSTR</th><td>0.44</td><td>0.33</td><td>263.69</td><td><strong>904.09</strong></td></tr><tr><th>C3SMALL</th><td>0.39</td><td>0.46</td><td>226.65</td><td><strong>823.85</strong></td></tr><tr><th>CLASSIC3</th><td>0.64</td><td>0.66</td><td>4106.90</td><td><strong>18104.75</strong></td></tr></tbody></table>

values for better projection, it was obtained the Table 9 which validates further the tensorization for these datasets after the previous improved error rates. The projections in Figure 4 induce a nearly spherical shape for a part of the clusters on the projection. It is able to separate better several clusters in comparison to the original projection from CA without tensorization. Visually, some cluster projections are reduced or more away but reducing the overlap may ask for further improvement. In summary, the numerical indicators for the clustering for several datasets suggest that in some cases, CA is able to reduce well the columns space before the (3-modes co-)clustering, leading to a better (co-)clustering than other methods without CA. Despite its rare study in machine learning, correspondence analysis leads to promising results to be improved as a perspective with for instance a removal

of outliers, a better clustering method or even the improvement of its definition.

# 6 Discussion and perspectives

This document presents a multiway generative model suitable for co-clustering and reduction of tensor data in a full probabilistic setting by a parameterization of a new extended latent block model, named WLBM. The need for a clustering along all the modes comes from the resulting quantization which partitions each mode in order to reduce its size accordingly into a low dimensional space. Parameterization of latent block model is very powerfull, as it is able to lead to visualization of large datasets while also able to help for a generic sparse model by selection of the blocks to have common distributions. The proposed model herein is related to classical tensor analysis in order to extend latent variables models for matrices with more generality. The main contributions are: a) Multiway block latent model (WLBM) for a tensor, b) Parameterization of WLBM for factorization or PCA of tensor, c) Parameterization of WLBM for a generic sparse model, d) Implementation of WLBM for gaussian or bernoullian or poissonian cells, e) Experiments with new models with tensors and tensorized matrices for comparison with state of art and with correspondence analysis. Perspectives may be alternative reductions or decompositions such that block-term one [62] for optimal results, and a further analysis of the tensorization with WLBM for deep CA [63] and deep co-clustering [64], with automatic selection of the best model.

# References

[1] A. K. Jain, “Data clustering: 50 years beyond K-means,” *Pattern Recognition Letters*, vol. 31, no. 8, pp. 651–666, Jun. 2010.

[2] J. A. Hartigan, *Clustering Algorithms*. New York: Wiley, 1975.

[3] G. Govaert, “Classification croisée,” Thèse d'État, Université Paris 6, France, 1983.

[4] ——, “Simultaneous clustering of rows and columns,” *Control and Cybernetics*, vol. 24, no. 4, pp. 437–458, 1995.

[5] I. S. Dhillon, S. Mallela, and D. S. Modha, “Information-theoretic co-clustering,” in *Proceedings of The Ninth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD-2003)*, 2003, pp. 89–98.

[6] A. Banerjee, I. Dhillon, J. Ghosh, S. Merugu, and D. S. Modha, “A generalized maximum entropy approach to bregman co-clustering and matrix approximation,” *J. Mach. Learn. Res.*, vol. 8, pp. 1919–1986, Dec. 2007.

[7] I. V. Mechelen, H. H. Bock, and P. D. Boeck, “Two-mode clustering methods: a structured overview,” *Statistical methods in medical research*, vol. 13, no. 5, pp. 363–394, Oct. 2004.

[8] D. Agarwal and S. Merugu, “Predictive discrete latent factor models for large scale dyadic data,” in *KDD*. ACM, 2007, pp. 26–35.

[9] T. Hofmann, “Probabilistic latent semantic analysis,” *SIGIR’99*, pp. 50–57, 1999.

[10] C. Ding, T. Li, and W. Peng, “On the equivalence between non-negative matrix factorization and probabilistic latent semantic indexing,” *Computational Statistics and Data Analysis*, vol. 52, pp. 3913–3927, 2008.
[11] D. Cai, X. Wang, and X. He, "Probabilistic dyadic data analysis with local and global consistency," in *ICML*, 2009.

[12] G. Govaert and M. Nadif, "Clustering with block mixture models," *Pattern Recognition*, vol. 36, no. 2, pp. 463-473, 2003.

[13] —, "Latent block model for contingency table," *Communications in Statistics-theory and Methods*, vol. 39, pp. 416-425, 2010.

[14] M. Nadif and G. Govaert, "Model-based co-clustering for continuous data," in *ICMLA, IEEE Computer Society*, 2010, pp. 175-180.

[15] G. Govaert and M. Nadif, "An EM algorithm for the block mixture model," *IEEE Trans. Pattern Anal. Mach. Intell.*, vol. 27, no. 4, pp. 643-647, 2005.

[16] R. Priam, M. Nadif, and G. Govaert, "Topographic bernoulli block mixture mapping for binary tables," *Pattern Analysis and Applications*, vol. 17, pp. 839-847, 2014.

[17] R. Priam and M. Nadif, "Data visualization via latent variables and mixture models: a brief survey," *Pattern Analysis and Applications*, vol. 19, no. 3, pp. 807-819, Aug 2016.

[18] R. Priam, M. Nadif, and G. Govaert, "Generalized topographic block model," *Neurocomputing*, vol. 173, pp. 442-449, 2016.

[19] E. Battaglia and R. G. Pensa, "A parameter-less algorithm for tensor co-clustering," *Machine Learning*, vol. 112, no. 2, pp. 385-427, Feb 2023.

[20] J. Kossaifi, Y. Panagakis, A. Anandkumar, and M. Pantic, "Tensorly: Tensor learning in python," *Journal of Machine Learning Research*, vol. 20, no. 26, pp. 1-6, 2019.

[21] F. Nattino, O. Ku, M. W. Grootes, E. Izquierdo-Verdiguier, S. Girgin, and R. Zurita-Mill, "CGC: a Scalable Python Package for Co- and Tri-Clustering of Geodata Cubes," *Journal of Open Source Software*, vol. 7, no. 72, 2022.

[22] M. Bhattarai, N. kharat, I. Boureima, E. Skau, B. Nebgen, H. Djidjev, S. Rajopadhye, J. P. Smith, and B. Alexandrov, "Distributed non-negative rescal with automatic model selection for exascale data," *Journal of Parallel and Distributed Computing*, vol. 179, 2023.

[23] K. Tsuyuzaki and I. Nikaido, "ntensor: An r package for non-negative matrix/tensor decomposition," *Journal of Open Source Software*, vol. 8, no. 84, 2023.

[24] A. P. Singh. and G. J. Gordon, "A unified view of matrix factorization models," in *ECML PKDD, LNAI 5212*, 2008, pp. 358-373.

[25] I. Jolliffe, *Principal Component Analysis*. Springer Verlag, 2002.

[26] L. Lebart, A. Morineau, and K. Warwick, *Multivariate Descriptive Statistical Analysis*. J. Wiley, 1984.

[27] M. E. Tipping and C. M. Bishop, "Probabilistic principal component analysis," *Journal of the Royal Statistical Society: Series B (Statistical Methodology)*, vol. 3, no. 61, pp. 611-622, 1999.

[28] G. Hinton and R. Salakhutdinov, "Reducing the dimensionality of data with neural networks," *Science (New York, N.Y.)*, vol. 313, pp. 504-7, 08 2006.

[29] F. L. Hitchcock, "The expression of a tensor or a polyadic as a sum of products," *J. Math. Phys*, vol. 6, no. 1, pp. 164-189, 1927.

[30] L. R. Tucker, "Implications of factor analysis of three-way matrices for measurement of change," *Problems in measuring change*, pp. 122-137, 1963.

[31] —, "The extension of factor analysis to three-dimensional matrices," *Contributions to mathematical psychology*, pp. 109-127, 1964.

[32] —, "Some mathematical notes on three-mode factor analysis," *Psychometrika*, vol. 31, pp. 279-311, 1966.

[33] T. G. Kolda and B. W. Bader, "Tensor decompositions and applications," *SIAM Review*, vol. 51, no. 3, pp. 455-500, 2009.

[34] A. Cichocki, D. Mandic, L. De Lathauwer, G. Zhou, Q. Zhao, C. Caiafa, and H. A. Phan, "Tensor decompositions for signal processing applications: From two-way to multiway component analysis," *IEEE Signal Processing Magazine*, vol. 32, no. 2, pp. 145-163, 2015.

[35] Q. Song, H. Ge, J. Caverlee, and X. Hu, "Tensor completion algorithms in big data analytics," *ACM Trans. Knowl. Discov. Data*, vol. 13, no. 1, Jan. 2019.

[36] G. Favier and D. S. Rocha, "Overview of tensor-based cooperative mimo communication systems - part 1: Tensor modeling," *Entropy*, vol. 25, no. 8, 2023.

[37] —, "Overview of tensor-based cooperative mimo communication systems - part 2: Semi-blind receivers," *Entropy*, vol. 26, no. 11, 2024.

[38] Z.-Y. Zhang, T. Li, and C. Ding, "Non-negative tri-factor tensor decomposition with applications," *Knowledge and Information Systems*, vol. 34, no. 2, pp. 243-265, Feb 2013.

[39] R. Henriques and S. Madeira, "Triclustering algorithms for three-dimensional data analysis: A comprehensive survey," *ACM Computing Surveys*, vol. 51, pp. 1-43, 09 2018.

[40] X. Wu, R. Zurita-Milla, E. I. Verdiguier, and M.-J. Kraak, "Tri-clustering georeferenced time series for analyzing patterns of intra-annual variability in temperature," *Annals of the American Association of Geographers*, vol. 108, no. 1, pp. 71-87, 2018.

[41] M. Selosse, A. Gourru, J. Jacques, and J. Velcin, "Tri-clustering pour données de comptage," in *Journées des Statistiques*, Nancy, France, Jun. 2019.

[42] R. Boutalbi, L. Labiod, and M. Nadif, "Tensor latent block model for co-clustering," *International Journal of Data Science and Analytics*, vol. 10, no. 2, pp. 161-175, Aug 2020.

[43] X. Wu, C. Cheng, R. Zurita-Milla, and C. Song, "An overview of clustering methods for geo-referenced time series: from one-way clustering to co- and tri-clustering," *International Journal of Geographical Information Science*, vol. 34, no. 9, pp. 1822-1848, 2020.

[44] E. Battaglia, F. Peiretti, and R. Pensa, "Co-clustering: A survey of the main methods, recent trends, and open problems," *ACM Computing Surveys*, vol. 57, 10 2024.

[45] K. Song, X. Yao, F. Nie, X. Li, and M. Xu, "Weighted bilateral k-means algorithm for fast co-clustering and fast spectral clustering," *Pattern Recognition*, vol. 109, p. 107560, 2021.

[46] A. Dempster, N. Laird, and D. Rubin, "Maximum-likelihood from incomplete data via the EM algorithm," *J. Royal Statist. Soc. Ser. B.*, 39, pp. 1-38, 1977.

[47] J. Tobin, M. Black, J. Ng, D. Rankin, J. Wallace, C. Hughes, L. Hoey, A. Moore, J. Wang, G. Horigan, P. Carlin, H. McNulty, A. M. Molloy, and M. Zhang, "Co-clustering multi-view data using the latent block model," 2024.

[48] M. Ailem, F. Role, and M. Nadif, "Model-based co-clustering for the effective handling of sparse data," *Pattern Recognition*, vol. 72, pp. 108-122, 2017.

[49] M. Selosse, J. Jacques, and C. Biernacki, "Textual data summarization using the self-organized co-clustering model," *Pattern Recognition*, vol. 103, 2020.

[50] D. R. Hunter and K. Lange, "A tutorial on mm algorithms," *Amer. Statist.*, no. 58, pp. 30-37, 2004.

[51] S. Lee, J. Z. Huang, and J. Hu, "Sparse logistic principal components analysis for binary data," *The Annals of Applied Statistics*, vol. 4, no. 3, pp. 579-1601, 2010.

[52] J. P. Benzecri, *L'analyse des données tome 1 et 2 : l'analyse des correspondances*. Paris:Dunod, 1980.
[53] M. Greenacre, *Theory and Applications of Correspondence Analysis*. London: Academic Press, 1983.

[54] W. Fithian and J. Josse, "Multiple correspondence analysis and the multilogit bilinear model," *Journal of Multivariate Analysis*, vol. 157, pp. 87–102, 2017.

[55] E. J. Beh and R. Lombardo, "Multiple and multiway correspondence analysis," *WIREs Computational Statistics*, vol. 11, no. 5, p. e1464, 2019.

[56] D. Hong, T. G. Kolda, and J. A. Duersch, "Generalized canonical polyadic tensor decomposition," *SIAM Review*, vol. 62, no. 1, pp. 133–163, 2020.

[57] M. van de Velden, A. I. D'Enza, and F. Palumbo, "Cluster correspondence analysis," *Psychometrika*, vol. 82, no. 1, pp. 158–185, Mar 2017.

[58] M. Ailem, F. Role, and M. Nadif, "Sparse poisson latent block model for document clustering," *IEEE Transactions on Knowledge and Data Engineering*, vol. 29, no. 7, pp. 1563–1576, 2017.

[59] E. Battaglia, F. Peiretti, and R. G. Pensa, "Fast parameterless prototype-based co-clustering," *Machine Learning*, vol. 113, no. 4, pp. 2153–2181, Apr 2024.

[60] T. Caliński and J. Harabasz, "A dendrite method for cluster analysis," *Communications in Statistics*, vol. 3, no. 1, pp. 1–27, 1974.

[61] P. Rousseeuw, "Silhouettes: a graphical aid to the interpretation and validation of cluster analysis," *J. Comput. Appl. Math.*, vol. 20, pp. 53–65, November 1987.

[62] A. A. Rontogiannis, E. Kofidis, and P. V. Giampouras, "Block-term tensor decomposition: Model selection and computation," in *2020 28th European Signal Processing Conference (EUSIPCO)*, 2021, pp. 1976–1980.

[63] H. Hsu, S. Salamatian, and F. P. Calmon, "Correspondence analysis using neural networks," in *Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics*, vol. 89. PMLR, 2019, pp. 2671–2680.

[64] D. Xu, W. Cheng, B. Zong, J. Ni, D. Song, W. Yu, Y. Chen, H. Chen, and X. Zhang, "Deep co-clustering," in *Proceedings of the 2019 SIAM International Conference on Data Mining*. SIAM, 2019.