code: def _half_log_det(self, M):
          
          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)
docs: 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.