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.)
 A: 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.
A: 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.
