Achieving Optimal Clustering in Gaussian Mixture Models with Anisotropic Covariance Structures

We thank you for your valuable comments and remarks.

Originality: A related work is

Thank you for pointing this out. Our hard-EM algorithm indeed shares similarities with the soft one. Specifically, our algorithm modifies the E-step of the soft-EM by implementing a maximization step that hard assigns data points to clusters instead of calculating probabilities. We will make this clearer in the final version of the paper.

Suggestions:

Thank you for pointing these out. We will surely address them in the final version.

What about individual weights (mixing proportions) of the Gaussian components?

In our numerical study, we initially opted for equal cluster sizes for simplicity, but our model and analysis are fully capable of accommodating variable cluster sizes. To better demonstrate the flexibility and applicability of our approach, we will have a variety of cluster sizes in the numerical settings in the final version of the paper.

"decent" initial guess required

Thank you for the question. In our manuscript, the term 'decent initial guess' refers to the need for initial cluster centers to be sufficiently close to the ground truth so that our algorithm achieves the rate-optimal result. It is due to that our theoretical analysis requires the initialization being within a specific proximity to the true parameters. In the final version of our paper, following your suggestion, we will explicitly detail the dependencies to provide a clearer understanding of when and how our algorithm can be expected to perform optimally.

We thank you for your valuable comments and remarks.

The paper misses some high-level intuition, and several proofs are tough to parse.

Thanks for pointing this out. In the final version, we will enhance our discussion to better articulate the derivation and significance of key terms. The "ideal error" is defined as the error remaining after one algorithm iteration with known ground truth $z^*$. It represents the minimum error under ideal conditions. The actual error emerges when iterating from an estimation $z$, with the difference between actual and ideal errors expressed through the terms $F$, $G$, and $H$: $F$ includes terms related to noise $\epsilon_j$, illustrating the impact of measurement noise. $G$ covers estimation errors for cluster centers ($\hat{\theta}_a(z) - \hat{\theta}_a(z^*)$) and covariance matrices ($\hat{\Sigma}(z) - \hat{\Sigma}(z^*)$), showing the effect of parameter estimation inaccuracies. $H$ contains all other terms from additional error sources.

The interpretation of the minimax rate for model 2 is quite complex.

Thanks for the suggestion. We will include connections with Chernoff information in the final version to clarify the understanding of the minimax rate. Regarding SNR', although it diverges from the traditional signal-to-noise ratio, we use this term as it extends the classical definition $\frac{||\theta_1^* - \theta_2^*||^2}{\sigma^2}$ in the context of isotropic Gaussian noise.

The assumption on the loss of the initialization is quite strong

We acknowledge that this assumption appears strong; it is primarily driven by technical challenges encountered in our analytical framework.

In large dimensions (say d >> n), the lower bound derived here cannot be attained

You are correct that our paper does not cover the large dimensional case. In the final version, we will add more comments and references to make it clearer.

Minor comments

Thank you and we will correct them in the final version.

Is the assumption SNR / log k \gg 1 needed?

Thank you for the question. On one hand, this condition helps simplify the complexity for establishing the lower bound. As noted in the literature you referenced, a more refined analysis might allow us to relax this assumption to simply $\text{SNR} \gg 1$. On the other hand, this condition is essential to establish a matching upper bound. The ideal error analysis does not explicitly require this assumption; however, it becomes necessary when considering the aggregate impact of errors $F$, $G$, and $H$. These components introduce multiple factors of $k$ that appear before the desired exponential term, and the assumption ensures that these factors remain negligible.

Line 42 the sentence:

