Subspace Clustering with the Multivariate-t Distribution

TL;DR

This paper introduces a robust subspace clustering method based on the multivariate-t distribution, improving upon Gaussian mixture models for high-dimensional data analysis, with applications demonstrated on simulated and real satellite imagery data.

Contribution

It extends the HDDC framework by incorporating the multivariate-t distribution, creating 28 models that enhance robustness and address Gaussian model limitations.

Findings

tHDDC outperforms HDDC in robustness to outliers

The method effectively clusters high-dimensional data in lower-dimensional subspaces

Application to satellite imagery demonstrates practical utility

Abstract

Clustering procedures suitable for the analysis of very high-dimensional data are needed for many modern data sets. In model-based clustering, a method called high-dimensional data clustering (HDDC) uses a family of Gaussian mixture models for clustering. HDDC is based on the idea that high-dimensional data usually exists in lower-dimensional subspaces; as such, an intrinsic dimension for each sub-population of the observed data can be estimated and cluster analysis can be performed in this lower-dimensional subspace. As a result, only a fraction of the total number of parameters need to be estimated and a computationally efficient parameter estimation scheme based on the EM algorithm was developed. This family of models has gained attention due to its superior classification performance compared to other families of mixture models; however, it still suffers from the usual limitations of Gaussian mixture model-based approaches. In this paper, a robust analogue of the HDDC approach is proposed. This approach, which extends the HDDC procedure to include the mulitvariate-t distribution, encompasses 28 models that rectify the aforementioned shortcomings of the HDDC procedure. Our tHDDC procedure is fitted to both simulated and real data sets and is compared to the HDDC procedure using an image reconstruction problem that arose from satellite imagery of Mars' surface.

Figures

Tables

Table 1. Table 1: Nomenclature, covariance structure, and number of free covariance parameters for each member of the GPCM family.

Model	Volume	Shape	Orientation	${𝚺}_g$	Free covariance parameters
EII	Equal	Spherical	–	$\lambda 𝐈$	1
VII	Variable	Spherical	–	${\lambda }_g𝐈$	$G$
EEI	Equal	Equal	Axis-Aligned	$\lambda 𝐀$	$p$
VEI	Variable	Equal	Axis-Aligned	${\lambda }_g𝐀$	$p+G-1$
EVI	Equal	Variable	Axis-Aligned	$\lambda {𝐀}_g$	$pG-G+1$
VVI	Variable	Variable	Axis-Aligned	${\lambda }_g{𝐀}_g$	$pG$
EEE	Equal	Equal	Equal	$\lambda {\mathrm{𝐃𝐀𝐃}}^{\prime }$	$p\left(p+1\right)/2$
EEV	Equal	Equal	Variable	$\lambda {𝐃}_g{\mathrm{𝐀𝐃}}_g^{\prime }$	$Gp\left(p+1\right)/2-\left(G-1\right)p$
VEV	Variable	Equal	Variable	${\lambda }_g{𝐃}_g{\mathrm{𝐀𝐃}}_g^{\prime }$	$Gp\left(p+1\right)/2-\left(G-1\right)\left(p-1\right)$
VVV	Variable	Variable	Variable	${\lambda }_g{𝐃}_g{𝐀}_g{𝐃}_g^{\prime }$	$Gp\left(p+1\right)/2$
EVE	Equal	Variable	Equal	$\lambda {\mathrm{𝐃𝐀}}_g{𝐃}^{\prime }$	$p\left(p+1\right)/2+\left(G-1\right)\left(p-1\right)$
VVE	Variable	Variable	Equal	${\lambda }_g{\mathrm{𝐃𝐀}}_g{𝐃}^{\prime }$	$p\left(p+1\right)/2+\left(G-1\right)p$
VEE	Variable	Equal	Equal	${\lambda }_g{\mathrm{𝐃𝐀𝐃}}^{\prime }$	$p\left(p+1\right)/2+\left(G-1\right)$
EVV	Equal	Variable	Variable	$\lambda {𝐃}_g{𝐀}_g{𝐃}_g^{\prime }$	$Gp\left(p+1\right)/2-\left(G-1\right)$

Table 2. Table 2: Nomenclature, covariance decomposition and number of free covariance parameters for the t 𝑡 t HDDC models. For constaints on a g ​ i subscript 𝑎 𝑔 𝑖 a_{gi} , U represents unconstrained, D represents constrained accross dimension, G represents constrained across groups and C represents constrained accross both dimension and group. For all other components, U and C are is unconstrained and constrained across groups, respectively. For the number of free parameters, ρ = G ​ p + G + 1 𝜌 𝐺 𝑝 𝐺 1 \rho=Gp+G+1 is the number of parameters required to estimate the mean and proportions. The number of parameters required to estimate D ~ g subscript ~ 𝐷 𝑔 \tilde{D}_{g} , D ~ ~ 𝐷 \tilde{D} and s = ∑ g = 1 G d g 𝑠 superscript subscript 𝑔 1 𝐺 subscript 𝑑 𝑔 s=\sum_{g=1}^{G}d_{g} are τ = d ​ [ p − ( d + 1 ) / 2 ] 𝜏 𝑑 delimited-[] 𝑝 𝑑 1 2 \tau=d[p-(d+1)/2] and τ ¯ = d g ​ [ p − ( d g + 1 ) / 2 ] ¯ 𝜏 subscript 𝑑 𝑔 delimited-[] 𝑝 subscript 𝑑 𝑔 1 2 \bar{\tau}=d_{g}[p-(d_{g}+1)/2]

