Efficient online linear regression

Question

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.

$$\mathbf y = \mathbf A\mathbf x + \mathbf b + \mathbf e$$

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

Ideally:

• I'd like an algorithm that is most $\mathcal O(N\cdot M)$ space and time complexity per update, where $N$ is the dimensionality of the independent variable ($\mathbf x$) and $M$ is the dimensionality of the dependent variable ($\mathbf 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.

Answer 1

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.

Answer 2

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.

Answer 3

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.

Answer 4

Surprised no one else touched on this so far. Linear regression has a quadratic objective function. So, a newton Raphson step from any starting point leads you straight to the optima. Now, let's say you already did your linear regression. The objective function is:

$$ L(\beta) = (y - X \beta)^t (y - X \beta) $$ The gradient becomes $$ \nabla L (\beta) = -2 X^t (y - X \beta)$$ And the hessian: $$ \nabla^2 L (\beta) = X^t X$$

Now, you got some past data and did a linear regression and are sitting with your parameters ($\beta$). The gradient at this point is zero by definition. The hessian is as given above. A new data point ($x_{new}, y_{new}$) arrives. You just calculate the gradient for the new point via:

$$\nabla L_{new}(\beta) = -2 x_{new} (y_{new}-x_{new}^T \beta)$$ and that will become your overall gradient (since the gradient from the existing data was zero). The hessian for the new data point is:

$$\nabla^2 L_{new} = x_{new}x_{new}^T $$.

Add this to the old hessian given above. Then, just take a Newton Raphson step.

$$\beta_{new} = \beta_{old} + (\nabla^2L)^{-1} \nabla L_{new}$$

And you're done.

Answer 5

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.

Answer 6

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.