Understanding the math behind regression from .fit method

Question

Just to preface, I am very new to machine learning and linear regression.

I have this simple data set with attributes: area, bedrooms, and age...trying to predict the price of a house.

Here is a screenshot of my Jupyter notebook

I am having trouble understanding how these coefficients were calculated (137.25, -26025, -6825). How were these values calculated? I have tried the least-squares solution. This was my work:

As you can see I used the formula $(A^T A)^{-1} A^T b$. The matrix I calculated does not match the coefficients calculated in reg.coef_. What is the math behind these calculations?

Answer 1

Here is the math:

$$ X = \begin{pmatrix} 1 & 2600 & 3 & 20 \\ 1 & 3000 & 4 & 15 \\ 1 & 3200 & 3 & 18 \\ 1 & 3600 & 3 & 30 \\ 1 & 4000 & 5 & 8 \\ \end{pmatrix}, \quad Y = \begin{pmatrix}550000\\ 565000\\ 610000\\ 595000\\ 760000\end{pmatrix} $$

import numpy as np
X = np.array([[1, 2600, 3, 20], 
              [1, 3000, 4, 15],
              [1, 3200, 3, 18], 
              [1, 3600, 3, 30], 
              [1, 4000, 5, 8 ]])
Y = np.array([550000, 565000, 610000, 595000, 760000])

np.linalg.inv(X.T @ X) @ X.T @ Y
array([ 3.83725e+05,  1.37250e+02, -2.60250e+04, -6.82500e+03])

These are the results given in your linear regression model above.

Note that in some python versions ie 3.6.11 when you do $X^TY$ then you post multiply it to the matrix $(X^TX)^{-1}$ you get a funky result. ie:

 np.linalg.inv(X.T @ X) @ (X.T @ Y)
 array([-9446381.39872008,   -17633.1771872 , 13819875.82048003, 997123.60544   ])

or even

 np.linalg.solve(X.T @ X, X.T @ Y)
 array([-9446381.39872008,   -17633.1771872 , 13819875.82048003, 997123.60544   ])

Further check in python:

 np.linalg.inv(X.T.dot(X)).dot(X.T.dot(Y))
 array([-9446381.39872007,   -17633.1771872 , 13819875.82048003, 997123.60544   ])
 
 np.linalg.inv(X.T.dot(X)).dot(X.T).dot(Y) # Look at the paranthesis:
 array([ 3.83725e+05,  1.37250e+02, -2.60250e+04, -6.82500e+03])

This has been fixed in later versions, eg 3.8. Therefore you might not be getting the correct results due to a bug in the python version you are using.

Using R, we end up with the same results.

X = array( c(1, 1, 1, 1, 1, 2600, 3000, 3200, 3600, 4000, 3, 4, 3, 3, 5, 
         20, 15, 18, 30, 8), dim = c(5, 4))
Y = c(550000, 565000, 610000, 595000, 760000)

solve(crossprod(X), crossprod(X, Y)) # solve(t(X) %*% X, t(X) %*% Y)
          [,1]
[1,] 383725.00
[2,]    137.25
[3,] -26025.00
[4,]  -6825.00

solve(t(X) %*% X) %*% t(X) %*% Y 
          [,1]
[1,] 383725.00
[2,]    137.25
[3,] -26025.00
[4,]  -6825.00

Although this is the basic notion for linear regression, note that all the regression platforms do not try to compute the inverse of the matrix directly. This is due to the the process being so complex and expensive. Usually QR decomposition is employed.

Answer 2

Whenever I want to learn how software works, I read the documentation.

In this case, when reading sklearn.linear_model.LinearRegression, I notice that the class includes an intercept by default.

fit_intercept - bool, default=True Whether to calculate the intercept for this model. If set to False, no intercept will be used in calculations (i.e. data is expected to be centered).

The design matrix you used to compute $\text{coef}=(A^\top A)^{-1}A^\top b$ does not include an intercept term. We would not expect these two models to match in general because you're comparing two different models.

To make the results match, you'll have to do either of two things

1. Add an intercept to your computations by hand. This means appending a column of 1s to $A$ and then compute $\text{coef}=(A^\top A)^{-1}A^\top b$.
2. Remove the intercept from the sklearn model (use fit_intercept=False). (But this is generally inadvisable. See: When is it ok to remove the intercept in a linear regression model?)

Notes that floating point arithmetic is susceptible to roundoff error, so your hand-cranked solutions may not match to a large number of digits (or, in the worst case, you may lose all digits of precision). QR factorization is generally used to compute regression estimates because it is less susceptible to roundoff error (see: Understanding QR Decomposition).