7
$\begingroup$

Suppose I have an online stream of data points $x_i,y_i$, where $i=1,2,\dots$. I want to compute the Pearson correlation coefficient between the vectors $\vec x$ and $\vec y$.

But here is the catch. I receive the points one by one, and computing the correlation from scratch with each new point would be too slow (at some point I cannot even store all the points at once).

So let $\rho_N$ be the Pearson correlation up to the $N$'th data point. Is there a way to efficiently update this to $\rho_{N+1}$ when I receive the next data point? (Probably I have to store some additional intermediate quantities as I receive more points).

$\endgroup$
1

4 Answers 4

9
$\begingroup$

Recall the formula for the sample Pearson correlation between two vectors $x\in\mathbb{R}^n$ and $y\in\mathbb{R}^n$ (Eq. 3 in Wikipedia):

$$ r = \frac{\sum_{i=1}^n(x_i-\overline{x})(y_i-\overline{y})}{\sqrt{\sum_{i=1}^n(x_i-\overline{x})^2}\sqrt{\sum_{i=1}^n(y_i-\overline{y})^2}} $$

We simply have to store and update the relevant quantities in this fraction:

  • $\overline{x}_{n+1}$ will contain the sample mean of $x_1, \dots, x_{n+1}$ (this is easily calculated online)
  • ditto for $\overline{y}_{n+1}$
  • $N_{n+1}=\sum_{i=1}^{n+1}(x_i-\overline{x})(y_i-\overline{y})$ will contain the numerator of $r$
  • $D_{n+1}=\sum_{i=1}^{n+1}(x_i-\overline{x})^2$ and $E_{n+1}=\sum_{i=1}^{n+1}(y_i-\overline{y})^2$ will contain the two components for the denominator.

Initialize:

$$ \overline{x}_0:=\overline{y}_0:=N_0:=D_0:=E_0:=0 $$

In updating, assume that $\overline{x}_n, \overline{y}_n, N_n, D_n, E_n$ are known, and that a new data pair $(x_{n+1}, y_{n+1})$ arrives. We update:

$$ \begin{array} \;\;\; \overline{x}_{n+1}:=& \frac{1}{n+1}(n\overline{x}_n+x_n) \\ \overline{y}_{n+1}:=& \frac{1}{n+1}(n\overline{y}_n+y_n) \\ N_{n+1}:=& N_n + (x_{n+1}-\overline{x}_{n+1})(y_{n+1}-\overline{y}_{n+1}) \\ D_{n+1}:=& D_n + (x_{n+1}-\overline{x}_{n+1})^2 \\ E_{n+1}:=& E_n + (y_{n+1}-\overline{y}_{n+1})^2. \end{array} $$

Then the correlation is

$$ r = \frac{N_{n+1}}{\sqrt{D_{n+1}}\sqrt{E_{n+1}}}. $$

