On the Use of the Matrix-Variate Tail-Inflated Normal Distribution for Parsimonious Mixture Modeling

Recent advances in the matrix-variate model-based clustering literature have shown the growing interest for this kind of data modelization. In this framework, finite mixture models constitute a powerful clustering technique, despite the fact that they tend to suffer from overparameterization problems because of the high number of parameters to be estimated. To cope with this issue, parsimonious matrix-variate normal mixtures have been recently proposed in the literature. However, for many real phenomena, the tails of the mixture components of such models are lighter than required, with a direct effect on the corresponding fitting results. Thus, in this paper we introduce a family of 196 parsimonious mixture models based on the matrix-variate tail-inflated normal distribution, an elliptical heavy-tailed generalization of the matrix-variate normal distribution. Parsimony is reached by applying the well-known eigen-decomposition of the component scale matrices, as well as by allowing the tailedness parameters of the mixture components to be tied across groups. An AECM algorithm for parameter estimation is presented. The proposed models are then fitted to simulated and real data. Comparisons with parsimonious matrix-variate normal mixtures are also provided.

Authors and Affiliations

1. Università degli Studi di Catania, Dipartimento di Economia e Impresa, Catania, Italia

   Salvatore D. Tomarchio & Antonio Punzo

2. Università Cattolica del Sacro Cuore, Dipartimento di Scienze Economiche e Sociali, Piacenza, Italia

   Luca Bagnato

Corresponding author

Correspondence to Salvatore D. Tomarchio .

Editors and Affiliations

1. Department of Economics and Management, University of Pisa, Pisa, Italy

   Nicola Salvati

2. Department of Economics and Statistics, University of Salerno, Fisciano, Salerno, Italy

   Cira Perna

3. Department of Economics and Management, University of Pisa, Pisa, Italy

   Stefano Marchetti

4. School of Mathematics and Applied Statistics, University of Wollongong, Wollongong, NSW, Australia

   Raymond Chambers

Appendix A

