# Principal component analysis with known uncertainty varying across both samples and variables

I want to do dimensionality reduction on a data set $X_{ij}$. In this case, $i$ indexes samples and $j$ indexes a large number of variables (densities at different locations in space). The units of all elements are the same, but the uncertainties (which I can estimate) vary quite a lot, mostly across samples. Each $X_{ij}$ is Gaussian, with standard deviation $\sigma_{ij}$ varying across both samples and variables.

My problem is that doing standard PCA without taking into account the uncertainties will tend to produce results which are dominated by variance due to measurement error. What I would like to do is estimate some similar quantities to the principal components, but which will be more robust to outliers with large $\sigma_{ij}$. The variables $X_{ij}$ are drawn from $\mathcal{N}(\mu_{ij},\sigma_{ij})$. I have one such sample for each $i$ and $j$, and I would like to calculate a transformation involving $X_{ij}$ and $\sigma_{ij}$ which will estimate principle components of the matrix $\mu_{ij}$, for example.

It is difficult for me to make progress analytically because I cannot calculate how uncertainty in $X$ propagates to its principle components. I can generate surrogate data to confirm that uncertainty dominates the principle component analysis, and I can also do this to validate any alternative techniques. My main goal is to look for independent modes across the variables, so my concerns are with dimensionality reduction for presentation, description and simplification, rather than physical interpretation.

## 1 Answer

An astrophysicist colleague of mine answered this problem in this article. Take a look at sections 2 and 3. It describes an iterative method on how to compute the principal components by criss-crossing between vectors and values, when the error bars as you describe are known. It is very simple to implement and cites some interesting related work.

• Thanks for the link, the technique presented in the article may be useful. I am giving it a try today. One thing I do not understand is their choice of normalizing by $\sigma_{ij}$. Is this because the signal to noise ratio is of fundamental interest? In that case, Why not just do PCA on $X_{ij}/\sigma_{ij}$? – OskarM May 20 '14 at 15:21
• My interpretation is that the objective is to take into account the variance variation from one image to the other. So it is a matter of S/N, but in relative terms, across images. – pedrofigueira May 20 '14 at 15:49
• I would like to add that the paper cited in the paper you suggested is very nice: it generalizes PCA to fitting a data matrix to the product of two matrices of lower rank than the original data, with a weight per data point on the fit. Lower Rank Approximation of Matrices by Least Squares With Any Choice of Weights, Gabriel and Zamir, 1979 – OskarM May 20 '14 at 22:05
• @user45826: Note that you can avoid iterative solution (with multiple local minima etc.) if you are happy to restrict your weights $1/\sigma_{ij}^2$ to vary only with samples. Then you can simply multiply each row of your data matrix $X$ with the square root of the corresponding weight (i.e. with $1/\sigma_i$), perform standard PCA of that matrix and you are done. If the weights vary with samples and variables then this trick does not work, and you have to resort to an iterative minimization. – amoeba May 21 '14 at 10:00
• This is one of the first things I tried. This is like the transpose of normalizing by per-variable deviance, which is often used to compare quantities with different units. I suspect some groups have used this transposed method, although I am not sure where. Have you seen any references? – OskarM May 22 '14 at 0:01