Algorithms for an Efficient Tensor Biclustering

TL;DR

This paper introduces spectral decomposition algorithms for tensor biclustering, identifying subsets of individuals and features with similar signal trajectories over time, validated on synthetic and real datasets.

Contribution

It proposes novel spectral decomposition-based algorithms for tensor biclustering that effectively identify low-dimensional biclusters in multi-dimensional data.

Findings

Algorithms successfully identify meaningful biclusters

Effective on both synthetic and real datasets

Outperforms existing methods in accuracy

Abstract

Consider a data set collected by (individuals-features) pairs in different times. It can be represented as a tensor of three dimensions (Individuals, features and times). The tensor biclustering problem computes a subset of individuals and a subset of features whose signal trajectories over time lie in a low-dimensional subspace, modeling similarity among the signal trajectories while allowing different scalings across different individuals or different features. This approach are based on spectral decomposition in order to build the desired biclusters. We evaluate the quality of the results from each algorithms with both synthetic and real data set.

Figures

Equations

$$\mathcal{T}=\mathcal{X}+\mathcal{Z}=\sum _{r=1}^q{\sigma }_r(u_r^{J_1^{(r)}}\otimes w_r^{J_2^{(r)}}\otimes v_r)+\mathcal{Z},$$

$$T_{(j_1,1)}=\mathcal{T}(j_1,:,:)\text{ and }{T}_{(j_2,2)}=\mathcal{T}(:,j_2,:)$$

$$C_1=\sum _{j_2=1}^n{T}_{(j_2,2)}^t{T}_{(j_2,2)}\text{ and }{C}_2=\sum _{j_1=1}^n{T}_{(j_1,1)}^t{T}_{(j_1,1)}$$

