Building Linear Regression From Scratch in Python: From Matrix Math to Real Predictions

When learning machine learning, it's always easy to rely on scikit-learn to build models, but importing these algorithms and estimators makes it difficult for one to explain what is happening with these numbers and data under the hood. Building a linear regression from scratch using just numpy forces you to meet matrix dimension, cost functions, learning rate directly.

What is Linear Regression

Linear Regression is a supervised learning machine learning algorithm that models a relationship between the features/input ( $\text{X}$ ) and the target/output ( $\text{y}$ ) by choosing optimal weights/parameters. So we will understand in simple terms how scikit-learn handles this by implementing it from scratch using the Ordinary Least Square (Normal Equation) and iterative Gradient Descent Approach through vectorized matrix operations.

I am not going to go through the whole mathematics behind it in this article, if you are interested, you can find my note here on github.

Designing a learning algorithm

I am gonna define symbols i will be using through out the article.

$\text{X}$ = input/features
$\text{y}$ = prediction/target/output
${\theta}$ = weight/parameter.
$\text{n}$ = number of features
$\text{m}$ = number of rows
$J(\theta)$ = cost function/error function

You have your training data -> you pass it through a learning algorithm (Linear Regresssor) -> you get an hypothesis -> pass an input through the hypothesis -> get an output

This hypothesis is just a function with weights. The learning algorithm learns the weights for each feature and by this it will use to make predictions.

Since our weights in our hypothesis determines the prediction we get

How do we choose a good weight/parameter for our algorithm?

• We choose ${\theta}$ s that when used, we get predictions $\hat{y}$ close to our actual values $\text{y}$ .
• We calculate our error = $\hat{y}$ - $\text{y}$ and try to the least possible value of the error. If the actaul value is 100 and our model's output is 80. The error is -20 and we can know that our model predicted 20 less than the actual.

Instead of just looking at the error for a single row, we need a way to measure the total error across all rows in our dataset. To do this, we square the differences so that positive and negative errors don't cancel each other out, and we sum them up using the Mean Squared Error (MSE) formula:

$$J(\theta) = \frac{1}{2} \sum_{i=1}^{m} (\hat{y}^{(i)} - y^{(i)})^2$$

(Note: The $\frac{1}{2}$ is just a math convenience that makes the calculus cleaner when we take derivatives later).

Preparing the Data and the Intercept Column

Before feeding data into our algorithm, we load our dataset (a small CSV with roughly 30 samples). Because our hypothesis requires a y-intercept ( $\theta_0$ ), we have to append a column of ones to our feature matrix $X$ using np.hstack so matrix multiplication can compute the intercept for every row simultaneously:

import numpy as np

# Suppose X_raw is our feature matrix from our ~30 sample CSV
m = X_raw.shape[0]
X = np.hstack((np.ones((m, 1)), X_raw))  # Appending column of 1s

How Gradient Descent Finds the Best Weights

Now that we know how to calculate our error, how do we actually find the optimal weights $(\theta)$ that give us the smallest possible error? This is where Gradient Descent comes in.

Imagine you are standing at the very top of a foggy mountain, and your goal is to walk down to the lowest possible valley (the lowest elevation). That lowest elevation represents the minimum of our cost function—where our error is as small as it can possibly get.

Since it's foggy and you can't see the whole map, what do you do? You look around your feet, feel the steepness of the slope, and take a step downward in the direction that goes down the steepest. You repeat this process over and over until you can't step down any further.That is exactly how Gradient Descent works:

1. Start Randomly: We initialize our weights $(\theta)$ with random values (or zeros). At this starting point, our model's predictions are terrible, and our cost (the elevation on our mountain) is massive.

2. Calculate the Slope (Gradient): We calculate the partial derivative of the cost function with respect to each feature's weight. This tells us the slope of the hill at our current position.

3. Take a Step: We update our weights by moving in the opposite direction of the slope to reduce the cost.

The Learning Rate ( $\alpha$ )

*How big of a step do we take as we walk down the mountain? *