Model	$a_{jg}=a_g/a_j$	$b_g=b$	${𝐃}_g=𝐃$	$d_g=d$	${\nu }_g=\nu$	Number of
						Covariance Parameters
UUUUU	U	U	U	U	U	$\rho +\overline{\tau }+3G+s$
UCUUU	U	C	U	U	U	$\rho +\overline{\tau }+2G+s+1$
DUUUU	D	U	U	U	U	$\rho +\overline{\tau }+4G$
CUUUU	C	U	U	U	U	$\rho +\overline{\tau }+3G+1$
DCUUU	D	C	U	U	U	$\rho +\overline{\tau }+3G+1$
CCUUU	C	C	U	U	U	$\rho +\overline{\tau }+2G+2$
UUUCU	U	U	U	C	U	$\rho +G\left(\tau +d+2\right)+1$
UCUCU	U	C	U	C	U	$\rho +G\left(\tau +d+1\right)+2$
DUUCU	D	U	U	C	U	$\rho +G\left(\tau +2+1\right)+1$
CUUCU	C	U	U	C	U	$\rho +G\left(\tau +2\right)+2$
DCUCU	D	C	U	C	U	$\rho +G\left(\tau +2\right)+2$
CCUCU	C	C	U	C	U	$\rho +G\left(\tau +1\right)+3$
GCCCU	G	C	C	C	U	$\rho +\tau +d+G+2$
CCCCU	C	C	C	C	U	$\rho +\tau +G+3$
UUUUC	U	U	U	U	C	$\rho +\overline{\tau }+2G+s+1$
UCUUC	U	C	U	U	C	$\rho +\overline{\tau }+G+s+2$
DUUUC	D	U	U	U	C	$\rho +\overline{\tau }+3G+1$
CUUUC	C	U	U	U	C	$\rho +\overline{\tau }+2G+2$
DCUUC	D	C	U	U	C	$\rho +\overline{\tau }+2G+2$
CCUUC	C	C	U	U	C	$\rho +\overline{\tau }+G+3$
UUUCC	U	U	U	C	C	$\rho +G\left(\tau +d+1\right)+2$
UCUCC	U	C	U	C	C	$\rho +G\left(\tau +d\right)+3$
DUUCC	D	U	U	C	C	$\rho +G\left(\tau +2\right)+2$
CUUCC	C	U	U	C	C	$\rho +G\left(\tau +1\right)+3$
DCUCC	D	C	U	C	C	$\rho +G\left(\tau +1\right)+3$
CCUCC	C	C	U	C	C	$\rho +G\tau +4$
GCCCC	G	C	C	C	C	$\rho +\tau +d+3$
CCCCC	C	C	C	C	C	$\rho +\tau +4$

Table 3. Table 3: Mean and standard deviation of ARI values returned by the best fitting HDDC and t 𝑡 t HDDC models found by the BIC for the simulated multivariate- t 𝑡 t data sets.

	Mean of ARI	Standard Deviation of ARI
HDDC	0.021	0.012
$t$HDDC	0.995	0.005

Table 4. Table 4: Model decomposition, number of components, BIC and ARI values for the best fitting tHDDC and HDDC models found for the Iris data.

	Model	$G$	BIC	ARI
$t$HDDC	UUUCC	$3$	$-646.327$	0.904
HDDC	${\text{A}}_{kj}{\text{B}}_k{\text{Q}}_k\text{D}$	3	$-588.01$	0.868

Table 5. Table 5: A classification table showing the results for the selected three-component UUUCC t 𝑡 t HDDC model for the iris data.

	A	B	C
Setosa	50	0	0
Versicolor	0	45	5
Virginica	0	0	50

Table 6. Table 6: Model decomposition, number of components, BIC and ARI values for the best fitting HDDC and t 𝑡 t HDDC models when fitted for g = 1 , … , 10 𝑔 1 … 10 g=1,\ldots,10 and g = 1 , … , 5 𝑔 1 … 5 g=1,\ldots,5 components for the Italian wines data.

	Model	$G$	BIC	ARI
HDDC	${\text{A}}_j\text{B}\text{Q}\text{D}$	$8$	$-12,071.63$	0.658
$t$HDDC	GCCCC	$4$	$-11,965.26$	0.758

Table 7. Table 7: Model decomposition, BIC and ARI values for the selected eight-component HDDC and t 𝑡 t HDDC models for the Martian surface data.

	$G$	Model	BIC	ARI
HDDC	$8$	$A{B}_k{Q}_k{D}_k$	$62,460,591$	0.319
$t$HDDC	$8$	CCCCC	$64,249,591$	0.351

Table 8. Table 8: Model decomposition, BIC and ARI values for the selected five-component HDDC and t 𝑡 t HDDC models for the Martian surface data.

	$G$	Best Model	BIC	ARI
HDDC	$5$	$A_{kj}{B}_k{Q}_k{D}_k$	$61,956,344$	0.472
$t$HDDC	$5$	UCUUC	$70,120,085$	0.645

