# Scikit learn - GMM log-likelihood. Why use Cholesky's precision matrix instead of covariance matrix?

This is my first post, please let me know if I am not being clear.

I am trying to understand the sklearn.mixture.GaussianMixture.score(X).

As I understand that the following is the equation for log likelihood:

$\ln&space;p(x|\mu_k,\Sigma)&space;=&space;-\frac{1}{2}\ln(2\pi)&space;-&space;\frac{1}{2}\ln\left&space;|&space;\Sigma&space;\right&space;|&space;-\frac{1}{2}(x-\mu_k)^T\Sigma^{-1}_k(x&space;-&space;\mu_k)$

In my case, I want to calculate log likelihood for a sklearn.mixture.GaussianMixture(covariances_type = 'diag'). Below is the code that is calculating $\ln&space;\left|&space;\Sigma\right|$.

def _compute_log_det_cholesky(matrix_chol, covariance_type, n_features):
[...]
elif covariance_type == 'diag':
log_det_chol = (np.sum(np.log(matrix_chol), axis=1))
[...]
return log_det_chol


This post explains how to use Cholesky's precision matrix. My question is:

• why are are not using sklearn.mixture.GaussianMixture.fit(X).covariances_.?

Since we are using sklearn.mixture.GaussianMixture(covariances_type = 'diag'), the output of sklearn.mixture.GaussianMixture.fit(X).covariances_ is k x m. Where k is number of components and m is the number of features.

Working with Cholesky factors is cheaper because determinants of triangular matrices are given by the product of the elements on the diagonal. This massively reduces the computation time for the determinants because we effectiely get to skip the computations which are needed to compute the determinant in the general case.

Working with Cholesky factors also economizes the computations because we are also computing Cholesky factors when working with the $$\Sigma^{-1}$$ portion of the log-likelihood, so we're re-using the same intermediate steps instead of computing two different, expensive things (the inverse of a matrix and the determinant of a matrix).

The code is just collapsing some operations to avoid redundant computations. In the full code, the determinant for covariance_type=="diag" has the -0.5 factor combined in log_det.

def _estimate_log_gaussian_prob(X, means, precisions_chol, covariance_type):
log_det = _compute_log_det_cholesky(
precisions_chol, covariance_type, n_features)

[...]

elif covariance_type == 'diag':
precisions = precisions_chol ** 2
log_prob = (np.sum((means ** 2 * precisions), 1) -
2. * np.dot(X, (means * precisions).T) +
np.dot(X ** 2, precisions.T))

[...]
return -.5 * (n_features * np.log(2 * np.pi) + log_prob) + log_det


Since the Gaussian log-likelihood doesn't ever use the covariance matrix $$\Sigma$$, but instead only uses its inverse $$\Sigma^{-1}$$ and a negative multiple of its determinant $$-\frac{1}{2} \log | \Sigma|$$, we actually only need to work with Cholesky factors of the precision matrix.

You can prove that this is correct using the properties of Cholesky factors, determinants, and matrix inverses. Our goal is to find a way to write $$-\frac{1}{2} \ln |\Sigma|$$ in terms of the Cholesky factor of the precision matrix: $$\Sigma^{-1} = A^\top A$$

\begin{align} \left| \Sigma^{-1} \right| &= \left| A^\top A \right| \\ &= \left| A \right| \left| A^{\top} \right| \\ &= \left| A \right|^2 \\ \ln \left| \Sigma^{-1} \right| &= \ln \left| A \right|^2 \\ \ln \left| \Sigma \right| &= \ln \left| A \right|^{-2} \\ &= -2 \ln \left| A \right| \\ -\frac{1}{2}\ln \left| \Sigma \right| &= \ln \left| A \right| \end{align}

This proves that the code is algebraically correct for all Cholesky factors of nonsingular $$\Sigma$$, including the special case when $$\Sigma$$ is diagonal.