Linear Regression: Part 1

Here we will be exploring:

• Definition and equation of Linear Regression
• Solving Linear Regression by Normal Equation Approach
• Solving Linear Regression by Gradient Descent
• Code both approaches from scratch by subclassing scikit-learn base class
• Understand why solve it by Gradient Descent even if Normal Equation Approach is available

Let’s Begin!!

To start with let’s see the equation of a line:

y = mx + c (Simple!!)

m = slope(change in y for 1 unit change in x)

c = point where line crosses the y-axis

Now, let’s expand this idea to multiple dimensions.

y = m₁x₁ + m₂x₂ + …….. + mₙxₙ + c

Machine learning likes the below notation:

y = β₁x₁ + β₂x₂ + ……. + β₀

They look the same right!!

Now this longer equation is of a Hyperplane. It’s like a line but for n-dimensions.

To be correct, we should think of it as line is a limiting condition of Hyperplane and that should be evident from their equations as well.

Why are we concerned about lines (Hyperplanes in n-dimension)?

As the topic is Linear Regression it is evident from the Linear part that line plays a crucial role here.

Now coming to the Regression part:

It means predicting values that are continuous with continuous numbers. Like weights of people, salary, height etc. which can take any value on a number line.

Now we are given a tabular data where we have multiple features (columns). What we will try to do here is using the equation of a hyperplane we will try to predict the value of a continuous feature (target variable) w.r.t. to the rest of the feature variables (predictor features).

If we take each row of the dataset and try to understand it w.r.t. the hyperplane equation:

y = β₁x₁ + β₂x₂ + ……. + β₀

we will see that the x₁, x₂,….xₙ are the values of data in each column.

y is the continuous feature (a number) which we are trying to predict.

Caution: One should not assume that y is perfectly linearly related to the other feature variables. What I am saying is y is somewhat linearly related to each of the feature variable. And this assumption forms the basis of linear regression. If we do observe that y is nowhere linearly related to the predictors, using Linear Regression is not going to help us to predict the target variable.

Then where are the betas?

The value for β are for us to find out.

We will be using various mathematical tools to achieve the same.

Before that let’s try to understand how to write the same equation in the for of Linear Algebra. (Here I assume reader has some prior knowledge of vectors, matrices and mathematical operations related to them)

Lets say in our dataset we have n-rows and m predictor features. Then the equation for each row of data can be written as below:

Figure-1

Now, this looks long and tedious even by notation. Lets convert them to Linear Algebra form as show below:

Figure-2

y_hat becomes a vector of dimension nX1

β becomes a vector of dimension (m+1) X 1

X becomes a matrix of dimension n X (m+1)

NOTE: In m+1, the +1 is to there to include the intercept. Since the data set only has predictor features (m columns)

Now using the concept of matrix multiplication the equation can be re-written as:

Figure-3

We see how small and intuitive it is to read. Not only that computers can perform the computations much faster in parallel as they are vectorised.

So out agenda here is to find the values of β vector here.

Now let’s allow intuition to kick in

We know the values of two out of the three variables in the equation above

y_hat and X RIGHT!!

It’s the data itself.

So to find β just take the product and y_hat and X inverse and there we have the values of β. If you too think so, yes we can.

BUT

There is an issue. If you are familiar with inverse of matrices, we know that finding X_inverse is difficult. Lets see why.

Below are the conditions a matrix should fulfil to be inverted:

• The matrix should be a square matrix
• The matrix should be singular. Which means each column in the matrix should be linearly independent of the other columns.

We know that X is nothing but our dataset. And in real world we hardly expect the number of rows and number of columns to be equal in our dataset. Also there might be columns in our data that can be linearly dependent on one or a combination of our features in our dataset.

So in all aspects X is violating the rules of being inverted.

Get Saurav Naik’s stories in your inbox

Join Medium for free to get updates from this writer.

Remember me for faster sign in