Table 9. Table 9: A classification table comparing the best fitting five-component t 𝑡 t HDDC and HDDC models for the Martian surface data.

		HDDC
		A	B	C	D	E
	1	10744	22	0	0	2973
	2	1019	1598	0	99	785
$t$HDDC	3	0	39	7807	4372	4
	4	2	319	1	4871	716
	5	605	914	4	283	1223

Equations

$$f(\mathbf{x}\,|\,\mbox{\boldmath$\vartheta$})=\sum_{g=1}^{G}\pi_{g}\,p_{g}(\mathbf{x}\,|\,\mbox{\boldmath$\theta$}_{g}),$$

$$\displaystyle p_{g}\left(\mathbf{x}\mid\mbox{\boldmath$\theta$}_{g}\right)=f_{t}(\mathbf{x}\mid\mbox{\boldmath$\mu$}_{g},\mbox{\boldmath$\Sigma$}_{g},\nu_{g})=\frac{\Gamma\left[\left(\nu_{g}+p\right)/2\right]\lvert\mbox{\boldmath$\Sigma$}_{g}\rvert^{-1/2}(\pi\nu_{g})^{-p/2}}{\Gamma\left[\nu_{g}/2\right]\left[1+\delta\left(\mathbf{x},\mbox{\boldmath$\mu$}_{g}\mid\mbox{\boldmath$\Sigma$}_{g}\right)/\nu_{g}\right]^{\left(\nu_{g}+p\right)/2}},$$

$$\mbox{\boldmath$\Delta$}_{g}=\left(\begin{array}[]{ccc|ccccc}a_{1g}&&0&&&&&\\ &\ddots&&&&\textbf{0}&&\\ 0&&a_{d_{g}g}&&&&&\\ \hline\cr&&&b_{g}&&&&0\\ &&&&\ddots&&&\\ &\textbf{0}&&&&\ddots&&\\ &&&&&&\ddots&\\ &&&0&&&&b_{g}\\ \end{array}\right),$$

$$P_{g}(\mathbf{x})=\tilde{\mathbf{D}}_{g}\tilde{\mathbf{D}}_{g}^{\prime}(\mathbf{x}-\mbox{\boldmath$\mu$}_{g})+\mbox{\boldmath$\mu$}_{g}$$

$$P_{g}^{\perp}(x)=\bar{\mathbf{D}}_{g}\bar{\mathbf{D}}_{g}^{\prime}(\mathbf{x}-\mbox{\boldmath$\mu$}_{g})+\mbox{\boldmath$\mu$}_{g}$$

$$z_{ig}=\left\{\begin{array}[]{ll}1&\mbox{if observation $\mathbf{x}_{i}$ belongs to component $g$ }\\ 0&\mbox{otherwise. }\end{array}\right.$$

$$\mathbb{E}[Z_{ig}\mid\mathbf{x}_{i}]=\frac{\pi_{g}p_{g}(\mathbf{x}\mid\mbox{\boldmath$\vartheta$}_{g})}{\sum_{h=1}^{G}\pi_{h}p_{h}(\mathbf{x}\mid\mbox{\boldmath$\vartheta$}_{h})}=:\hat{z}_{ig},$$

$$\displaystyle-2\log{f_{t}(\mathbf{x}\mid\mbox{\boldmath$\mu$},\mbox{\boldmath$\Sigma$},\nu)}$$

$$\displaystyle+\left(\nu+p\right)\log{\left[1+\frac{1}{\nu}\left(\lvert\lvert\mbox{\boldmath$\mu$}-P\left(\mathbf{x}\right)\rvert\rvert^{2}_{\mathbf{A}}+\frac{1}{b}\lvert\lvert\mathbf{x}-P^{\prime}\left(\mathbf{x}\right)\rvert\rvert^{2}\right)\right]}$$

$$+\sum _{j=1}^d\log a_j+(p-d)\log b,$$

$${\hat{z}}_{ig}=\frac{1}{{\sum }_{h=1}^G\exp \left[\frac{1}{2}\left(K_g\left({\mathbf{x}}_i\right)-K_h\left({\mathbf{x}}_i\right)\right)\right]},$$

$$K_g\left({\mathbf{x}}_i\right)$$

$$\displaystyle~{}~{}~{}+\left(\nu_{g}+p\right)\log{\left[1+\frac{1}{\nu_{g}}\left({\color[rgb]{0,0,0}\lvert\lvert\mbox{\boldmath$\mu$}_{g}-P_{g}\left(\mathbf{x}\right)\rvert\rvert^{2}_{\mathbf{A}_{g}}}+{\color[rgb]{0,0,0}\frac{1}{b_{g}}\lvert\lvert\mathbf{x}-P_{g}^{\prime}\left(\mathbf{x}\right)\rvert\rvert^{2}}\right)\right]}$$

$$\displaystyle~{}~{}~{}+{\color[rgb]{0,0,0}\sum_{j=1}^{d_{g}}{\log{a_{jg}}}}-2\log{\pi_{g}},$$

$$\mathbb{E}[U_{ig}\mid {\mathbf{x}}_i,z_{ig}=1]$$

$$\hat{\pi}_{g}=\frac{n_{g}}{n}\quad\text{and}\quad\hat{\mbox{\boldmath$\mu$}}_{g}=\frac{\sum_{i=1}^{n}\hat{z}_{ig}\hat{u}_{ig}\mathbf{x}_{i}}{\sum_{i=1}^{n}\hat{z}_{ig}\hat{u}_{ig}},$$

$${\hat{\nu }}_g\approx \frac{-\exp (k)+2\exp \left(\phi \left(\frac{{\hat{\nu }}_g^{\text{old}}}{2}\right)+\frac{1-{\hat{\nu }}_g^{\text{old}}}{2}\right)\exp (k)}{1-\exp (k)}$$

$$k=-1$$

$$\displaystyle l\left(\hat{\mbox{\boldmath$\vartheta$}}\right)$$

$$+\log \pi -2\log \Gamma \left[\left({\nu }_g+p\right)/2\right]+2\log \Gamma \left({\nu }_g/2\right))$$

