how to implement linear regression from scratch

Introduction to Linear Regression

Imagine you're a real estate agent trying to predict a house's price. You know one key fact: its size. You might suspect bigger houses cost more. If you plot size on the x-axis against price on the y-axis, you'd likely see points trending upward. Linear regression is simply the process of drawing the straight line that best follows that trend.

That line has an equation you might remember from school: y = mx + b.

y is the Price (what we want to predict).
x is the Size (the data we have).
m is the slope (how much price changes per square foot).
b is the intercept (the base price when size is zero).

The "best" line is the one where, on average, its predictions are closest to the actual data points. We measure "closeness" with a loss function, typically Mean Squared Error (MSE), which penalizes larger errors more heavily.

Visualizing the "Best Fit"

Professor's Tip: The blue dots are real data. The red line is our model's guess. When you click "Run Training," the model adjusts the line to minimize the distance (error) between the line and the dots.

Common Misconception: Linear models are always accurate

A straight line is simple. That's its power and its limit. Real-world data—like house prices—is messy. Location, age, and market conditions also matter. A single straight line using just size will be a rough, foundational approximation, not a perfect predictor.

Thinking of it as your first, honest attempt at finding a pattern is key. Its value is in interpretability and as a baseline, not in capturing complex, non-linear relationships.

Why Learn Linear Regression from Scratch?

You could call a library function in one line. But doing it yourself—defining the loss, calculating gradients, and updating parameters—is how you internalize the core loop of machine learning.

The "Guess, Measure, Adjust" Cycle

This process via gradient descent is the engine for nearly every model you'll build later.

From User to Builder

When you see a neural network with thousands of parameters, you'll recognize the same fundamental process at work, just scaled up.

Setting Up Your Python Workshop

Before we write a single line of code, we need to prepare our workshop. Think of Python as the carpenter, but it needs a clean, organized space to work. This setup isn't just busywork; it ensures your code runs smoothly and stays consistent for anyone else (or your future self) who tries to reproduce it.

The Essential Toolkit

We need NumPy for the heavy math (fast arrays and matrix operations) and Matplotlib for drawing our charts.

pip install numpy matplotlib

If you try to run code without these, you'll hit the dreaded ModuleNotFoundError. This simply means the library isn't installed in your current environment. The command above downloads them from the Python Package Index (PyPI).

The "Dependency Conflict" Problem

Why we isolate our projects (Virtual Environments)

Project A Needs Library v2.0

Project B Needs Library v1.0

Global Python Environment

v1.0

v2.0

Status: Both projects are fighting for control of the single Global Environment. One library version will break the other.

The Solution: Virtual Environments

A virtual environment creates a private, project-specific toolbox. When you activate it, pip install only affects that one project, keeping its dependencies neatly contained.

Terminal

➜ python -m venv venv

# Creates a folder named 'venv' for your project

➜ source venv/bin/activate # (Mac/Linux)

➜ venv\Scripts\activate # (Windows)

# Your terminal prompt will now show (venv)

Why Bother?

This is a non-negotiable best practice for any serious work. It ensures that if you share your project, others can create their own identical environment from your requirements.txt file, guaranteeing they run the exact same code with the exact same tools.

Understanding Gradient Descent

You have your line y = mx + b. Right now, m and b are just guesses. You calculate the Mean Squared Error (MSE) for those guesses—that's your loss, a single number telling you how wrong you are.

Your goal is simple: find the m and b that make this loss as small as possible.

The "Blindfolded Hiker" Analogy

Imagine the loss is a landscape of hills and valleys. The height of the land at any point represents your error. The bottom of the deepest valley is your perfect m and b.

You are standing on this landscape blindfolded. How do you get to the bottom? You feel the ground under your feet. The steepest downhill direction is the negative gradient. You take a small step in that direction, feel the ground again, and repeat.

This "guess, measure slope, step" loop is Gradient Descent. The math behind the "step" is governed by a parameter called the Learning Rate.

Visualizing the Descent

Watch how the "Learning Rate" affects our path to the minimum.