That is controlled by a hyperparameter called the learning rate ($\alpha$).

• If your learning rate is too big, you take giant leaps down the mountain. Instead of landing smoothly in the valley, you overshoot, bounce back and forth, or launch your weights straight into infinity.

• If your learning rate is too small, you take microscopic steps. Your model will survive, but it might take tens of thousands of epochs just to inch its way down the slope.Each step of the update rule looks like this:

$$\theta_j := \theta_j - \alpha \frac{\partial}{\partial \theta_j} J(\theta)$$

Computing Cost and Gradient Descent in Code

Now that we have our cost function, we can write the Python functions using NumPy to calculate the cost and update our weights iteratively without using any slow for loops.

Here is how the cost function looks in code:

def compute_cost(X, y, theta):
  predictions = X @ theta
  errors = predictions - y
  return 0.5 * np.sum(errors ** 2)

And here is the gradient descent step that calculates how to shift our weights $(\theta)$ to lower that cost:

def gradient_descent_step(theta, X, y, alpha=0.001):
  predictions = X @ theta # @ stands for dot product operation
  errors = predictions - y
  # X.T @ errors calculates the sum across all rows for every feature simultaneously
  gradient = X.T @ errors
  return theta - alpha * gradient

Running the Training Loop

With our functions ready, we can initialize our hyperparameters, set up our training loop to run across thousands of epochs(how much iteration), and track our best cost and weights:

EPOCHS = 200000
LEARNING_RATE = 1e-3 #0.001
np.random.seed(42) # for reproducibilty and make the random initializion of theta always same

best_theta = None
best_epoch = 0
best_cost = float("inf")

theta = np.random.randn(n, 1)
cost_history = []

for epoch in range(EPOCHS):
  theta = gradient_descent_step(theta, X_train, y_train, LEARNING_RATE)
  cost = compute_cost(X_train, y_train, theta)
  cost_history.append(cost)

  # Check for improvement
  if cost < best_cost:
    best_cost = cost
    best_epoch = epoch
    best_theta = theta

print("Initial Cost:", cost_history[0])
print(
    f"Final Cost: {cost_history[-1]} (Best Cost: {best_cost} at epoch"
    f" {best_epoch})"
)

When I ran this loop, my initial cost started massive at over 1.5 billion, but as gradient descent iterated, it dropped down to 337 million, finding its optimal weights right around epoch 6,128.

(Note: Alternatively, for linear regression, you can also solve for weights in a single closed-form step using the Normal Equation: $\theta = (X^T X)^{-1} X^T y$ , but iterative gradient descent gives us a clear look at how models learn dynamically!)

def normal_equation(X,y):
    return np.linalg.inv(X.T@X) @ X.T @ y
one_shot_weight = normal_equation(X_train,y_train)
print(one_shot_weight)

Where It Broke (Mistakes and Fixes)

• The Learning Rate Trap: When I first tried a learning rate of 0.01, my model didn't gently roll down the mountain, it violently overshot and threw overflow warnings and my predictions were NaN. Dropping it down to 0.001 was what finally kept the weights stable.

• Forgetting the Intercept Column: Early on, my matrix multiplications kept throwing shape mismatch errors until I remembered to explicitly use np.hstack to add a column of ones so that $\theta_0$ had a feature to multiply against.

What I'd Tell My Past Self

• Always check your shapes: Print out X.shape and theta.shape before debugging matrix multiplication errors.
• Watch your initial cost: If your first epoch's cost is in the trillions or throws NaN, stop the script immediately, your learning rate is too aggressive.
• Vectorization over loops: Let NumPy handle the dot products (@) instead of writing slow Python for loops over your rows.

If you want to dive deeper into the full mathematical derivations and proofs behind these equations, you can check out my complete notes stored in my public Obsidian vault on GitHub.

Resources I Used

• Stanford CS229 Lecture Notes (Andrew Ng & Tengyu Ma)
• Gemini Notebook (pka Notebook LM) + Gemini
• NumPy Official Documentation for Matrix Multiplication (@) and random methods (default_rng)