$$\hat{\mbox{\boldmath$\Sigma$}}_{g}=\frac{1}{n_{g}}\sum_{i=1}^{n}\hat{z}_{ig}\hat{u}_{ig}(\mathbf{x}_{i}-\hat{\mbox{\boldmath$\mu$}}_{g})(\mathbf{x}_{i}-\hat{\mbox{\boldmath$\mu$}}_{g})^{\prime}.$$

$$\hat{b}_{g}=\frac{1}{(p-d_{g})}\left(\text{tr}\left(\hat{\mbox{\boldmath$\Sigma$}}_{g}\right)-\sum_{j=1}^{d_{g}}a_{jg}\right).$$

$$l_{\infty }^{(k+1)}=l^{(k)}+\frac{1}{1-a^{(k)}}\left(l^{(k+1)}-l^{(k)}\right),$$

$$a^{(k)}=\left(l^{(k+1)}-l^{(k)}\right)/\left(l^{(k)}-l^{(k-1)}\right).$$

$$\text{ARI}=\frac{\text{number of pairwise agreements}}{\text{number of pairs}}.$$

Full text

1 Introduction

Cluster analysis refers to the practice of using statistical approaches to detect subgroups within a given data set. These subgroups can represent a physical attribute not described by the given explanatory variables, e.g., gender, income tax bracket or blood type, which can reveal important relationships among the observed data and may be a crucial component in the effective analysis of a given data set. Due to their construction, finite mixture models are very useful when modelling data that contain a finite collection of sub-populations because each component of the model can be used to represent one of these sub-populations. Reviews of the application of finite mixture models for clustering are given by Fraley and Raftery (2002); Bouveyron and Brunet-Saumard (2014) and McNicholas (2016b), and extensive details can be found in the monographs by McLachlan and Peel (2000) and McNicholas (2016a).

The density of a parametric finite mixture distribution is

[TABLE]

where $\pi_{g}>0$, such that $\sum_{g=1}^{G}\pi_{g}=1$, are called mixing proportions, p_{g}(\mathbf{x}\,|\,\mbox{\boldmath\theta}_{g}) are the component densities and \mbox{\boldmath\vartheta}=(\pi_{1},\ldots,\pi_{G},\mbox{\boldmath\theta}_{1},\ldots,\mbox{\boldmath\theta}_{G}) is a vector containing the model parameters. Herein, we follow convention and refer to the application of finite mixture models for clustering as model-based clustering.

The general mixture given in (1) can be specified to contain components of any univariate or multivariate probability distribution. Until the last decade or so, the majority of work on model-based clustering using multivariate component densities focused on the Gaussian mixture model. One of the first notable departures from Gaussianity was provided by Peel and McLachlan (2000), who utilized mixtures of multivariate-$t$ distributions for clustering. Despite rectifying a well known short-coming of the Gaussian mixture model by formulating a model that is robust to outliers, mixtures of multivariate-$t$ distributions have only gained popularity in the last few years (see McLachlan et al., 2007; Andrews and McNicholas, 2011a, b; Baek and McLachlan, 2011; Steane et al., 2012; Andrews and McNicholas, 2012; Lin et al., 2014, for examples). As the name suggests, mixtures of multivariate-$t$ distributions assume that each sub-population of the observed data follow the multivariate-$t$ distribution. As such, the density of a mixture of multivariate-$t$ distribution is formulated by writing the component density in (1) as

[TABLE]

where $\Gamma\left(\cdot\right)$ is the Gamma function, $p$ is the number of dimensions in the observed data set, \mbox{\boldmath\mu}_{g} is the component location parameter, \mbox{\boldmath\Sigma}_{g} is the component covariance matrix, $\nu_{g}$ parameterizes the degrees of freedom in each component, and \delta(\mathbf{x},\mbox{\boldmath\mu}_{g}\mid\mbox{\boldmath\Sigma}_{g})=\left(\mathbf{x}-\mbox{\boldmath\mu}_{g}\right)^{\prime}\mbox{\boldmath\Sigma}_{g}^{-1}\left(\mathbf{x}-\mbox{\boldmath\mu}_{g}\right) is the Mahalanobis distance between $\mathbf{x}$ and \mbox{\boldmath\mu}_{g} for $g=1,\ldots,G$.