Thank you for the question. To address it, consider a simpler case with two Gaussian distributions: whether the SNR calculated from $N(\theta^*_1, \Sigma^*_1)$ and $N(\theta^*_2, \Sigma^*_2)$ is always larger than that from $N(\theta^*_1, \Sigma^*_1)$ and $N(\theta^*_2, \Sigma^*_1)$. In this case, the answer to your question is not necessarily yes, as demonstrated by the following counterexample. When $\theta^*_1=(0,0)$, $\theta^*_2=(1,0)$, $\Sigma^*_1 = I_2$ and $\Sigma^*_2$ is a diagonal matrix with 2 in its (1,1) entry and 1 in its (2,2) entry, one can verify SNR is not larger. In fact, the effect of replacing $\Sigma^*_2$ with $\Sigma^*_1$ on SNR depends on the shapes of $\Sigma^*_1$ and $\Sigma^*_2$ and the direction of $\theta^*_2 - \theta^*_1$. This is different from the sub-Gaussian setting. The rationale why $N(0, \sigma^2 I_d)$ leads to the smallest SNR among sub-Gaussian distributions with variance proxy $\sigma^2$ is that it is flatter in all directions compared to any other sub-Gaussian distribution.

I am not sure how to obtain the second to last inequality of the inequations starting line 473.

The inequality does not result directly from Cauchy-Schwarz. Instead, it first formulates a quadratic form and then upper bounds it by the operator norm of the matrix. Let $x$ be a vector. The inequality is essentially about showing $\sum_j |\epsilon_j^Tx|^2$ can be upper bounded by $||x||^2 ||\sum_j \epsilon_j \epsilon_j^T||$. This can be proved by $\sum_j |\epsilon_j^Tx|^2 = \sum_j (x^T\epsilon_j)(\epsilon_j^T x) = \sum_j x^T(\epsilon_j \epsilon_j^T)x = x^T(\sum_j \epsilon_j \epsilon_j^T)x \leq ||x||^2 ||\sum_j \epsilon_j \epsilon_j^T||$. In the final version, we will add this intermediate argument.

I do not understand the last statement

Thanks for the feedback. We agree it needs clarification. Our intent was to indicate a potential relaxation of the well-conditioned assumption on the covariance matrices' condition numbers. In the final version we will make it clearer.

Albeit this is already a long and dense paper, I would have enjoyed a section

Thank you for your suggestion. In the final version of the manuscript, we will incorporate a new section dedicated to evaluating the performance of the algorithms on real datasets. This section will also include practical guidelines for practitioners on selecting the most suitable algorithm.

While I went quite far down the proof, I found them hard to read.

Thanks for the comment. In the final version, we will add more details to make the proof more accessible. For example, in Appendix C, we will give an overview of lemmas to be proved.

We thank you for your valuable comments and remarks.

My major concern is about its experiments.

Thank you for your comment. We plan to include a new experiment in the final version of our paper using a real dataset from the Fashion-MNIST collection. This experiment will focus on the clustering of two classes, T-shirt/top and Trouser, each comprising 6000 images. Numerical results show that Algorithm 2 achieves a misclustering error rate of 5.7%, which outperforms the vanilla Lloyd’s algorithm (8%). This addition will provide a more comprehensive evaluation of our methods, highlighting their effectiveness in practical scenarios.

In Algorithm 2, it needs to compute the inverse and determinant of the covariance matrices

Thanks for the comment. You are correct in noting that it incurs additional computational overhead due to the necessity of computing the inverse and determinant of covariance matrices. The time complexity of Algorithm 2 is $O(nkd^3T)$, where $n$ is the number of points, $d$ is the dimensionality of each data point, $k$ is the number of clusters, and $T$ is the number of iterations. This contrasts with the vanilla Lloyd's algorithm, which has a lower time complexity of $O(nkdT)$. The increased complexity is primarily due to the matrix operations in $d$ dimensions, which scale as $O(d^3)$ for matrix inversion and determinant computation. To provide a clearer perspective on the performance impact, our experiments show that at a dimensionality of 5, Algorithm 2 requires approximately twice the computation time of the vanilla Lloyd’s algorithm. This ratio increases to approximately 14 when the dimensionality is increased to 100. In the final version of the paper, we will include a detailed time complexity analysis and experimental results to illustrate how the overhead scales with increased dimensionality.