\section{Introduction}

Superconductivity occurs as a result of the second-order phase
transition, accompanied by a sudden appearance of the
superconducting order parameter below the critical temperature
$T_c$ and the upper critical field $H_{c2}$ \cite{SaintJames}.
This should lead to an abrupt vanishing of resistance. However, in
reality resistive transitions are always broadened, especially in
magnetic field. This is usually ascribed to flux-flow phenomenon
caused by motion of Abrikosov vortices \cite{BardeenStephen_1965}.
Broadening can also be caused by spatial inhomogeneity (e.g.
variation of $T_c$), or superconducting fluctuations
\cite{Varlamov_book,Finkelstein_2012,Glatz_2014}, particularly for
high temperature superconductors. Finally, surface
superconductivity (SSC) may survive up to a significantly higher
field $H_{c3}\simeq 1.69 H_{c2}$ than bulk superconductivity
\cite{SaintJames}, which can also smear the superconducting
transition. Although SSC is quite profound in polished clean
superconductors
\cite{Bellau_1966,Park_2004,Casalbuoni_2005,Das_2008}, it is
usually ignored for disordered, polycrysstaline films because SSC
is considered to be very sensitive to the quality of the surface
(e.g., surface passivation \cite{Casalbuoni_2005} and order
parameter suppression \cite{SaintJames}), surface roughness
\cite{Bellau_1966,Casalbuoni_2005} and surface scattering
\cite{Agterberg_1996}.

The broadening is detrimental for superconducting devices such as
transition-edge sensors \cite{TES1,TES2} and resonators
\cite{Casalbuoni_2005}. Presence of several mechanisms makes the
interpretation of broadening ambiguous. The lack of understanding
does not allow confident extraction of fundamental parameters of
superconductors, such as $H_{c2}$, because it is unclear which
point at the transition curve corresponds to $H=H_{c2}$. Arbitrary
criteria, such as 10, 50 or 90 $\%$ of the normal state resistance
$R_n$, are commonly applied which apparently does not work for
high-$T_c$ superconductors with very broad transitions
\cite{Zavaritsky}. Therefore, clarification of the mechanism of
broadening is important both for fundamental and applied research
on superconductors.

In this work we study the origin of broadening of superconducting
transitions in sputtered Nb films. We perform simultaneous
tunneling spectroscopy and transport measurements, which allow
unambiguous ascription of $H_{c2}$ to the bottom of resistive
transition $R(H_{c2})/R_n \sim 0$. The top of transition
corresponds to $\sim 1.7$ times higher fields, which is close to
the third critical field $H_{c3}$ for destruction of surface
superconductivity. The two-dimensional (2D) nature of SSC at
$H_{c2}<H<H_{c3}$ is confirmed by observation of cusp-like angular
dependence of magnetoresistance. Thus we conclude that surface
superconductivity, rather than flux-flow, inhomogeneity or
fluctuations, is the primary cause of broadening of
superconducting transitions in magnetic field. Our data indicates
that surface superconductivity is remarkably robust even in
disordered polycrystalline films and, surprisingly, even in
perpendicular magnetic fields.

\begin{figure*}[t]
    \centering
    \includegraphics[width=0.9\textwidth]{Fig1_RvsTvsH.pdf}
    \caption{(Color online). Resistive transitions of a 150 nm thick Nb film. (a) Temperature dependence of
    the resistance in zero magnetic field. (b) and (c) Field dependencies of
    resistances at $T=1.8$ K in fields (b) perpendicular and (c)
    parallel to the film. Black and red circles mark the upper
    critical field $H_{c2}$ and the field of the onset of resistive
    transition, which coincides with the critical field of surface
    superconductivity $H_{on} \sim H_{c3}$.
    }
    \label{fig:fig1}
\end{figure*}

\section{Experimental}

The studied sample contains several Nb/Al-AlOx/Nb tunnel junctions
with sputtered Nb electrodes of thicknesses $d=150$ and 50 nm.
Junction characteristics in perpendicular fields were reported in
Ref. \cite{MR}. Due to different thicknesses, electrodes have
slightly different $T_c$ of 9.2 and 8.8 K. Parameters extracted
from tunneling characteristics are determined by the thinner
electrode, while transport measurements are made at the thicker
electrode. This explains a minor difference in $H_{c2}$ values
obtained by those techniques. Measurements are performed in a
gas-flow $^4$He cryostat with a superconducting solenoid.
Samples are mounted on a rotatable holder with the alignment
accuracy $\sim 0.02^{\circ}$. Details of the setup can be found
elsewhere \cite{MR}.

\section{Results}

In Figures \ref{fig:fig1} (a-c) we show superconducting
transitions of a 150 nm thick Nb film: (a) $R(T)$ in zero field
and $R(H)$ at $T=1.8$ K for field perpendicular (b) and parallel
(c) to the film. It is seen that at zero field the $R(T)$
transition is very sharp and does not show any extended
fluctuation region or inhomogeneity. However, $R(H)$ transitions
are quite broad. Interestingly, $R(H)$ is broader when the field
is parallel to the film. This confronts interpretation of
broadening in terms of vortex motion because the driving Lorentz
force is most effective in perpendicular and vanishes in parallel
field. Therefore, this broadening is not consistent with either
flux-flow, inhomogeneity or fluctuation mechanisms.

\begin{figure*}[th]
    \centering
    \includegraphics[width=1.0\textwidth]{Fig2_Hperp_T196ExpTheor.pdf}
    \caption{(Color online). Comparison of theoretically calculated (a-c) and experimentally measured (d-f)
    characteristics of Nb/AlAlOx/Nb tunnel junctions. (a) and (d) Temperature dependence of $I$-$V$ characteristics at zero field.
    (b) and (e) Field dependence of $I$-$V$ characteristics for field perpendicular to the junction/films at $T\simeq 2$ K.
    (c) and (f) The corresponding differential conductances for curves from panels (b) and (e).
    The field scale in (f) is normalized by the upper critical
    field $H_{c2}^{\perp}=1.61$ T, which is obtained as a single
    scaling factor for all the curves at different fields. Data
    from Ref. \cite{MR}.
    }
    \label{fig:fig2}
\end{figure*}

\subsection{Determination of $H_{c2}$ from tunneling spectroscopy}

In order to analyze surface superconductivity scenario, first of
all, it is necessary to determine bulk $H_{c2}$. For this we
perform magneto-tunneling spectroscopy. Figure \ref{fig:fig2}
represents a comparison of theoretically calculated (a-c) and
experimentally measured tunneling characteristics of our
Nb/AlAlOx/Nb junction (data from Ref. \cite{MR}). Details of
calculations are described in Ref. \cite{MR}). Panels (a) and (d)
show temperature dependencies of $I$-$V$ characteristics at zero
field. Panels (b) and (e) show field dependence of $I$-$V$
characteristics for field perpendicular to the junction/films at
$T\simeq 2$ K. Panels (c) and (f) show the corresponding
differential conductances for $I$-$V$ curves from panels (b) and
(e). There is a good quantitative agreement between theoretical
and experimental characteristics. The main spectroscopic features
are the sharp sum-gap peak at $V_p=2\Delta/e$, where $\Delta$ is
the superconducting gap and the suppressed quasiparticle current
and conductance below the sum-gap voltage. With increasing field
the peak is decreasing in height and is moving to lower voltages.
Simultaneously the sub-gap conductance is increased. All this is
due to suppression of the superconducting gap by magnetic field.
The extent of suppression depends solely on $H/H_{c2}$. Above
$H_{c2}$ the superconducting gap vanishes and the $I$-$V$ becomes
linear (Ohmic). Thus, the ratio $H/H_{c2}$ uniquely determines the
shape of tunneling characteristics in magnetic field. Therefore,
as discussed in Ref. \cite{MR}, the ratio $H/H_{c2}$ can be
uniquely determined from analysis of the shape of tunneling
characteristics. In Fig. \ref{fig:fig2} (f) the field is
normalized by thus obtained $H_{c2}^{\perp}=1.61$ T. We emphasize
that this value is obtained as a single fitting parameter for the
whole set of d$I/$d$V(V)$ characteristics at different $H$. This
removes the uncertainty in determination of $H_{c2}$.

Fig. \ref{fig:fig3} (a) shows a set of tunneling d$I$/d$V(V)$
characteristics of a Nb/AlOx/Nb junction at $T=1.8$ K in fields
parallel to Nb films. From comparison of Figs. \ref{fig:fig2} and
\ref{fig:fig3} (a), it can be seen that the influence of magnetic
field is qualitative similar both for parallel and perpendicular
field orientations. In Ref. \cite{MR} it was shown that the peak
height and the peak voltage exhibit universal almost
$T$-independent quasi-linear scaling as a function of
$H/H_{c2}(T)$. Fig. \ref{fig:fig3} (b) and (c) demonstrate such a
scaling at different temperatures for field parallel to Nb films.
Dashed and dashed-dotted lines in Fig. \ref{fig:fig3} (b)
represent theoretical results from Ref. \cite{MR} for $T=1.96$ and
4.7 K, correspondingly. The overall quality of scaling is quite
good, which allows confident extraction of $H_{c2}(T)$.
Thus obtained $H_{c2}$ is unambiguous because it is deduced as a
single fitting parameter for the whole set of d$I/$d$V(V)$
characteristics at fixed $T$ for different $H$.

\begin{figure*}[t]
    \centering
    \includegraphics[width=0.9\textwidth]{Fig3_Hll_Nb_And_Theory.pdf}
    \caption{(Color online). (a) Differential conductances of a Nb/AlOx/Nb junction
    (in a semi-logarithmic scale) at parallel to the films magnetic fields and $T=1.8$ K.
(b) and (c) Scaling of the sum-gap peak resistance (b) and voltage
(c) as a function of $H/H_{c2}(T)$ at different temperatures and
parallel fields. Dashed and dashed-dotted lines in (b) represent
theoretical curves at $T=1.96$ K and 4.7 K, respectively. (d) and
(e) Black squares represent upper critical
    fields perpendicular (d) and parallel (e) to the films,
    obtained from the scaling of magneto-tunnelling
    characteristics. Dashed lines represent the expected third critical field for surface
    superconductivity $H_{c3}=1.69 H_{c2}$. Red circles mark the top onset of the resistive
    transition.(f) Angular dependence of resistance of Nb electrodes at fields slightly
    below and above $\mu_0 H_{c2}(\parallel)=2.52$ T. A cusp-like feature at $H>H_{c2}$ indicates occurrence of the
    two-dimensional surface superconductivity.
    }
    \label{fig:fig3}
\end{figure*}

Squares in Figs. \ref{fig:fig3} (d) and (e) represent obtained
$H_{c2}(T)$ dependencies for perpendicular and parallel field
orientations, respectively. Using the relation
$H_{c2}^{\perp}=\Phi_0/2\pi\xi^2$ we calculate the coherence
length $\xi_0\simeq 14$ nm. This small value indicates that the
film is in the dirty limit with a very short scattering length due
to a disordered film structure with nm-scale crystallites. Thus,
the studied Nb film $d=150$ nm is an order of magnitude thicker
than $\xi_0$.
This leads to an important for a further discussion conclusion
that our films are bulk three-dimensional (3D) superconductors
practically in the whole temperature range $T<T_c$. Red circles in
Figs. \ref{fig:fig1} (b) and (c) represent top onsets
$R(H_{on})\simeq R_n$ of resistive transitions.
Dashed blue lines correspond to $H_{c3} =1.69 H_{c2}$ expected for
surface superconductivity, which is close to the onset field.
Remarkably this is true even for the perpendicular field
orientation when SSC in the uniform case is not expected
\cite{SaintJames}.

\begin{figure*}[t]
    \centering
    \includegraphics[width=0.9\textwidth]{Fig4_BiasDep.pdf}
    \caption{(Color online). (a) Bias dependence
    of $R(\Theta)$ at $T=1.8$ K and $H=2.81$ T. (b) Deduced
    current-voltage characteristics at different angles from the data in panel
    (a). Dashed line represents the normal state $I$-$V$.
    (c) Simulated angular dependence of $R(\Theta)$ in the
    surface superconductivity model, taking into account
    non-linearity of $I$-$V$ characteristics in the vicinity of
    critical current. Solid lines represent calculated of $R(\Theta)$
    at different bias. Dashed curve represents the standard
    flux-flow dependence $H/H_{c3}$.
    }
    \label{fig:fig4}
\end{figure*}

\subsection{Angular dependence of magnetoresistance}

The 2D nature of SSC should be reflected in a cusp-like angular
dependence of $H_{c3}$, given by the equation
\cite{Yamafuji_1966}:
\begin{equation}\label{Hc3_Q}
\left[\frac{H_{c3}(\Theta)}{H_{c3}^\parallel}\cos\Theta\right]^2A(\Theta)
+ \left|\frac{H_{c3}(\Theta)}{H_{c3}^\perp}\sin\Theta\right| =1,
\end{equation}
where $A(\Theta)=1+(1-\sin \Theta)\tan\Theta$. It is only slightly
different from Tinkham's 2D result with $A(\Theta)=1$.

Fig. \ref{fig:fig3} (f) shows angular-dependencies of
magnetoresistance at $T=1.8$ K and at different fields below and
above $\mu_0 H_{c2}^{\parallel} = 2.52$ T. Zero angle $\Theta=0$
corresponds to field parallel to the film. It is seen that at 2.5
T very slightly below $H_{c2}$, $R(\Theta)$ is flat at $\Theta=0$,
which is characteristic for 3D bulk Nb. However, at
$H>H_{c2}^{\parallel}$ angular dependence acquires a 2D cusp.
Since the film thickness is much larger than $\xi$, the observed
2D behavior at low $T$ may originate solely from SSC.