A family of mixture models emerges when we introduce constraints on the component densities. Some families of Gaussian mixture models are well established and widely used, e.g., the Gaussian parsimonious clustering models (GPCM; Celeux and Govaert, 1995) which arise from constraints being imposed on the eigen-decomposed covariance structure in a Gaussian mixture model. This eigen-decomposition is \mbox{\boldmath\Sigma}_{g}=\lambda_{g}\mathbf{D}_{g}\mathbf{A}_{g}\mathbf{D}_{g}^{{}^{\prime}}, where $\lambda_{g}$ is a constant, $\mathbf{D}_{g}$ is a matrix of eigenvectors, and $\mathbf{A}_{g}$ is a diagonal matrix with $|\mathbf{A}_{g}|=1$ and entries proportional to the eigenvalues of \mbox{\boldmath\Sigma}_{g}. Applying a combination of the constraints: $\lambda_{g}=\lambda,\mathbf{A}_{g}=\mathbf{A},\mathbf{D}_{g}=\mathbf{D},\mathbf{D}_{g}=\mathbf{I}$, $\mathbf{A}_{g}=\mathbf{I}$, where $\mathbf{I}$ is the identity matrix of the appropriate dimension, across the groups in the data creates a family of fourteen models and allows for various shapes and sizes of clusters (Table 1). In more than half of these fourteen models there are $\mathcal{O}(p^{2})$ free parameters to be estimated; hence, with higher dimensions, it can be very computationally inefficient to use these models. The GPCMs are supported by the R packages mclust (Fraley and Raftery, 2002) and mixture (Browne and McNicholas, 2014; Browne et al., 2015).

The multivariate-$t$ analog of the GPCM family for the mixtures of multivariate-$t$ distributions is the $t$EIGEN family (Andrews and McNicholas, 2012). These models use the same eigen-decomposition as the GPCM family and therefore the same constraints mentioned above can be applied, in addition to $\nu_{g}=\nu$. By combining these constraints, a total of 28 different models are derived. In R, all 28 models are supported by the teigen package (Andrews and McNicholas, 2014; Andrews et al., 2017).

Another popular family of mixture models are the Parsimonious Gaussian Mixture Models (PGMM; McNicholas and Murphy, 2008). These models are an extension of the mixture of factor analyzers (Ghahramani and Hinton, 1997) whose component covariance matrices are written as \mbox{\boldmath\Sigma}_{g}=\mbox{\boldmath\Lambda}_{g}\mbox{\boldmath\Lambda}_{g}^{\prime}+\mbox{\boldmath\Psi}_{g}, where \mbox{\boldmath\Lambda}_{g} is a $p\times q$ loading matrix with $q<p$, and \mbox{\boldmath\Psi}_{g} is a diagonal $p\times p$ matrix with positive entries for $g=1,\ldots,G$. By imposing constraints on \mbox{\boldmath\Lambda}_{g} and \mbox{\boldmath\Psi}_{g} across the components, McNicholas and Murphy (2008) introduced eight parsimonious models in which the number of free parameters is $\mathcal{O}(p)$ so that the number of covariance parameters grows linearly with dimension. For this reason, these models are more appropriate than the GPCMs for high-dimensional data. These models can be implemented via the pgmm package for R (McNicholas et al., 2015). Note: McNicholas and Murphy (2010) extended the PGMMs to include four new models by setting \mbox{\boldmath\Psi}_{g}=\omega_{g}\mbox{\boldmath\Delta}_{g}, where $\omega_{g}\in\mathbb{R}^{+}$ and \mbox{\boldmath\Delta}_{g}=\,\mbox{diag}{\left\{\delta_{1},\delta_{2},\ldots,\delta_{p}\right\}} is a noise matrix, such that \lvert\mbox{\boldmath\Delta}_{g}\rvert=1.

The multivariate-$t$ analogue of the PGMMs, known as the mixture of multivariate $t$-factor analyzers (MM$t$FA) were introduced by McLachlan et al. (2007) and extended by Andrews and McNicholas (2011a, b). In the MM$t$FAs, the component covariance stucture is also parameterized as \mbox{\boldmath\Sigma}_{g}=\mbox{\boldmath\Lambda}_{g}\mbox{\boldmath\Lambda}_{g}^{\prime}+\mbox{\boldmath\Psi}_{g}. By applying the constraints: \mbox{\boldmath\Psi}_{g}=\psi_{g}\mathbf{I}, \mbox{\boldmath\Lambda}_{g}=\mbox{\boldmath\Lambda}, and $\nu_{g}=\nu$, Andrews and McNicholas (2011a) created a family of six models, whose covariance parameters grow linearly with $p$, and Andrews and McNicholas (2011b) extended this to a family of 24 models. It is worth noting that the probabilistic principal $t$-component analyzer model MPP$t$CA model is a special case of the MM$t$FA model, where \mbox{\boldmath\Psi}_{g}=\psi_{g}\mathbf{I}. This family of 24 models is supported by the mmtfa package (Andrews et al., 2015).

Bouveyron et al. (2007) proposed a high-dimensional data clustering (HDDC) technique that is also based on an eigen-decomposition of the covariance structure of the Gaussian mixture model. This technique projects the data into a lower-dimensional subspace spanned by a subset of the eigenvectors of \mbox{\boldmath\Sigma}_{g}. Formally, given a data set $\{\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\}$ of $n$ data points in p with $G$ sub-populations, this method assumes that high-dimensional data mostly rests in lower-dimensional subspaces. This assumption can drastically reduce the number of covariance parameters that require estimation and result in an efficient parameter estimation scheme. As with the GPCMs, Bouveyron et al. (2007) lets $\mathbf{D}_{g}$ be the orthogonal matrix of eigenvectors of \mbox{\boldmath\Sigma}_{g}, but instead considers a block-diagonal matrix, \mbox{\boldmath\Delta}_{g}, which contains the eigenvalues of \mbox{\boldmath\Sigma}_{g}. Formally, \mbox{\boldmath\Delta}_{g} has the following form:

