LinearRegression predictions differ depending on the type of input

Viewed 170

I am working on a model which predicts car prices by their attributes. I have noticed that the predictions of LinearRegression model differ depending on the type of input (numpy.ndarray, scipy.sparse.csr.csr_matrix).

My data consists of a few numerical and categorical attributes, there are no NaNs.

This is a simple data preparation code (it is common for every case I describe later):

from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.linear_model import LinearRegression

# Splitting to test and train
X = data_orig.drop("price", axis=1)
y = data_orig[["price"]]
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Numerical attributes pipeline
num_pipeline = Pipeline([ ("scaler", StandardScaler()) ])

# Categorical attributes pipeline
cat_pipeline = Pipeline([ ("encoder", OneHotEncoder(handle_unknown="ignore")) ])

# Complete pipeline
full_pipeline = ColumnTransformer([
    ("cat", cat_pipeline, ["model", "transmission", "fuelType"]),
    ("num", num_pipeline, ["year", "mileage", "tax", "mpg", "engineSize"]),
])

Let's build a LinearRegression model (X_train and X_test will be instances of scipy.sparse.csr.csr_matrix):

...   
X_train = full_pipeline.fit_transform(X_train)
X_test = full_pipeline.transform(X_test)

from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression().fit(X_train, y_train)
pred = lin_reg.predict(X_test)

r2_score(y_test, pred) # 0.896044623680753 OK

If I convert X_test and X_train to the numpy.ndarray, the predictions of the model are completely incorrect:

...   
X_train = full_pipeline.fit_transform(X_train).toarray() # Here
X_test = full_pipeline.transform(X_test).toarray() # And here

from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression().fit(X_train, y_train)
pred = lin_reg.predict(X_test)
    
r2_score(y_test, pred) # -7.919935999010152e+19 Something is wrong

I also tested DecisionTreeRegressor, RandomForestRegressor and SVR but the problem occurs only with LinearRegression.

2 Answers

In the source code, you can see if your input is a sparse matrix, it does some centering then calls on the sparse version of linear least square. If the array is dense, it calls on the numpy version of linear least square.

However the larger issue with this example is that before you perform onehot encoding, you should check whether any of the categorical values have only 1 entry:

data_orig.select_dtypes(['object']).apply(lambda x:min(pd.Series.value_counts(x)))

model              1
transmission    2708
fuelType          28

And if we check model:

data_orig['model'].value_counts().tail()

 SQ7    8
 S8     4
 S5     3
 RS7    1
 A2     1

So if RS7 and A2 are in your test but not in your training, then this coefficient will be totally rubbish because all it's values are zeros. If we try using another seed to split the data, you can see both fits are similar:

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=22)

A function to fit using sparse/dense data differently:

import matplotlib.pyplot as plt

def fit(X_tr,X_ts,y_tr,y_ts,is_sparse=True):
    data_train = full_pipeline.fit_transform(X_tr)
    data_test = full_pipeline.transform(X_ts)
    
    if not is_sparse:
        data_train = data_train.toarray()
        data_test = data_test.toarray()
    
    lin_reg = LinearRegression().fit(data_train, y_tr)
    pred = lin_reg.predict(data_test)
    
    plt.scatter(y_ts,pred)
    return {'r2_train':r2_score(y_tr, lin_reg.predict(data_train)),
            'r2_test':r2_score(y_ts, pred),
            'pred':pred}

We can see the r2 for training and test:

sparse_pred = fit(X_train,X_test,y_train,y_test,is_sparse=True)
[sparse_pred['r2_train'],sparse_pred['r2_test']]

[0.8896333645670668, 0.898030271986993]

enter image description here

sparse_pred = fit(X_train,X_test,y_train,y_test,is_sparse=False)
[sparse_pred['r2_train'],sparse_pred['r2_test']]

[0.8896302753422759, 0.8980115229388697]

enter image description here

You can try the above with the seed (42) in your example and you will see the r^2 for training is similar. It's the prediction that goes haywire.

So if you use a sparse matrix, the least square would most likely return a less nonsense coefficient for an all zero column (most likely what @piterbarg was pointing to).

Still, I think what makes sense is to check the data for such missing factors between test and train, before plunging the pipeline. For this dataset, it is most likely not over-determined so sparse vs dense should not be the difference.

Linear regression involves numerical inversion of certain matrices or, more precisely, solving certain linear equations, see here. In some cases these matrices are either singular or nearly so ("poorly conditioned"), which is often the case when the "features" of your model are either not independent or nearly not independent.

in this case straightforward inversion can lead to catastrophic error accumulation with the resulting inversion/solution basically blowing up. This is more prevalent for larger number of dimensions as you have here (36 as far as I can tell)

When a matrix is nearly singular, its inversion depends strongly on subtle differences in representations of numbers involved (their base), roundoff errors, the precise order of calculations, etc. In your case, I believe, this is the reason why the two answers are dramatically different. It seems the calculation path involving sparse-matrix representation can handle nearly-singular matrices somewhat better than the numpy representation. Why is it so I cannot say for sure beyond what I have written above.

It is easy enough to test my assertion by regularizing the problem a bit. For example, one can use Tikhonov regularization, aka Ridge regression, instead of the plain Linear Regression. Tikhonov regression adds a (small) term to the matrix that makes nearly-singular matrices behave much better under inversion.

It is a one-line change in your code:

from sklearn.linear_model import Ridge
lin_reg = Ridge(alpha=1e-4).fit(X_train, y_train)

Now both cases should produce (nearly) identical answers. Note alpha parameter that specifies the strength of regularization, and the small value I used there.

Related