4
$\begingroup$

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

$\endgroup$
6
  • $\begingroup$ I'm voting to close this question as off-topic because this site doesn't deal with timing issues for algorithms. $\endgroup$ – Michael R. Chernick Mar 12 '18 at 14:03
  • 2
    $\begingroup$ @MichaelChernick There are currently 135 questions with the tag "online" which explicitely refers to algorithms. $\endgroup$ – Tom Cornebize Mar 12 '18 at 14:18
  • $\begingroup$ We don't deal with algorithm timing issues. Whether there are online issues that we might deal with I don't know. $\endgroup$ – Michael R. Chernick Mar 12 '18 at 14:45
  • 3
    $\begingroup$ @MichaelChernick: "online" does not refer to "being online". It refers to algorithms that work iteratively, as new data become available. This is not a question of algorithm timing. It's a question about updating an already-calculated MSE with newly available data, without having to touch all the original data again (or even needing the original data to be available). This is clearly on-topic (see the 135 questions tagged online). I have upvoted the question and voted to leave open. $\endgroup$ – Stephan Kolassa Mar 12 '18 at 14:48
  • $\begingroup$ I didn't say that online means being online. The question deals with the existence of an algorithm that does an iteration at a specified speed. I still see this as off topic. $\endgroup$ – Michael R. Chernick Mar 12 '18 at 16:04
0
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ Thank you for the answer, I will have a look at the gradient descent. Regarding the computation of the MSE, I think my question was ambiguous. I need the current value (as if the whole regression was done from scratch), so adding the squared error of the current observation does not work. Re-calculating everything obviously works, but is too long. I finally managed to do this, see my answer. $\endgroup$ – Tom Cornebize Mar 14 '18 at 12:29
0
$\begingroup$

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

$\endgroup$
4
  • $\begingroup$ This algorithm will only work if you have a small amount of data, otherwise the numerical issue will quickly make your result meaningless, because you are trying to get a small number as the difference of two larger number. A true useful online algorithm usually has the form MSE_n = MSE_{n-1} + some small correction term. See the online algorithm of variance calculation. $\endgroup$ – DiveIntoML Jun 28 '18 at 20:19
  • $\begingroup$ @DiveIntoML Thanks for your comment. I have implemented the algorithm of my answer and indeed, I have a bad numerical accuracy in comparison with the off-line algorithm. $\endgroup$ – Tom Cornebize Jun 29 '18 at 6:15
  • $\begingroup$ @DiveIntoML If you have an online algorithm which is numerically stable, please post an answer, I will accept it. $\endgroup$ – Tom Cornebize Jun 29 '18 at 6:16
  • $\begingroup$ I do not know any existing online algorithm for MSE, so I derived my own algorithm in my answer. I am not an expert in this field, so I will appreciate it if you can try it out and inform me if it is indeed better than what you have. $\endgroup$ – DiveIntoML Jun 30 '18 at 15:33
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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