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

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$ 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 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 $d$ orthonormal vectors). $$ \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 !

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 !

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 \log(\sigma_i^2) + \frac1N \sum_{n=1}^N\tilde{x}^\ast \Sigma^{-1} \tilde{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(\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 !

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 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$ 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 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 $d$ orthonormal vectors). $$ \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 !

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 \log(\sigma_i^2) + \frac1N \sum_{n=1}^N\tilde{x}^\ast \Sigma^{-1} \tilde{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(\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 !

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 \log(\sigma_i^2) + \frac1N \sum_{n=1}^N\tilde{x}^\ast \Sigma^{-1} \tilde{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(\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 !