How can I deal with numerical errors in a large-scale linear regression?

Question

I am currently conducting a linear regression on a large-scale data set which has many sparse features ($\simeq 10^5$) and many observations ($\simeq 10^6$) by using scikit-learn. (Most of the features are categorical variables, so my data set is very sparse and of size linear in the number of the observations.) While the solver outputs a solution in some minutes, the solution appears to contain a significant error (presumably due to numerical errors).

Specifically, I encountered two unreasonable behaviors:

1. When I added some new features, the R-squared score of the train data set decreased significantly.
2. When I normalized a feature, the R-squared score increased significantly.

I am not sure, but I suspect that these unreasonable behaviors are due to numerical errors.

My question is:

1. Are there more numerically stable linear regression libraries (compared to sklearn.linear_model.LinearRegression, or equivalently, scipy.linalg.lstsq)?
2. Is there any other way to fix these two problems?
3. Why do these problems occur?

Answer 1

I would like to know how and what features did you add that made your R squared decrease. To answer your questions:

1. If the model is sparse, I would suggest using lasso regression which helps get rid of most features, by L1 penalization. You would need to find the optimal parameter of regularizer though, I would suggest cross validation.

2. (and 3) It could be numerical errors, I do not know, but they happen mostly when the magnitude of your features are too big or too small, or they are, relatively, not on the same scale. I would suggest you used either normalize=True in your Linear_regression, or just call scale:

   from sklearn.preprocessing import scale
   scaled_X = scale(X, axis=0)
   

to make the variance of all your features go to 1. This is a good practice because partially you will remedy the relative difference in your feature values, and partially, you will help your gradient descent converge faster. The latter might be a reason you are getting a better R squared when normalizing features. Imagine a bowl, when features are not normalized, the bowl will have an ellipsoid shape. While when the features are normalized the bowl will have a circular shape. When the shape is circular, given a reasonable step size, your gradient descent or stochastic gradient descent will converge faster, because they would jump back and forth less often as is the case when your features are not normalized.