LinearModelFit with millions of rows

Question

I'd like to do a simple regression with millions of rows of (x,y) data. It seems to be too long for Mathematica. Is there a way to speed it up? It shouldn't be that hard to do, I am not sure what Mathematica is doing.

For example here's how the timing grows with the size of the input:

    First@Timing[LinearModelFit[RandomReal[NormalDistribution[], {#, 2}], x, x];] & 
        /@ {10, 100, 1000, 10000, 100000}

generating on my machine:

    {0., 0.015600, 0.046800, 0.670804, 26.722971}

It doesn't seem like it should be taking so much longer for longer data, and yet it does.

Answer 1

LinearModelFit does too much--it computes residuals, etc., etc. When working with large problems, just compute what you need when you need it. It all begins with the fit itself, which should be done with LeastSquares:

First[Timing[
LeastSquares[{ConstantArray[1, #], RandomReal[NormalDistribution[], #]}\[Transpose], 
 RandomReal[NormalDistribution[], #]]]] & /@ (10^Range[7])

{0., 0., 0., 0., 0.047, 0.437, 4.508}

More than half that time is just the random number generation. In brief, computing the fit only takes about one second per four million points. The timing with more independent variables will be proportionately longer (it should scale linearly with the number of points and cubically with the number of variables).

Once datasets get this large, regression (with just a single independent variable) is usually pointless: you can do just fine by randomly subsampling the data, sometimes with as few as a hundred points or so, depending on the sizes of the residuals and how precise the estimates need to be.