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

  • 1
    $\begingroup$ It will be virtually impossible to provide a correct answer without knowing if $\epsilon_i$ and $\alpha_i$ are independent for each $i.$ $\endgroup$ Jan 27, 2015 at 17:36
  • $\begingroup$ We can assume they are independent $\endgroup$ Jan 27, 2015 at 21:35
  • 1
    $\begingroup$ 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. $\endgroup$
    – whuber
    Jan 28, 2015 at 17:25
  • 1
    $\begingroup$ 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. $\endgroup$
    – whuber
    Jan 28, 2015 at 17:38
  • $\begingroup$ @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. $\endgroup$ Jan 28, 2015 at 17:38

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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.