How to combine covariance matrices?

I have a data set of points in three dimensions.

I'm calculating the barycenter (mean) and $3\times3$ covariance matrix from this data set. I store the average, the $3\times3$ matrix (where really only 6 elements are unique) and the number of points that went into this calculation. Let's call this a "measure."

Now, I do the same for a second set of points. After doing this, I want to combine the two "measures" into one measure that's as close as possible to what I would get if I were to calculate one big "measure" based on the union of all the data points.

What's the right (or least wrong) way of doing this? I know how to accurately update the mean by calculating a weighted average based on the number of points in each of the measures, but it feels wrong to do the same thing with the covariance matrix.

$\newcommand{\cov}{\operatorname{cov}}$

$$\cov(X,Y) = \mathbb E(\cov(X,Y)\mid Z) + \cov(\mathbb E(X\mid Z), \mathbb E(Y\mid Z)).$$

Let $Z$ be $0$ or $1$ according as the data point in question belongs to one or the other data set. Then the two covariance matrices you've already got, which you want to combine, are the two conditional covariance matrices. The probabilities of belonging to one data set or the other are proportional to their sizes (if you're giving every point equal weight).

So the first term above, the conditional expected value, is the appropriate weighted average of the two matrices.

PS: OK, now I'm realizing that this assumes that for each data point you can say what it's value is when it's in one of the two data sets and what it's value is when it's in the other, so that they're paired. Maybe your data are not paired.

Still later edit: Forget pairing. This would have been better if answered more patiently. You've got a cloud of points $\{(x_i,y_i) : i\in I\}$, and another cloud $\{(x_i,y_i)$, $i\in J\}$, $I\cap J=\varnothing$. You pick one member of $I\cup J$ at random, all equally probable. That gives you an $x$ and a $y$. So there's a covariance. And it obeys the rule above. However: Letting $n=I\cup J$, don't divide by $n-1$ in finding the covariances, but by $n$. Statistics textbooks often tell students that a sample variance is $$\frac{1}{n-1}\sum_{i=1}^n (x_i - \bar x)^2,\text{ where }\bar x =\frac{x_1+\cdots+x_n}{n}.$$ And maybe they don't emphasize enough that that is to be used ONLY when one is estimating a population variance by using a small sample or an i.i.d. sample to estimate the variance of a large population from which the sample was taken. The principal reason why standard deviation is used as a measure of dispersion rather than seemingly simpler things like mean absolute deviation is that the variance of a sum of independent random variables is the sum of their variances, BUT that doesn't work if variances are defined using $n-1$ rather than $n$.

• At the point where I have the two measures, I have no knowledge of the actual data values. However, the sources of the different measures are similar -- specifically, the points are RGB values of pixels around a particular neighborhood in a digital image. If I understand /pairing/ correctly, if the two images are taken with the same camera in the same conditions, they could be considered to be paired. – Jon Watte May 1 '13 at 1:43
• The question is whether each point in the first set corresponds to one point in the second set. – Michael Hardy May 1 '13 at 22:00
• each of the sets are subsamplings of the same very big data set, but each sample set contains independent data point samples, if that's what you're asking. – Jon Watte May 2 '13 at 6:08
• OK, look at the revised comments above. – Michael Hardy May 3 '13 at 1:58
• Thanks for the later edit! – Jon Watte May 4 '13 at 6:46