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

  • $\begingroup$ (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. $\endgroup$
    – Glen_b
    Commented Mar 15, 2014 at 19:53
  • $\begingroup$ 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$. $\endgroup$
    – TooTone
    Commented Mar 15, 2014 at 21:53
  • $\begingroup$ @TooTone yes, this is all correct. $\endgroup$
    – greg
    Commented Mar 16, 2014 at 6:18
  • $\begingroup$ @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. $\endgroup$
    – greg
    Commented Mar 16, 2014 at 7:16

2 Answers 2


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
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
    sum(allwgt^2 * 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)
predoptgr = X %*% alloptgrwgt
predoptgrerr = X %*% alloptgrwgt - x

# Plot results
plot(prederr^2 ~ x, ylab="Squared error")
points(predoptgrerr^2 ~ x, col="green")
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.

enter image description here

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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