# Slow mixture fitting when having ill-conditioned covariance matrices. How to make it faster?

I'm using Mathematica 9.0 to fit multidimensional observations with a gaussian mixture. To do so I first adapted these steps to ensure that the mixture parameters are well defined (see https://mathematica.stackexchange.com/a/16833/12672). Runing time was unacceptable (I waited some hours and aborted it afterwards), so I tried the fitting in Matlab. Matlab instantly retrieved an error message from the clustering algorithm that states "ill-conditioned covariance created at iteration 3". To overcome this I tried 2 solutions concerning the covariance matrices definitions following the suggestions in https://stackoverflow.com/a/16263971 and http://www.mathworks.com/matlabcentral/answers/100210.

To reproduce my problem and illustrate my experiments here I enclose some data, part of my original data, this being 7 samples each 5-dimensional, which I try to fit in a 2 gaussian mixture:

data = {{0, 2.76, 10.68, 12.4133, 9.26667}, {0, 0, 0, 0, 1.13333}, {6.06667, 8.32, 6.86667, 6.6, 4.49333}, {1.25333, 5.96, 7.49333, 5.81333, 5.05333}, {0, 0, 1.13333, 6.42667, 7.41333}, {0, 0, 0, 3.25333, 3.82667}, {0.76, 6.28, 10.28, 9.05333, 5.12}};
data // MatrixForm
d = Dimensions[data]
nDim = Last[d];
nMixtures = 2;


One solution is to use a single covariance matrix for all the mixtures. What follows is the implementation in Mathematica:

For[k = 1, k < nMixtures + 1, k++,
m[k] = Table[mm[k, i], {i, nDim}]];(*means*)
For[k = 1, k < nMixtures, k++,
w[k] = Exp[ww[k]]/((1 + Exp[ww[k]])*(nMixtures - 1))]; (*mixture weights*)
c = Table[If[i >= j, cc[i, j], 0], {i, nDim}, {j, nDim}];(*triangular matrix*)
s = c.Transpose[c];(*covariance matrix*)
M = Append[Array[w, nMixtures - 1], 1 - Total[Array[w, nMixtures - 1]]];(*all weights, which sum up to 1*)
Q = {};
For[i = 1, i < nMixtures + 1, i++,
Q = Insert[Q, MultinormalDistribution[m[i], s], i]];(*all mixtures, which share the same covariance matrix*)


The corresponding fitting takes 85 seconds in my PC using the following command:

edf_shared = EstimatedDistribution[data, MixtureDistribution[M, Q]]; // Timing


The same implementation in Matlab takes 0.0047 seconds:

data = [0,0,6.0667,1.2533,0,0,0.76; 2.76,0,8.32,5.96,0,0,6.28; 10.68,0,6.8667,7.4933,1.1333,0,10.28; 12.4133,0,6.6,5.8133,6.4267,3.2533,9.0533; 9.2667,1.1333,4.4933,5.0533,7.4133,3.8267,5.1200]';
edf_shared = gmdistribution.fit(data,nMixtures,'SharedCov',true);


Another way to go around this problem is, instead of using the same covariance matrix for all gaussians, add a small epsilon to the diagonal of the covariance matrices. In Mathematica it takes 29,247 seconds (about 8 hours):

For[k = 1, k < nMixtures + 1, k++,
m[k] = Table[mm[k, i], {i, nDim}]];(*means*)
For[k = 1, k < nMixtures, k++,
w[k] = Exp[ww[k]]/((1 + Exp[ww[k]])*(nMixtures - 1))];(*mixture weigths*)
For[k = 1, k < nMixtures + 1, k++,
c[k] = Table[If[i >= j, cc[k, i, j], 0], {i, nDim}, {j, nDim}]];(*triangular matrices*)
epsilon = 0.00005;
epsilonMatrix = DiagonalMatrix[ConstantArray[epsilon, nDim]];(*diagonal matrix*)
For[k = 1, k < nMixtures + 1, k++,
s[k] = c[k].Transpose[c[k]] + epsilonMatrix];(*covariance matrices*)
M= Append[Array[w, nMixtures - 1], 1 - Total[Array[w, nMixtures - 1]]];
Q = {};
For[i = 1, i < nMixtures + 1, i++,
Q = Insert[Q, MultinormalDistribution[m[i], s[i]], i]];(*all mixtures, with regularized covariance matrices*)
edf_epsilon = EstimatedDistribution[data,MixtureDistribution[M, Q]]; // Timing


The corresponding implementation in Matlab takes 0.0055 seconds:

edf_epsilon = gmdistribution.fit(data,nMixtures,'Regularize', 1e-5);


Neadless to say, when evaluating the estimated distribution function on a point, execution time is also an issue. The implementation using shared covariances in Mathematica requires 276 seconds, while in Matlab is just 0.41 seconds. The first approach involving a shared matrix, when applied to this new point takes 107 seconds in Mathematica and 0.20 seconds in Matlab.

I also tried playing with the EstimatedDistribution function parameters such as WorkingPrecision, AccuracyGoal, ParameterEstimator, etc., I didn't see any substantial improvement. At most, when I reduced the WorkingPrecision to 5 I managed to cut down the estimation time to just 40 seconds.

My question is whether there is a better way to do this, I mean to estimate the parameters of a gaussian distribution in Mathematica while coping for ill-conditioned covariance matrices. By better I mean faster while keeping a relatively good accuracy in the estimations, in the Matlab direction I enclosed here. Perhaps there is some issue with my code, or the approaches themselves are intrinsically inconvenient in Mathematica.

Any suggestion will be welcome, thanks in advance!