I was given a dataset (a mat file) of $100\: 000$ observations, each with $50$ dimensions (coordinates). Denote matrix $X$ a $50\times 100\:000$ matrix in which each column was generated according to:
$$\mathbf x_i = a_i \mathbf u + \boldsymbol \epsilon_i,$$
where $\mathbf u \in \mathbb R^{50}$ is some fixed 50-dimensional vector, $\mathbf x_i$ is a $50$-dimensional vector with $i$ indexing observations ($i=1...100\:000$). Each $a_i$ is a scalar, $a_i \sim \mathcal N(0, \sigma_a)$ i.i.d Gaussian with zero mean and unknown finite variance. Noise term $\boldsymbol \epsilon_i$ is a 50-dimensional vector, with each coordinate $\epsilon_j \sim \mathcal N(0, \sigma_\epsilon)$ i.i.d. Gaussian with zero mean and unknown finite variance. We can assume that $a_i$ and $\boldsymbol \epsilon_i$ are independent.
I need to estimate $\sigma_a$ (a scalar), $\sigma_\epsilon$ (a scalar) and $\mathbf u$ (a vector).
My attempts so far
Denote $\mathbf z_i=a_i \mathbf u$. Since $E[a_i^2]=\sigma_a^2$, the $p\times p$ covariance matrix of $\mathbf z$ is given by:
$$\Sigma=E[z_iz_i^T]=\mathbf u\mathbf u^T\sigma_a^2.$$
Then, the covariance matrix of $\mathbf x$ is given by $$S=\mathbf u\mathbf u^T\sigma_a^2+\sigma_\epsilon^2\mathbf I_{p\times p}.$$
Since $n\gg p$ in my case, the sample covariance matrix $S_n$ converges to $S$. Hence, $\sigma_a^2$ (the variance of $a_i$) and the noise variance $\sigma_\epsilon^2$ can be estimated from the sample covariance matrix.
So, we computed $S$ and we have an equation of $S$ with $\mathbf u$, $\sigma_a$, and $\sigma_\epsilon$. But I couldn't find a way to recover $\mathbf u$ and I feel like I'm missing something.
I applied SVD and found the spectral decomposition of $X$, but couldn't figure out how the eigenvalues/eigenvectors can help me to estimate the variances $\sigma_a$ and $\sigma_\epsilon$ or to find $\mathbf u$.