# Efficient online linear regression

I'm analysing some data where I would like to perform ordinary linear regression, however this is not possible as I am dealing with an on-line setting with a continuous stream of input data (which will quickly get too large for memory) and need to update parameter estimates while this is being consumed. i.e. I cannot just load it all into memory and perform linear regression on the entire data set.

I'm assuming a simple linear multivariate regression model, i.e.

y = Ax + b + e

What's the best algorithm for creating a continuously updating estimate of the linear regression parameters A and b?

Ideally:

• I'd like an algorithm that is most O(N*M) space and time complexity per update, where N is the dimensionality of the independent variable (x) and M is the dimensionality of the dependent variable (y).
• I'd like to be able to specify some parameter to determine how much the parameters are updated by each new sample, e.g. 0.000001 would mean that the next sample would provide one millionth of the parameter estimate. This would give some kind of exponential decay for the effect of samples in the distant past.
-
Look up (1) Flexible linear regression, (2) Kalman Filters. –  Jase Aug 10 at 7:44

Maindonald describes a sequential method based on Givens rotations. (A Givens rotation is an orthogonal transformation of two vectors that zeros out a given entry in one of the vectors.) At the previous step you have decomposed the design matrix $\mathbf{X}$ into a triangular matrix $\mathbf{T}$ via an orthogonal transformation $\mathbf{Q}$ so that $\mathbf{Q}\mathbf{X} = (\mathbf{T}, \mathbf{0})'$. (It's fast and easy to get the regression results from a triangular matrix.) Upon adjoining a new row $v$ below $\mathbf{X}$, you effectively extend $(\mathbf{T}, \mathbf{0})'$ by a nonzero row, too, say $t$. The task is to zero out this row while keeping the entries in the position of $\mathbf{T}$ diagonal. A sequence of Givens rotations does this: the rotation with the first row of $\mathbf{T}$ zeros the first element of $t$; then the rotation with the second row of $\mathbf{T}$ zeros the second element, and so on. The effect is to premultiply $\mathbf{Q}$ by a series of rotations, which does not change its orthogonality.

When the design matrix has $p+1$ columns (which is the case when regressing on $p$ variables plus a constant), the number of rotations needed does not exceed $p+1$ and each rotation changes two $p+1$-vectors. The storage needed for $\mathbf{T}$ is $O((p+1)^2)$. Thus this algorithm has a computational cost of $O((p+1)^2)$ in both time and space.

A similar approach lets you determine the effect on regression of deleting a row. Maindonald gives formulas; so do Belsley, Kuh, & Welsh. Thus, if you are looking for a moving window for regression, you can retain data for the window within a circular buffer, adjoining the new datum and dropping the old one with each update. This doubles the update time and requires additional $O(k (p+1))$ storage for a window of width $k$. It appears that $1/k$ would be the analog of the influence parameter.

For exponential decay, I think (speculatively) that you could adapt this approach to weighted least squares, giving each new value a weight greater than 1. There shouldn't be any need to maintain a buffer of previous values or delete any old data.

### References

J. H. Maindonald, Statistical Computation. J. Wiley & Sons, 1984. Chapter 4.

D. A. Belsley, E. Kuh, R. E. Welsch, Regression Diagnostics: Identifying Influential Data and Sources of Collinearity. J. Wiley & Sons, 1980.

-
In that case see also the extensions by Alan Miller jstor.org/stable/2347583. An archive of his Fortran software site is now at jblevins.org/mirror/amiller –  onestop Feb 5 '11 at 22:44
An explicit algorithm appears at the bottom of p. 4 of saba.kntu.ac.ir/eecd/people/aliyari/NN%20%20files/rls.pdf . This can be found by Googling "recursive least squares." It does not look like an improvement on the Gentleman/Maindonald approach, but at least it is clearly and explicitly described. –  whuber Feb 5 '11 at 23:20
The last link looks like the method I was going to suggest. The matrix identity they use is known in other places as the Sherman--Morrison--Woodbury identity. It is also quite numerically efficient to implement, but may not be as stable as a Givens rotation. –  cardinal Feb 6 '11 at 0:28
@cardinal Thank you for the reference and for the advice. I was mainly concerned about stability. In a long sequence of updates, there may be circumstances where the data are nearly collinear and there could also be a lengthy accumulation of numerical errors, so stability seems to be of paramount concern. That is the justification for the Givens rotations (instead of Householder reflections, which take less computation). –  whuber Feb 6 '11 at 19:16
@suncoolsu Hmm... Maindonald's book was newly published when I started using it :-). –  whuber Feb 8 '11 at 6:02

