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, it takes about one second per four million points. 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.

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.