Gaussian Mixture Models and Expectation-Maximization › Disadvantages

The likelihood surface is difficult to optimize. With multiple components the log-likelihood is not concave, so EM can settle near different stationary points. Convergence of parameters requires additional regularity; likelihoods can also approach singular boundary configurations. Different initializations can produce visibly different answers with different likelihoods, and there is no certificate telling us that the best one was found. Multiple restarts and a structured initialization such as k-means are useful safeguards, although neither guarantees the best fit.

For an unrestricted mixture with at least two components, the sample likelihood is unbounded above. Put one component mean on an observation and shrink its covariance to zero while another nondegenerate component gives the remaining observations positive density. The first contribution diverges, even though this is a degenerate density fit. This is not a numerical quirk that better code removes; it is a property of the objective. The standard remedy is to add a small constant $\epsilon$ to the diagonal of every covariance matrix, that is, replace $\boldsymbol{\Sigma}_k$ with $\boldsymbol{\Sigma}_k + \epsilon\mathbf{I}$, where $\epsilon>0$ is a positive regularization strength (often called reg_covar in code) and $\mathbf{I}$ is the $d \times d$ identity matrix. For a positive-semidefinite scatter matrix, this imposes a positive covariance eigenvalue floor and caps the density. It prevents singular factorization but does not establish global optimality or exact EM monotonicity. The size of $\epsilon$ is therefore a real modelling decision, not a rounding detail.

There are further costs. The method assumes each component is Gaussian. Within a fixed feature space, approximating strongly skewed or heavy-tailed groups may require extra Gaussian components, which can make interpretation harder. Components are not the same thing as real-world groups: a single natural group with a non-Gaussian shape may be split across two fitted components, and two natural groups may share one. Component labels are arbitrary: permuting the component indices gives exactly the same likelihood, so comparisons across runs require matching components before comparing parameters. Finally, the computational cost grows with the square of the dimension: evaluating the full-covariance densities costs $O(nKd^2)$ after $O(Kd^3)$ factorizations, and forming the weighted covariance scatters costs $O(nKd^2)$. Here $n$ is the observation count, $d$ the feature count, and $K$ the number of components. For a thousand well-populated clusters in high dimensions, k-means remains far cheaper.