\subsection{Non-linear bias dependence}

The sheet surface critical current (in A/cm) is
\cite{Abrikosov_1965} :
\begin{equation}\label{JcH}
I_c \simeq \frac{5}{2\sqrt{3\pi}}\frac{H_c}{\kappa}
\left(1-\frac{H}{H_{c3}}\right)^{3/2}.
\end{equation}
Here $H_c$ is the thermodynamic critical field (in Oe) and
$\kappa$ is the Ginzburg-Landau parameter. Typically such $I_c$ is
in the range from few to few tens of A/m. For our films with
$\kappa \gg 1$ and the width of few microns the $I_c$ is in the
$\mu$A range, comparable to the probe current. Therefore, the
results do depend on the bias, as illustrated in Fig.
\ref{fig:fig4} (a). This is due to strong non-linearity of
current-voltage characteristics at $I\sim I_c$, as demonstrated in
Fig. \ref{fig:fig4} (b). In order to demonstrate how such the
non-linearity affects experimental characteristics we consider a
standard shape of $I$-$V$:
\begin{equation}\label{IV}
V=R_n\sqrt{I^2-I_c^2}.
\end{equation}
Together with Eq. (\ref{JcH}) it yields
\begin{equation}\label{R_Q}
1-\left[\frac{R(\Theta)}{R_n}\right]^2=\frac{I_c^2}{I^2}\left[1-\frac{H}{H_{c3}(\Theta)}\right]^{\nu}.
\end{equation}
The exponent $\nu$ depends on the shape of the $I$-$V$ and the
quality of the surface \cite{Bellau_1966}. For the case of Eq.
(\ref{IV}) it is $\nu=3$.

\begin{figure*}[t]
    \centering
    \includegraphics[width=0.6\textwidth]{Fig5_NbFluct.pdf}
    \caption{(Color online). (a) and (b) Excess conductance vs. magnetic field at $T=1.8$ K plotted in (a) a double-logarithmic and (b)
    a semi-logarithmic scale.
    Quasi-exponential decay $\Delta S(H)$ is seen. (c) The high-field part of excess
    conductance. It reveals fluctuation contribution, which is
    positive for parallel, but negative for perpendicular field
    orientation. (d) Analysis of power-law scaling expected for surface superconductivity
    according to Eq. (\ref{R_Q}). Dashed line represents a power-law with the exponent $\nu =4.8$.
    }
    \label{fig:fig5}
\end{figure*}

In Fig. \ref{fig:fig4} (c) we show $R(\Theta)$ curves for the SSC
model calculated from Eqs. (\ref{Hc3_Q}) and (\ref{R_Q}) for
$H=1.1H_{c2}$ at different bias. Calculations are made for
$H_{c3}^{\perp}=H_{c2}$ and $H_{c3}^{\parallel}=1.69 H_{c2}$ and
for $I_c(\Theta)=$const. For comparison we also show flux-flow
type dependence $R/R_n = H/H_{c3}$. It is seen that the cusp in
the SSC model is much sharper, primarily due to non-linearity of
the $I$-$V$. Overall behavior is similar to experimental data from
Fig. \ref{fig:fig4} (a), even though in experiment a very sharp
cusp at $\Theta=0$ survives up to much higher current. The
difference is due to an oversimplified assumption of
angular-independent $I_c(\Theta)=$const, used in calculations. In
reality $I_c(\Theta)$ has a sharp maximum at $\Theta=0$ because
the Lorenz force vanishes as $\sin(\Theta)$. It is possible to get
a better fit using a more realistic $I_c(\Theta)$, but we don't
want to go in to more complicated modelling because the main
purpose of calculations was to demonstrate how non-linearity of
$I$-$V$'s leads to a much sharper (compared to a simple 2D
flux-flow model) cusp in $R(\Theta)$.

\subsection{Analysis of fluctuation contribution}

Finally we discuss fluctuation contribution to the broadening of
resistive transition. In Fig. \ref{fig:fig5} (a) we show
normalized excess conductance $\Delta S= (1/R-1/R_n)R_n$ for the
data from Figs. \ref{fig:fig1} (b) and (c) in a double-logarithmic
scale. Such graph is usually used for analysis of fluctuation
contribution to conductivity. Dashed line shows $\Delta S \propto
(H-H_{c2})^{-1}$ dependence expected for 2D fluctuations
\cite{Varlamov_book}. It is seen that although there is a narrow
range of fields close to $H_{c2}$ with similar behavior, the
overall agreement is poor. In Fig. \ref{fig:fig5} (b) we replot
the same data in a semi-logarithmic scale. It is seen that $\Delta
S$ decays quasi-exponentially with increasing field at
approximately the same rate for both field orientations. A similar
exponentially decay versus both $T$ and $H$ has been reported for
high-$T_c$ cuprates \cite{MR,Katterwe_2014,Alloul2011,Luo2014}.
Here we demonstrate that it is generic also for conventional
superconductors. Such behavior is not expected for fluctuations
\cite{Varlamov_book,Finkelstein_2012,Glatz_2014} and we argue that
it is rather a signature of SSC.

It is possible to distinguish fluctuation contribution from
non-fluctuating SSC. SSC always leads to excess conductance, but
fluctuation contribution to magnetoresistance can be both positive
and negative \cite{Varlamov_book,Finkelstein_2012,Glatz_2014}. In
particular, at low $T$ and at field perpendicular to the film the
density-of-state contribution to fluctuations leads to excess
resistance, rather than excess conductance \cite{Glatz_2014}. We
clearly see such a contribution in our data. In Fig.
\ref{fig:fig5} (c) we show high-field part of excess conductance
in the linear scale. It is seen that in parallel field there is
always an excess conductance $\Delta S >0$, which rapidly
decreases upon approaching the surface critical field
$H_{c3}\simeq 1.7 H_{c2}^{\parallel}$, but never really vanishes.
The remaining tail is a signature of fluctuations that persist at
any field. For perpendicular field, $\Delta S$ at high fields
becomes {\it negative}, which is consistent with theoretical
expectations for fluctuation contribution at $T\ll T_c$
\cite{Glatz_2014}.

In Fig. \ref{fig:fig5} (d) we check the power-law scaling
suggested by Eq. (\ref{R_Q}) for SSC. It is seen that there is a
good scaling in a broad field range, although extraction of the
exponent $\nu$ is not very confident because it depends on the
chose of $H_{c3}$. The dashed line corresponds to $\nu=4.8$. Upon
approaching $H_{c3}$, deviations with opposite signs for parallel
and perpendicular field orientations appear, signalizing
fluctuation contributions. This indicates that SSC makes a
dominant contribution to excess conductivity at $H_{c2}<H\lesssim
0.8 H_{c3}$, while fluctuation contribution starts to become
significant only upon weakening of SSC at $H>0.8H_{c3}$ and takes
over completely at $H>H_{c3}$.

\section{Discussion}

Our results suggest that surface superconductivity is the primary
cause of broadening of superconducting transition in magnetic
field. As indicated in Figs. \ref{fig:fig1} (b) and (c),
$H=H_{c2}$ corresponds to the bottom of transition, consistent
with earlier studies \cite{Burger_1965}, and $H_{c3}$ to the top
of the resistive transition. Thus the full width of the transition
is dominated by SSC. Although SSC is well known for carefully
polished single crystals \cite{Park_2004,Das_2008}, it is usually
considered to be insignificant for disordered, rough or
inhomogeneous superconducting films because of its assumed
fragility and sensitivity to surface conditions
\cite{SaintJames,Casalbuoni_2005,Bellau_1966,Agterberg_1996}.
Therefore, observation of a very robust SSC in our strongly
disordered polycrystalline films is rather surprising, especially
for field perpendicular to the film. In perfectly uniform films
SSC should not occur at perpendicular field orientation
\cite{SaintJames,Abrikosov_1965}. Yet, SSC in perpendicular fields
has been directly visualized by scanning laser microscopy for
similar films \cite{Werner_2013} and also reported for some
layered supercondcutors \cite{LevyBertrand_2016} and sintered
polycrystalline MgB$_2$ samples \cite{Tsindlekht_2006}. Presumably
it is the polycrystallinity of our films that allows SSC at grain
boundaries even in perpendicular fields. Thus we conclude that
surface superconductivity is a robust phenomenon that should be
carefully considered in analysis of data close to superconducting
transition.

\begin{acknowledgements}

The work was supported by the Swedish Research Council Grant No.
621-2014-4314 and the Swedish Foundation for International
Cooperation in Research and Higher Education Grant No.
IG2013-5453.
\end{acknowledgements}


\section{Introduction}
Let $\mathbf{x}$ be a $J$-variate random variable following a Gaussian mixture model (GMM) with density
\begin{equation}
\label{md}
h(\mathbf{x};\theta)=\sum_{g=1}^G p_g \phi(\mathbf{x};\boldsymbol{\mu}_g,\boldsymbol{\Sigma}_g),
\end{equation}
where the $p_g$'s are the mixing proportions, with $p_g > 0$ $\forall g$ and $\sum_g p_g=1$, and
the component $\phi(\mathbf{x};\boldsymbol{\mu}_g,\boldsymbol{\Sigma}_g)$
represents the density of a $J$-variate normal distribution
with mean vector $\boldsymbol{\mu}_g$ and covariance matrix $\boldsymbol{\Sigma}_g$; furthermore, let us indicate the set of model parameters with $\theta=\left\{p_g, \boldsymbol{\mu}_g,\boldsymbol{\Sigma}_g \right\}_G=\left\{p_1,\dots p_G, \boldsymbol{\mu}_1,\dots,\boldsymbol{\mu}_G,\boldsymbol{\Sigma}_1,\dots,\boldsymbol{\Sigma}_G \right\} \in \mathbb{R}^{G(1+J+J^{2})},$
and the parameter space with
\begin{equation}
\label{ps}
\Theta=\left\{\theta \in \mathbb{R}^{G(1+J+J^{2})} :\sum_g p_g=1, p_g > 0, \boldsymbol{\Sigma}_g \succ \boldsymbol{0},
g=1, \dots, G\right\},
\end{equation}
where the symbol $\succ$ refers to L\"{o}wner ordering on symmetric matrices and, in this case, is equivalent to requiring that $\boldsymbol{\Sigma}_g$ be positive definite.
The GMM is frequently used to classify a sample of observations. The idea is to consider the sample as drawn from a heterogeneous population where each sub-population is described by one component of the mixture. In other terms, each observation is assumed to come from one of the $G$ different groups characterized by the mixture components.
The observations are classified into the groups by computing the posterior probabilities
\begin{equation}
p(g|\mathbf{x}) = \frac{p_g \phi(\mathbf{x};\boldsymbol{\mu}_g,\boldsymbol{\Sigma}_g)}
{\sum_h p_h \phi(\mathbf{x};\boldsymbol{\mu}_h,\boldsymbol{\Sigma}_h)},
\end{equation}
and assigning each observation to the group with the largest posterior probability.

The parameters of the GMM are generally unknown and estimated from the data. Given a sample of i.i.d. observations $\left\lbrace \mathbf{x}_i \right\rbrace_n= \left\{ \mathbf{x}_1, \mathbf{x}_2, \dots \mathbf{x}_n \right\}$, the estimation is usually done by maximizing the likelihood
\begin{equation}
\label{eq:like}
L\left(\theta; \left\lbrace \mathbf{x}_i \right\rbrace_n \right)=\prod_{i=1}^n  \left(\sum_{g=1}^G p_g \phi(\mathbf{x}_i;\boldsymbol{\mu}_g,\boldsymbol{\Sigma}_g)\right).
\end{equation}
The likelihood in Equation \eqref{eq:like} is known to be unbounded and it is cursed by the presence of several local maxima. As a consequence, the EM algorithm may fail to converge, leading to such degenerate solutions. To face degeneracy, several methods have been proposed by the literature in which constraints or penalties are added to the log-likelihood. Their main objective is to keep the eigenvalues of the class conditional covariance matrices bounded away from zero.

This paper considers the sufficient condition formulated by Ingrassia (2004) such that Hathaway's (1985) constraints hold: we propose a generalization that enforces the equivariance with respect to linear affine transformations of the data. The idea is to shrink the class conditional covariance matrices towards a pre-specified matrix $\boldsymbol{\Psi}.$ We investigate possible data-driven methods for choosing the matrix $\boldsymbol{\Psi},$ when {\em a priori} information on the group-specific covariance structure is not available, and we let the data determine the optimal amount of shrinkage. The equivariance property the method possesses is a key feature for twofold reasons. First, it means that irrespective of the kind of standardization performed on the data, the final clustering will be the same - provided that $\boldsymbol{\Psi}$ is transformed accordingly. Second, whatever the scale of the data is as they come in, there will be no \emph{best} pre-processing of the data ensuring a \emph{better} result, as the final clustering is not affected by changes in scale.