[TABLE]

where the upper left block is of size $d_{g}\times d_{g}$, where $d_{g}\in\{1,p-1\}$ is the intrinsic dimension in each component, or cluster, and the lower right block is of size $(p-d_{g})\times(p-d_{g})$, with $a_{jg}>b_{g}$, for $j=1,\ldots,d_{g}$ for $g=1,\ldots,G$. Bouveyron et al. (2007) proposed two methods for estimating the intrinsic dimension in each component of this eigen-decomposed GMM. The first approach utilizes the scree-test of Cattell (1966) and the second approaches utilizes the probabilistic Bayesian Information Criterion (BIC; Schwarz, 1978), which is given by \text{BIC}=2l(\hat{\mbox{\boldmath\vartheta}})-\rho\log n, where $\rho$ is the number of parameters in the model, $n$ is the number of observations, and l(\hat{\mbox{\boldmath\vartheta}}) is the maximized log-likelihood value. The eigenvectors associated with the eigenvalues $a_{jg}$, for $j=1,\ldots,d_{g}$, span a subspace $\mathbb{E}_{g}\in{}^{d_{g}}$ for each cluster, such that \mbox{\boldmath\mu}_{g}\in\mathbb{E}_{g}. The affine subspace $\mathbb{E}_{g}^{\perp}$ is defined such that $\mathbb{E}_{g}\otimes\mathbb{E}_{g}^{\perp}={}^{p}$ and \mbox{\boldmath\mu}_{g}\in\mathbb{E}_{g}^{\perp}. Each observation $\mathbf{x}_{i}$ is then projected onto the subspace $\mathbb{E}_{g}$, which is called the specific subspace of the $g$th group since most of the data are assumed to live on or near this subspace. This decomposition leads to 28 possible models by constraining the parameters $[a_{jg},b_{g},\mathbf{D}_{g},d_{g}]$ across the $G$ components. Of these 28 models, 14 have been implemented in the R package HDclassif (Berge et al., 2012).

Like the other Gaussian mixture model based approaches, this technique will suffer from the usual limitations, i.e., its parameter estimation scheme will not be robust to outliers. Herein, we discuss the derivation of a multivariate-$t$ analogue of the HDDC approach. This paper proceeds as follows: in Section 2, we outline the derivation of a multivariate-$t$ high dimensional data clustering (tHDDC) approach and present a computationally efficient parameter estimation scheme to fit the resulting models, in Section 3 we assess the classification performance of this novel family of models using a simulation study and three real data sets, and in Section 4 we conclude with a discussion and suggestions for future work.

2 Methodology

We now lay out some groundwork for the newly proposed tHDDC approach. As previously mentioned, it is the $t$-analogue of HDDC method described in Section 1. As with all the families introduced earlier, the goal is to cluster a given data set $\{\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\}$ in p into $G$ homogeneous groups. The data are modelled by the general density in (1), with the multivariate component density given in (2). The general multivariate-$t$ mixture model requires the estimation of the full covariance structure, so the number of parameters to estimate is $\mathcal{O}(p^{2})$. As Bouveyron et al. (2007) describe, via the empty space phenomenon (Scott and Thompson, 1983), we can assume that most of the data live around lower-dimensional subspaces. By performing clustering in these lower-dimensional subspaces, the number of parameters to be estimated is reduced significantly.

2.1 The General $t$HDDC model

Analogous to the HDDC approach, we specify $\mathbf{D}_{g}$ to be the orthogonal matrix of eigenvectors and \mbox{\boldmath\Delta}_{g}=\mathbf{D}_{g}^{\prime}\mbox{\boldmath\Sigma}_{g}\mathbf{D}_{g}, where \mbox{\boldmath\Delta}_{g} is a class specific matrix of the form given in (3). Following Bouveyron et al. (2007) we define

[TABLE]

as the projection of $\mathbf{x}$ on $\mathbb{E}_{g}$ and

[TABLE]

as the projection of $\mathbf{x}$ on $\mathbb{E}_{g}^{\perp}$, where $\tilde{\mathbf{D}}_{g}$ consists of the first $d_{g}$ columns of $\mathbf{D}_{g}$, concatenated with $p-d_{g}$ zero columns and $\bar{\mathbf{D}}_{g}=\mathbf{D}_{g}-\tilde{\mathbf{D}}_{g}$. Each $t$HDDC model has parameters $a_{jg},b_{g},\mathbf{D}_{g},d_{g},\nu_{g}$ for $j=1,\ldots,d_{g}$ and $g=1,\ldots,G$. Applying group-wide constraints to these parameters can lead to a total of 56 possible models. (The $t$HDDC analogues of the HDDC models available in HDclassif are listed in Table 2).

2.2 The ECM for the General $t$HDDC Model