$\endgroup$
3
  • $\begingroup$ The LHS of the first line of the update should be $\bar{x}_{n+1}$, shouldn't it? $\endgroup$
    – Glen_b
    Commented Oct 1, 2019 at 13:10
  • $\begingroup$ @Glen_b: you are right, it should. And it is. At least in the source code of my answer. I have no idea why MathJax decided to eat the overline. (I believe I already saw this when I first wrote the answer, but didn't see anything to do about it.) I'm looking forward to the next typo of mine you unearth! $\endgroup$ Commented Oct 1, 2019 at 13:33
  • $\begingroup$ See @Marcos answer below. $\endgroup$
    – a06e
    Commented Sep 13, 2021 at 14:18
5
$\begingroup$

A few observations on Stephan Kolassa's answer:

Computing the incremental averages as follows is more numerically robust in practice: $$ \bar{x}_{n+1} = \bar{x}_{n} + \frac{x_{n+1} - \bar{x}_{n}}{n+1} $$ $$ \bar{y}_{n+1} = \bar{y}_{n} + \frac{y_{n+1} - \bar{y}_{n}}{n+1} $$

Note how this formulation avoids (potentially) large products such as $n\bar{x}_{n}$.

$$ \\ $$

Then there's a small error in the first factors of $N_{n+1}$, $D_{n+1}$ and $E_{n+1}$: $$ N_{n+1} = N_{n} + (x_{n+1} - \bar{x}_{n})(y_{n+1} - \bar{y}_{n+1}) $$ $$ D_{n+1} = D_{n} + (x_{n+1} - \bar{x}_{n})(x_{n+1} - \bar{x}_{n+1}) $$ $$ E_{n+1} = E_{n} + (y_{n+1} - \bar{y}_{n})(y_{n+1} - \bar{y}_{n+1}) $$

Note the $\bar{x}_n$ (and $\bar{y}_{n}$) in the first factor, instead of $\bar{x}_{n+1}$ (and $\bar{y}_{n+1}$). When there's not much variance in the data series, $\bar{x}_n$ and $\bar{x}_{n+1}$ tend to be nearly identical, so it wouldn't matter much. However, for all other scenarios, the difference can be substantial.

Some useful references:
https://jonisalonen.com/2013/deriving-welfords-method-for-computing-variance/
https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Online

$\endgroup$
1
  • $\begingroup$ +1 for the corrections. But "avoid[ing] ... large products" is a chimera, and "the difference can be substantial" is a mere theoretical circumstance, because the products become large only when the $x_i$ exceed $10^{290}$ or so in absolute value. After all, surely you're not doing this calculation for $n\gg 10^{18}$ and IEEE doubles will handle products up to $10^{290+18}$ (and a little greater). $\endgroup$
    – whuber
    Commented yesterday
3
$\begingroup$

It's worth noting that the other answers here do not calculate the same thing as calculating $r$ at the end. They use the current value of $\bar{x}$ for each step of accumulation, rather than the true value for all steps. They are an approximation, and are likely to be similar, but not exactly the same.

On Wikipedia, there is a re-arrangement of the formula that allows "a convenient single-pass algorithm" which allows you to calculate the correct value at any time-step.

$$ r_{xy} = \frac{\left(\sum_i x_iy_i\right) - n\bar{x}\bar{y}} {\sqrt{\left(\sum_i x_i^2\right) - n\bar{x}^2}\sqrt{\left(\sum_i y_i^2\right) - n\bar{y}^2}} $$

So you can accumulate 5 quantities:

$$ \hat{x} := \sum^n_{i=1} x_i \;\;\;\;\;\;\; \hat{y} := \sum^n_{i=1} y_i $$ $$ \hat{a} := \sum^n_{i=1} x_i^2 \;\;\;\;\;\;\; \hat{b} := \sum^n_{i=1} y_i^2 \;\;\;\;\;\;\; \hat{c} := \sum^n_{i=1} x_iy_i $$

And plug those into the formula whenever you want $r$. For completeness here it is with my terminology for accumulated values:

$$ r = \frac{\hat{c} - \frac{\hat{x}\hat{y}}{n}}{\sqrt{\hat{a} - \frac{\hat{x}^2}{n}}\sqrt{\hat{b} - \frac{\hat{y}^2}{n}}} $$

This will give you identical results to storing all the points and calculating $r$ at the end.

Note: The $\hat{x}^2$ terms can get very large, so to use this method, you might need to handle overflow.

$\endgroup$
1
$\begingroup$

Multihunter's answer is great, and can be modified slightly to mitigate the problem that they point out of the large cumulative squared terms.

Each accumulator variable can be divided by $n$...

$$ \hat{x} = \frac{\sum_{i=1}^{n}x_i}{n}=\bar{x} \qquad \hat{y} = \frac{\sum_{i=1}^{n}y_i}{n}=\bar{y} $$ $$ \hat{a} = \frac{\sum_{i=1}^{n}x_i^2}{n} \qquad \hat{b} = \frac{\sum_{i=1}^{n}y_i^2}{n} \qquad \hat{c} = \frac{\sum_{i=1}^{n}x_iy_i}{n} $$

...which slightly modifies the on-demand correlation function to...

$$ r=\frac{\hat{c}-\hat{x}\hat{y}}{\sqrt{\hat{a}-\hat{x}^2}\sqrt{\hat{b}-\hat{y}^2}} $$

For an additional $k$ new data samples (resulting in a new total of $n+k$ samples), the accumulator variables can be updated as follows...

$$ \hat{x}_{n+k} = \hat{x}_n+\frac{\sum_{i=n+1}^{n+k}{x_i}-k\hat{x}_n}{n+k} \qquad \hat{y}_{n+k} = \hat{y}_n+\frac{\sum_{i=n+1}^{n+k}{y_i}-k\hat{x}_n}{n+k} $$ $$ \hat{a}_{n+k} = \hat{a}_n+\frac{\sum_{i=n+1}^{n+k}{x_i^2}-k\hat{a}_n}{n+k} \qquad \hat{b}_{n+k} = \hat{b}_n+\frac{\sum_{i=n+1}^{n+k}{y_i^2}-k\hat{b}_n}{n+k} \qquad \hat{c}_{n+k} = \hat{c}_n+\frac{\sum_{i=n+1}^{n+k}{x_iy_i}-k\hat{c}_n}{n+k} $$

Those update functions can be simplified somewhat by multiplying the leading term by $\frac{n+k}{n+k}$ and combining everything into a single fraction, but that reintroduces the issue of potentially having to deal with very large intermediate terms.

New contributor
cg49me is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
$\endgroup$

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.