Current Status:

Select a rate to start

Professor's Note:

Low: Safe, but very slow.
Medium: The sweet spot.
High: Might overshoot and bounce around!

The Code: How We Update

Here is the core update rule derived from the math. The gradient tells you the direction (uphill), and the learning rate controls the distance.

gradient_descent.py

# 1. Calculate the gradient (slope of the loss)

grad_m

= (2 / n) * sum(error * X)

grad_b

= (2 / n) * sum(error)

# 2. Update parameters (Move opposite to gradient)

# m = m - (learning_rate * slope)

= m - (lr * grad_m)

= b - (lr * grad_b)

Notice the minus sign in the update rule (m = m - lr * grad). Since the gradient points uphill (towards higher error), we subtract it to move downhill towards the minimum.

Does it always work?

Gradient descent finds a local minimum. In our simple Linear Regression case, the loss landscape is a smooth "bowl" (convex), so any local minimum is actually the global best. However, in complex Deep Learning models, the landscape is bumpy, and we might get stuck in a small dip that isn't the absolute best solution.

Implementing Linear Regression from Scratch

You already know the goal: find m and b that minimize MSE. You also know the gradient descent update rule. Now, let's turn that rule into a repeating process—a loop that gradually improves your guesses.

The Core Mental Model: The Training Loop

Start with random guesses for m and b. Then, repeatedly:

Predict using current m and b.
Measure the error (difference from true y).
Calculate how to adjust m and b to reduce error (the gradients).
Nudge m and b a small step in the right direction (using the learning rate).

This cycle is your training loop. Each full pass through all your data points is one iteration (or epoch). With each iteration, your line should fit the data a little better, and the total loss should decrease.

Visualizing Convergence

Watch the Loss (Error) decrease with every iteration.

Current Iteration

Current Loss (MSE)

High

Status

Idle

Professor's Note: Notice how the curve is steep at first (rapid improvement) and then flattens out. That flat part is where we hit our tolerance and stop training.

Common Pitfall: Incorrect Convergence Criteria

How do you know when to stop? A common beginner mistake is to pick an arbitrary fixed number of iterations (e.g., "run 1000 times") and hope it's enough. This is fragile: too few and you stop before reaching a good solution; too many and you waste time computing after you've already converged.

Better Approach: Dynamic Stopping

Stop when the improvement becomes negligible. Track the loss after each iteration. If the change in loss is smaller than a tiny tolerance (e.g., 1e-6) for several consecutive iterations, declare convergence.

# Pseudocode for convergence check
if abs(previous_loss - current_loss) < tolerance:
    convergence_counter += 1
else:
    convergence_counter = 0 # reset if we make significant progress

if convergence_counter >= patience:
    break # stop early

This prevents both under-fitting (stopping too soon) and wasted computation (running long after convergence).

Coding the Gradient Descent Loop

1. Initializing Parameters

Before the loop, you need starting values for m (slope) and b (intercept). For this simple, convex problem, initializing both to 0 is perfectly safe and common.

setup.py

# 1. Initialize parameters

= 0.0

# Hyperparameters

= 0.01 # learning rate: step size

tolerance

= 1e-6 # convergence tolerance

patience

= 10 # iterations with tiny improvement needed to stop

Why `0.0`?