But we still have a way out of it. We can find the Moore-Penrose inverse of X also known simply as pseudo-inverse. I will not be going into the details of pseudo inverse here but intuitively we can say that pseudo-inverse will give us the inverse of matrices even if it violates the fundamental conditions of matrix to be inverted.

So we can write β as below:

Figure-4

Numpy as a method np.linalg.pinv() which can calculate the pseudo-inverse.

Let’s see an example:

# importing libraries
from sklearn.datasets import load_diabetes
import numpy as np
import pandas as pd
from sklearn.utils import check_X_y

# loading a toy dataset
data = load_diabetes()
X, y = data['data'], data['target']
col_names = [f"f_{idx}" for idx,_ in enumerate(range(len(X[0])),start=1)]

# reshaping to keep the data as shown in equation
X = pd.DataFrame(X, columns=col_names)
y = y.reshape(-1,1)

# checking the shapes of data
print(f"Number of rows: {X.shape[0]}")
print(f"Number of columns: {X.shape[1]}")
print(f"Shape of the target variable: {y.shape}")

Number of rows: 442

Number of columns: 10

Shape of the target variable: (442, 1)

Let’s make sure our data violates both conditions for it to be inverted.

From the number of rows it is already evident that the matrix X is not a square matrix by any means.

For the independence part, it will not be possible for me here to show that a single feature is a near resultant of a combination of features, but just showing that it is one feature is sufficiently correlated to another should be enough of a proof to show that the features are not entirely independent.

import seaborn as sns
import matplotlib.pyplot as plt

plt.figure(figsize = (10,10))
plt.title('Correlation between Predictor Features')
sns.heatmap(X.corr(), square= True, annot=True, cmap="inferno")
plt.show()

Figure-5

We see from the above image that f_7 and f_8 are somewhat highly correlated showing a linear relation between them.

So I believe I am able to convince you that matrix X cannot be inverted.

Before going forward, remember in Figure-1 we added a columns of 1 to the data matrix to include the intercept lets do the same with our dataset here.

X['const'] = 1

Our dataset now looks like this:

Figure-6

Now our data set has 11 columns and 442 rows

Let’s apply the equation we saw in Figure-4 to find our weight vector β

# checks data formatting like non-null and so on
# you can read more on it in the scikit-learn documentation
x,y = check_X_y(X,y)
beta = np.linalg.pinv(x)@y
beta

β₁ → -10.00987

β₂ → -239.81564

β₃→ 519.84592

β₄ → 324.38465

β₅ → -792.17564

β₆ → 476.73902

β ₇ → 101.04327

β₈ → 177.06324

β₉ → 751.2737

β₁₀ → 67.62669

β₀ → 152.13348

We get the weight vectors Simple!!

Now let me tell you an Interesting thing.

This is the implementation that scikit-learn uses internally for Linear Regression. Please feel free to search the web on the same.

Let me show you whether we get similar results with scikit-learn as well. That should be able to convince us to a certain extent that what we did is indeed a correct approach.

from sklearn.linear_model import LinearRegression
# using fit_intercept as false as we have already included a
# 'const' column in our data with all values 1 to calculate the intercept
reg_model = LinearRegression(fit_intercept=False)
x,y = check_X_y(X,y)
reg_model.fit(x,y)
params = reg_model.coef_.tolist()
params

β₁ → -10.00987

β₂ → -239.81564

β₃→ 519.84592

β₄ → 324.38465

β₅ → -792.17564

β₆ → 476.73902

β ₇ → 101.04327

β₈ → 177.06324

β₉ → 751.2737

β₁₀ → 67.62669

β₀ → 152.13348

We get the same values for the beta vector.

Well all seems good right….Well not everything

What I personally feel is that we calculated a pseudo inverse for a matrix which cannot be inverted. Intuitionally I feel that since its pseudo, the inverse is not a proper inverse so that values we get may not be the best one.

Also, the time complexity for np.linalg.pinv() is approx. O(n²)

This means that as our data doubles, our computation time is multiplied roughly by 4 times.

So, there must be other…lets say better ways to find the weights β.

We will continue on this journey in Part 2……..