Variance estimation in a one-factor linear model

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

• It will be virtually impossible to provide a correct answer without knowing if $\epsilon_i$ and $\alpha_i$ are independent for each $i.$ – StatsStudent Jan 27 '15 at 17:36
• We can assume they are independent – Serendipity Jan 27 '15 at 21:35
• It would be best to tell us what your data really are and what you are trying to learn about them: it looks like essential information has been lost or corrupted in your effort to state your question abstractly. – whuber Jan 28 '15 at 17:25
• Much better, thank you! Notice that your model does not change if you simultaneously change $u$ to $\lambda u$ and $\sigma_a$ to $\sigma_a/|\lambda|$ for some nonzero number $\lambda$, because that would not alter the distribution of $x_i$ at all. Therefore $\sigma_a$ is not identifiable and only the direction of $u$ is identifiable, but not its magnitude. If you were to constrain $u$--insist it's a unit vector, for instance--then almost everything would be identifiable; the only remaining ambiguity would be that both $u$ and $-u$ would work as solutions. – whuber Jan 28 '15 at 17:38
• @whuber thank you for the help and good will! I re-added the post with all I know about the given data and problem definition. – Serendipity Jan 28 '15 at 17:38

First, @whuber is right: your $\sigma_a^2$ and $\newcommand{\u}{\mathbf u}\u$ are underdetermined: as you computed yourself, they enter the covariance matrix as a $\u\u^\top\sigma^2_a$ term, and so if you e.g. multiply $\u$ by two and divide $\sigma_a$ by two, this term will not change. Therefore you should either fix the length of $\u$ or the value of $\sigma_a$. Let's fix $$\sigma_a=1.$$
Second, let's rewrite your formula in a somewhat more standard notation. We have: $$a \sim \mathcal N(0,1) \\\mathbf x =a\u+\boldsymbol \epsilon \sim \mathcal N (a\u, \sigma^2_\epsilon \mathbf I).$$ It is a probabilistic model that is almost factor analysis with a single factor $a$. It would be exactly factor analysis, if you had an arbitrary diagonal covariance of $\boldsymbol \epsilon$, i.e. if different components of $\boldsymbol \epsilon$ were allowed to have different variances. Factor analysis model with one factor can be fit e.g. via expectation-maximization (and there are standard routines for that). Vector $\u$ is called a vector of loadings.
If your noise covariance matrix is not only diagonal, but also isotropic, i.e. $\sigma^2_\epsilon \mathbf I$, then factor analysis reduces to probabilistic PCA (pPCA). It can also be fit via expectation-maximization, however it turns out that there is an analytical maximum likelihood solution (unlike for FA). The solution is as follows:
where $\mathbf w_1$ and $\lambda_1$ are the first eigenvector of the observed covariance matrix $\mathbf S$ and its eigenvalue, and $\lambda_i$ -- other eigenvalues.
I think this is what you need. Note that $\u$ is pointing in the direction of the first principal axis and is almost equal to the PCA loading vector $\lambda_1 \mathbf w_1$, but compared to it is slightly scaled down in length to "make room" for the noise covariance.