Starting at zero is neutral. Since the loss landscape is a smooth bowl (convex), any starting point leads to the same bottom. Zero is simple and reproducible. (For more complex, non-convex models like neural networks, small random initializations help break symmetry, but that's not needed here.)

2. The Training Loop

Now, the loop itself. You'll compute predictions, error, gradients, and updates in vectorized form using NumPy—this is fast and avoids slow Python for loops over individual data points.

training_loop.py

# 2. Training loop

convergence_counter

= 0

prev_loss

= float('inf') # start with infinitely large previous loss

# Safe upper bound, though we usually break early
for iteration in range(10000):
# Forward pass: predict

y_pred

= m * X + b # vectorized: all predictions at once

    # Compute error and loss

error

= y_pred - y

loss

= np.mean(error ** 2) # MSE

    # Check convergence
    if abs(prev_loss - loss) < tolerance:
        convergence_counter += 1
    else:
        convergence_counter = 0

    if convergence_counter >= patience:
        print(f"Converged after {iteration} iterations.")
        break

    prev_loss = loss

    # Compute gradients

grad_m

= (2 / len(X)) * np.sum(error * X)

grad_b

= (2 / len(X)) * np.sum(error)

    # Update parameters (move downhill)
    m = m - lr * grad_m
    b = b - lr * grad_b

Key Takeaways from the Code

Vectorization: y_pred = m * X + b uses broadcasting. m (scalar) multiplies every element in X, then b (scalar) is added to every result. This is one of NumPy's most powerful features—it replaces explicit loops.
Loss Calculation: np.mean(error ** 2) is the MSE. Squaring the error array element-wise, then averaging.
The Update: The update subtracts lr * gradient from the current parameter. Remember: gradient points uphill; subtracting moves downhill.
Robustness: The convergence check uses abs(prev_loss - loss) < tolerance. The patience counter prevents stopping on a single lucky small improvement; we want sustained near-zero change.

What You've Built

You have just built a complete, from-scratch linear regression engine. It learns m and b from data using gradient descent. The loop runs until the loss stabilizes, leaving you with a fitted line y = m*x + b ready for prediction.

Data Preparation for Predictive Modeling

Before your gradient descent loop can even start, you need data it can actually learn from. This means collecting relevant features (like house size) and the true target values (like sale price), then ensuring that data is trustworthy.

The Golden Rule: Garbage In, Garbage Out

A single NaN or an outlier (like a house listed for $1) will propagate through your calculations. Your predictions y_pred = m*X + b will be NaN if X contains NaN. Your gradients np.sum(error * X) will be wildly distorted by an outlier, pulling your learned line in the wrong direction.

Collecting and Cleaning Data

The most critical step here is handling missing or erroneous values. Your first task is to curate a clean, consistent dataset. This might mean:

Removing rows with critical missing values.
Imputing reasonable values for missing entries (e.g., median size for a missing house size).
Identifying and investigating outliers—deciding if they are true rare events or data entry mistakes.

A smaller, pristine dataset is a far better teacher than a large, noisy one.

Common Misconception: More data always improves model

Quantity cannot compensate for poor quality. Adding 10,000 records with 30% missing values or systematic measurement errors will degrade performance.

Your gradient descent algorithm is a precise, mathematical procedure. It needs a clean signal to find the pattern.

Feature Scaling and Normalization

You have your clean X (feature) and y (target) arrays. But there's a hidden performance factor: the scale of your features.

Imagine your X is house size in square feet, ranging from 500 to 5000. Now imagine a second feature, "number of bedrooms," ranging from 1 to 5. These features exist on vastly different scales.

Visualizing the Loss Landscape

How feature scaling changes the shape of the "valley" we are walking down.

Unscaled Data: The loss landscape is a narrow, steep "V". Gradient descent struggles here, often zig-zagging across the valley walls, requiring tiny steps and many iterations.

The Solution: Min-Max Normalization

To fix this, we squish every feature value into a fixed range, typically [0, 1]. This makes the loss landscape "rounder" and easier to navigate.

scaling.py

# Before scaling

# X

ranges from 500 to 5000 (sq ft)

# After Min-Max scaling

X_min

= X.min()

X_max

= X.max()

X_scaled

= (X - X_min) / (X_max - X_min)

# Now X_scaled ranges from 0.0 to 1.0

Why this works: Now all features contribute to the loss in a comparable numerical range. The gradient components for all parameters will be on a similar scale, leading to balanced updates.

Crucial Note on Prediction

You must apply the same scaling transformation to any new data you want to predict on. If you scaled your training data with X_min and X_max, a new house size x_new must be transformed as (x_new - X_min) / (X_max - X_min) before plugging it into your final equation.

Linear Regression in the Machine Learning Workflow

You've built a working linear regressor from scratch. You know it finds the best straight line through your data. Now, where does this tool live in a real machine learning workflow?

Think of linear regression as your first, honest question to the data. Before you stack neural networks or ensemble complex models, you ask: "Is there even a linear signal here?" You fit your simple y = mx + b model and check its performance.

The "Start Simple" Philosophy

If your clean, scaled data achieves decent predictive accuracy (say, R² = 0.80), that tells you something profound: the core relationship is fundamentally linear. Your feature X has a strong, consistent influence on y.

This is a win. A simple, interpretable model that works well is better than a complex black box. In industry, any complex model you try next must beat this baseline to justify its added cost. If a neural network doesn't significantly outperform your simple line, you should just use the line.

The Baseline Benchmark

Is a simple line enough, or do we need complexity?

Click "Check R² Score" to evaluate how well this simple line explains the data.

Common Misconception: Linear models capture everything

This is the critical limitation. A straight line is by definition a constraint. It assumes the effect of X on y is constant everywhere—the slope m doesn't change. It cannot capture curves, thresholds, or interactions unless you manually engineer new features (like adding X²).

High Bias (Underfitting)

If your data shows a clear parabolic trend, your linear model will systematically fail in specific regions. Its errors will be biased—it will have a persistent, structured pattern in its residuals.

This isn't a flaw in your gradient descent code; it's a model assumption violation. The math guarantees the best possible straight line, but if the true pattern isn't straight, "best possible" still means "wrong in a predictable way."

Visualizing High Bias (Underfitting)

When a straight line tries to fit a curve.

Professor's Note: Notice how the red line (the linear model) cuts through the middle of the curve. It misses the data at the extremes. This is underfitting. The model is too simple to capture the underlying complexity.

Integrating with Machine Learning Pipelines

Your from-scratch implementation is a learning tool. In practice, you'd use sklearn.linear_model.LinearRegression for production. But the concepts you implemented—loss, gradient descent, convergence, scaling—are the same. Here's how it slots into a standard pipeline:

pipeline_concept.py

# 1. Load & clean data → (X_clean, y_clean)

# 2. Split data → (X_train, X_test, y_train, y_test)

# 3. Scale features → fit scaler on train, transform train & test

# 4. Train model → model.fit(X_train_scaled, y_train)

# 5. Evaluate → model.score(X_test_scaled, y_test) # R² score

# 6. Deploy → model.predict(new_scaled_data)

Your implementation already handles steps 3 and 4 (scaling and the training loop). Step 5—evaluation—is crucial. The primary metric for regression is R-squared (R²), which tells you the proportion of variance in y your model explains. An R² of 1.0 is perfect; 0.0 means your model is no better than just predicting the mean of y.

metrics.py

# Compute R² manually

y_pred

= m * X_test_scaled + b

ss_res

= np.sum((y_test - y_pred) ** 2)

ss_tot

= np.sum((y_test - np.mean(y_test)) ** 2)

= 1 - (ss_res / ss_tot)

When to Transition to Advanced Models

You move beyond linear regression when one or both of these conditions are true:

High Bias (Underfitting): Your baseline model's performance is poor (low R², high MSE), and residual plots show clear structure (e.g., a curve). This means your linear assumption is too restrictive. You need a more flexible model: Polynomial Regression (adding X², X³... features), decision trees, or kernel methods.
Complex, High-Dimensional Data: You have many features (dozens or thousands), possibly with interactions. While you could manually engineer interaction terms (X1 * X2), it becomes impractical. Models like Regularized Linear Regression (Ridge/Lasso) automatically handle feature selection and multicollinearity, or you move to non-linear models like random forests.

The Practical Workflow

Always start with your clean, scaled linear regression baseline. If it works well, stop. If it fails due to bias, try a slightly more complex model (e.g., polynomial features). If it fails due to noise or high dimensionality, try regularized linear models. Only when those fail do you reach for the heavy, less interpretable tools like deep neural networks. This disciplined approach saves immense time and produces more robust, understandable solutions.

Your from-scratch experience gives you the intuition to diagnose why a model is failing. Is the loss landscape ill-conditioned? Is it stuck in a poor local minimum? That diagnostic sense is what separates a practitioner from a mere library caller.

Evaluation Metrics and Model Assessment

You have a fitted line y = mx + b. But how do you know if it's actually good? You need a way to measure performance objectively.

Measuring Performance with R-squared (R²)

The most common metric for regression is R-squared. Think of it as a score from 0.0 to 1.0 (though it can be negative for very bad models).

The Intuition

R-squared answers: "What percentage of the variation in my target variable (y) is explained by my model?"

0.0: Your model is no better than just guessing the average value of y every time.
1.0: Your model perfectly predicts every data point.

Technically, it compares your model's error to the error of a simple baseline (predicting the mean).

metrics.py

# 1. Calculate Total Sum of Squares (Variance in y)

ss_tot

= np.sum((y_test - np.mean(y_test)) ** 2)

# 2. Calculate Residual Sum of Squares (Error of your model)

ss_res

= np.sum((y_test - y_pred) ** 2)

# 3. R-squared Formula

= 1 - (ss_res / ss_tot)

The Trap: R-squared always increases

If you add more features (variables) to your model, R-squared always goes up, even if those features are garbage. It's a "lazy" metric.

A high R² is good, but it doesn't mean your model is correct. You need to look deeper.

Visualizing the Leftovers (Residuals)

Residuals = Actual Value - Predicted Value

Select a model type to see what the residuals look like.

Residual Analysis: The Diagnostic Tool

While R-squared gives you a score, residual plots tell you why you got that score. A residual is simply the error for each data point: error = y_true - y_pred.

We plot these errors against the predicted values. If your linear model is doing its job correctly, the residuals should look like random static noise centered around zero.

Random Scatter (Good): The points are scattered all over the place with no shape. This means the model captured the trend, and the remaining error is just natural noise.
Curved Pattern (Bad): If you see a U-shape or an inverted U, your model is underfitting. It's trying to fit a straight line to a curved reality. You need a non-linear model (like Polynomial Regression).
Funnel Shape (Bad): If the error gets bigger as the predictions get bigger (a cone shape), this is Heteroscedasticity. Your model is less confident for larger values.

Professor's Checklist

Never trust an R² score blindly. Always check the residuals. If you see a pattern, your model has failed to capture that pattern, and it is systematically wrong in specific ways.

Real-World Constraints and Limitations

You've built a gradient descent loop that works perfectly on a small, tidy dataset. But what happens when you point it at a real-world problem—like predicting prices from millions of records with hundreds of features? The first shock is memory and computation time.

The Normal Equation's Breaking Point

Earlier, you saw the elegant closed-form solution: w = (Xᵀ X)⁻¹ Xᵀ y. This is fast and exact for small data. But its cost is dominated by inverting Xᵀ X, a matrix of size (n_features + 1) x (n_features + 1).

The Math Trap

The inversion is roughly O(n_features³).
With 100 features, that's ~1 million operations—trivial.
With 10,000 features, that's ~1 trillion operations—prohibitively expensive.
Worse, you must hold Xᵀ X in memory. For 10,000 features, that's an 800 MB matrix. Scale to 100,000 features, and you need 80 GB of RAM just for that matrix.

The normal equation is fundamentally memory-bound and computationally cubic in feature count. It does not scale to modern "big data" sizes. Your gradient descent loop, in contrast, only ever needs to hold a few vectors in memory at once.

Scaling Gradient Descent for Big Data

The gradient descent loop you wrote processes the entire dataset on every iteration—this is Batch Gradient Descent. For 1 million samples, each iteration means 1 million predictions, 1 million error calculations, and summing over all of them. That's a lot of work per update.

Mental Model Shift

Instead of asking "What is the true gradient over all data?" before each step, ask: "What is a good gradient estimate from a small, random subset?" You trade the exact, smooth path of batch GD for a noisier but vastly more frequent and memory-light update.

The Solution: Mini-Batch Gradient Descent

Mini-Batch Gradient Descent picks a small random batch (e.g., 32, 64, 128 samples) per iteration.

Pros: Best of both worlds. Much faster per iteration than batch GD, and the gradient estimate is less noisy than Stochastic GD (SGD), leading to smoother convergence. Leverages vectorized NumPy operations efficiently on the batch.
Cons: Requires tuning an extra hyperparameter: batch_size.

Visualizing the Update Paths

Compare the "True Gradient" (Batch) vs. "Estimated Gradient" (Mini-Batch).

Select a method to see how the model learns.

Implementing Mini-Batch GD from Scratch

Here is how you implement the core logic. Instead of iterating over the whole dataset X, you iterate over batches generated from a shuffled index.

mini_batch_gd.py

# 1. Helper to generate batches

def

generate_batches(X, y, batch_size):
    n_samples = X.shape[0]
    indices = np.random.permutation(n_samples) # Shuffle indices
    for start in range(0, n_samples, batch_size):
        end = min(start + batch_size, n_samples)
        batch_idx = indices[start:end]
        yield X[batch_idx], y[batch_idx]

# 2. The Training Loop

batch_size

= 64
for epoch in range(100):
    batches = generate_batches(X, y, batch_size)
    for X_batch, y_batch in batches:
        y_pred = m * X_batch + b
        error = y_pred - y_batch
        # Gradients are averaged over the BATCH size
        grad_m = (2 / batch_size) * np.sum(error * X_batch)
        grad_b = (2 / batch_size) * np.sum(error)
        # Update parameters
        m -= lr * grad_m
        b -= lr * grad_b

Key Implementation Details

Shuffling: We shuffle the data at the start of every epoch to prevent the algorithm from learning order artifacts.
Batch Size: We divide by batch_size instead of len(X) because we are averaging the gradient over just that small slice.
Convergence: Because the loss is noisy, we usually check convergence by evaluating the full dataset loss once per epoch, rather than every batch update.

Practical Considerations

Learning Rate Scheduling: With mini-batches, start with a higher learning rate and decay it over time. This lets you make rapid progress initially, then fine-tune precisely as you near the minimum.
Feature Scaling is Non-Negotiable: If one feature has values in the thousands and another in decimals, the gradient components will have wildly different scales. This makes choosing a single learning rate impossible. Always scale features before using mini-batch GD.
When to Abandon GD: If you have billions of samples, consider using sklearn.linear_model.SGDRegressor (highly optimized) or streaming algorithms where you never store all data.

FAQ: Troubleshooting & Best Practices

You've built the model, you've trained it, and you've evaluated it. But when you run your own code, things might not go perfectly. Here are the most common hurdles students face when implementing Linear Regression from scratch, and how to clear them.

Why is my Loss not going down?

Visualizing the effect of the Learning Rate (LR) on convergence.

Learning Rate 0.01 is usually a safe starting point for scaled data. The loss decreases steadily.

Why does my gradient descent not converge?

The most common culprits are:

Learning rate too large: You are overshooting the minimum, causing the loss to oscillate or increase.
Unscaled features: If features have vastly different ranges, the loss landscape becomes a narrow, steep valley that is hard to navigate.
Loose convergence criteria: You might be stopping too early or checking the wrong metric.

Fix: Check your learning rate first—try 0.001 or 0.0001. Ensure you scaled your features to similar ranges (e.g., [0,1]).

Monitor the full training loss periodically:

monitor.py

# Example: Safe convergence monitoring per epoch

iteration % 100 == 0:
full_loss = np.mean((y - (m*X + b))**2)
print(f"Iter {iteration}: Loss = {full_loss:.6f}")

# If loss is increasing, reduce learning rate immediately!

How do I choose the right learning rate?

There is no single "correct" value—it depends on your data's scale. Start with 0.01 after scaling features to [0,1].

If loss increases, halve it (0.005, 0.001).
If loss decreases very slowly, try 0.05.
For mini-batch GD, you may need a slightly higher rate since batch gradients are noisier.

Pro Tip: Run a quick test with [0.1, 0.01, 0.001] for just 100 iterations. Pick the one where loss decreases steadily without exploding.

Can I use this implementation for real-world predictive modeling?

For learning and prototyping—yes. For production—no, use sklearn.

Your from-scratch version correctly implements the core algorithm and will work on clean, scaled, small-to-medium datasets. However, real-world modeling requires robust handling of edge cases: automatic feature scaling, regularization (L1/L2), efficient mini-batching, and thorough validation.

sklearn.linear_model.LinearRegression (for normal equations) or SGDRegressor (for large data) are battle-tested, optimized, and integrate seamlessly into pipelines. Build your own to understand; use libraries to deliver.

What are the main assumptions of linear regression?

Linear regression assumes:

Linearity: The true relationship between each feature and target is linear.
Independence: Data points are independent (no autocorrelation, e.g., time-series).
Homoscedasticity: The error variance is constant across all prediction levels.
Normality of errors: Residuals are normally distributed (important for statistical inference, not prediction).
No multicollinearity: Features aren't highly correlated (causes unstable coefficient estimates).

Why they matter

Violating assumptions doesn't break prediction entirely, but it invalidates statistical guarantees (like p-values, confidence intervals) and often indicates a poorly specified model. For example, non-linearity (violation #1) means your model is systematically wrong—residuals will show a clear curve. Diagnose with residual plots (as shown in the Evaluation section).

How does data preprocessing affect model performance?

Preprocessing is non-optional for gradient descent. Its effects are dramatic:

Cleaning: Missing values or outliers distort gradients. A single NaN makes predictions NaN; an outlier can pull your learned line toward itself.
Scaling: Without scaling, features on different scales create an elongated loss valley. Gradient descent takes tiny, inefficient steps along shallow dimensions and may diverge along steep ones. Scaling to [0,1] makes the valley roughly spherical, allowing a single learning rate to work well.
Adding intercept: For multiple regression, ensure your design matrix X includes a column of 1s for the bias term b.

Bottom line: Preprocessing isn't "extra work"—it's part of the model. The parameters m, b you learn are only valid for the preprocessed feature space.

Is linear regression suitable for non-linear relationships?

No—by design it is unsuitable. A straight line cannot capture curves, thresholds, or interactions. If your residual plot (predicted vs. residuals) shows a clear pattern (e.g., U-shaped), that's a signature of high bias—your model is underfitting because the true pattern is non-linear.

What to do instead:

Feature engineering: Add polynomial terms (e.g., X², X³) or interaction terms (X1*X2). This is still linear regression (in the parameters), but now in a transformed feature space.
Switch models: Use decision trees, random forests, or kernel methods (like SVM with RBF kernel) that naturally capture non-linearities.
Accept the baseline: If a simple line gives decent R² (e.g., >0.7) and residuals look random, the relationship might be approximately linear enough for your needs. Interpretability is a virtue.

What are common pitfalls when implementing from scratch?

Forgetting to scale features: Gradient descent fails to converge or is painfully slow.
Using batch loss for convergence checks: With mini-batch GD, batch loss is noisy; you'll stop too early or never stop. Monitor full training loss every epoch.
Not shuffling batches: If you always iterate in the same order, the model may learn artificial sequences. Shuffle data each epoch.
Integer division in Python 2: Ensure you use float division: 2/len(X) is fine in Python 3, but 2/float(len(X)) is safer.
Mismatched array shapes: X should be 1D for simple regression; for multiple regression, ensure X is 2D (n_samples, n_features) and use matrix multiplication X @ w.
Ignoring multicollinearity: In multiple regression, correlated features cause unstable weight estimates. Check correlations and consider regularization (L2) or feature removal.

The safest approach: Test on a tiny synthetic dataset where you know the true m and b (e.g., y = 2.5*x + 1.0 + noise). Your implementation should recover parameters very close to 2.5 and 1.0. If it doesn't, debug on this controlled case before touching real data.