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 bases.

Setup

We consider a dataset $\{x_i \in \mathcal{H}\}_{i=1}^N$ with $(\mathcal{H},\langle\cdot,\cdot\rangle)$ a Hilbert space of finite dimension $d$. For simplicity however, we will consider the centered data $\{\tilde{x}_i\}_{i=1}^N$ with $\tilde{x}_i = x_i - \frac1N\sum_{j=1}^N x_j$ for all $i = 1,\ldots,N$. The adjoint of any $\tilde{x}$ is defined as $\tilde{x}^\ast = \langle \tilde{x},\cdot\rangle$. If we define a (centered) Gaussian as $$ p(x) = \frac1Z\exp\left(-\frac12\tilde{x}^\ast \Sigma^{-1} \tilde{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$ and orthonormal base 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 deduct that the normalization term is equal to $$ Z = (2\pi)^{d/2}\prod_{i=1}^d | \sigma_i |. $$

Maximum Likelihood Estimator (MLE)

We can now compute the MLE based on our dataset. The likelihood function is given by $$ \begin{eqnarray} \mathcal{L} &=& -\frac N2\left\{d\log(2\pi) + \sum_{i=1}^d |\sigma_i| + \frac1N \sum_{n=1}^N\tilde{x}^\ast \Sigma^{-1} \tilde{x} \right\}, \\ &=& -\frac N2\left\{d\log(2\pi) + \sum_{i=1}^d |\sigma_i| + \frac1N \sum_{n=1}^N\sum_{i=1}^d\sigma^{-2}_i(\tilde{x}^\ast u_i)^2\right\}. \end{eqnarray} $$

Issue

When computing the stationnary points 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 \tilde{x}\tilde{x}^\ast \right) = 0, $$ which is absurd as it should be equal to $1$ in order to obtain the eigendecomposition of $\frac1N\sum_{n=1}^N \tilde{x}\tilde{x}^\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 !