Why does the normalization term disappear when computing the MLE of decomposed Gaussians

Question

Computing the Maximum Likelihood Estimator of Gaussians in arbitrary finite Hilbert spaces seems no easy task and I must admit to lamentably fail at it. The classical theory most often relies on matrix notations, but I cannot make it work when decomposing the operators into their basis.

Setup

We consider a set of observations $\{x_n \in \mathcal{H}\}_{n=1}^N$ with $(\mathcal{H},\langle\cdot,\cdot\rangle)$ a Hilbert space of finite dimension $d$. The adjoint of any $x$ is defined as $x^\ast = \langle x,\cdot\rangle \in \mathcal{H}^\ast \sim \mathcal{H}$ with $\mathcal{H}^\ast$ the Fréchet-Riesz dual space of $\mathcal{H}$. If we define the distribution of a (centered) Gaussian on $\mathcal{H}$ as $$ p(x) = \frac1Z\exp\left(-\frac12x^\ast \Sigma^{-1} x\right), $$ and decompose the variance as $\Sigma = \sum_{i=1}^d \sigma_i^2 u_i u^\ast_i$ with $\{u_i\}_{i=1}^N$ an orthonormal basis of $\mathcal{H}$ and the eigenvectors $\sigma^2_i \in \mathbb{R}_{>0}$. We thus assume the variance to be of full rank. By consequence, we can deduce that the normalization term is equal to $$ Z = (2\pi)^{d/2}\left(\prod_{i=1}^d\sigma_i^2\right)^{1/2}. $$

Maximum Likelihood Estimator (MLE)

We can now compute the MLE based on our observations. More formally, we consider the parameters to be optimized as $\theta = (\{\sigma_i^2\}_{i=1}^N,\left\{u_i\}_{i=1}^d\right) \in \Theta$ with $\Theta = \mathbb{R}_{\>0}^d \times V_d(\mathcal{H})$ and $V_d(\mathcal{H})$ the Stiefel manifold on $\mathcal{H}$ (the set of the sets of $d$ orthonormal vectors of $\mathcal{H}$). $$ \min_{\theta \in \Theta} \prod_{n=1}^N p(x_n; \theta) = \min_{\theta \in \Theta} \log\left(\prod_{n=1}^N p(x_n; \theta) \right) = \min_{\theta \in \Theta} \sum_{n=1}^N\log\left( p(x_n; \theta) \right) = \min_{\theta \in \Theta} \mathcal{L}(\theta). $$

We refer to $\mathcal{L}(\theta)$ as the likelihood function and it is given by:

$$ \begin{eqnarray} \mathcal{L}(\theta) &=& -\frac N2\left\{d\log(2\pi) + \sum_{i=1}^d \log(\sigma_i^2) + \frac1N \sum_{n=1}^Nx^\ast \Sigma^{-1} x \right\}, \\ &=& -\frac N2\left\{d\log(2\pi) + \sum_{i=1}^d \log(\sigma_i^2) + \frac1N \sum_{n=1}^N\sum_{i=1}^d\sigma^{-2}_i(x^\ast u_i)^2\right\}. \end{eqnarray} $$

Issue

When computing the stationnary points of the likelihood in function of $u_i$, we get $$ \frac{\partial \mathcal{L}}{\partial u_i} = 0 \Longleftrightarrow \sigma_i^{-2} u_d^\ast \left(\frac1N\sum_{n=1}^N x\tilde{x}^\ast \right) = 0, $$ which is absurd as we should get $\sigma_i^{-2} u_d^\ast \left(\frac1N\sum_{n=1}^N x\tilde{x}^\ast \right) = 1$ in order to get $\{(\sigma_i^2,u_i)\}_{i=1}^d$ to be the eigenpairs of $\frac1N\sum_{n=1}^N xx^\ast$. The problem according to me is that the normalization now is independent of the basis $\{u_i\}_{i=1}^N$, which was not the case with the matrix notation. I don't understand however where the error is in an algebraic point of view.

I thank you in advance for your help and insights !

Answer 1

$\newcommand\th\theta\newcommand\si\sigma\newcommand\Si\Sigma$(i) Your parameter $\th=((\si_i^2)_{i=1}^d,(u_i)_{i=1}^d)$ (not $(\si_i^2)_{i=1}^N$) is not even identifiable. One reason for the non-identifiability is that the corresponding (centered multivariate Gaussian) distribution $P_\th$ will not change if you replace any $u_i$ by $-u_i$. The other reason is that the eigenvalues $\si_i^2$ of $\Si$ may be not distinct, and then even the two-sets $\{u_i,-u_i\}$ cannot be determined even if the distribution $P_\th$ is completely known. So, your parameter $\th$ cannot possibly be consistently estimated by any method. In particular, $\th$ cannot possibly be consistently estimated by the maximum likelihood (ML) method.

(ii) Your differentiation of the likelihood function with respect to $u_i$ is incorrect. It appears you are doing it without taking into account the restriction that the $u_i$'s form an orthonormal system. So, you need to use here Lagrange multipliers or something like that.

(iii) You do not need any parametrization like yours anyway, because the ML estimator of $\Si$ is well known and easy to get: In your case, the ML estimator of $\Si$ is just $$\tilde\Si:=S_N:=\frac1N\sum_{n=1}^N x_nx_n^\top.$$