Online update of Pearson coefficient

Question

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

Answer 1

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}}}. $$

Answer 2

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

Answer 3

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.

Answer 4

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.