$$u_{1}(j_{1})=\left\{\begin{array}[]{rcr}1/\sqrt{k}&&\quad\text{ for }j_{1}\in J_{1}^{(1)}\\ 0&&\quad\text{ if not }\\ \end{array}\right.,\qquad w_{1}(j_{2})=\left\{\begin{array}[]{rcr}1/\sqrt{k}&&\quad\text{ for }j_{2}\in J_{2}^{(1)}\\ 0&&\quad\text{ if not }\\ \end{array},\right.$$

$$u_{2}(j_{1})=\left\{\begin{array}[]{rcr}1/\sqrt{k}&&\quad\text{ for }j_{1}\in J_{1}^{(2)}\\ 0&&\quad\text{ if not }\\ \end{array}\right.,\qquad w_{2}(j_{2})=\left\{\begin{array}[]{rcr}1/\sqrt{k}&&\quad\text{ for }j_{2}\in J_{2}^{(2)}\\ 0&&\quad\text{ if not }\\ \end{array}\right.$$

$$0\le \frac{|{\widehat{J_1}}^{(1)}\cap J_1^{(1)}|}{4k}+\frac{|{\widehat{J_1}}^{(2)}\cap J_1^{(2)}|}{4k}+\frac{|{\widehat{J_2}}^{(1)}\cap J_2^{(1)}|}{4k}+\frac{|{\widehat{J_2}}^{(2)}\cap J_2^{(2)}|}{4k}\le 1.$$

Full text

Index terms— Multilinear Algebra, Tensor Decomposition, Principal Component Analysis.

1 Introduction

Clustering analysis has become a fundamental tool in statistics and machine learning. Many clustering algorithms have been developed with the general idea of seeking groups among different individuals in all space of features. Biclustering consists of simultaneous partitioning of a set of observations and a set of their features into subsets often called bicluster. Consequently, a subset of rows exhibiting significant coherence within a subset of columns in the matrix can be extracted, which corresponds to a specific coherent pattern [2, 8]. Nowadays, there is a new type of data collection, in which we may collect data by individual-feature pair at multiple times. The variation of a couple (individual-feature) at different instants is called trajectory. This data can be represented as a three dimensional object called tensor $\mathcal{T}\in\mathbb{R}^{n_{1}\times n_{2}\times m}$, where $n_{1}$ and $n_{2}$ are respectively the size of observations and features at $m$ different times. Tensor biclustering selects a subset of individual indices and a subset of features indices whose trajectories are highly correlated. Grouping those trajectories according to the correlation or similarity behaviour between them is useful in different area such as decision making, but it is still a very challenging topic in research.

In [7], the authors proposed different methods based on the spectral decomposition of matrix and the length of trajectory, although they provide a unique bicluster. Many tools on tensor manipulation already exist in literature to solve this tensor biclustering problem [4, 1, 6, 5, 3]. Our algorithms are based on a spectral decomposition as proposed in [7]. This article is structured as follows. In Section 2, we start by a brief summary of problem formulation. Section 3 introduces our algorithm extensions. Section 4 is related to the experiments. We make some concluding remarks in Section 5.

2 Problem Formulation

We use the common notation where $\mathcal{T},\;\mathcal{X}$ and $\mathcal{Z}$ are used respectively to denote input, signal and noise tensors. For any set $J$, $|\bar{J}|$ denotes its cardinality. $[n]$ denotes the set $\{1,2,\dots,n\}$. $|\bar{J}|=[n]-J$. $\|x\|_{2}=(x^{t}x)^{1/2}$ is the second norm of the vector $x$. $x\otimes y$ is the Kronecker product of two vectors $x$ and $y$. We also use Matlab notation to denote the elements in tensor. Specifically, $\mathcal{T}(:,:,i)$, $\mathcal{T}(:,i,:)$ and $\mathcal{T}(i,:,:)$ are respectively the $i-th$ frontal, lateral and horizontal slice. $\mathcal{T}(:,i,j)$, $\mathcal{T}(i,:,j)$ and $\mathcal{T}(i,j,:)$ denote respectively the $\mathrm{mode-1}$, $\mathrm{mode-2}$ and $\mathrm{mode-3}$ fiber. Let $\mathcal{T}\in\mathbb{R}^{n_{1}\times n_{2}\times m}$ a third-order tensor, $\mathcal{T}=\mathcal{X}+\mathcal{Z}$ where $\mathcal{X}$ is the signal tensor and $\mathcal{Z}$ is the noise tensor. Consider

[TABLE]

where $J_{1}^{(i)}$ and $J_{2}^{(i)}$ are respectively the sets of observations indices and features indices in the $i-th$ bicluster and $u_{r}\in\mathbb{R}^{n_{1}},\;w_{r}\in\mathbb{R}^{n_{2}}$ and $v_{r}\in\mathbb{R}^{m}$ are unit vectors. We assume that $u_{i}^{J_{1}^{(i)}}$ and $w_{i}^{J_{2}^{(i)}}$ have zero entries outside of ${J_{1}^{(i)}}$ and ${J_{2}^{(i)}}$ respectively for $i\in\{1,2\cdots,q\}$ and $\sigma_{1}\geq\sigma_{2}\geq\dots\geq\sigma_{q}>0$. we define $J_{1}~{}=~{}\bigcup_{i}J_{1}^{(i)}$ and $J_{2}~{}=~{}\bigcup_{i}J_{2}^{(i)}$. Under this model, trajectories $\mathcal{X}(J_{1},J_{2},:)$ form at most $q$ dimensional subspace.

Concerning the noise model, if $(j_{1},j_{2})\notin J_{1}\times J_{2}$, we assume that entries of the noise trajectory $\mathcal{Z}(j_{1},j_{2},:)$ are independent and identically distributed (i.i.d) and each entry has a standard normal distribution. If $(j_{1},j_{2})\in J_{1}\times J_{2}$, we assume that entries of $\mathcal{Z}(j_{1},j_{2},:)$ are i.i.d and each entry has a Gaussian distribution with zero means and $\sigma_{z}^{2}$ variance. We analyse tensor biclustering problem under two variances models of the noise trajectory:

Noise Model I: in this model, we assume $\sigma_{z}^{2}~{}=~{}1$, i.e., the variance of the noise within and outside of the clustering is assumed to be the same. Although this model simplifies the analysis, it has the following drawback: under this noise model, for every value of $\sigma_{1}$, the average trajectory lengths in the bicluster is larger than the average trajectory lengths outside the bicluster.

Indeed, let $T_{1}\in\mathbb{R}^{m\times k^{2}}$ be a matrix whose columns include trajectories $\mathcal{T}(j_{1},j_{2},:)$ for $(j_{1},j_{2})\in J_{1}\times J_{2}$ (i.e $T_{1}$ is the unfolded $\mathcal{T}(j_{1},j_{2},:)$). We can write $T_{1}~{}=~{}X_{1}~{}+~{}Z_{1}$ where $X_{1}$ and $Z_{1}$ are unfolded $\mathcal{X}(j_{1},j_{2},:)$ and $\mathcal{Z}(j_{1},j_{2},:)$, respectively. The squared Frobinius norm of $X_{1}$ is equal to $\|X_{1}\|_{F}^{2}~{}=~{}\sigma_{1}^{2}$. Morever, the squared Frobenius norm of $Z_{1}$ has a Chi-squared distribution with $mk^{2}$ degrees of freedom i.e $\chi^{2}(mk^{2})$. Thus, the average squared Frobenius norm of $T_{1}$ is equal to $\sigma_{1}^{2}+\sigma_{z}^{2}mk^{2}$. Let $T_{2}\in\mathbb{R}^{m\times k^{2}}$ be a matrix whose columns include only noise trajectories. Using a similar argument, we have $\mathbb{E}[\|T_{2}\|_{F}^{2}]~{}=~{}mk^{2}$, which is smaller than $\sigma_{1}^{2}~{}+~{}\sigma_{z}^{2}mk^{2}$. 2. -

Noise Model II: in this model, we assume \sigma_{z}^{2}~{}=~{}\max\big{(}0,1-(\sigma_{1}^{2}/mk^{2})\big{)}, i.e., $\sigma^{2}_{z}$ is modeled to minimize the difference between average trajectory lengths within and outside the bicluster.

Indeed, if $\sigma_{1}^{2}<mk^{2}$, without noise, the average trajectory in the bicluster is smaller than the one outside the bicluster. In this regime, having $\sigma_{z}^{2}~{}=~{}1~{}-~{}\sigma_{1}^{2}/mk^{2}$ makes the average trajectory lengths within and outside the bicluster comparable. This regime is called the low-SNR (signal noise ratio) regime. If $\sigma_{1}^{2}>mk^{2}$, the average trajectory lengths in the bicluster is larger than the one outside the bicluster. This regime is called high-SNR regime. In this regime, adding noise to signal trajectories increases their lengths and makes solving the tensor biclustering problem easier. Therefore, in this regime we assume $\sigma_{z}^{2}~{}=~{}0$ to minimize the difference between average trajectory lengths within and outside of the bicluster.

2.1 Tensor Folding and Spectral (FS)

The algorithm and the asymptotic behaviour of this method are available in [7]. Under the assumption $q=1$ and $n=|n_{1}|=|n_{2}|$, we drop the subscript $(1)$ from $J_{1}^{(1)}$ and $J_{2}^{(1)}$. We assume that $|J_{1}|=|J_{2}|=k$. The author propose to provide only one bicluster. This method separates the selection of the two sets $J_{1}$ and $J_{2}$ using lateral slice and horizontal slice of the tensor respectively.

[TABLE]

The aim is to select the row and column indices whose trajectories are highly correlated. The elements of $J_{1}$ and $J_{2}$ are the indices of the top $k$ elements of the top eigenvector of the matrix $C_{1}$ and $C_{2}$ respectively (algorithm 1). We denoted by $\hat{J}_{1}$ and $\hat{J_{2}}$ the subset of individuals and features respectively in the bicluster given from the algorithm.

Tensor FS method have the best performance in both noise models compared to the three another methods (tensor unfolding+spectral, thresholding sum of squared and individual trajectory lengths) proposed by Soheil Feizi, Hamid Javadi, David Tse [7].

3 Extension of Tensor Folding and Spectral

In this section, we aim to extract many biclusters in the tensor data and improve the quality of the result. We propose several methods in order to do this task. However instead of seeking only one bicluster, we assume that in equation (1) $q=r\geq 2$ where $r$ is defined by the number of gap in the eigenvalues of the covariance matrix $C_{1}$ and $C_{2}$ (equation (3)).

3.1 Recursive extension

The classical method extracts one rank of the low dimensional subspace which is not very interesting because it neglect the majority of the data sets. So, Direct improvement of this method is to compute recursively according to the number of gap shown in the eigenvalues (algorithm 2). In this method, there is no intersection in two different blocks of tensor biclustering.

3.2 Multiple biclusters

This method extract simultaneously the $r$ biclusters in our tensor by using the idea of top $r$ principal component analysis (PCA). The orthogonality of the principal components favor the quality of the result (algorithm 3).

The illustration step of tensor FS method is showed in the Fig.1. For each fix individual, we have a horizontal slice of the tensor represented by $m\times n_{2}$ matrix (equation (2)). Then, we compute the covariance matrix for each horizontal slice and their sum give us only one squared matrix of order $n_{2}$ ($C_{2}$ in equation (3)). We apply singular value decomposition (SVD) in $C_{2}$, the top $r$ eigenvectors in the matrix $C_{2}$ ensure the selection of the elements of the features index set $(J_{2}^{(i)})_{i\in[r]}$ (algorithm 3). A similar step is applied to each lateral slice of the tensor to find all the element of the index set $(J_{1}^{(i)})_{i\in[r]}$.

Since $k$ is a fix parameter, multiple bicluster method allow some trajectory belong to many blocks of tensor biclustering. We call them a boundary of bicluster. Those boundaries are very important as they belong to the intersection of all the biclusters. Thus they have all their properties.

4 Experimentation

4.1 Synthetic data

We build synthetic data to evaluate the implementation of our methods. In this dataset, we have two biclusters with signal strength $\sigma_{1}$ and $\sigma_{2}$ such that $\sigma_{1}~{}>~{}\sigma_{2}$. We assume that $v_{1}$ and $v_{2}$ are fixed unit vectors in $\mathbb{R}^{m}$ and $v_{1}~{}=~{}v2$. We assume also that $J_{1}^{(1)}\bigcap J_{1}^{(2)}=\emptyset$ and $J_{2}^{(1)}\bigcap J_{2}^{(2)}=\emptyset$. We have $n~{}=~{}n_{1}~{}=~{}n_{2}~{}=~{}150$, $m~{}=~{}40$ and $k=|J_{1}^{(1)}|=|J_{1}^{(2)}|=|J_{2}^{(1)}|=|J_{2}^{(2)}|=30$, we assume:

[TABLE]

[TABLE]

We apply the assumption above to generate the input tensor $\mathcal{T}$ with the noise model II define in section 2. Let $\hat{J_{1}}^{(1)}\times\hat{J_{2}}^{(1)}$ and $\hat{J_{1}}^{(2)}\times\hat{J_{2}}^{(2)}$ be the two estimated biclusters indices of $J_{1}^{(1)}\times J_{2}^{(1)}$ and $J_{1}^{(2)}\times J_{2}^{(2)}$ respectively where $|J_{1}^{(1)}|=|J_{2}^{(1)}|=|J_{1}^{(2)}|=|J_{2}^{(2)}|=k$. We fix the signal strength $\sigma_{2}=2\sigma_{1}/3$, if the value of $\sigma_{1}>90$, the bar plot of the top five eigenvalues of both covariance matrices tell us that there is two block of tensor biclustering in the data (see Fig.2).

In this case, we know the value of parameter $k=30$. To evaluate the inference quality of the result given from the algorithm, we compute the recovery rate:

[TABLE]

Recovery rate return zero if the algorithm do not find any of the element of the two biclusters and return one if the algorithm find all the elements of the two biclusters.

We did the experiment with different value of signal strength ($\sigma_{1}$) and for each value of $\sigma_{1}$ we repeat 10 times. Then we compute the average of the recovery rate (Figure 2(c)).

4.2 Real Data

We apply the both contribution algorithms to an electricity load diagrams 111http://archive.ics.uci.edu/ml/datasets/ElectricityLoadDiagrams20112014data set during four years (2011-2014). This data set contains electricity consumption of 370 clients for each 15 minutes during four years. After the data prepossessing, we have a tensor $\mathcal{T}\in\mathbb{R}^{n_{1}\times n_{2}\times m}$ where $n_{1}~{}=~{}365$ is the number of day in one year, $n_{2}~{}=~{}161$ is the number of clients and $m~{}=~{}4$ is the number of years.

As illustrated in Fig.(3(a),3(d)), the gap on the eigenvalues shows the existence of two tensor biclustering (section 3) in the data set. The parameter $k$ cardinality of each index sets are defined from the multiple bicluster method, we choose $k$ with few intersection of two blocks of bicluster $|J_{1}^{(1)}|~{}=~{}|J_{1}^{(2)}|~{}=~{}50$ and $|J_{2}^{(1)}|~{}=~{}|J_{2}^{(2)}|~{}=~{}25$. After compilation, we note that the two individuals sets $J_{1}^{(1)}$ and $J_{1}^{(2)}$ are disjoint in both methods. Besides, the two features sets have 22 intersection elements for multiple biclusters method and 19 intersection elements for the recursive method. So, we have two distinct blocks with two distinct subsets of individuals and one subset of feature.

To evaluate the quality of the bicluster given for each algorithm, we compute the total absolute pairwise correlations of the trajectories among each bicluster. With the recursive method, the the trajectories in first bicluster is highly correlated but the quality of the second bicluster is a little bit low as seen in Figure.(3(b), 3(c)). Besides, with multiple bicluster method, the trajectories on both biclusters are highly correlated as seen in Figure.(3(e), 3(f)).

5 Conclusion

In this article, we introduced two methods to increase the number of bicluster selected in the tensor data set based on [7], which depends on the number of rank of the low dimensional subspace. The goal is to extract $r$ subsets of tensor ($r\geq 2$) rows and columns such that each block of the trajectories form a low dimensional subspace. We proposed two algorithms to solve this problem, tensor recursive and multiple bicluster. The performance of both algorithms depends on the parameter $k$, one way to choose this parameter is in the multiple bicluster method. If the parameter chosen gives a lot of index intersections, decreasing the value of $k$ is a good idea to improve the quality of the results.