The plan of the paper is the following. Section \ref{degen} gives insights on the notion of degeneracy for multivariate GMM, and Section \ref{remdegen} reviews some of the workarounds proposed by the existing literature. In Section \ref{invar} we state the property of equivariance of GMM and we show, in Section \ref{constr1}, that the property holds in the constrained approach of Hathaway (1985), whereas it does not hold in the sufficient condition provided by Ingrassia (2004). In Section \ref{constr2} we illustrate how these constraints can be generalized, by introducing the matrix $\boldsymbol{\Psi},$ to become equivariant under linear affine transformations of the data, provided that $\boldsymbol{\Psi}$ is transformed accordingly, and how their configuration can be tuned from the data (Section \ref{crossvalid}). Section \ref{alg} summarizes the algorithm. The proposal is evaluated through a simulation study (Section \ref{simulation}) and an empirical application (Section \ref{wineapp}). Section \ref{concl} concludes with a final discussion.
\section{The issue of degeneracy in GMM}\label{degen}
In the univariate case, the likelihood function increases without bound if some variances tend to zero and the corresponding component's mean coincides with a sample observation (Kiefer and Wolfowitz, 1956; Day, 1969). Biernacki and Chr\'etien (2003) showed that if mixture parameters are close to a degenerate solution, then the EM is attracted by it and the divergence is extremely fast. Although Kiefer (1978) proved that maximum likelihood does not fail, as there exists a local maximizer strongly consistent and asymptotically efficient, several local maximizers can exist for a given sample. That is, some local maximizers are spurious, i.e. with a high likelihood but of little practical use because highly biased. They are characterized by some component variances and mixing proportions very small relative to the others (Day, 1969; McLachlan and Peel, 2000). Detecting the desired solution, among the many available, can therefore be a complicated task.