In model-based clustering, the EM algorithm (McLachlan and Krishnan, 2008) is the usual choice for parameter estimation. It is an iterative procedure that alternates between two steps: an E-step and a M-step, until convergence is reached. On the E-step, the expected value of the complete-data log-likelihood is updated given the current estimates of the parameters. In the M-step, the same complete-data log-likelihood is maximized in terms of the model parameters. We use a variation of the EM algorithm called the expectation conditional-maximization (ECM) algorithm (Meng and Rubin, 1993), which replaces each M-step with multiple CM-steps. For each $t$HDDC model the complete-data is made up of the observed $\mathbf{x}_{i}$, the latent $u_{ig}$, and the missing $z_{ig}$ for $i=1,\ldots,n$ and $g=1,\ldots,G$. Note that the $u_{ig}$ is a realization of a gamma distributed random variable, $U_{ig}$ that arises because we exploit the fact that the multivariate-$t$ distribution is a normal-variance mean mixture (Barndorff-Nielsen et al., ; Peel and McLachlan, 2000), whereas the $z_{ig}$ are introduced to represent component membership. Formally, we write that

[TABLE]

2.2.1 The E-step

For the general finite mixture model, the component indicator variables are usually replaced by their expected values,

[TABLE]

for $i=1,\ldots,n$ and $g=1,\ldots,G$. Unfortunately, this usually requires the computation of both the determinant and inverse of a $p\times p$ covariance matrix. To avoid these potentially cumbersome calculations, we follow Bouveyron et al. (2007) and derive a cost function that utilizes the projection functions: $P_{g}(x)$ and $P_{g}^{\perp}(x)$, defined in Section 2.1. The derivation of the cost function is as follows: first, note that we can write

[TABLE]

where $\lvert\lvert\mathbf{x}\rvert\rvert^{2}_{\mathbf{A}}=\mathbf{x}\mathbf{A}\mathbf{x}^{\prime}$ with \mathbf{A}=\tilde{\mathbf{D}}_{g}\mbox{\boldmath\Delta}_{g}\tilde{\mathbf{D}}_{g}^{\prime}, and all other values are as previously defined. So, on the $E$-step of the proposed ECM algorithm we replace each $z_{ig}$ with

[TABLE]

where we refer to

[TABLE]

for $i=1,\ldots,n$ and $g=1,\ldots,G$, as the cost function. Each $u_{ig}$ is then replaced by their expected values

[TABLE]

for $i=1,\ldots,n$ and $g=1,\ldots,G$ (cf. Peel and McLachlan, 2000; Andrews and McNicholas, 2012).

2.2.2 The CM-steps

On the first CM-step we update the mixing proportions and component location parameter using

[TABLE]

respectively, where $n_{g}=\sum_{i=1}^{n}\hat{z}_{ig}$. The degrees of freedom parameter, ${\nu}_{g}$, is updated using the closed form approximation given in Andrews et al. (2017). Formally, we let

[TABLE]

with

[TABLE]

where $\hat{\nu}_{g}^{\text{old}}$ is the estimate of $\nu_{g}$ from the previous iteration of this ECM algorithm, and $\varphi\left(\cdot\right)$ is the digamma function.

For each tHDDC model, the updates on the second CM step are analogous to those given in Bouveyron et al. (2007). For illustrative purposes, we outline how to update each covariance parameter for the UUUUU model, i.e., the model where $a_{jg},b_{g},\mathbf{D}_{g},d_{g},\nu_{g}$ are free to vary across all $g=1,\ldots,G$ and $j=1,\ldots,d_{g}$.

First, we calculate the intrinsic dimension, $d_{g}$. For each value of $j\in\left\{1,p-1\right\}$ we compute

[TABLE]

and set $d_{g}$ equal to the value of $d_{j}$ that maximizes the BIC values associated with the log-likelihood values found using (2.2.2). Then we let $\mathbf{D}_{g}$ be the eigenvectors of \hat{\mbox{\boldmath\Sigma}}_{g}, where

[TABLE]

Each $a_{jg}$, for $j=1,\ldots,d_{g}$ and $g=1,\ldots,G$ is then replaced with the first $d_{g}$ eigenvalues of \hat{\mbox{\boldmath\Sigma}}_{g} and we estimate $b_{g}$ using

[TABLE]

2.3 Computational Considerations

For the proposed ECM algorithm, we initialize each model using either k-means clustering or random starting values and use the Aitkin’s acceleration (Aitken, 1926) procedure to determine if the algorithm has converged. That is, we consider this ECM algorithm to have converged when $l_{\infty}^{(k+1)}-l^{(k)}<\epsilon$, where $\epsilon=10^{-2}$ (see Lindsay, 1995; McNicholas et al., 2010). In this criterion, $l^{(k)}$ is the log-likelihood value at iteration $(k)$ and $l_{\infty}^{(k+1)}$ is the asymptotic estimate of the log-likelihood at iteration $(k+1)$. Formally,

[TABLE]

where

[TABLE]

All analyses are performed in R version 3.3.2 (R Core Team, 2016) for Linux 6.5111Using a 32-core Intel Xeon E5 server with 256GB RAM running 64-bit CentOS. The HDDC models are fit using the default number of groups, whereas the $t$HDDC models are fit with $G=1,\ldots,4$. It is important to note that the only HDDC models considered are ones with a monotonic likelihood.

