Use of Moore-Penrose Inverse in Likelihood Computation I've been running MLE on a mixed-Gaussian model (let's just focus on the multivariate Gaussian case for now). In order to take the inverse of the covariance matrix, I've been using MATLAB's pinv() command, which is the Moore-Penrose pseudo-inverse, $\Sigma^+$. The likelihood is then. 
$l(\mu,\Sigma;Y) = (2\pi)^{-n/2} |\Sigma|^{-n/2} \exp\big(-\frac{1}{2} \sum\limits_{t=1}^T (Y_t - \mu)^\prime \Sigma^{+} (Y_t - \mu)\big)$
What I've been observing are cases where the optimization chooses a very small error variance (or close to singular in the multivariate case), even when no variable is close to being perfectly predicted by the model. In other words, the variances of $(Y - \mu)$ are all of the same order of magnitude, yet the optimization algorithm is still able to blow up the likelihood via $|\Sigma|^{-n/2}$, because $\Sigma^+ << \Sigma^{-1}$.
Is it known that using the MP-Inv is problematic in the computation of the likelihood? Is the what I am saying generally making sense? How would I test this?
 A: I think this are three more fundamental issues here. 


*

*We do not need the Moore-Penrose pseudo-invense (to solve y' * inv(S) * y).

*We do not need the determinant of $\Sigma$ (to compute the term det(S)).

*We do not need to work with the likelihood directly (to get the MLE -  Log-likelihood is fine too). So instead of using the likelihood as in the question's main body we can use the logarithm of it and  get $\log L = -0.5 y^T S^{-1} y - 0.5 \log(|S|) - 0.5\log(2\pi)n$ (see the section on the Likelihood function of the MVN in the Wikipedia article here for more details)


In particular, the first two points are due to the fact that $\Sigma$ is positive definite and it can be re-expressed as $\Sigma = LL^T$, i.e. we can work with the Cholesky decomposition of it. The third point allows us to simplify our lives even further by allowing us to avoid any exponentiation procedures.. As such, assuming that our covariance matrix is S,  we can:


*

*We can make the solution of the system inv(S)* y faster by using the Cholesky decomposition - this step can be done implicitly in MATLAB by using the backslash operator \, i.e. the first term becomes  -(1/2)*y'*(S\y). 

*Avoid the need to exponentiate these results. 

*Avoid computing the determinant of $\Sigma$ by using the fact that the sum of the logarithms of the diagonal entries of the Cholesky decomposition of a matrix S equals half of the logarithm of determinant of matrix S. i.e. 0.5 * log(det(S)) and sum(log(diag(chol(S)))) are equal up to numerical precision (or $0.5\log(|S|)  = \sum_i \log(L_{i,i})$). 


A simplified version of the code in MATLAB would be something like:
  L = chol(S)';     
  alpha =  (S)\Y; % solves L'*alpha=y 
  logLikelihood = -(1/2)*Y'*alpha - sum(log(diag(L))) - (N/2)*log(2*pi);

I hope it is clear that while the Moore-Penrose pseudo-inverse itself is not problematic, the need to exponentiate any results as well as compute the determinant of a full matrix, make this procedure suboptimal. As a first reference on this matters, I have found Rasmussen and William's  Gaussian Processes for Machine Learning invaluable. Especially Chapter 5 on Model Selection and Adaptation of Hyperparameters has a wealth of information on the computations around Gaussian likelihoods.
As noted by @whuber, there are cases that using the Moore-Penrose pseudo-inverse might result into an sub-optimal estimator because the influence of certain data-point might be ignored. This can happen if we are dealing with a semi-positive definite covariance matrix $S$. 
A proper fix in that case would be to include a minimal amount of additive noise variation $\sigma_{noise}^2$ such that $S = K_{data} + \sigma_{noise}^2I$. Here, $K_{data}$ is the data-derived original covariance. This is actually how R&W formulate the computation in their book where the noise variance $\sigma_{noise}^2$ is explicitly included as a (hyper-)parameter to optimise.
