I recall to have look at this formula sometimes ago (#21481) and reformulated the mathematical section: https://scikit-learn.org/stable/modules/linear_model.html#lasso-lars-ic

The n constant is coming from making the assumption to have a Gaussian model and plug it in the BIC formula. I don't recall that it was an issue there.

The current formula is:
return -2 * self.score(X) * X.shape[0] + self._n_parameters() * np.log(X.shape[0])

The proposed correction is:
return self._n_parameters() * np.log(X.shape[0]) -2 * self.score(X)

The general formula of the BIC is
BIC = k*ln(n)-2ln(L)

Now let

k = self._n_parameters()
n = X.shape[0]
ln(L) = self.score(X)

Then

BIC = k*ln(n)-2ln(L)
    = self._n_parameters() * np.log(X.shape[0]) -2 * self.score(X)

Which in turn is the proposed corretion.

Example

When one computes the silhouette scores (to evaluate the best number of clusters) in contrast to the BIC, one will retrieve

For n_clusters = 2 The average silhouette_score is : 0.7669465622459366
For n_clusters = 3 The average silhouette_score is : 0.41538251879888394
For n_clusters = 4 The average silhouette_score is : 0.36208909622898705
For n_clusters = 5 The average silhouette_score is : 0.06036439944496741
For n_clusters = 6 The average silhouette_score is : 0.08875521208115827
For n_clusters = 7 The average silhouette_score is : 0.13079963966127484

So in conclusion, silhoutte indices that 2 clusters are optimal.
As a reminder from above: GMM with BIC suggest 5 clusters.
Now have a look at the clustering for n=2 and n=5.

For convenience here are the BIC for 2,...,7 clusters again. The score for n=2 clusters is the best here and consistent with the visual clustering and silhouette.

[57.180072607292026,
86.96960303077948,
116.92208468419805,
145.6026996264713,
175.62029248920572,
206.73427255638083]

My current understanding is that since GaussianMixture.score is defined to:

Compute the per-sample average log-likelihood of the given data X.

then the total log-likelihood (the actual factor BIC uses) is then attained as self.score(X) * X.shape[0] as is currently the case.

@Micky774 agreed. Then I have absolutely no idea what's going wrong here. In any case, here is the colab notebook for reproducibility.

@Micky774 agreed. Then I have absolutely no idea what's going wrong here. In any case, here is the colab notebook for reproducibility.

I'm not sure there even is anything wrong here -- AIC/BIC and the silhouette score are all heuristics and it isn't altogether surprising to me that they don't agree in this case. These are good rule-of-thumb heuristics of course, but none are nearly robust enough where this divergence in opinion raises a red flag for me.

I am going to close this one, since it seems there is no issue in scikit-learn and this turned into a more generic machine learning discussion.

Thanks @Micky774 for the input!