# Why doScikit Learn's regression coefficients only match formal equation with standardized data?

I was playing around with LinearRegression in Scikit Learn and I found a peculiarity that I'm trying to make sense out of.

If you compare values from the coef_ attribute, they only match the values given from the equation:

if you standardize the X feature matrix.

As an example, let's look at the results from the boston housing data set:

import pandas as pd
import numpy as np
from sklearn.linear_model import LinearRegression

X             = pd.DataFrame(boston.data)
X.columns     = boston.feature_names
y             = boston.target


With non-standardized data the results between the two methods differ:

# Comparison of Coefficients Without Standardized Data
np_coeffs = (np.linalg.inv(X.T.dot(X))).dot(X.T.dot(y))

lreg      = LinearRegression()
sk_coeffs = lreg.fit(X, y).coef_

# Put results into a dataframe for easy comparison
coeff_comparison1 = pd.DataFrame({
'Variable' : X.columns,
'np_coeffs': np_coeffs,
'sk_coeffs': sk_coeffs
})


Which yields the following results:

If we standardize the data though, the results match exactly:

# comparison for coefficients with standardized data
X_std = (X - X.mean()) / X.std()

np_std_coeffs = (np.linalg.inv(X_std.T.dot(X_std))).dot(X_std.T.dot(y))
sk_std_coeffs = lreg.fit(X_std, y).coef_

coeff_comparison2 = pd.DataFrame({
'Variable' : X.columns,
'np_std_coeffs': np_std_coeffs,
'sk_std_coeffs': sk_std_coeffs
})


Which gives the following dataframe:

And as you can see the results match exactly.

I'm curious what Scikit Learn does that causes one set of results to match those of the matrix equation, and the others to diverge.

Thank you.

• Just a heads up. Instead of this np_coeffs = (np.linalg.inv(X.T.dot(X))).dot(X.T.dot(y)), it's always better to do this: np_coeffs = np.linalg.solve(X.T.dot(X), X.T.dot(y)). The equation is written mathematically with an inverse, but it's best to avoid explicitly computing an inverse wherever possible. As an aside, there's now a better way to write that out in python: np_coeffs = np.linalg.solve(X.T @ X, X.T @ y). The @ for matrix multiplication makes the whole thing more readable. Commented Feb 26, 2019 at 17:06

When you wrote

np_coeffs = (np.linalg.inv(X.T.dot(X))).dot(X.T.dot(y))


you were solving for the coefficients for a linear model with no intercept term since you did not prepend a column of 1's to the feature matrix. When you standardize your data beforehand, the resulting linear model does not have an intercept term, so your results agree in that case.

X1 = np.c_[np.ones(len(X)), X]  # Append a column for the intercept term
np_coeffs = (np.linalg.inv(X1.T.dot(X1))).dot(X1.T.dot(y))


and then compare:

lreg = LinearRegression().fit(X, y)  # Fit with original X matrix
sk_coeffs = [lreg.intercept_] + list(lreg.coef_)

coeff_comparison1 = pd.DataFrame({
'np_coeffs': np_coeffs,
'sk_coeffs': sk_coeffs},
index=['Intercept'] + list(X.columns))

print(coeff_comparison1)


The printout is

            np_coeffs  sk_coeffs
Intercept  36.459488  36.459488
CRIM       -0.108011  -0.108011
ZN          0.046420   0.046420
INDUS       0.020559   0.020559
CHAS        2.686734   2.686734
NOX       -17.766611 -17.766611
RM          3.809865   3.809865
AGE         0.000692   0.000692
DIS        -1.475567  -1.475567

• I should mention that, as a comment on the question already mentioned, explicitly inverting the matrix $X^T X$ is generally considered bad since it introduces a lot of numerical instability. Consider using numpy.linalg.lstsq or another approach for solving least squares systems like QR factorization Commented Feb 26, 2019 at 17:37