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 !