The same problems hold in the multivariate case (as an example, see Ingrassia and Rocci, 2011, for an extension of Biernacki and Chr\'etien, 2003), with additional complications. To notice how unboundedness is caused, first of all let us express the density of the $i$-th observation on the $g$-th component as follows
\begin{equation}\label{Lispec}
\phi\left(\mathbf{x}_i;\boldsymbol{\mu}_g, \boldsymbol{\Sigma}_g \right) = \frac{1}{\sqrt{2\pi \prod_{j=1}^{J}\lambda_{jg}}} \exp\left\{ -\frac{1}{2} (\mathbf{x}_{i} - \boldsymbol{\mu}_{g})'\mathbf{Q}_{g}\mathbf{L}_{g}^{-1}\mathbf{Q}_{g}'(\mathbf{x}_{i} - \boldsymbol{\mu}_{g}) \right\},
\end{equation}
where $\mathbf{Q}_{g}$ is the square $J \times J$ matrix whose $j-$th column is the eigenvector $\mathbf{q}_{jg}$ of $\mathbf{\Sigma}_{g},$ and $\mathbf{L}_{g}$ is the diagonal matrix whose diagonal elements are the corresponding eigenvalues $\{\lambda_{jg}\}_{J},$ ordered such that $\lambda_{1g} \geq \dots \geq \lambda_{Jg}.$ Equation \eqref{Lispec} can be rewritten as
\begin{align*}
\phi\left(\mathbf{x}_i;\boldsymbol{\mu}_g, \boldsymbol{\Sigma}_g \right) &= \frac{1}{\sqrt{2\pi \prod_{j=1}^{J}\lambda_{jg}}} \exp\left\{ -\frac{1}{2} (\mathbf{x}_{i} - \boldsymbol{\mu}_{g})'(\sum_{j=1}^{J} \lambda_{jg}^{-1} \mathbf{q}_{jg} \mathbf{q}_{jg}')(\mathbf{x}_{i} - \boldsymbol{\mu}_{g}) \right\} \\
&= \frac{1}{\sqrt{2\pi \prod_{j=1}^{J}\lambda_{jg}}} \exp\left\{ -\frac{1}{2} \sum_{j=1}^{J}\lambda_{jg}^{-1}[(\mathbf{x}_{i} - \boldsymbol{\mu}_{g})'\mathbf{q}_{jg}][\mathbf{q}_{jg}'(\mathbf{x}_{i} - \boldsymbol{\mu}_{g})] \right\} \\
&= \frac{1}{\sqrt{2\pi \prod_{j=1}^{J}\lambda_{jg}}} \exp\left\{ -\frac{1}{2} \sum_{j=1}^{J}\lambda_{jg}^{-1}[(\mathbf{x}_{i} - \boldsymbol{\mu}_{g})'\mathbf{q}_{jg}]^{2} \right\}. \numberthis \label{eq:Li2parts}
\end{align*}
As Policiello (1981) argued, the likelihood in Equation \eqref{eq:like} can be written as the sum of non negative terms.
Among them, it is possible to isolate the product of the density of the $i$-th observation on
the $g$-th component - Equation \eqref{eq:Li2parts} - and the densities of the other observations on the other components and the corresponding mixing proportions. If observation $i$ is such that $\mathbf{x}_{i}'\mathbf{q}_{Jg} - \boldsymbol{\mu}_{g}'\mathbf{q}_{Jg} = 0,$ then, as $\lambda_{Jg} \rightarrow 0,$ there would be no exponential term involving $\lambda_{Jg}$ who can attenuate the effect of $\frac{1}{\sqrt{2\pi \prod_{j=1}^{J}\lambda_{jg}}} \rightarrow \infty.$ In words, the sample likelihood diverges when in one component the covariance matrix is close to singularity and the projection of the component's mean on the eigenvector corresponding to the smallest eigenvalue coincides with the projection of one of the observations on the same eigenvector.

\section{Remedies to degeneracy}\label{remdegen}
The easiest way to handle degeneracy is to initialize the EM algorithm from several starting points until a local maximum is found (Biernacki and Chr\'etien, 2003). McLachlan and Peel (2000) proposed monitoring the local maximizers by inspecting the relative size of the estimated mixing proportions and component variances. This leads, in practice, to performing maximum likelihood estimation by looking for the correct local maximum and discarding those that seem to be spurious.

Further methods exploit constraints on the covariance matrices. This approach is based on the seminal work of Hathaway (1985), where he studied how to avoid the divergence of the likelihood in the univariate case by imposing a lower bound, say $c$, to the ratios of the scale parameters. In this way the variances cannot be arbitrarily different. Hathaway proved the boundedness of the likelihood  and the consistency of the ML estimator under such constraints. In the multivariate case, the lower bound is imposed on the generalized eigenvalues of each pair of covariance matrices and the ML estimator results to be equivariant under linear affine transformations of the data. This implies that, as in the unconstrained case, if the data are linearly transformed, the estimated posterior probabilities do not change and the clustering remains unaltered (see Sections \ref{invar} and \ref{constr1}).

An important issue is the choice of the constant $c,$ which controls the strength of the constraints. In the context of univariate mixtures of Gaussians or linear regression models, some authors have shown that the maximum likelihood constrained estimator is consistent if $c$ decreases to zero at a certain rate as the sample size increases to infinity (e.g. Tanaka and Takemura (2006), Tan et al. (2007), Xu et al. (2010)). Nevertheless, finite-sample sensible choice of $c$ is still an open issue.

Hathaway's constraints are very difficult to apply within iterative procedures like the EM algorithm. To solve this problem, Ingrassia (2004) proposed to simplify the constraints by putting bounds on the eigenvalues of the covariance matrices. Although putting lower bounds on the group conditional covariance matrices was already common practice, Ingrassia (2004) found a way to reconcile Hathaway's contribution with the common practice: his bounds on the eigenvalues give a sufficient condition such that Hathaway's constraints are satisfied. The simplification is such that the constraints can be easily implemented within the EM algorithm, preserving its monotonicity property (as shown in Ingrassia and Rocci, 2007).

Several authors extended the constrained setup of Ingrassia (2004). Greselin and Ingrassia (2013) applied this setup to mixtures of factor analyzers. They proposed a tuning procedure for selecting the bounds for the eigenvalues of the covariance matrices, based on the final likelihood over a set of runs. Ingrassia and Rocci (2011) modified the constrained algorithm, allowing for stringent constraints which are lifted during the iterations. Browne et al. (2013) combined the ideas in Ingrassia and Rocci (2007, 2011), constraining dynamically the smallest eigenvalue, the largest eigenvalue and both the smallest and the largest ones. All of these proposals share the drawback of not being affine equivariant.

Gallegos and Ritter (2009a; 2009b), and Ritter (2014) applied Hathaway's constraints to robust clustering. They proposed to obtain all local maxima of the trimmed likelihood and, for each solution, investigate the value of $c$ such that it fulfills the constraints. The idea is to choose, \emph{a posteriori}, the solution  with the highest trade-off between scale balance ($c$) and fit (log-likelihood). This approach can be viewed as a refined version of what was proposed in McLachlan and Peel (2000). Garcia-Escudero et al. (2008), from the same strand of literature, introduced the TCLUST algorithm, based on controlling the relative sizes  of the eigenvalues of the cluster scatter matrices. The TCLUST algorithm implies solving several complex optimization problems. Fritz et al. (2013) and Garcia-Escudero et al. (2014) proposed further improvements to the algorithm in order to make it more efficient. The constraints considered therein are not affine equivariant.

Seo and Kim (2012) pointed out that singular and spurious solutions overfit random localized patterns composed of few observations in the dataset. Such observations have a strong influence on the formation of the likelihood-based solutions. Their proposal was to take out such, say, $k$ observations with the highest likelihood (likelihood-based $k$-deleted method), or with the highest value for a score-based statistic (score-based $k$-deleted method). In this way the likelihood of the reduced samples is evaluated at each local maximizer previously found: the root they suggested to select is the one with the highest $k$-deleted likelihood. Kim and Seo (2014) show that their score-based method can be fairly well approximated with the computationally more efficient gradient-based version of the $k$-deleted method.

The degeneracy problem may also be addressed by  adding a penalty to the log-likelihood (penalized approach). Ciuperca et al. (2003) have shown the consistency of the penalized likelihood estimators proposed in Ridolfi and Idier (1999, 2000) for univariate GMM. Chen and Tan (2009) extended the consistency result for the multivariate case. In this framework, the penalty term on the component covariance is added to the log-likelihood (Snoussi and Djafari, 2001; Chen et al, 2008). This penalty can be interpreted as the log of the prior distribution in a Maximum-A-Posteriori estimation setup. Yet, the penalized methods are not affine equivariant, unless the prior's hyperparameters are suitably transformed. MAP estimation, with an \emph{a priori} distribution for the covariance matrices, is what Fraley and Raftery (2007) suggested to use, instead of Maximum-Likelihood, to circumvent the issues of degeneracy and spurious solutions.

\section{Equivariance in the Gaussian Mixture model}\label{invar}

The maximum likelihood estimators (MLE) of Equation (\ref{md}) are equivariant with respect to linear affine transformations of the data. That is, if the data are linearly transformed, the MLE are transformed accordingly. This property is particularly important
in classification because it implies that linear affine transformations of the data
do not change the posterior estimates (Kleinberg, 2002; Ritter, 2014).

The equivariance property can be shown in the following way.
Let us define a linear affine transformation $\mathbf{x}^*=\mathbf{Ax}+\mathbf{b}$, where $\mathbf{A}$ is non singular.
It is well known that
\begin{align*}
\phi(\mathbf{x};\boldsymbol{\mu},\boldsymbol{\Sigma}) & = |\mathbf{A}|\phi(\mathbf{Ax}+\mathbf{b};
\mathbf{A}\boldsymbol{\mu}+\mathbf{b},\mathbf{A}\boldsymbol{\Sigma}\mathbf{A}')\\
& =  |\mathbf{A}|\phi(\mathbf{x}^*;\boldsymbol{\mu}^*,
\boldsymbol{\Sigma}^*) \numberthis  \label{eq:equiGauss}
\end{align*}
where $\boldsymbol{\mu}^*=\mathbf{A}\boldsymbol{\mu}+\mathbf{b}$ and $\boldsymbol{\Sigma}^*=\mathbf{A}\boldsymbol{\Sigma}\mathbf{A}'$.
 This implies that, denoting the likelihood of the original data with $\mathcal{F}$ and the likelihood of the transformed data with $\mathcal{F}^*$, we have, with obvious notation
\begin{align*}
\mathcal{F} & =  L\left(\left\lbrace p_g,\boldsymbol{\mu}_g, \boldsymbol{\Sigma}_g \right\rbrace_G; \left\lbrace \mathbf{x}_i \right\rbrace_n  \right) = \prod_{i=1}^n \sum_{g=1}^G  p_g \phi(\mathbf{x}_i;\boldsymbol{\mu}_g,\boldsymbol{\Sigma}_g) \\
& =   \prod_{i=1}^n \sum_{g=1}^G  p_g |\mathbf{A}| \phi(\mathbf{x}_i^*;\boldsymbol{\mu}_g^*,\boldsymbol{\Sigma}_g^*) =
|\mathbf{A}|^n L\left(\left\lbrace p_g,\boldsymbol{\mu}_g^*, \boldsymbol{\Sigma}_g^* \right\rbrace_G; \left\lbrace \mathbf{x}_i^* \right\rbrace_n \right)\\
& =  |\mathbf{A}|^n \mathcal{F}^*. \numberthis \label{eq:equiGaussmix}
\end{align*}
It follows that there exists a one to one correspondence among the local maxima of $\mathcal{F}$ and $\mathcal{F}^*$.
In particular, if
$\left\{\hat{p}_g,\hat{\boldsymbol{\mu}}_g, \hat{\boldsymbol{\Sigma}}_g \right\}_G$ is a local maximizer for $\mathcal{F},$ then $\left\{\hat{p}_g,\mathbf{A}\hat{\boldsymbol{\mu}}_g+\mathbf{b}, \mathbf{A}\hat{\boldsymbol{\Sigma}}_g \mathbf{A}' \right\}_G$ will be a local maximizer for $\mathcal{F}^*$. Analogously, if $\left\{\hat{p}_g^*,\hat{\boldsymbol{\mu}}_g^*, \hat{\boldsymbol{\Sigma}}_g^* \right\}_G$ is a local maximizer for $\mathcal{F}^*$, then
$\left\{\hat{p}_g^*,\mathbf{A}^{-1}(\hat{\boldsymbol{\mu}}_g^*-\mathbf{b}), \mathbf{A}^{-1}\hat{\boldsymbol{\Sigma}}_g^* (\mathbf{A}')^{-1} \right\}_G$
will be a local maximizer for $\mathcal{F}$.
It is interesting to note that every pair of local maximizers produces the same estimates of the posterior probabilities, that is
\begin{eqnarray*}
\hat{p}(g|\mathbf{x}_i) & = &\frac{\hat{p}_g \phi(\mathbf{x}_i;\hat{\boldsymbol{\mu}}_g,\hat{\boldsymbol{\Sigma}}_g)}{\sum_h \hat{p}_h \phi(\mathbf{x}_i;\hat{\boldsymbol{\mu}}_h,\hat{\boldsymbol{\Sigma}}_h)}\\
& = & \frac{\hat{p}_g |\mathbf{A}| \phi(\mathbf{x}_i^*;\hat{\boldsymbol{\mu}}_g^*,\hat{\boldsymbol{\Sigma}}_g^*)}{\sum_h \hat{p}_h |\mathbf{A}| \phi(\mathbf{x}_i^*;\hat{\boldsymbol{\mu}}_h^*,\hat{\boldsymbol{\Sigma}}_h^*)}=
\hat{p}^*(g|\mathbf{x}^*_i).
\end{eqnarray*}

The above equality proves that the classification obtained via the GMM model is invariant under the group of linear affine transformations on the data $\left\lbrace \mathbf{x}_i\right\rbrace_n$.
This property is crucial when dealing with practical applications as it implies that the clustering does not depend on the choice of a particular method of data standardization - which could instead affect the inference.

\section{Constraints on covariance eigenvalues}\label{constr1}
Hathaway (1985) proposed to impose the following restrictions on the covariance matrices
\begin{equation}
\label{ch}
\lambda_j (\boldsymbol{\Sigma}_g\boldsymbol{\Sigma}_h^{-1}) \geq c,
\hspace{1cm}j=1,\dots,J;\hspace{0.2cm} g,h=1 \dots G
\end{equation}
where $\lambda_j(\mathbf{A})$ is the $j$-th eigenvalue of $\mathbf{A}$ and $0<c\leq 1$.
This prevents the likelihood from diverging and reduces the number of spurious maximizers. However, the method is difficult to implement and a correct choice of $c$ is not simple in practice.  A value of $c$ close to $1$ could exclude the correct solution, whereas a value too close to $0$ is likely to increase the chance of converging to a spurious maximizer.

Ingrassia (2004) simplified Hathaway's constraints as
\begin{equation}
\label{nc}
\sqrt{c} \leq \lambda_j (\boldsymbol{\Sigma}_g)\leq \frac{1}{\sqrt{c}}, \hspace{1cm}j=1,\dots,J;\hspace{0.2cm} g=1 \dots G.
\end{equation}
It is easy to show that (\ref{nc}) implies Hathaway's constraints (\ref{ch}) while the reverse is not necessarily true (Ingrassia, 2004). This ensures a bounded likelihood, and a reduction in the number of spurious maximizers. The constraints are easy to implement, as shown in  Rocci and Ingrassia (2007); however, choosing an optimal \textit{c} is still an issue.

It is important to check if the above constrained approaches offer equivariant estimators under linear affine transformations.
The property can be shown to hold for Hathaway's approach as follows.

Let $\left\lbrace \mathbf{x}_i \right\rbrace_n= \left\{ \mathbf{x}_1, \mathbf{x}_2, \dots \mathbf{x}_n \right\}$ be a sample of i.i.d. observations. The estimates are computed as the solution of the optimization problem
\begin{equation}\label{eq:maxcon1}
\begin{cases}
 \max\limits_{\left\lbrace p_g,\boldsymbol{\mu}_g, \boldsymbol{\Sigma}_g \right\rbrace_G}& L\left(\left\lbrace p_g,\boldsymbol{\mu}_g, \boldsymbol{\Sigma}_g \right\rbrace_G; \left\lbrace \mathbf{x}_i \right\rbrace_n  \right) \\
 \qquad \text{s.t.} & \lambda_j (\boldsymbol{\Sigma}_g\boldsymbol{\Sigma}_h^{-1}) \geq c,\hspace{1cm}j=1,\dots,J;\hspace{0.2cm} g,h=1 \dots G.
\end{cases}
\end{equation}
Given the transformation $\mathbf{x}^*=\mathbf{Ax}+\mathbf{b},$ the maximand in Equation (\ref{eq:maxcon1}) can be rewritten (see Section \ref{invar}) as
\[
L\left(\left\lbrace p_g,\boldsymbol{\mu}_g, \boldsymbol{\Sigma}_g \right\rbrace_G; \left\lbrace \mathbf{x}_i \right\rbrace_n  \right)
=|\mathbf{A}|^n L\left(\left\lbrace p_g,\boldsymbol{\mu}_g^*, \boldsymbol{\Sigma}_g^* \right\rbrace_G; \left\lbrace \mathbf{x}_i^* \right\rbrace_n \right),
\]
where $\boldsymbol{\mu}^*=\mathbf{A}\boldsymbol{\mu}+\mathbf{b}$ and $\boldsymbol{\Sigma}^*=\mathbf{A}\boldsymbol{\Sigma}\mathbf{A}'$.

Noting that
\begin{eqnarray*}
\lambda_j (\boldsymbol{\Sigma}_g \boldsymbol{\Sigma}_h^{-1})
&=&
\lambda_j(\mathbf{A}^{-1}\mathbf{A}\boldsymbol{\Sigma}_g \mathbf{A}'(\mathbf{A}')^{-1} \boldsymbol{\Sigma}_h^{-1})\\
&=&
\lambda_j(\mathbf{A}\boldsymbol{\Sigma}_g \mathbf{A}'(\mathbf{A}')^{-1} \boldsymbol{\Sigma}_h^{-1}\mathbf{A}^{-1})\\
&=&
\lambda_j (\boldsymbol{\Sigma}_g^* (\boldsymbol{\Sigma}_h^*)^{-1}),
\end{eqnarray*}
we can equivalently write the optimization problem in Equation (\ref{eq:maxcon1}) as
\begin{equation}\label{eq:maxcon12}
\begin{cases}
 \max\limits_{\left\lbrace p_g,\boldsymbol{\mu}_g^*, \boldsymbol{\Sigma}_g^* \right\rbrace_G}& L\left(\left\lbrace p_g,\boldsymbol{\mu}_g^*, \boldsymbol{\Sigma}_g^* \right\rbrace_G; \left\lbrace \mathbf{x}_i^* \right\rbrace_n \right) \\
 \qquad \text{s.t.} & \lambda_j (\boldsymbol{\Sigma}_g^* (\boldsymbol{\Sigma}_h^*)^{-1}) \geq c,\hspace{1cm}j=1,\dots,J;\hspace{0.2cm} g,h=1 \dots G.
\end{cases}
\end{equation}
It follows that if $\left\{\hat{p}_g,\hat{\boldsymbol{\mu}}_g, \hat{\boldsymbol{\Sigma}}_g \right\}_G$ is a maximizer for (\ref{eq:maxcon1}), then $\left\{\hat{p}_g,\mathbf{A}\hat{\boldsymbol{\mu}}_g+\mathbf{b}, \mathbf{A}\hat{\boldsymbol{\Sigma}}_g \mathbf{A}' \right\}_G = \left\{\hat{p}_g,\hat{\boldsymbol{\mu}}^*_g, \hat{\boldsymbol{\Sigma}}^*_g \right\}_G$ is a maximizer for (\ref{eq:maxcon12}) and vice-versa, and the two maximization problems are equivalent. As in the unconstrained case, every pair of local maximizers produces the same estimates of the posterior probabilities.
This property does not hold for the constraints given in (\ref{nc}). That is, if $\left\lbrace \hat{\boldsymbol{\Sigma}}_g \right\rbrace_G$ is a constrained local maximizer for $\mathcal{F}$ subject to (\ref{nc}), $\left\lbrace\mathbf{A}\hat{\boldsymbol{\Sigma}}_g \mathbf{A}'\right\rbrace_{G}$ does not necessarily satisfy (\ref{nc}). As an example, let us suppose that $s_{max}(\mathbf{A})< \sqrt{c},$ where $s_{max}(\mathbf{A})$ is the largest singular value of $\mathbf{A}$. In this case, for a given $g,$
\begin{eqnarray*}
\lambda_j(\mathbf{A}\hat{\boldsymbol{\Sigma}}_g \mathbf{A}') &\leq &  \lambda_{max}(\mathbf{A}\hat{\boldsymbol{\Sigma}}_g \mathbf{A}')
\leq s_{max}(\mathbf{A})^2 \lambda_{max}(\hat{\boldsymbol{\Sigma}}_g)\\
& \leq &s_{max}(\mathbf{A})^2 \frac{1}{\sqrt{c}} < c \frac{1}{\sqrt{c}}=\sqrt{c}.
\end{eqnarray*}
We conclude that $\mathbf{A}\hat{\boldsymbol{\Sigma}}_g \mathbf{A}'$ does not satisfy the constraints in (\ref{nc}) because $\lambda_j(\mathbf{A}\hat{\boldsymbol{\Sigma}}_g \mathbf{A}') < \sqrt{c}$, and then it cannot be a constrained local maximizer for $\mathcal{F}^*$.

Constrains in (\ref{nc}) are such that there is no one to one correspondence between the set of local maximizers of $\mathcal{F}$ and $\mathcal{F}^*$.
Thus, the method suffers the disadvantage that the clustering depends on the choice of matrix $\mathbf{A}$. To fix this, data standardization is not the best way to go for two main reasons. First, the standardization requires a choice for the matrix $\mathbf{A}$ and, second, there is no single best approach to data standardization (Milligan and Cooper, 1988; Doherty \textit{et al}, 2007).

It is now clear that affine equivariance is not just a desirable property. It is one of the basic requirements of any clustering method, which should not be sensitive to the changes in the units of measurement of the data. In the next section, our goal will be that of deriving a new set of constraints that are affine equivariant. With an affine equivariant clustering method, researchers and practitioners shall not be concerned anymore with choosing what method to adopt to standardize their data.
 \section{Equivariant constraints} \label{constr2}
 Our proposal is to generalize the constraints (\ref{nc}) by
\begin{equation}
\label{eq:gnc}
\sqrt{c} \leq \lambda_j (\boldsymbol{\Sigma}_g\boldsymbol{\Psi}^{-1})\leq \frac{1}{\sqrt{c}},
\hspace{1cm}j=1,\dots,J;\hspace{0.2cm} g=1 \dots G
\end{equation}
where $\boldsymbol{\Psi}$ is a symmetric positive definite matrix representing our \textit{prior} information about the covariance structure. Clearly, \eqref{eq:gnc} is equal to (\ref{nc}) when $\boldsymbol{\Psi}=\mathbf{I}$.

It can be shown that the above constraints imply Hathaway's constraints.
It is known that (Anderson and Gupta, 1963)
\begin{equation}
\lambda_{min}(\mathbf{AB}^{-1})\geq \lambda_{min}(\mathbf{AC}^{-1})
\lambda_{min}(\mathbf{CB}^{-1}),
\end{equation}
where $\mathbf{A}$ is a positive semi-definite matrix and $\mathbf{B}$ and $\mathbf{C}$ are positive definite matrices.
Now, if \eqref{eq:gnc} holds, then
\begin{equation}
\label{re}
\lambda_{min}(\boldsymbol{\Sigma}_g \boldsymbol{\Sigma}_h^{-1}) \geq
\lambda_{min}(\boldsymbol{\Sigma}_g\boldsymbol{\Psi}^{-1})
\lambda_{min}(\boldsymbol{\Psi} \boldsymbol{\Sigma}_h^{-1})=
\frac{\lambda_{min}(\boldsymbol{\Sigma}_g\boldsymbol{\Psi}^{-1})}
{\lambda_{max}(\boldsymbol{\Sigma}_h \boldsymbol{\Psi}^{-1})}\geq \frac{\sqrt{c}}{\frac{1}{\sqrt{c}}}= c.
\end{equation}
Thus, \eqref{eq:gnc} implies (\ref{ch}).

Furthermore, it can be shown that \eqref{eq:gnc} is invariant under linear and affine transformations provided that $\boldsymbol{\Psi}$ is transformed accordingly, i.e. it is replaced by $\boldsymbol{\Psi}^*=\mathbf{A}\boldsymbol{\Psi} \mathbf{A}'$.
If $\left\lbrace \hat{\boldsymbol{\Sigma}}_g \right\rbrace_G$ is a constrained local maximizer for $\mathcal{F}$ subject to \eqref{eq:gnc}, then $\left\lbrace \hat{\boldsymbol{\Sigma}}^*_g = \mathbf{A}\hat{\boldsymbol{\Sigma}}_g \mathbf{A}' \right\rbrace_G$ is a local maximizer for $\mathcal{F}^*$ subject to \eqref{eq:gnc} for $g=1 \dots G$. We have that
\begin{eqnarray*}
\label{eiglt}
\lambda_j (\hat{\boldsymbol{\Sigma}}_g \boldsymbol{\Psi}^{-1}) &=&
\lambda_j  (\hat{\boldsymbol{\Sigma}}_g\mathbf{A}'(\mathbf{A}')^{-1} \boldsymbol{\Psi}^{-1}\mathbf{A}^{-1}\mathbf{A})\\
&=&\lambda_j(\mathbf{A}\hat{\boldsymbol{\Sigma}}_g \mathbf{A}'(\mathbf{A}')^{-1} \boldsymbol{\Psi}^{-1}\mathbf{A}^{-1}) \\
&=& \lambda_j (\hat{\boldsymbol{\Sigma}}^*_g \boldsymbol{(\Psi^*)}^{-1}).
\end{eqnarray*}

In words, if a linear affine transformation is performed on the data, $\boldsymbol{\Psi}$ must be changed accordingly. This scheme of transforming $\boldsymbol{\Psi}$ ensures the equivariance of the method. By contrast, holding $\boldsymbol{\Psi}$ fixed breaks the equivariance property.

The constraints \eqref{eq:gnc} have the effect of shrinking the covariance matrices to $\boldsymbol{\Psi},$ and the level of shrinkage is given by the value of $c.$ Note that for
$c=1$, $\hat{\boldsymbol{\Sigma}}_g = \boldsymbol{\Psi},$ whereas for $c\rightarrow 0,$ $\hat{\boldsymbol{\Sigma}}_g$ equals the unconstrained ML estimate. Furthermore we can show that the Stein's discrepancy - known as Stein's loss (James and Stein, 1961) - between the matrices $\hat{\boldsymbol{\Sigma}}_g$ and $\boldsymbol{\Psi}$ goes to zero as $c$ approaches one. The Stein's discrepancy between the matrices $\hat{\boldsymbol{\Sigma}}_g$ and $\boldsymbol{\Psi}$ is
\begin{equation}\label{eq:Lstein1}
\text{L}(\hat{\boldsymbol{\Sigma}}_g,\boldsymbol{\Psi}) = \text{tr}(\hat{\boldsymbol{\Sigma}}_g\boldsymbol{\Psi}^{-1}) - \log |\hat{\boldsymbol{\Sigma}}_g\boldsymbol{\Psi}^{-1}| - J \geq 0.
\end{equation}
Let us rewrite Equation \eqref{eq:Lstein1} as follows.
\begin{equation}\label{eq:Lstein2}
\text{L}(\hat{\boldsymbol{\Sigma}}_g,\boldsymbol{\Psi}) = \sum_{j=1}^{J} \lambda_{j} (\hat{\boldsymbol{\Sigma}}_g\boldsymbol{\Psi}^{-1}) - \sum_{j=1}^{J} \log(\lambda_{j} (\hat{\boldsymbol{\Sigma}}_g\boldsymbol{\Psi}^{-1})) - J
\end{equation}
Using the constraints in \eqref{eq:gnc}, we can derive the following majorizing function
\begin{equation}\label{eq:minorz}
\text{L}(\hat{\boldsymbol{\Sigma}}_g,\boldsymbol{\Psi}) \leq \frac{J}{\sqrt{c}} - J\log(\sqrt{c}) - J,
\end{equation}
which is decreasing in $c.$ This can be shown by noting that the first derivative of the right-hand side of \eqref{eq:minorz} with respect to $c$ is equal to $-\frac{J}{2c\sqrt{c}} - \frac{J}{2c},$ and is negative when $0<c\leq 1$. This implies that the function is decreasing when $c$ increases within the interval $(0,1]$.

Intuitively, the constraints \eqref{eq:gnc} provide with a way to obtain a model in between a too restrictive model, the homoscedastic, and an ill-conditioned model, the heteroscedastic.

\section{Data-driven choice of $\boldsymbol{\Psi}$ and $c$} \label{crossvalid}

Issues arise when \emph{a priori} information about the structure of the class conditional covariance matrices is not available. In that case, $\boldsymbol{\Psi}$ and $c$ have to be selected from the data. From the previous discussion, for a given $c,$ every $\hat{\boldsymbol{\Sigma}}_g$ cannot be too far from $\boldsymbol{\Psi}$ in terms of Stein's discrepancy. Thus $\boldsymbol{\Psi}$ can be seen as the barycenter of the cloud of the $\hat{\boldsymbol{\Sigma}}_g$'s: the \emph{average} conditional covariance matrix. Therefore, the most natural choice is to estimate such \emph{average} as the within covariance matrix of the homoscedastic Gaussian model. How close the final clustering will be to the homoscedastic model will depend on the value of the tuning constant: for values of $c$ close to $0$, the resulting clustering will be close to that of the heteroscedastic mixture model, whereas $c\rightarrow 1$ implies a clustering close to that of the homoscedastic mixture model.

Other possible choices of $\boldsymbol{\Psi},$ which guarantee the equivariance of the constraints, are available: the sample covariance matrix, which is computationally faster and is frequently used as hyperparameter in Bayesian Gaussian mixtures (for instance, see Fraley and Raftery, 2007), or the within covariance matrix of a homoscedastic mixture of Student-\textit{t}. To motivate this, let us recall that a random vector conditionally distributed as a multivariate Gaussian, given Wishart inverse covariance matrix, has a multivariate Student-$t$ distribution (Dawid, 1981; Dickey, 1967). Using similar arguments as in Peel and McLachlan (2000), if $\mathbf{x}| \boldsymbol{\Sigma}_1, \dots, \boldsymbol{\Sigma}_G$ is a GMM, and $\boldsymbol{\Sigma}_1^{-1}, \dots, \boldsymbol{\Sigma}_G^{-1}$ are i.i.d. Wishart random variables, the marginal distribution of $\mathbf{x}$ is a homoscedastic mixture of Student-$t$'s.

The choice of $c$ is crucial. A value of $c$ too large could exclude the right solution, whereas a too small value of $c$ is likely to increase the chance to converge to spurious local maxima: such solutions overfit random localized pattern composed of few data points being almost co-planar (Ritter, 2014; Seo and Kim, 2012). Hence, selecting $c$ jointly with the mixture parameters by maximizing the likelihood on the entire sample would trivially yield a scale balance approaching zero.

A practical alternative would be to split the data into a training set, where model parameters are estimated, and a test set, where the log-likelihood is evaluated for a given value of $c.$ The optimal tuning parameter $c$ would then be selected such that the test set log-likelihood is maximized.

The use of the test set log-likelihood as a model selection tool is advocated by Smyth (1996; 2000), in the context of estimating the number of mixture components. The motivation behind its use is that it can be showed to be an unbiased estimator (within a constant) of the Kullback-Leibler divergence between the \emph{truth} and the model under consideration (Smyth, 2000). This means that, even under a misspecified model, the procedure renders a $c$ such that the Kullback-Leibler divergence is minimized.

In spite of the usual unavailability of large independent test sets, a valid alternative is to use the \emph{cross-validated} log-likelihood in order to estimate the test set log-likelihood. This consists in repeatedly partitioning the data into training and test sets and, for a given $c,$ estimate the mixture parameters on the training sets. The model fit is then measured summing the log-likelihoods of the test sets evaluated at the parameters computed on the training sets, obtaining the so-called \emph{cross-validated} log-likelihood. The constant $c$ is chosen such that the \emph{cross-validated} log-likelihood is maximized. This can be viewed as a function of $c$ only (Smyth, 1996), and would solve the issue of overfitting as training and test sets are independent (Arlot and Celisse, 2010).

In details, let us partition $K$ times the full data set $\left\lbrace \mathbf{x}_i; i \in N \right\rbrace_n$ into two parts, a training set $\mathbf{x}_{S} = \left\lbrace \mathbf{x}_i ;i \in S \right\rbrace_{n_S},$ and a test set $\mathbf{x}_{\bar{S}} = \left\lbrace \mathbf{x}_i ; i \in \bar{S} \right\rbrace_{n_{\bar{S}}}$ with $S \cup \bar{S} = N$ and $n_S+n_{\bar{S}}=n.$ For the $k$-th partition, let $\hat{\theta}(c,S_{k})$ be the constrained maximum likelihood estimator based on the training set $\mathbf{x}_{S_{k}}.$ Furthermore, let $l_{\bar{S}_{k}} [\hat{\theta}(c,S_k)]$ be the log-likelihood function evaluated at the test set $\mathbf{x}_{\bar{S}_k}.$ The \emph{cross-validated} log-likelihood is defined as the sum of the contributions of each test set to the log-likelihood
\begin{equation}
\text{CV}(c)=\sum_{k=1}^K l_{\bar{S}_k} [\hat{\theta}(c,S_k)].
\end{equation}
The best $c$ is chosen as the maximizer of $\text{CV}(c)$.

Further details on the choice of the number of random partitions $K$ and of the sizes of training and test sets are given in Section \ref{simulation}.

\section{Algorithm} \label{alg}
The objective is to maximize \eqref{eq:like} under the constraints \eqref{eq:gnc}.
Thanks to the equivariance property of the constraints, we can act any linear affine transformation to the data. This is useful since it will suffice to transform the data so to have $\boldsymbol{\Psi}=\mathbf{I}_J$ and the existing algorithm of Ingrassia and Rocci (2007) can be applied on the transformed data.

The transformation is $\mathbf{x}^*=\mathbf{L}^{-\frac{1}{2}}\mathbf{Q}^{'}\mathbf{x},$ where $\boldsymbol{\Psi}=\mathbf{QLQ}'$ is the singular value decomposition of $\boldsymbol{\Psi}.$ This leads to $\boldsymbol{\Psi}^*=\mathbf{L}^{-\frac{1}{2}}\mathbf{Q}' \boldsymbol{\Psi}\mathbf{QL}^{-\frac{1}{2}}=\mathbf{I}_J$.

For sake of completeness, we recall briefly the updates of the algorithm proposed by Ingrassia and Rocci (2007).\\
\\
\textbf{Update} $u_{ig}$, $p_g,$ $\boldsymbol{\mu}_g$ \\
As in the case of a normal mixture, the updates are
\begin{equation}
u_{ig}=\frac{p_g  \phi(\mathbf{x}_i;\boldsymbol{\mu}_g,\boldsymbol{\Sigma}_g)}{\sum_{h=1}^G  \phi(\mathbf{x}_i;\boldsymbol{\mu}_h,\boldsymbol{\Sigma}_h)};
\end{equation}
\begin{equation}
p_g=\frac{1}{n}\sum_{i=1}^n u_{ig};
\end{equation}
\begin{equation}
\mu_g=\frac{\sum_i u_{ig}\mathbf{x}_i}{\sum_i u_{ig}}.
\end{equation}
\textbf{Update} $\boldsymbol{\Sigma}_g$\\
Compute
\begin{equation}
\mathbf{S}_g=\frac{1}{\sum_i u_{ig}}\sum_i u_{ig}(\mathbf{x}_i-\mu_g)(\mathbf{x}_i-\mu_g)',
\end{equation}
and set
\begin{equation}
\lambda_{qg}=\min\left( \frac{1}{\sqrt{c}},\max\left(\sqrt{c},l_{qg} \right) \right),
\end{equation}
where $\mathbf{L}_g=\textnormal{diag}\left(l_{1g},\dots l_{Jg} \right)$ is the diagonal matrix of the eigenvalues in non decreasing order of $\mathbf{S}_g,$ and $\mathbf{S}_g=\mathbf{Q}_g \mathbf{L}_g \mathbf{Q}_g'$ its singular value decomposition.
Letting $\boldsymbol{\Lambda}_g=\textnormal{diag}\left(\lambda_{1g},\dots \lambda_{jg} \right)$, the update of $\boldsymbol{\Sigma}_g$ is given by
\begin{equation}
\boldsymbol{\Sigma}_g=\mathbf{Q}_g\boldsymbol{\Lambda}_g \mathbf{Q}_g'.
\end{equation}

\section{Simulation study}\label{simulation}
\subsection{Design}
In this section we perform a simulation experiment in order to compare the performance of the proposed methods with respect to some existing approaches in the literature. In particular we consider the following seven algorithms:
\begin{enumerate}
\item Unconstrained
\begin{enumerate}
\item homoscedastic normal (homN), within covariance matrix $\boldsymbol{\Sigma};$
\item heteroscedastic normal (hetN), $0.0000001 \leq \lambda_j(\boldsymbol{\Sigma}_g) \leq 10000000$ to prevent degeneracy and numerical instability;
\item homoscedastic Student \textit{t} (hom\textit{t}), scale matrix $\boldsymbol{\Xi},$ $\beta=4$ (McLachlan and Peel, 1998).
\end{enumerate}
\item Constrained
\begin{enumerate}
\item sample covariance (con$\mathbf{S}$), $\boldsymbol{\Psi}=\mathbf{S};$
\item normal (conN), $\boldsymbol{\Psi}=\boldsymbol{\Sigma};$
\item Student \textit{t} (con\textit{t}), $\boldsymbol{\Psi}=\frac{\beta \boldsymbol{\Xi}}{(\beta -2)}.$
\end{enumerate}
\end{enumerate}
For each sample, we randomly split the data $K = 25$ times into a training set $\mathbf{x}_{S},$ and a test set $\mathbf{x}_{\bar{S}}.$ Choosing how many times to partition the full data set is a trade-off between variability of the estimates and computational burden. As Smyth (2000) argues in the context of model selection for probabilistic clustering using cross-validation, the larger the value of $K,$ the less the variability in the log-likelihood estimates. In practice - the Author argues - values of $K$ between 20 and 50 appear adequate for most applications.

The choice of the size of the test set must be such that the training set has all components represented. If one component is not represented in the test set, but the parameters are correctly estimated using the training set, the test set log-likelihood will correctly display the fit of the model. By contrast, if one component is not represented in the training set, although estimation of the other components’ parameters can be correct, the fit displayed by the test set log-likelihood will be poor. Van der Laan, Dudoit, and Keles (2004) found, in their simulation study, that the likelihood-based cross-validation procedure is performing equally well with any choice of the relative size of the test set between 0.1 and 0.5. As argued in Kearns (1996), the importance of choosing an optimal size for the training set increases as the target function becomes more complex relative to the sample size. Bearing this in mind, we choose to consider a training set of size $n_S=n - \frac{n}{10},$ and a test set $\mathbf{x}_{\bar{S}}$ of size $n_{\bar{S}}=\frac{n}{10}.$

Then the cross-validation scheme, as described in Section \ref{crossvalid}, is applied and the optimal $c$ is chosen by using a line search with six function evaluations.
The sample data have been generated from $G$-class mixtures of heteroscedastic $J$-variate normal distributions with:
\begin{itemize}
\item $n=50,$ $100,$ $200;$
\item $J=5,$ $8;$
\item prior membership probabilities $\mathbf{p}=(0.2, 0.3, 0.5)',$ $(0.1, 0.4, 0.5)',$ $(0.1, 0.1, 0.2, 0.3, 0.3)'.$
\end{itemize}
This yields a total of $2 \times 3 \times 3$ simulation conditions.

For each simulation condition, we generate 250 data sets, each with different means and covariance matrices, where:
\begin{itemize}
\item component means $\mu_{jg} \sim N(0,1.5^2)$, independent;
\item eigenvalues of the covariance matrices $\lambda_{jg} \sim U(0,\frac{g}{\text{sep}})$, independent with $\text{sep}=2$;
\item eigenvectors of the covariance matrices generated by orthonormalizing matrices generated independently from a standard normal distribution.
\end{itemize}

It is well known that the EM for GMM is sensitive to the initial position, especially in the multivariate context (among others, McLachlan and Peel, 2000). We choose to adopt the standard \emph{multiple random starts} strategy. That is, for each data set, 10 random initial partitions are generated: these are used as starting values for the M-step of all the seven algorithms under analysis. For conN, conS, and cont, a constrained algorithm with arbitrary lower and upper bounds of respectively  $0.5$ and $2$ is run in order to exclude degenerate (and some spurious) solutions, and the estimated clustering is used to initialize the cross-validation scheme. The alternative option of directly generating 10 different starts for each training set - within the cross-validation scheme - would have added little in terms of accuracy of the final estimates.

Concerning the root selection criterion, for the unconstrained algorithm, we select the roots yielding the highest likelihood, whereas for the constrained algorithms we select the roots based on the cross-validated likelihood.

The performance of the different techniques has been analyzed in terms of:
\begin{itemize}
\item MAD (Mean absolute deviation): $\sum_{g=1}^G \sum_{i=1}^n |p(g|x_i)-\hat{p}(g|x_i)|;$
\item ARand (Adjusted Rand index; Hubert and Arabie, 1985);
\item computational time needed to analyze a single data set;
\item the value of the calibrated constant $c$ (for the constrained approach only).
\end{itemize}
The MAD is computed evaluating the above expression for all possible permutations of the estimated classes. The final MAD reported refers to the permutation which yields the lowest difference, and measures inaccuracy of estimated \emph{fuzzy} classification - whereas ARand measures accuracy of estimated \emph{crisp} classification.

In addition, we tested the robustness of the results with respect to changes in 1) the cross-validation settings, and 2) the level of class separation. In order to test robustness with respect to cross-validation settings, we considered a subset of the above simulation conditions as follows. 250 samples, of 50, 100, and 200 observations, were generated from a 3-group 8-variate heteroscedastic Gaussian mixture model, with prior class membership probabilities of 0.1, 0.4, and 0.5.

The same setting was used in order to count how many different local maxima each algorithm converged to over the 10 random initializations considered. This serves the purpose of providing some information on the likelihood surface.

Class separation has been manipulated by controlling the dispersion of the group conditional covariance matrices' eigenvalues (through the above sep value): higher dispersion levels correspond to overlap between the classes. Considering the above full simulation as corresponding to fixed moderate separation ($\text{sep}=2$), this final setup compares results for low, moderate, and high separation levels - respectively $\text{sep}=1,$ $\text{sep}=2,$ and $\text{sep}=3$. The subset of simulation conditions considered is as follows. 250 samples, of 50 observations each, were generated from a 3-group and 5-group heteroscedastic Gaussian mixture model, with prior class membership probabilities of respectively 0.2, 0.3, and 0.5; 0.1, 0.4, and 0.5; 0.1, 0.1, 0.2, 0.3, and 0.3. Table \ref{tabsimcon} summarizes the conditions explored in all testing setups.

\FloatBarrier
\begin{table}[h!]
	\centering
		\begin{tabular}{lcccccccc}
		\hline \hline
	&	Full simulation	&	Cross-val settings	&	Class-sep	&	N. local max	\\
	&		&		&		&		\\
$J=5$, $p=(0.2,0.3,0.5)'$, $n=50$	&	\checkmark	&	$\times$	&	\checkmark	&	$\times$	\\
$J=5$, $p=(0.2,0.3,0.5)'$, $n=100$	&	\checkmark	&	$\times$	&	$\times$	&	$\times$	\\
$J=5$, $p=(0.2,0.3,0.5)'$, $n=200$	&	\checkmark	&	$\times$	&	$\times$	&	$\times$	\\
$J=5$, $p=(0.1,0.4,0.5)'$, $n=50$	&	\checkmark	&	$\times$	&	\checkmark	&	$\times$	\\
$J=5$, $p=(0.1,0.4,0.5)'$, $n=100$	&	\checkmark	&	$\times$	&	$\times$	&	$\times$	\\
$J=5$, $p=(0.1,0.4,0.5)'$, $n=200$	&	\checkmark	&	$\times$	&	$\times$	&	$\times$	\\
$J=5$, $p=(0.1,0.1,0.2,0.3,0.3)'$, $n=50$	&	\checkmark	&	$\times$	&	\checkmark	&	$\times$	\\
$J=5$, $p=(0.1,0.1,0.2,0.3,0.3)'$, $n=100$	&	\checkmark	&	$\times$	&	$\times$	&	$\times$	\\
$J=5$, $p=(0.1,0.1,0.2,0.3,0.3)'$, $n=200$	&	\checkmark	&	$\times$	&	$\times$	&	$\times$	\\
$J=8$, $p=(0.2,0.3,0.5)'$, $n=50$	&	\checkmark	&	$\times$	&	\checkmark	&	$\times$	\\
$J=8$, $p=(0.2,0.3,0.5)'$, $n=100$	&	\checkmark	&	$\times$	&	$\times$	&	$\times$	\\
$J=8$, $p=(0.2,0.3,0.5)'$, $n=200$	&	\checkmark	&	$\times$	&	$\times$	&	$\times$	\\
$J=8$, $p=(0.1,0.4,0.5)'$, $n=50$	&	\checkmark	&	\checkmark	&	\checkmark	&	\checkmark	\\
$J=8$, $p=(0.1,0.4,0.5)'$, $n=100$	&	\checkmark	&	\checkmark	&	$\times$	&	\checkmark	\\
$J=8$, $p=(0.1,0.4,0.5)'$, $n=200$	&	\checkmark	&	\checkmark	&	$\times$	&	\checkmark	\\
$J=8$, $p=(0.1,0.1,0.2,0.3,0.3)'$, $n=50$	&	\checkmark	&	$\times$	&	\checkmark	&	$\times$	\\
$J=8$, $p=(0.1,0.1,0.2,0.3,0.3)'$, $n=100$	&	\checkmark	&	$\times$	&	$\times$	&	$\times$	\\
$J=8$, $p=(0.1,0.1,0.2,0.3,0.3)'$, $n=200$	&	\checkmark	&	$\times$	&	$\times$	&	$\times$	\\

		\hline \hline
		\end{tabular}
	\caption{Cross-table of simulation condition and simulation type.}
	\label{tabsimcon}
\end{table}

 \FloatBarrier

\subsection{Simulation results}
Tables \ref{tabfullj5} and \ref{tabfullj8} present the results obtained with, respectively, $J=5$ and $J=8.$

Among the unconstrained approaches, as expected, for small samples ($n = 50$), the heteroscedastic
normal (hetN) performs poorly, while the homoscedastic Student-t (homt) works nicely. However, the constrained approach (cont) is able to cope with such a small sample size and to improve the performance of the unconstrained approach. The heteroscedastic performs poorly, especially with higher model complexity in terms of number of components and variables relative to the sample size. A similar pattern is, indeed, observed for $n=100$ when $J=5$ with $G=5$ and when $J=8$.
As the sample size gets larger ($n = 200$), the homoscedastic models are the worst, while the unconstrained heteroscedastic increases in quality of classification. Interestingly, even for a large sample size, conS, conN and cont yield higher or equal quality estimation compared to the unconstrained approach. However, on average, cont seems to be the best, especially for small sample sizes. When $G=3$, the gains observed by the constrained approaches, in terms of cluster recovery, are more pronounced in presence of a component with a small weight ($p=(0.1,0.4,0.5)'$). In general we observe that, whereas increasing the sample size improves the performance of all methods, an increasing number of components lowers the quality of the clustering results.

The results point out that hetN and conS suffer higher values of $J.$ This is not surprising, as a growing number of variables causes, all else equal, parameter proliferation and a consequent loss in the accuracy of the estimation of the class conditional covariance matrices. Although parameter proliferation is limited with conS, this is constructed based on the choice of the sample covariance matrix as target, which is also very sensible to an increasing $J$ (holding $n$ fixed). In addition, the sample covariance matrix is the sum of the within and the between variance: as such it does not seem to be the best choice for $\boldsymbol{\Psi}.$ On the other hand, homN, homt, conN and cont process the increase in $J$ from 5 to 8 as additional information that, all else equal, improves the quality of the estimation (see Table \ref{tabfullj8}, and Figures \ref{fig235}, \ref{fig145}, and \ref{fig11233}. This is typically the case in finite mixtures with discrete variables (among others, Vermunt, 2010; Di Mari, Oberski, and Vermunt, 2016).

Overall the results show that the constant $c$ decreases in all constrained approaches as the sample size increases, coherently with the results of Xu et al. (2010) in the univariate case.

In terms of computational time, even if the sample covariance matrix is faster to compute, simulation results show that conS converges slower than conN and cont.
\FloatBarrier
\begin{table}[h!]
	\centering
		\begin{tabular}{cclcccccc}
		\hline \hline
	&		&		&	homN	&	hetN	&	homt	&	conS	&	conN	&	cont	\\ \hline
	&		&		&		 	&			&		   &		&		&				\\
p=(0.2,0.3,0.5)'	&	n=50	&	MAD	&	0.11	&	0.25	&	0.11	&	0.16	&	0.11	&	0.09	\\
	&		&	ARand	&	0.78	&	0.52	&	0.79	&	0.68	&	0.79	&	0.82	\\
	&		&	time	&	0.10	&	0.08	&	0.09	&	1.17	&	0.67	&	0.73	\\
	&		&	c	&		&		&		&	0.32	&	0.93	&	0.79			\\
	&		&		&		&		&		&		&		&		\\
	&	n=100	&	MAD	&	0.08	&	0.07	&	0.06	&	0.06	&	0.06	&	0.05	\\
	&		&	ARand	&	0.84	&	0.86	&	0.87	&	0.87	&	0.88	&	0.91	\\
	&		&	time	&	0.19	&	0.14	&	0.14	&	1.22	&	0.81	&	0.85	\\
	&		&	c	&		&		&		&	0.12	&	0.77	&	0.53			\\
	&		&		&		&		&		&		&		&		\\
	&	n=200	&	MAD	&	0.06	&	0.02	&	0.05	&	0.02	&	0.02	&	0.02	\\
	&		&	ARand	&	0.88	&	0.96	&	0.90	&	0.96	&	0.96	&	0.96	\\
	&		&	time	&	0.42	&	0.21	&	0.21	&	1.59	&	1.11	&	1.16	\\
	&		&	c	&		&		&		&	0.07	&	0.39	&	0.28	\\ \hline
	&		&		&		&		&		&		&		&		\\
p=(0.1,0.4,0.5)'	&	n=50	&	MAD	&	0.12	&	0.23	&	0.17	&	0.19	&	0.12	&	0.12	\\
	&		&	ARand	&	0.77	&	0.54	&	0.72	&	0.65	&	0.78	&	0.79	\\
	&		&	time	&	0.11	&	0.08	&	0.10	&	1.2	&	0.70	&	0.77	\\
	&		&	c	&		&		&		&	0.34	&	0.93	&	0.78			\\
	&		&		&		&		&		&		&		&		\\
	&	n=100	&	MAD	&	0.09	&	0.09	&	0.15	&	0.09	&	0.07	&	0.06	\\
	&		&	ARand	&	0.82	&	0.83	&	0.76	&	0.83	&	0.86	&	0.89	\\
	&		&	time	&	0.19	&	0.15	&	0.17	&	1.28	&	0.84	&	0.91	\\
	&		&	c	&		&		&		&	0.14	&	0.79	&	0.52			\\
	&		&		&		&		&		&		&		&		\\
	&	n=200	&	MAD	&	0.07	&	0.03	&	0.11	&	0.03	&	0.03	&	0.03	\\
	&		&	ARand	&	0.87	&	0.95	&	0.82	&	0.93	&	0.94	&	0.94	\\
	&		&	time	&	0.47	&	0.25	&	0.35	&	1.78	&	1.28	&	1.38	\\
	&		&	c	&		&		&		&	0.07	&	0.45	&	0.31	\\ \hline
	&		&		&		&		&		&		&		&		\\
p=(0.1,0.1,0.2,0.3,0.3)'	&	n=50	&	MAD	&	0.27	&	0.46	&	0.28	&	0.31	&	0.27	&	0.26	\\
	&		&	ARand	&	0.58	&	0.29	&	0.57	&	0.52	&	0.58	&	0.58	\\
	&		&	time	&	0.12	&	0.09	&	0.15	&	1.49	&	0.99	&	1.05	\\
	&		&	c	&		&		&		&	0.42	&	0.95	&	0.92	\\
	&		&		&		&		&		&		&		&		\\
	&	n=100	&	MAD	&	0.23	&	0.35	&	0.21	&	0.26	&	0.21	&	0.19	\\
	&		&	ARand	&	0.65	&	0.48	&	0.68	&	0.60	&	0.67	&	0.70	\\
	&		&	time	&	0.26	&	0.24	&	0.28	&	2.03	&	1.49	&	1.64	\\
	&		&	c	&		&		&		&	0.27	&	0.90	&	0.66	\\
	&		&		&		&		&		&		&		&		\\
	&	n=200	&	MAD	&	0.21	&	0.17	&	0.18	&	0.16	&	0.13	&	0.14	\\
	&		&	ARand	&	0.68	&	0.73	&	0.72	&	0.74	&	0.79	&	0.78	\\
	&		&	time	&	0.71	&	0.62	&	0.65	&	3.87	&	2.75	&	2.8	\\
	&		&	c	&		&		&		&	0.13	&	0.53	&	0.40	\\ \hline \hline
\end{tabular}
	\caption{Average results over 250 simulated data sets of five-variate ($J=5$) GMM estimation, with three and five components. Initialization from ten random starts.}
	\label{tabfullj5}
\end{table}
\FloatBarrier

\begin{table}[h!]
	\centering
		\begin{tabular}{cclccccccc}
		\hline \hline
	&		&		&	homN	&	hetN	&	homt	&	conS	&	conN	&	cont	\\ \hline
	&		&		&		&		&		&		&		&		\\
p=(0.2,0.3,0.5)'	&	n=50	&	MAD	&	0.13	&	0.39	&	0.05	&	0.25	&	0.13	&	0.04	\\
	&		&	ARand	&	0.74	&	0.27	&	0.91	&	0.50	&	0.75	&	0.91	\\
	&		&	time	&	0.08	&	0.06	&	0.09	&	1.51	&	0.66	&	0.70	\\
	&		&	c	&		&		&		&	0.58	&	0.94	&	0.84	\\
	&		&		&		&		&		&		&		&		\\
	&	n=100	&	MAD	&	0.03	&	0.12	&	0.02	&	0.07	&	0.02	&	0.01	\\
	&		&	ARand	&	0.95	&	0.75	&	0.96	&	0.84	&	0.96	&	0.97	\\
	&		&	time	&	0.17	&	0.14	&	0.13	&	1.51	&	0.81	&	0.85	\\
	&		&	c	&		&		&		&	0.20	&	0.94	&	0.60	\\
	&		&		&		&		&		&		&		&		\\
	&	n=200	&	MAD	&	0.01	&	0.00	&	0.01	&	0.00	&	0.00	&	0.00	\\
	&		&	ARand	&	0.97	&	0.99	&	0.98	&	0.99	&	0.99	&	0.99	\\
	&		&	time	&	0.41	&	0.21	&	0.24	&	1.71	&	1.14	&	1.18	\\
	&		&	c	&		&		&		&	0.07	&	0.50	&	0.34	\\ \hline
	&		&		&		&		&		&		&		&		\\
p=(0.1,0.4,0.5)'	&	n=50	&	MAD	&	0.15	&	0.37	&	0.14	&	0.27	&	0.14	&	0.10	\\
	&		&	ARand	&	0.72	&	0.28	&	0.77	&	0.49	&	0.73	&	0.83	\\
	&		&	time	&	0.09	&	0.06	&	0.11	&	1.66	&	0.72	&	0.76	\\
	&		&	c	&		&		&		&	0.61	&	0.95	&	0.81	\\
	&		&		&		&		&		&		&		&		\\
	&	n=100	&	MAD	&	0.05	&	0.12	&	0.09	&	0.10	&	0.04	&	0.03	\\
	&		&	ARand	&	0.90	&	0.77	&	0.86	&	0.81	&	0.93	&	0.94	\\
	&		&	time	&	0.18	&	0.15	&	0.17	&	1.62	&	0.88	&	0.95	\\
	&		&	c	&		&		&		&	0.24	&	0.89	&	0.58	\\
	&		&		&		&		&		&		&		&		\\
	&	n=200	&	MAD	&	0.02	&	0.02	&	0.05	&	0.02	&	0.01	&	0.01	\\
	&		&	ARand	&	0.96	&	0.97	&	0.92	&	0.97	&	0.98	&	0.98	\\
	&		&	time	&	0.44	&	0.26	&	0.35	&	1.77	&	1.27	&	1.34	\\
	&		&	c	&		&		&		&	0.08	&	0.56	&	0.36	\\ \hline
	&		&		&		&		&		&		&		&		\\
p=(0.1,0.1,0.2,0.3,0.3)'	&	n=50	&	MAD	&	0.22	&	0.52	&	0.19	&	0.33	&	0.21	&	0.18	\\
	&		&	ARand	&	0.65	&	0.19	&	0.70	&	0.49	&	0.66	&	0.72	\\
	&		&	time	&	0.09	&	0.05	&	0.16	&	1.98	&	1.05	&	1.10	\\
	&		&	c	&		&		&		&	0.60	&	0.95	&	0.93	\\
	&		&		&		&		&		&		&		&		\\
	&	n=100	&	MAD	&	0.13	&	0.40	&	0.10	&	0.25	&	0.13	&	0.10	\\
	&		&	ARand	&	0.80	&	0.38	&	0.84	&	0.61	&	0.80	&	0.84	\\
	&		&	time	&	0.22	&	0.18	&	0.26	&	2.11	&	1.49	&	1.63	\\
	&		&	c	&		&		&		&	0.37	&	0.93	&	0.69	\\
	&		&		&		&		&		&		&		&		\\
	&	n=200	&	MAD	&	0.09	&	0.18	&	0.07	&	0.10	&	0.05	&	0.06	\\
	&		&	ARand	&	0.86	&	0.74	&	0.90	&	0.84	&	0.91	&	0.91	\\
	&		&	time	&	0.57	&	0.58	&	0.52	&	3.46	&	2.44	&	2.40	\\
	&		&	c	&		&		&		&	0.18	&	0.68	&	0.46	\\ \hline \hline
\end{tabular}
	\caption{Average results over 250 simulated data sets of eight-variate ($J=8$) GMM estimation, with three and five components. Initialization from ten random starts.}
	\label{tabfullj8}
\end{table}
\FloatBarrier
\begin{figure}
\centering
\includegraphics[clip, trim=0cm 8.75cm 0cm 2.5cm,width=\linewidth]{020305.pdf}
\caption{Boxplot of the ARand values observed in $250$ simulated data sets, $n=50,$ $n=100,$ and $n=200,$ $p=(0.2,0.3,0.5)',$ and $J=8.$}
\label{fig235}
\end{figure}
\begin{figure}
\centering
\includegraphics[clip, trim=0cm 8.75cm 0cm 2.5cm,width=\linewidth]{010405.pdf}
\caption{Boxplot of the ARand values observed in $250$ simulated data sets, $n=50,$ $n=100,$ and $n=200,$ $p=(0.1,0.4,0.5)',$ and $J=8.$}
\label{fig145}
\end{figure}
\begin{figure}
\centering
\includegraphics[clip, trim=0cm 8.75cm 0cm 2.5cm,width=\linewidth]{0101020303.pdf}
\caption{Boxplot of the ARand values observed in $250$ simulated data sets, $n=50,$ $n=100,$ and $n=200,$ $p=(0.1,0.1,0.2,0.3,0.3)',$ and $J=8.$}
\label{fig11233}
\end{figure}
 \FloatBarrier

Table \ref{tablocmax} displays the means and the medians of the number of local maxima the methods under comparison found in the reduced simulation setup (see Table \ref{tabsimcon}). We observe that conN and cont, on average, yield the lowest number of local maxima in all the three sample size conditions. Increasing the sample size yields a more stable behavior also for the other methods.

\FloatBarrier
\begin{table}[h!]
	\centering
		\begin{tabular}{clllllll}
		\hline \hline
	&		&	homN	&	hetN	&	homt	&	conS	&	conN	&	cont		\\ \hline
n=50	& 	mean	&	9.64  	&	9.99	&	7.01	&	9.86	&	5.31	&	4.99	\\
	& 	median	&	10  	&	10	&	7	&	10	&	5	&	5	\\
	&		&		&		&		&		&		&		 \\
n=100	& 	mean	&	7.08  	&	9.94	&	4.71	&	8.50	&	3.96	&	4.64	\\
	& 	median	&	7  	&	10	&	5	&	9	&	4	&	4	\\
	&		&		&		&		&		&		&		\\
n=200	& 	mean	&	3.82  	&	7.04	&	3.16	&	3.78	&	3.37	&	3.59	\\
	& 	median	&	4  	&	7	&	3	&	3	&	3	&	3	\\ \hline
\end{tabular}
	\caption{Mean and median number of local maxima, over $250$ simulated data sets of eight-variate ($J=8$) GMM estimation, $G=3$ and $p=(0.1,0.4,0.5)'.$ Initialization from ten random starts.}
	\label{tablocmax}
\end{table}

\FloatBarrier
In Table \ref{tabcrossval} we display results for different cross-validation settings. Perhaps surprisingly, this points out that the results are not sensitive to the choice of the cross-validations settings in terms of classification. Concerning computational time, this is maximal for $n_{\bar{S}}=n/2$ and $K=n,$ whereas it is minimal for $n_{\bar{S}}=n/10$ and $K=n/10.$ Interestingly however, we observe a systematic decrease in $c$ for lower levels of $n_{\bar{S}}:$ this could be due to the fact that accuracy in parameter estimation, with a larger training set, is higher. Hence, the prior information incorporated in the target matrix becomes less important - and the optimal $c$ relatively smaller.

\FloatBarrier
\begin{table}[h!]
	\centering
		\begin{tabular}{clccccccccccc}
		\hline \hline
			&		&	\multicolumn{3}{c}{$K=n/10$}		&&       \multicolumn{3}{c}{$K=n/5$}		&&	\multicolumn{3}{c}{$K=n$}		\\
\cmidrule{3-5} \cmidrule{7-9} \cmidrule{11-13}
	$n_{\bar{S}}=n/2$	&		&	conS	&	conN	&	cont	&&	conS	&	conN	&	cont	&&	conS	&	conN	&	cont	\\
			&		&		&		&		&&		&		&		&&		&		&		\\
			&	MAD	&	0.10	&	0.04	&	0.03	&&	0.10	&	0.04	&	0.03	&&	0.10	&	0.04	&	0.03	\\
			&	ARand	&	0.80	&	0.93	&	0.94	&&	0.81	&	0.93	&	0.94	&&	0.81	&	0.93	&	0.94	\\
			&	time	&	0.94	&	0.45	&	0.52	&&	1.29	&	0.68	&	0.73	&&	4.04	&	2.39	&	2.49	\\
			&	c	&	0.25	&	0.89	&	0.58	&&	0.25	&	0.89	&	0.58	&&	0.25	&	0.90	&	0.58	\\
		 	&		&		&		&		&&		&		&		&&		&		&		\\
		 	&		&		&		&		&&		&		&		&&		&		&		\\ \hline
			&		&	\multicolumn{3}{c}{$K=n/10$}		&&       \multicolumn{3}{c}{$K=n/5$}		&&	\multicolumn{3}{c}{$K=n$}		\\
\cmidrule{3-5} \cmidrule{7-9} \cmidrule{11-13}
	$n_{\bar{S}}=n/5$	&		&	conS	&	conN	&	cont	&&	conS	&	conN	&	cont	&&	conS	&	conN	&	cont	\\
			&		&		&		&		&&		&		&		&&		&		&		\\
			&	MAD	&	0.10	&	0.03	&	0.04	&&	0.09	&	0.03	&	0.04	&&	0.09	&	0.03	&	0.03	\\
			&	ARand	&	0.81	&	0.94	&	0.93	&&	0.82	&	0.94	&	0.93	&&	0.83	&	0.95	&	0.94	\\
			&	time	&	0.87	&	0.44	&	0.48	&&	1.11	&	0.65	&	0.69	&&	3.15	&	2.33	&	2.39	\\
			&	c	&	0.10	&	0.69	&	0.47	&&	0.10	&	0.70	&	0.47	&&	0.10	&	0.70	&	0.47	\\
		 	&		&		&		&		&&		&		&		&&		&		&		\\
		 	&		&		&		&		&&		&		&		&&		&		&		\\ \hline
			&		&	\multicolumn{3}{c}{$K=n/10$}		&&       \multicolumn{3}{c}{$K=n/5$}		&&	\multicolumn{3}{c}{$K=n$}		\\
\cmidrule{3-5} \cmidrule{7-9} \cmidrule{11-13}
	$n_{\bar{S}}=n/10$	&		&	conS	&	conN	&	cont	&&	conS	&	conN	&	cont	&&	conS	&	conN	&	cont	\\
			&		&		&		&		&&		&		&		&&		&		&		\\
			&	MAD	&	0.10	&	0.03	&	0.04	&&	0.09	&	0.03	&	0.04	&&	0.09	&	0.03	&	0.03	\\
			&	ARand	&	0.80	&	0.94	&	0.93	&&	0.82	&	0.95	&	0.94	&&	0.82	&	0.95	&	0.94	\\
			&	time	&	0.82	&	0.43	&	0.48	&&	1.01	&	0.63	&	0.67	&&	2.84	&	2.32	&	2.34	\\
			&	c	&	0.09	&	0.63	&	0.43	&&	0.09	&	0.63	&	0.43	&&	0.09	&	0.63	&	0.43	\\ \hline \hline
\end{tabular}
	\caption{Average results over 250 simulated data sets, each of sample size 100, of eight-variate ($J=8$) GMM constrained estimation, $G=3$ and $p=(0.1,0.4,0.5)',$ for different cross-validation settings. Initialization from ten random starts.}
	\label{tabcrossval}
\end{table}

\FloatBarrier

Table \ref{tabsepar} gives results on three different levels of class separation, for three and five components with small sample size ($n=50$). Both the unconstrained and the constrained methods have a relatively stronger effect on performance when class separation increases from $\text{sep}=1$ to $\text{sep}=2$ than when it increases from $\text{sep}=2$ to $\text{sep}=3$. Among all the approaches considered, whatever the class separation, cont yields on average the most accurate clustering.

\FloatBarrier
\begin{table}[h!]
	\centering
		\resizebox{0.75\hsize}{!}{\begin{tabular}{lllcccccc}
		\hline \hline
$\text{sep}=1$	&		&		&	homN	&	hetN	&	homt	&	conS	&	conN	&	cont	\\ \cmidrule{3-9}
	&	$p=(0.2,0.3,0.5)'$	&		&		&		&		&		&		&		\\
	&		&	MAD	&	0.22	&	0.40	&	0.13	&	0.32	&	0.22	&	0.13	\\
	&		&	Arand	&	0.57	&	0.23	&	0.73	&	0.42	&	0.57	&	0.74	\\
	&		&	time	&	0.09	&	0.06	&	0.11	&	1.75	&	0.72	&	0.79	\\
	&		&	$c$	&		&		&		&	0.66	&	0.95	&	0.84	\\
	&		&		&		&		&		&		&		&		\\
	&	$p=(0.1,0.4,0.5)'$	&		&		&		&		&		&		&		\\
	&		&	MAD	&	0.23	&	0.41	&	0.22	&	0.33	&	0.23	&	0.20	\\
	&		&	Arand	&	0.55	&	0.22	&	0.62	&	0.41	&	0.55	&	0.64	\\
	&		&	time	&	0.09	&	0.06	&	0.12	&	1.83	&	0.76	&	0.84	\\
	&		&	$c$	&		&		&		&	0.67	&	0.95	&	0.81	\\
	&		&		&		&		&		&		&		&		\\
	&	$p=(0.1,0.1,0.2,0.3,0.3)'$	&		&		&		&		&		&		&		\\
	&		&	MAD	&	0.38	&	0.55	&	0.34	&	0.43	&	0.37	&	0.34	\\
	&		&	Arand	&	0.40	&	0.15	&	0.45	&	0.34	&	0.41	&	0.46	\\
	&		&	time	&	0.11	&	0.05	&	0.19	&	2.37	&	1.13	&	1.22	\\
	&		&	$c$	&		&		&		&	0.69	&	0.95	&	0.94	\\
	&		&		&		&		&		&		&		&		\\
\cmidrule{3-9}
$\text{sep}=2$	&		&		&		&		&		&		&		&		\\
	&	$p=(0.2,0.3,0.5)'$	&		&		&		&		&		&		&				\\
	&		&	MAD	&	0.13	&	0.39	&	0.05	&	0.25	&	0.13	&	0.04	\\
	&		&	Arand	&	0.74	&	0.27	&	0.91	&	0.50	&	0.75	&	0.91	\\
	&		&	time	&	0.08	&	0.06	&	0.09	&	1.51	&	0.66	&	0.70	\\
	&		&	$c$	&		&		&		&	0.58	&	0.94	&	0.84	\\
	&		&		&		&		&		&		&		&		\\
	&	$p=(0.1,0.4,0.5)'$	&		&		&		&		&		&		&				\\
	&		&	MAD	&	0.15	&	0.37	&	0.14	&	0.27	&	0.14	&	0.10	\\
	&		&	Arand	&	0.72	&	0.28	&	0.77	&	0.49	&	0.73	&	0.83	\\
	&		&	time	&	0.09	&	0.06	&	0.11	&	1.66	&	0.72	&	0.76	\\
	&		&	$c$	&		&		&		&	0.61	&	0.95	&	0.81	\\
	&		&		&		&		&		&		&		&				\\
	&	$p=(0.1,0.1,0.2,0.3,0.3)'$	&		&		&		&		&		&		&		\\
	&		&	MAD	&	0.22	&	0.52	&	0.19	&	0.33	&	0.21	&	0.18	\\
	&		&	Arand	&	0.65	&	0.19	&	0.70	&	0.49	&	0.66	&	0.72	\\
	&		&	time	&	0.09	&	0.05	&	0.16	&	1.98	&	1.05	&	1.10	\\
	&		&	$c$	&		&		&		&	0.60	&	0.95	&	0.93	\\
	&		&		&		&		&		&		&		&		\\
\cmidrule{3-9}
$\text{sep}=3$	&		&		&		&		&		&		&		&		\\
	&	$p=(0.2,0.3,0.5)'$	&		&		&		&		&		&		&			\\
	&		&	MAD	&	0.10	&	0.37	&	0.03	&	0.22	&	0.09	&	0.03	\\
	&		&	Arand	&	0.81	&	0.31	&	0.94	&	0.56	&	0.81	&	0.95	\\
	&		&	time	&	0.08	&	0.06	&	0.09	&	1.45	&	0.66	&	0.67	\\
	&		&	$c$	&		&		&		&	0.55	&	0.94	&	0.83	\\
	&		&		&		&		&		&		&		&		\\
	&	$p=(0.1,0.4,0.5)'$	&		&		&		&		&		&		&			\\
	&		&	MAD	&	0.12	&	0.36	&	0.13	&	0.24	&	0.11	&	0.06	\\
	&		&	Arand	&	0.77	&	0.29	&	0.79	&	0.54	&	0.78	&	0.89	\\
	&		&	time	&	0.09	&	0.06	&	0.11	&	1.61	&	0.72	&	0.74	\\
	&		&	$c$	&		&		&		&	0.58	&	0.94	&	0.80	\\
	&		&		&		&		&		&		&		&		\\
	&	$p=(0.1,0.1,0.2,0.3,0.3)'$	&		&		&		&		&		&		&		\\
	&		&	MAD	&	0.15	&	0.49	&	0.14	&	0.28	&	0.14	&	0.12	\\
	&		&	Arand	&	0.77	&	0.23	&	0.79	&	0.55	&	0.78	&	0.81	\\
	&		&	time	&	0.09	&	0.05	&	0.15	&	1.81	&	1.04	&	1.08	\\
	&		&	$c$	&		&		&		&	0.55	&	0.95	&	0.93	\\ \hline

\end{tabular}}
	\caption{Average results over 250 simulated data sets, each of sample size 50, of eight-variate ($J=8$) GMM estimation, three and five components, for different class separation levels. Initialization from ten random starts.}
	\label{tabsepar}
\end{table}

\FloatBarrier

\section{Empirical application: the wine data set}\label{wineapp}
In the present Section we evaluate the seven algorithms on the basis of a data set available at \url{http://www.ics.uci.edu/~mlearn/MLRepository.html}. These data are the results of a chemical analysis of three types of wine - \emph{Barolo, Grignolino} and \emph{Barbera} - grown in the same region in Italy. The analysis determined the quantities of 13 constituents found in each of the three types of wines: alcohol, malic acid, ash, alcalinity of ash, magnesium, total phenols, flavanoids, nonflavanoid phenols, proanthocyanins, color intensity, hue, OD280/OD315 of diluted wines, and proline.

All of the six algorithms have been initialized from the same 50 random starts, assuming $G=3.$ The selected solutions are the ones with the highest likelihood. The cross-validation scheme is the same used in the simulation study. Results are shown in Table \ref{wine}.\\

\begin{table}[h!]
	\centering
		\begin{tabular}{lcccccc}
		\hline
		\hline
	&	homN	&	hetN	&	homt	&	conS	&	conN	&	cont	\\ \hline
	&		&		&		&		&		&		\\
ARand	&	0.92	&	0.39	&	0.87	&	0.54	&	0.92	&	0.93	\\
time	&	2.22	&	1.89	&	1.68	&	23.34	&	11.32	&	10.71	\\
$c$	&		&		&		&	0.45	&	0.89	&	0.57	\\ \hline \hline
\end{tabular}
	\caption{Comparison of the 6 algorithms in terms of ARand, computational time and optimal $c.$}
	\label{wine}
\end{table}
\indent The homoscedasticity assumption seems to fit well the data. The constrained approach conN equals homN in terms of ARand, whereas cont yields an ARand of 0.93, compared to 0.87 of homt. Confirming the results obtained in the simulation study, cont seems to be the most accurate approach among the ones considered in this work.

Interestingly, however, all of the constrained approaches improve upon the unconstrained heteroscedastic approach.

\section{Discussion}\label{concl}
In this paper we have proposed affine equivariant constraints for the class conditional covariance matrices of multivariate GMM in order to circumvent the well-known issue of degenerate and spurious solutions in ML estimation. Our approach generalizes the sufficient condition for Hathaway (1985)'s constraints to hold as formulated by Ingrassia (2004). Previous constrained approaches lacked affine equivariance and suffered the choice of an optimal finite-sample scale balance ($c$). The setup we propose is such that the class specific covariance matrices are shrunk towards a pre-specified matrix $\boldsymbol{\Psi}.$ We have been able to show that this yields a clustering method which is equivariant with respect to linear affine transformations of the data, provided that $\boldsymbol{\Psi}$ is changed accordingly.

A natural choice for the shrinkage target matrix, whenever \emph{a priori} information on the covariance structure of the components is not available, seems to be the covariance matrix of a homoscedastic mixture of normals. For a given choice of the target matrix, we let the data decide, through the constant $c$,  how close to the target the final clustering will be. The tuning constant $c$ is chosen by cross-validation. We have also shown that, given a matrix $\boldsymbol{\Psi},$ our constrained ML estimate can be computed by applying the algorithm of Ingrassia and Rocci (2007) to the data appropriately linearly transformed. This allows us to interpret our proposal as a way to decide how to standardize the data before applying Ingrassia (2004)'s constraints.

The validity of the proposal has been assessed through a simulation study and an empirical example. All constrained approaches yield more accurate estimates than the unconstrained ones. More specifically, cont has been shown to be the best among the constrained approaches this work has been concerned with. This is not surprising, since a random vector conditionally distributed as a Gaussian mixture, given random inverse Wishart covariance matrices, has a marginal homoscedastic mixture of Student $t$'s distribution.

Given an affine transformation of the data, the equivariance property of the method is guaranteed if also $\boldsymbol{\Psi}$ is adapted accordingly. This requires that the methods used to estimate $\boldsymbol{\Psi}$ from the data be also equivariant. This is the case for the sample covariance matrix and the homoscedastic model, for which \ref{eq:equiGauss} and \ref{eq:equiGaussmix} apply. For a homoscedastic mixture of Student $t$'s, this can also be shown expressing each marginal as a combination of a multivariate Gaussian and Gamma random variables. Then affine equivariance results by applying \ref{eq:equiGauss} (Roth, 2013). All in all, different choices of $\boldsymbol{\Psi}$ can as well be considered, according to the data specificity: however, in order for the method to preserve equivariance, also the method used to estimate $\boldsymbol{\Psi}$ from the data has to be equivariant.

The equivariant method developed in Gallegos and Ritter (2009a; 2009b) and extended in Ritter (2014) requires to obtain all local maxima of the trimmed likelihood. Our method has the virtue of being easily implementable with a minimal extra computational effort, as we have shown in the simulation study and in the empirical example.

There are cases where the clustering model might assume a specific structure on the relationship between the variables, like local independence (within-cluster diagonal matrices). Such a model is not affine equivariant because some (non diagonal) affine transformations on the data might destroy the local independence. In cases like these, the affine equivariance property of the constraints is not required. Yet our approach can be applied using a diagonal matrix as target. This would prevent the likelihood from degenerating, still improving upon the unconstrained algorithm thanks to the cross-validation strategy we have proposed. Clearly, when all variables in a data set are measured in a common scale, non equivariant constraints are a competitive choice.

An additional issue, pointed out by both the simulation study and the empirical example, is the computational time cross-validation requires to select an optimal $c.$ Whether different cross-validation schemes can speed up the constrained routines can be a topic for future research.

\newpage