If your samples dimensionality is less than the vector space dimensionality, singular matrices may arise. If you have less samples than $d+1$ (when $d$ is your dimensionality), this situation will even necessarily arise: $k+1$ samples span at most a $d$ dimensional hyperplane. Given such a small sample, you obviously cannot compute a variance in the orthogonal space.
This is why it's common to not use literal PCA, but instead perform singular value decomposition, which can be used to compute the pseudoinverse of a matrix.
If the matrix is invertible, the pseudoinverse will be the inverse.
However, if you are seeing non-invertible matrixes, chances are that your distance from the cluster will be meaningless if the vector is outside of the hyperplane the cluster repesents, because you do not know the variance in the orthogonal space (you can think of this variance as 0!) SVD can compute the pseudoinverse, but the "variances" will still be not determined by your data.
In this case, you should probably have been doing global dimensionality reduction first. Increasing the sample size will only help when you actually have non-redundant dimensions: no matter how many samples you draw from a distributions with $y=x$, the matrix will always be non-invertible, and you will not be able to judge the deviation $x-y$ with respect to a standard deviation (which is 0).
Furthermore, depending on how you compute the covariance matrix, you might be running into numerical issues due to catastrophic cancellation. The simplest workaround is to always center the data first, to get zero mean.