I think recasting your linear regression model into a state-space model will give you what you are after. If you use R, you may want to use package dlm and have a look at the companion book by Petris et al.

-
maybe I'm confused but this appears to refer to a time series model? my model is actually simpler in that the samples are not a time series (effectively they are independent (x->y) samples, they are just accumulated in large volumes over time) –  mikera Feb 5 '11 at 19:26
Yes, in the general case this is used for time series with non independent observations; but you can always assume incorrelation between successive observations, which gives the special case of interest to you. –  F. Tusell Feb 5 '11 at 19:46

You can always just perform gradient descent on the sum of squares cost $E$ wrt the parameters of your model $W$. Just take the gradient of it but don't go for the closed form solution but only for the search direction instead.

Let $E(i; W)$ be the cost of the i'th training sample given the parameters $W$. Your update for the j'th parameter is then

$$W_{j} \leftarrow W_j + \alpha \frac{\partial{E(i; W)}}{\partial{W_j}}$$

where $\alpha$ is a step rate, which you should pick via cross validation or good measure.

This is very efficient and the way neural networks are typically trained. You can process even lots of samples in parallel (say, a 100 or so) efficiently.

Of course more sophisticated optimization algorithms (momentum, conjugate gradient, ...) can be applied.

-
Seems very similar to this paper eprints.pascal-network.org/archive/00002147/01/… . It's been implemented in an open source project called jubatus. –  saccharine Mar 7 at 0:24

The problem is more easily solved when you rewrite things a little bit:

Y = y

X = [x, 1 ]

then

Y = A*X

A one time-solution is found by calculating

V = X' * X

and

C = X' * Y

note the V should have size N-by-N and C a size of N-by-M. The parameters you're looking for are then given by:

A = inv(V) * C

Since both V and C are calculated by summing over your data, you can calculate A at every new sample. This has a time complexity of O(N^3), however.

Since V is square and semi-definite positive, a LU decomposition does exists, that makes inverting V numerically more stable. There are algorithms to perform rank-1 updates to the inverse of a matrix. Find those and you'll have the efficient implementation you're looking for.

The rank-1 update algorithms can be found in "Matrix computations" by Golub and van Loan. It's tough material, but it does have a comprehensive overview of such algorithms.

Note: The method above gives a least-square estimate at each step. You can easily adding weights to the updates to X and Y. When the values of X and Y grow too large, they can be scaled by a single scalar, without affecting the result.

-

The standard least-square fit gives regression coefficients

$\beta = ( X^T X )^{-1} X^T Y$

where X is a matrix of M values for each of N data points, and is NXM in size. Y is a NX1 matrix of outputs. $\beta$ of course is a MX1 matrix of coefficients. (If you want an intercept just make one set of x's equal always to 1.)

To make this online presumably you just need to keep track of $X^T X$ and $X^T Y$, so one MXM matrix and one MX1 matrix. Every time you get a new data point you update those $M^2+M$ elements, and then calculate $\beta$ again, which costs you an MXM matrix inversion and the multiplication of the MXM matrix and the MX1 matrix.

For example, if M=1, then the one coefficient is

$\beta = \frac{\sum_{i=1}^N{x_i y_i}}{\sum_{i=1}^N{x_i^2}}$

so every time you get a new data point you update both sums and calculate the ratio and you get the updated coefficient.

If you want to damp out the earlier estimates geometrically I suppose you could weight $X^T X$ and $X^T Y$ by $(1-\lambda)$ each time before adding the new term, where $\lambda$ is some small number.

-
It's nice to see this simple case explained. Did you notice, though, that the question specifically asks about multivariate regression? It's not so easy to update the denominator of $\beta$ in that case! –  whuber Apr 19 at 21:07
I think my answer works still: ie you need to keep track of the MxM matrix $X^T X$ and the Mx1 matrix $X^T Y$. Each element of those matrices is a sum like in the M=1 example. Or am I missing something? –  Mark Higgins Apr 24 at 13:17
Yes: in addition to computing a matrix product and applying a matrix to a vector, you now need to invert $X'X$ at each step. That's expensive. The whole point to online algorithms is to replace wholesale expensive steps by cheaper updating procedures. –  whuber Apr 24 at 13:24