3 Applications

3.1 Performance Assessment and Model Selection

The data analyses will be treated as genuine clustering problems, where the true classifications are not known. Since we do have the true class labels, the adjusted Rand index (ARI; Hubert and Arabie, 1985) will be used to assess class agreement between the true class labels and the predicted labels rendered by the clustering techniques. The ARI was introduced to correct the Rand Index (Rand, 1971) for chance since the expected value of the Rand Index is greater than 0 for a random classification, making it hard to interpret. So, the ARI has expected value equal to ‘0’, with a perfect classification being represented by a score of ‘1’. Formally, the ARI can be written as

[TABLE]

In all applications, the best fitting models will be chosen using the BIC.

3.2 Simulation Studies

We use a simulation study to highlight the aforementioned drawback of the considered mixture of multivariate Gaussian distributions. Ten data sets were simulated from a two-component multivariate-$t$ distribution with $\nu_{1}=2$ and $\nu_{2}=3$. Figure 1 provides an illustration of the first three dimensions of one of the simulated data sets. In each component, observations are scattered from the mean, with many outliers on far ends of the clusters.

Table 3 gives the classification results for the tHDDC and HDDC models when fitted to the simulated data set. As expected, the tHDDC approach outperforms the HDDC approach, achieving a near perfect classification. The relatively small standard deviation reveals that the selected $t$HDDC models are consistently returning an exceptional classification performance, whereas the selected HDDC models are using extra components to account for the increased variation in the simulated data sets.

3.3 Fisher’s Irises

In our first real data analysis, we consider Fisher’s famous iris data set, which is available in the R package datasets. It is composed of four explanatory variables: sepal length, sepal width, petal length, and petal width, measured in centimetres. There are three species of the plant: Setosa, Versicolour and Virginica. Table 4 gives the classification results.

Both the best fitting HDDC and tHDDC models return a very good classification of the irises, with the best fitting tHDDC model outperforming the corresponding HDDC model. Across the three selected components, the best fitting tHDDC model uses a varying number of eigenvalues and eigenvectors with a constant intrinsic dimension and degrees of freedom. In total, this model misclassifies only 5 irises (see Table 5).

3.4 Italian Wines

The Italian wines data set, available as wine in the R package pgmm (McNicholas et al., 2015), is composed of 178 Italian wines on which 27 measurements are taken (Forina et al., 1986). The wines come from three different cultivars and are classified based on which one they come from: Barolo, Grignolino and Barbera. The $t$HDDC models are fit using $G=1,\ldots,5$, since the BIC will select a four-component mixture model. Table 6 gives the classification results. Note: when fitting only three-component HDDC and tHDDC models, the selected $A_{j}BQD$ and GCCCC models gave the same classification result (ARI = 0.933).

We can see that the best fitting HDDC model overfits by selecting a model with eight components. Although the best $t$HDDC model did not have three groups, the four group solution gives a superior classification performance (See Table 6).

3.5 Martian Surface

This data set was retrieved by the OMEGA instrument (Mars Express, ESA; Bibring et al., 2004). The OMEGA instrument is used for characterization of the Martian surface based on physical and chemical composition. This can include classes of silicates, hydrated minerals, ices and more. The data used is based on one $300\times 128$ raw image. It contains 255 variables on 38,400 observations. With a physical model, eight groups were found and for the purpose of this analysis, these will be treated as true groups; however, the best determination of model performance here is based on efficacy for image reconstruction. For $G=8$ components, $t$HDDC does a little better than HDDC; however, neither performs well (Table 7).

Although the physical model suggests eight groups, experts in the field are interested in exploring a five group solution (Bouveyron et al., 2007). Both HDDC and $t$HDDC models are applied to this data with $G=5$, and the selected $t$HDDC model recovers the clusters better than the selected HDDC model (see Table 8).

In Table 9, we can see that the classification results returned by the selected $t$HDDC and HDDC models are quite different.

Furthermore, comparing the recovered image based on the predicted classes to the original image (see Figure 2), the utility of the model becomes clear, i.e., the physical details are generally recovered very well.

4 Discussion

A new family of multivariate-$t$ mixture models has been proposed. The $t$HDDC approach is an extension of the HDDC approach that incorporates the multivariate-$t$ distribution, allowing for a more robust clustering scheme. A total of 28 models have been developed and the need for these models was shown through a simulation study which demonstrated their flexibility in recognizing outliers. The models were tested on both simulated and real data sets and show superior results when compared to the HDDC family. In particular, the results on the high-dimensional Martian surface data show that image recovery can be greatly improved. Overall, the added degrees of freedom parameter allows for more flexible clusters and a more flexible modelling structure than HDDC. In future work, this method can be extended to include skewed mixture models. Examples include the mixture of multivariate skew-$t$ distributions (Lin, 2010; Murray et al., 2014, 2017), the mixture of shifted asymmetric Laplace distributions (Franczak et al., 2014), the mixture of variance-gamma distributions (McNicholas et al., 2017), and the mixture of generalized hyperbolic distributions (Browne and McNicholas, ).

Acknowledgements

This work was supported by a Discovery Grants from the Natural Sciences and Engineering Research Council of Canada (Franczak, McNicholas) and the Canada Research Chairs Program (McNicholas).