Online algorithm for the mean square error

Question

Given a dataset $\{(x_1, y_1), (x_2, y_2), \dots\}$, we can compute incrementally (or "online") the linear regression for those points.

In other words, given a new point $(x_i, y_i)$, we can recompute in time $O(1)$ the coefficients $\alpha$ and $\beta$ to get the best fit for the equation $y = \alpha x + \beta$. To do so, simply use the online algorithms to compute the variance and the covariance.

However, I don't know how to compute the mean square error incrementally in time $O(1)$ per new point:

$\text{MSE} = \frac{1}{N}\sum_i (y_i - \alpha x_i - \beta)^2$

Is there any known algorithm for that?

As a bonus, it would also be great if it was possible to remove any point from the dataset and recompute the MSE in constant time (i.e., doing the reverse operation).

Answer 1

Recall that simple linear regression can be described as

$$Y = m(X) + \frac{Cov(X,X)}{Var(X)}(X-m(X))$$

so to calculate it, you need means $m(X), m(Y)$, variance of $X$ and covariance between them. As you can learn from here mean and variance can be calculated using online algorithm in single pass. There is also online algorithm to calculate covariance as described in Wikipedia. So the only thing you need to do, is to substitute the regular estimators with their online counterparts. See also the Are there algorithms for computing "running" linear or logistic regression parameters? thread.

Taking this apart, what people would usually do in such cases is use the stochastic gradient descent as it is easy to implement (already implemented in many packages, e.g. Vowpal Wabbit, TensorFlow, Scikit-learn) and to adapt for more complicated cases. Also notice that while it would be possible to estimate it in single pass, then using more then one pass would give you more accurate results.

As about calculating mean squared error, this is very simple, you need only to store the sum of squared errors and at each iteration increase it by squared error for the current observation, and then just divide it by the number of the iterations made. Notice that MSE would change the same as your model changes, so how to calculate it is also a question what exactly does it have to measure. If you want to catch how good is your model currently but accounting for the past performance (previous iterations), then the above solution seems to work for this. If you want to measure only the current performance, then you'd need to re-calculate it at each iteration on the whole dataset. On another hand, if you want to measure the out-of-data performance, then probably the best you can do is to calculate it at each step for the hold-out sample, or for the $t+1,\dots,t+k$ observations ahead as described by Hyndman and Athanasopoulos.

Answer 2

Recall that the MSE of the linear regression would be:

$\frac{1}{N} \sum_i (y_i - (\alpha x_i + \beta))^2$

We sum over the $N$ pairs $(x_i, y_i)$ that we had so far.

To avoid having to recompute the whole sum each time a new pair $(x_i, y_i)$ is added, a naive idea could be to store the sum of squares somewhere:

sum = 0
N = 0
while True:
    x, y = get_new_pair()
    a, b = update_regression(x, y)
    N += 1
    sum += (y - (a*x + b))**2
    MSE = sum/N

This solution does not work to compute the current MSE since the values of $\alpha$ and $\beta$ will change for each new point (thus, the squared errors we added in previous iterations are not correct anymore).

A better solution can be obtained by developing the square inside the sum:

$\begin{align}\text{MSE} &= \frac{1}{N} \sum_i (y_i - (\alpha x_i + \beta))^2\\ &= \frac{1}{N} \left( \sum_i y_i^2 - 2\sum_i y_i(\alpha x_i + \beta) + \sum_i (\alpha x_i + \beta)^2 \right)\\ &= \frac{1}{N} \left( \sum_i y_i^2 -2\left(\alpha \sum_i x_i y_i + \beta\sum y_i\right) + \left(\alpha^2 \sum x_i^2 + 2\alpha\beta\sum x_i + N\beta^2\right) \right)\end{align}$

Note that, in the last equation, $\alpha$ and $\beta$ do not appear in any sum. This means that we can compute the updated value of the MSE without having to compute any sum from 0. We just have to update the values of each sub-sums (e.g. $\sum_i y_i^2$ and $\sum_i x_iy_i$).

Answer 3

The answer of @Tom Cornebize is mathematically correct, but leads to numerical problem. The numerical problem originates from the fact that we are calculating the difference of two larger number to get a small number, which is limited by the numerical precision. An extreme example is we want to calculate $A - B$ where $A$ and $B$ are two floats with 5 digits and their exact value is $A = 5 \times 10^{100}$ and $B = 5 \times 10^{100} - 1$. The representation of float value will totally ignore the 1 at the end and gives a zero.

To be more precise, each item in the expression given by @Tom Cornebize has order of magnitude $\sim N^2$, while their sum/difference is supposed to $\sim N$ for $N \times MSE$, and the large difference in their order of magnitudes will give incorrect result.

Instead, we want to calculate $N \times MSE$ by the following way: we already have an analytical expression for $N \times MSE$, and we would like to calculate the analytical expression where those large terms cancel each other, so we are only left with those terms up to order of magnitude $\sim N$ or even smaller.

We already have the online algorithm for variance $var_n(x)$ and $var_n(y)$ and covariance $cov_n(x, y)$ from Wikipedia , and each of them has order of magnitude $O(1)$, so we can calculate the MSE by the expression $$MSE_n = var_n(y) - \frac{cov_n(x, y)^2}{var_n(x)}$$ where $n$ is the index of sample the algorithm is up to.

An even better way is to cancel out the $O(1)$ items from the above expression, and calculate the incremental value as $$MSE_n = MSE_{n-1} + \Delta MSE$$ We make use of the result from online algorithm of variance and covariance in the following:

$$n MSE_n = (n-1) MSE_{n-1} + \Delta$$ $$n var_n(y) = (n-1) var_{n-1} + \Delta_y$$ $$n var_n(x) = (n-1) var_{n-1} + \Delta_x$$ $$n cov_n(x, y) = (n-1) cov_{n-1}(x, y) + \Delta_{x,y}$$

where $\Delta_y$, $\Delta_x$ and $\Delta_{x, y}$ is defined by the algorithm in the link above, then we get the ugly expression for $\Delta$ $$\Delta = \Delta_y + \frac{cov^2_{n-1}(x,y)}{var_n(x)var_{n-1}(x)}\Delta_x - \frac{cov_{n-1}(x,y) + cov_{n}(x,y)}{var_n(x)}\Delta_{x,y}$$ Since the term $n MSE_n$ and $(n-1) MSE_{n-1}$ are both $O(n)$ and $\Delta$ is $O(1)$, this should give better numerical accuracy.