Let \(\ddot{\textbf{V}}=\displaystyle \sum _{g=1}^G \ddot{\textbf{V}}_g\), where \(\ddot{\textbf{V}}_g = \displaystyle \sum _{i=1}^N \ddot{z}_{ig}\ddot{w}_{ig}\left( \textbf{X}_{i}-\ddot{\textbf{M}}_g\right) \dot{\boldsymbol{\Psi }}_g^{-1}\left( \textbf{X}_{i}-\ddot{\textbf{M}}_g\right) '\). Then, we have the following updates:

• Model EII

  $$\begin{aligned} {\ddot{\lambda }} = \frac{{{\,\text{ tr }}\left\{ \ddot{\textbf{V}}\right\} }}{prN}; \end{aligned}$$

• Model VII

  $$\begin{aligned} \ddot{\lambda }_g = \frac{{{\,\text{ tr }}\left\{ \ddot{\textbf{V}}_g\right\} }}{pr \displaystyle \sum _{i=1}^N \ddot{z}_{ig}}; \end{aligned}$$

• Model EEI

  $$\begin{aligned} \ddot{\boldsymbol{\Delta }} =\frac{\text {diag}\left( \ddot{\textbf{V}}\right) }{ \left| \text {diag}\left( \ddot{\textbf{V}}\right) \right| ^\frac{1}{p}} \quad \text {and} \quad \ddot{\lambda } = \frac{\left| \text {diag}\left( \ddot{\textbf{V}}\right) \right| ^\frac{1}{p}}{rN}; \end{aligned}$$

• Model VEI

  $$\begin{aligned} \ddot{\boldsymbol{\Delta }} = \frac{\text {diag}\left( \displaystyle \sum \limits _{g=1}^G \dot{\lambda }_g^{-1}\ddot{\textbf{V}}_g\right) }{\left| \text {diag}\left( \displaystyle \sum \limits _{g=1}^G \dot{\lambda }_g^{-1}\ddot{\textbf{V}}_g\right) \right| ^\frac{1}{p}} \quad \text {and} \quad \ddot{\lambda }_g = \frac{{{\,\text{ tr }}\left\{ \ddot{\boldsymbol{\Delta }}^{-1} \ddot{\textbf{V}}_g \right\} }}{pr\displaystyle \sum _{i=1}^N \ddot{z}_{ig}}; \end{aligned}$$

• Model EVI

  $$\begin{aligned} \ddot{\boldsymbol{\Delta }}_g = \frac{\text {diag}\left( \ddot{\textbf{V}}_g\right) }{\left| \text {diag}\left( \ddot{\textbf{V}}_g\right) \right| ^\frac{1}{p}} \quad \text {and} \quad \ddot{\lambda } = \frac{\displaystyle \sum \limits _{g=1}^G\left| \text {diag}\left( \ddot{\textbf{V}}_g\right) \right| ^\frac{1}{p}}{rN}; \end{aligned}$$

• Model VVI

  $$\begin{aligned} \ddot{\boldsymbol{\Delta }}_g = \frac{\text {diag}\left( \ddot{\textbf{V}}_g\right) }{\left| \text {diag}\left( \ddot{\textbf{V}}_g\right) \right| ^\frac{1}{p}} \quad \text {and} \quad \ddot{\lambda }_g = \frac{\left| \text {diag}\left( \ddot{\textbf{V}}_g\right) \right| ^\frac{1}{p}}{r\displaystyle \sum _{i=1}^N \ddot{z}_{ig}}; \end{aligned}$$

• Model EEE

  $$\begin{aligned} \ddot{\boldsymbol{\Sigma }}= \frac{\ddot{\textbf{V}}}{rN}; \end{aligned}$$

• Model VEE

  $$\begin{aligned} \ddot{\boldsymbol{\Gamma }}\ddot{\boldsymbol{\Delta }}\ddot{\boldsymbol{\Gamma }}' = \frac{\displaystyle \sum \limits _{g=1}^G \dot{\lambda }_g^{-1}\ddot{\textbf{V}}_g}{\left| \displaystyle \sum \limits _{g=1}^G \dot{\lambda }_g^{-1}\ddot{\textbf{V}}_g\right| ^{\frac{1}{p}}} \quad \text {and} \quad \ddot{\lambda }_g = \frac{{{\,\text{ tr }}\left\{ (\ddot{\boldsymbol{\Gamma }}\ddot{\boldsymbol{\Delta }}\ddot{\boldsymbol{\Gamma }}')^{-1} \ddot{\textbf{V}}_g \right\} }}{pr\displaystyle \sum _{i=1}^N \ddot{z}_{ig}}; \end{aligned}$$

