Implementation of Linear Regression Closed Form Solution

Viewed 413

I wonder if you all know if backend of sklearn's LinearRegression Module uses something different to calculate the optimal beta coefficients. I implemented my own using the closed form solution

    if self.solver == "Closed Form Solution":
        ### optimal beta = (XTX)^{-1}XTy
        XtX = np.transpose(X, axes=None) @ X
        XtX_inv = np.linalg.inv(XtX)
        Xty = np.transpose(X, axes=None) @ y_true
        self.optimal_beta = XtX_inv @ Xty

However, I do not get an exact match when I print the coefficients comparing with sklearn's one. I thought that having a closed form solution may guarantee similar results (my code doesn't handle not invertible).

I noticed the MSE are not that far apart, but it just makes me wonder if my implementation is wrong (i took in account of biases)

The deviations seem to appear in the intercept:

SKLEARN INTERCEPT 1.2490009027033011e-15
HN INTERCEPT 4.440892098500626e-16
First Value HN 183.22200150945497
First Value SKLEARN 183.22200150945548
HN MSE 3.1084228599404546e-27
SKLEARN MSE 2.126667544039961e-25
1 Answers

Figures are in the order of 10^-15 which means they are practically 0, and the same!

As to why there is a difference: you are solving the normal equations by directly inverting the X.T @ X matrix. On the other hand, scikit-learn uses scipy.linalg.lstsq under the hood, which uses for example an SVD-based approach. That is, the mechanism there does not invert the matrix and is therefore different than yours. Note that there are many ways to solve the linear least squares problem.

But we can say that the results you get are the same as their way, because these operations over matrices are subject to round-off errors anyway and yours "matches" theirs (or vice versa) for many significant digits.

Related