By reading your code i realized that you were using the sklearn lib. But according to what I found in the develop branch, you are not using the intersect parameter of the LinearRegression method.

It means that using LInearRegression all-together is useless as the optimisation problem have a closed solution given by

theta = (H^T*H)^-1*H^T*Y

I am proposing a change to the lib using only numpy to compute the coefficients. According to the tests I've performed on my local computer, I get the same results twice as fast (44s with sklearn, 23 without).

It will also reduce the number of dependencies which is alway important.