• Model EVE

  For this model, there is no analytical solution for \(\boldsymbol{\Gamma }\). Thus, an iterative minorization-maximization (MM) algorithm [2] is implemented. Specifically, the following surrogate function is defined

  $$\begin{aligned} f\left( \boldsymbol{\Gamma }\right) = \sum \limits _{g=1}^G \,\text{ tr }\left\{ \textbf{V}_g\boldsymbol{\Gamma }\boldsymbol{\Delta }_{k}^{-1}\boldsymbol{\Gamma }'\right\} \le S + \,\text{ tr }\left\{ \dot{\textbf{F}}\boldsymbol{\Gamma }\right\} , \end{aligned}$$

  where S is a constant and \(\dot{\textbf{F}} = \displaystyle \sum _{g=1}^G\left( \boldsymbol{\Delta }_{k}^{-1} \dot{\boldsymbol{\Gamma }}' \textbf{V}_g - e_g \boldsymbol{\Delta }_{k}^{-1} \dot{\boldsymbol{\Gamma }}'\right) \), with \(e_g\) being the largest eigenvalue of \(\textbf{V}_g\). The update of \(\boldsymbol{\Gamma }\) is given by \(\ddot{\boldsymbol{\Gamma }} = \dot{\textbf{G}} \dot{\textbf{H}} '\), where \(\dot{\textbf{G}}\) and \(\dot{\textbf{H}}\) are obtained from the singular value decomposition of \(\dot{\textbf{F}}\). This process is repeated until a specified convergence criterion is met and the estimate \(\ddot{\boldsymbol{\Gamma }}\) is obtained from the last iteration. Then, we obtain

  $$\begin{aligned} {\ddot{{\boldsymbol{\Delta }}}}_{g} = \frac{{\text {diag}}\left( {\ddot{{\boldsymbol{\Gamma }}}'} {\ddot{{\textbf{V}}}}_{g} {\ddot{{\boldsymbol{\Gamma }}}}\right) }{{\left| {\text {diag}}\left( {\ddot{{\boldsymbol{\Gamma }}}'}{\ddot{{\textbf{V}}}}_{g} {\ddot{{\boldsymbol{\Gamma }}}}\right) \right| ^{\frac{1}{p}}}} \quad \text {and} \quad {\ddot{\lambda }} = \frac{\sum \limits _{g=1}^{G}{\text {tr}}\left( \ddot{\boldsymbol{\Gamma }} \ddot{\boldsymbol{\Delta }}_{g}^{-1} \ddot{\boldsymbol{\Gamma }}'\ddot{\textbf{V}}_{g}\right) }{prN}; \end{aligned}$$

• Model VVE

  Similarly to the EVE case, there is no analytical solution for \(\boldsymbol{\Gamma }\), and the MM algorithm described above is implemented. Then, we have

  $$\begin{aligned} \ddot{\boldsymbol{\Delta }}_g = \frac{\text {diag}\left( \ddot{\boldsymbol{\Gamma }}' \ddot{\textbf{V}}_g \ddot{\boldsymbol{\Gamma }}\right) }{\left| \text {diag}\left( \ddot{\boldsymbol{\Gamma }}' \ddot{\textbf{V}}_g \ddot{\boldsymbol{\Gamma }}\right) \right| ^\frac{1}{p}} \quad \text {and} \quad \ddot{\lambda }_g = \frac{\left| \text {diag}\left( \ddot{\boldsymbol{\Gamma }}' \ddot{\textbf{V}}_g \ddot{\boldsymbol{\Gamma }}\right) \right| ^{\frac{1}{p}}}{r\displaystyle \sum _{i=1}^N \ddot{z}_{ig}}; \end{aligned}$$

• Model EEV

  Consider the eigen-decomposition \(\textbf{V}_g=\textbf{L}_g \textbf{D}_g \textbf{L}_g'\), with eigenvalues in the diagonal matrix \(\textbf{D}_g\) following descending order and orthogonal matrix \(\textbf{L}_g\) composed of the corresponding eigenvectors. Then, we obtain

  $$\begin{aligned} \ddot{\boldsymbol{\Gamma }}_g=\ddot{\textbf{L}}_g , \quad \ddot{\boldsymbol{\Delta }} = \frac{\displaystyle \sum \limits _{g=1}^G \ddot{\textbf{D}}_g}{\left| \displaystyle \sum \limits _{g=1}^G \ddot{\textbf{D}}_g\right| ^\frac{1}{p}} \quad \text {and} \quad \ddot{\lambda } = \frac{\left| \displaystyle \sum \limits _{g=1}^G \ddot{\textbf{D}}_g\right| ^\frac{1}{p}}{rN}; \end{aligned}$$

• Model VEV

  By using the same algorithm applied for the EEV model, we have

  $$\begin{aligned} \ddot{\boldsymbol{\Gamma }}_g=\ddot{\textbf{L}}_g , \quad \ddot{\boldsymbol{\Delta }} = \frac{\displaystyle \sum \limits _{g=1}^G \lambda _g^{-1} \ddot{\textbf{D}}_g}{\left| \displaystyle \sum \limits _{g=1}^G \lambda _g^{-1} \ddot{\textbf{D}}_g\right| ^\frac{1}{p}} \quad \text {and} \quad \ddot{\lambda }_g = \frac{{{\,\text{ tr }}\left\{ \ddot{\textbf{D}}_g \ddot{\boldsymbol{\Delta }}^{-1} \right\} }}{pr\displaystyle \sum _{i=1}^N \ddot{z}_{ig}}; \end{aligned}$$

• Model EVV

  $$\begin{aligned} \ddot{\boldsymbol{\Gamma }}_g\ddot{\boldsymbol{\Delta }}_g\ddot{\boldsymbol{\Gamma }}_g' = \frac{\ddot{\textbf{V}}_g}{\left| \ddot{\textbf{V}}_g\right| ^{\frac{1}{p}}} \quad \text {and} \quad \ddot{\lambda } = \frac{\displaystyle \sum \limits _{g=1}^G \left| \ddot{\textbf{V}}_g\right| ^{\frac{1}{p}}}{rN}; \end{aligned}$$

• Model VVV

  $$\begin{aligned} \ddot{\boldsymbol{\Sigma }}_g = \frac{\ddot{\textbf{V}}_g}{r\displaystyle \sum _{i=1}^N \ddot{z}_{ig}}. \end{aligned}$$

Appendix B

Let \(\ddot{\textbf{W}}=\displaystyle \sum _{g=1}^G \ddot{\textbf{W}}_g\), where \(\ddot{\textbf{W}}_g = \displaystyle \sum _{i=1}^N \ddot{z}_{ig}\ddot{w}_{ig}\left( \textbf{X}_{i}-\ddot{\textbf{M}}_{g}\right) '\ddot{\boldsymbol{\Sigma }}_g^{-1}\left( \textbf{X}_{i}-\ddot{\textbf{M}}_{g}\right) \). With the exclusion of the II model, for which no parameters need to be estimated, we have the following updates:

• Model EI

  $$\begin{aligned} \ddot{\boldsymbol{\Delta }} = \frac{\text {diag}\left( \ddot{\textbf{W}}\right) }{\left| \text {diag}\left( \ddot{\textbf{W}}\right) \right| ^\frac{1}{r}}; \end{aligned}$$

• Model VI

  $$\begin{aligned} \ddot{\boldsymbol{\Delta }}_g = \frac{\text {diag}\left( \ddot{\textbf{W}}_g\right) }{\left| \text {diag}\left( \ddot{\textbf{W}}_g\right) \right| ^\frac{1}{r}}; \end{aligned}$$

• Model EE

  $$\begin{aligned} \ddot{\boldsymbol{\Psi }} = \frac{\ddot{\textbf{W}}}{\left| \ddot{\textbf{W}} \right| ^\frac{1}{r}}; \end{aligned}$$

• Model VE

  As for the EVE and VVE models, there is no analytical solution for \(\boldsymbol{\Gamma }\) and a MM algorithm of the type described for the EVE model is implemented, after replacing \(\textbf{V}\) with \(\textbf{W}\). Then, we have

  $$\begin{aligned} \ddot{\boldsymbol{\Delta }}_g= \frac{\text {diag}\left( \ddot{\boldsymbol{\Gamma }}' \ddot{\textbf{W}}_g \ddot{\boldsymbol{\Gamma }}\right) }{\left| \text {diag}\left( \ddot{\boldsymbol{\Gamma }}' \ddot{\textbf{W}}_g \ddot{\boldsymbol{\Gamma }}\right) \right| ^\frac{1}{r}}; \end{aligned}$$

• Model EV

  By using the same approach of the EEV and VEV models, after replacing \(\ddot{\textbf{V}}\) with \(\ddot{\textbf{W}}\), we have

  $$\begin{aligned} \ddot{\boldsymbol{\Gamma }}_g=\ddot{\textbf{L}}_g \quad \text {and} \quad \ddot{\boldsymbol{\Delta }} = \frac{\displaystyle \sum \limits _{g=1}^G \ddot{\textbf{D}}_g}{\left| \displaystyle \sum \limits _{g=1}^G \ddot{\textbf{D}}_g\right| ^\frac{1}{r}}; \end{aligned}$$

• Model VV

  $$\begin{aligned} \ddot{\boldsymbol{\Psi }}_g = \frac{\ddot{\textbf{W}}_g}{\left| \ddot{\textbf{W}}_g\right| ^\frac{1}{r}}. \end{aligned}$$

Reprints and permissions

© 2022 The Author(s), under exclusive license to Springer Nature Switzerland AG

Policies and ethics