Recovering true data from many noisy samples with varying unknown amounts of noise

Input: $k$ vectors $x^1,\ldots,x^k \in \mathbb{R}^n$, where $x^i \sim \mathcal{N}(x,\mathbb{1} \cdot \sigma_i^2)$.

Goal: approximate the vector $x$ as well as possible.

The quality of approximation is measured in any reasonable metric such as $\mathbb{L}_1$ norm or Pearson correlation. Note that the values of $\sigma_i$ are not given as part of the input. In my input, $k \gg n$.

What is a good solution to this problem?

I'm pretty sure one can do better than simply returning the average of all the vectors $x^i$. One possibility I'm thinking of is performing a variant of canonical correlation analysis between $x^1$ and $x^2,\ldots,x^k$: the variant I'm thinking of would find the best linear combination only among those that have only positive coefficients. (I think this can be solved using semidefinite programming.)

• (1) Why would Pearson correlation be a reasonable metric for quality of approximation? (2) Please ask a definite question. So far you have only statements. Mar 15, 2014 at 19:53
• A couple of clarifications would be useful. First of all, just to check, by $x^i$ you mean the $i$th vector (often written $x_i$) not the vector $x$ to the $i$th power? Second, by $x^i\sim N(x,1\cdot\sigma_i^2)$ do you mean that $x^i_j\sim N(x_j,1\cdot\sigma_i^2)$, where $y_j$ means the $j$th component of vector $y$? I.e. the $j$th component of the $i$th vector is normally distributed with variance $\sigma_i^2$ and mean = the $j$th component of vector $x$. Mar 15, 2014 at 21:53
• @TooTone yes, this is all correct.
– greg
Mar 16, 2014 at 6:18
• @Glen_b (1) Well, if we managed to reconstruct $x$ perfectly, then the Pearson correlation is $1$. The worse we were at approximating $x$, the lower the Pearson correlation will be. (2) fixed.
– greg
Mar 16, 2014 at 7:16

Yes it is possible to do better than returning the average of all the vectors. Let $X=[\mathbf{x_1},\ldots,\mathbf{x_k}]$, where $\mathbf{x_i}$ is the $i$th vector, of length $n$ (note I've written $\mathbf{x_i}$ where you wrote $x^i$). The vector of weighted row means of $X$, i.e. $\sum_{i=1}^nw_i\mathbf{x_i} \text{ s.t. } \sum_{i=1}^nw_i=1$, is an unbiased estimator of $\mathbf{x}$. The variance of each element of this estimator vector is proportional to the sum of the variance of the $x_k$s, i.e. $\propto \sum_{i=1}^nw_i^2\sigma_i^2$. This sum can be minimized, as shown below, so the key part of the problem is estimating the $\sigma_i^2$s.

I have chosen a very simple method of estimating the $\sigma_i^2$s, namely first to estimate $\mathbf{x}$ based on equal weights $w_i=1/k$ and then subtract this from each $\mathbf{x_i}$ so that the estimate of the variance $\sigma^2_i$ is $s^2_i=\mathrm{Var}(\mathbf{x_i}-\overline{\mathbf{x_i}})$. As well as being simple, this method has the benefit that the estimated variances are trivially greater than zero! (I considered other methods to estimate $\mathrm{Var}(\mathbf{x})$ by taking differences between variances and ran into various problems. I wasn't convinced that a more complex estimation method is justified although I would love to hear otherwise.)

The variance of the weighted sum of columns $\sum_{i=1}^nw_i^2\sigma_i^2$ can be estimated with $\sum_{i=1}^nw_i^2s_i^2$. Given $w_k=1-\sum_i^{k-1}w_i$, the parameters to be minimized are $w_1,\ldots,w_{k-1}$ with a cost function $$C = \sum_i^{k-1}w_i^2s_i^2 + \left(1 - \sum_{i=1}^{n-1}w_i\right)^2s_k^2$$ The gradient is given by $$\frac{\partial C}{\partial w_i} = 2w_is_i^2 - 2\left(1 - \sum_{i=1}^{n-1}w_i\right)s^2_k$$

This can be incorporated into a simple program, as in the R code below. For simplicity and ease of demonstration, in this program $\mathbf{x} = (1,2,\ldots,n)$, and the $\sigma$s are also set to some increasing function.

# Generate random matrix X
set.seed(0)
k=100
n=10
x=1:n
sigma=(1:n) # interesting to try other functions, e.g. sqrt(1:n) or (1:n)^2
X=matrix(rnorm(k*n, rep(x,k), rep(sigma,each=k)), n,k)

# Estimate column variances
vars=apply(X-apply(X,1,mean), 2, var)

# Estimate means and errors with equal weights
allwgt= rep(1/k,k)
pred= X %*% allwgt
prederr = X %*% allwgt - x

# Estimate means and errors by minimizing weighted sum of column variances
wgtvar=function(wgt,vars)
{
allwgt=c(wgt,1-sum(wgt))
sum(allwgt^2 * vars)
}
wgtvargr=function(wgt,vars)
{
2*(wgt*vars[1:length(wgt)] - 2*(1-sum(wgt))*vars[length(wgt)+1])
}
wgt=rep(1/k,k-1) # initial weights are equal
optresgr=optim(wgt, wgtvar, vars=vars,method="BFGS",gr=wgtvargr)
alloptgrwgt=c(optresgr$par,1-sum(optresgr$par))
predoptgr = X %*% alloptgrwgt
predoptgrerr = X %*% alloptgrwgt - x

# Plot results
par(mfrow=c(2,1))
plot(prederr^2 ~ x, ylab="Squared error")
points(predoptgrerr^2 ~ x, col="green")
plot(allwgt,ylab="Weights",ylim=c(0,max(allwgt,alloptgrwgt)))
points(alloptgrwgt, col="green")


The errors and weights are compared in the plot below, with black for the equal weights and green for the optimized weights.

As you can see, the optimized weights produce a better distribution of errors, and, in general, higher weights $w_i$ correspond to lower variances $\sigma_i^2$. The reduction in errors with optimized weights is more pronounced as the $\sigma_i$s are increased, e.g. by setting sigma=(1:n)^2 in the above program.

Here is what might be a decent solution:

Given two vectors $x^i$ and $x^j$ as in the question, the difference $x_i-x_j$ is a noise vector distributed according to $\mathcal{N} \left( x,\mathbb{1} \cdot \left(\sigma_i^2 + \sigma_j^2 \right) \right)$. Note that the base vector $x$ cancelled out. Now we can get a decent estimate for $\sigma_i^2 + \sigma_j^2$ by just measuring the empirical standard deviation of $x^i-x^j$. We can thus write $\binom{k}{2}$ approximate equations in $k$ variables, get a decent solution (since if $n \gg k$ each of the equations separately is a pretty good approximation), and once approximate values for the $\sigma_i$'s are obtained, we can solve the problem in a standard way, probably by taking a weighted average of the $x^i$'s.