Variance estimation in a one-factor linear model

Question

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

Answer 1

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:

\begin{align}\sigma_\epsilon^2 &= \frac{1}{p-1}\sum_{i=2}^p \lambda_i \\ \u &= \mathbf w_1 (\lambda_1 - \sigma_\epsilon^2)^{1/2},\end{align}

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.