def _half_log_det(self, M):
        """ Return log(|M|)*0.5. For positive definite matrix M
            of more than 2 dimensions, calculate this for the
            last two dimension and return a value corresponding
            to each element in the first few dimensions.
        """
        chol = np.linalg.cholesky(M)
        if M.ndim == 2:
            return np.sum(np.log(np.abs(np.diag(chol))))
        else:
            return np.sum(np.log(np.abs(np.diagonal(
                chol, axis1=-2, axis2=-1))), axis=-1)