Linear regression draws the best-fitting straight line through your data. This linear regression Python guide takes you from the math behind that line, to gradient descent written from scratch, to scikit-learn’s LinearRegression, R-squared, residual analysis, and multiple regression on real data.

“All models are wrong, but some are useful.”

George Box, statistician

Last Updated: July 2026 | Tested on: Python 3.14.6, scikit-learn 1.9.0, NumPy 2.4.6, pandas 2.3.3 | Difficulty: Intermediate | Reading Time: 20 minutes

Here is the whole idea in one line. Linear regression finds the straight line that sits as close as possible to all your data points at once. “As close as possible” means the line that makes the total squared gap between itself and every point as small as it can be. That is it. The formulas, the gradient descent, the scikit-learn calls: those are just different ways to find that one line.

Think of stretching a single rubber band so it passes through a cloud of pins on a board. You cannot touch every pin, so you settle the band where the total tension is lowest. That resting position is your best-fit line. Linear regression does the same thing with numbers instead of tension.

It is the simplest model you will meet, and the easiest to explain to a non-technical person, which is exactly why it refuses to go away. Linear regression is the baseline that every fancier model gets measured against. If your gradient-boosted ensemble cannot beat a plain line, the line wins. It is also the seed for a lot of what comes later: logistic regression, neural networks, and regularized models all grow from this same core.

A quick story. Say an analyst named Aviraj, 27, was asked to predict monthly electricity bills from past usage. He started with linear regression and got R-squared = 0.94, so 94% of the swing in bills was explained by usage alone. His manager wanted to try a neural network. Aviraj built one. It scored 0.95, barely better, while running far slower and being impossible to explain to the finance team. The line won, and so did Aviraj.

Prerequisites

The diagram walks the whole pipeline. You start with data points, look for the line y = wx + b, measure how wrong that line is with a cost function (Mean Squared Error, MSE), then push the cost down until the line settles. You can find that minimum two ways: the normal equation solves it exactly in one shot, while gradient descent inches toward it step by step. Both land at the same trained model, which you then use to predict and evaluate.

The key picture to hold in your head is the residual: the vertical gap between a data point and the line. Every regression metric is just a different way of summarizing those gaps. MSE is the average squared gap, RMSE (Root Mean Squared Error) puts that back into real units (rupees, kilograms, whatever you measured), and R-squared tells you what slice of the variation the line managed to explain.

The Math, With Real Numbers

Think of a fruit vendor who, after weighing a few bags of apples, works out a rough price per kilo in his head. Linear regression does the same thing, only it writes the rule down as numbers. We want the line y = wx + b that best fits our data, where “best” means the smallest Mean Squared Error. Before reaching for any library, let us crank the handle by hand on a tiny dataset so you can see exactly where each number comes from. Five apartments, their area, and their monthly rent. That is the whole world for now.

📄 linear_regression_math.py: linear regression worked out by hand

import numpy as np

# Tiny dataset: apartment area (100s sqft) vs monthly rent (thousands of rupees)
X = np.array([4, 6, 8, 10, 12])    # area in 100 sqft
y = np.array([10, 15, 18, 22, 27]) # rent in thousands

n = len(X)
x_mean = X.mean()
y_mean = y.mean()

# Closed-form solution (the normal equation):
#   w = Σ(xi - x̄)(yi - ȳ) / Σ(xi - x̄)²
#   b = ȳ - w * x̄

numerator = np.sum((X - x_mean) * (y - y_mean))
denominator = np.sum((X - x_mean) ** 2)
w = numerator / denominator
b = y_mean - w * x_mean

print("Step-by-step calculation:")
print(f"  x̄ = {x_mean}, ȳ = {y_mean}")
print(f"  Σ(xi-x̄)(yi-ȳ) = {numerator}")
print(f"  Σ(xi-x̄)² = {denominator}")
print(f"  w (slope) = {numerator}/{denominator} = {w:.2f}")
print(f"  b (intercept) = {y_mean} - {w:.2f}×{x_mean} = {b:.2f}")
print(f"\nEquation: rent = {w:.2f} × area + {b:.2f}")
print(f"Interpretation: each 100 sqft adds about ₹{w*1000:.0f} in rent")

# Predict, then measure how wrong we were
predictions = w * X + b
residuals = y - predictions
mse = (residuals ** 2).mean()
ss_res = np.sum(residuals ** 2)
ss_tot = np.sum((y - y_mean) ** 2)
r_squared = 1 - ss_res / ss_tot

print(f"\nPredictions: {predictions}")
print(f"Residuals:   {residuals}")
print(f"MSE: {mse:.4f}")
print(f"R²: {r_squared:.4f} ({r_squared:.0%} of variance explained)")

▶ Output

Step-by-step calculation:
  x̄ = 8.0, ȳ = 18.4
  Σ(xi-x̄)(yi-ȳ) = 82.0
  Σ(xi-x̄)² = 40.0
  w (slope) = 82.0/40.0 = 2.05
  b (intercept) = 18.4 - 2.05×8.0 = 2.00

Equation: rent = 2.05 × area + 2.00
Interpretation: each 100 sqft adds about ₹2050 in rent

Predictions: [10.2 14.3 18.4 22.5 26.6]
Residuals:   [-0.2  0.7 -0.4 -0.5  0.4]
MSE: 0.2200
R²: 0.9935 (99% of variance explained)

What happened here: No library, nothing hidden, just the two formulas from the comment. The slope w came out to 2.05, so each extra 100 sqft adds about ₹2,050 to the monthly rent, on top of a ₹2,000 base (the intercept b). The residuals are the leftover errors at each point, and notice they are tiny and mixed (some above the line, some below), which is exactly what a good fit looks like. R-squared of 0.9935 says the line explains 99.35% of why the rents differ. Square the residuals, average them, and you get the MSE; that single number is what the next two methods will work to shrink.

From Scratch: Gradient Descent

The normal equation solves for the line in one step, but it is worth seeing the other way: gradient descent. Picture standing on a foggy hillside, wanting the lowest point. You cannot see the valley, so you feel which way the ground slopes down and take a small step that way. Repeat. That is gradient descent. The “hill” is the MSE, and each step nudges the slope and intercept downhill until the error stops dropping. We will fit salary from years of experience.

📄 linear_regression_scratch.py: gradient descent, step by step

import numpy as np

rng = np.random.default_rng(42)

# Make up 100 employees: salary ≈ 5000*experience + 30000, plus noise
n = 100
experience = rng.uniform(1, 15, n)
salary = 5000 * experience + 30000 + rng.normal(0, 5000, n)

# Scale both columns so gradient descent takes even, stable steps
X = (experience - experience.mean()) / experience.std()
y_norm = (salary - salary.mean()) / salary.std()

# Gradient descent
w, b = 0.0, 0.0
lr = 0.1
history = []

for epoch in range(100):
    y_pred = w * X + b
    mse = ((y_norm - y_pred) ** 2).mean()
    history.append(mse)

    dw = (-2 / n) * np.sum(X * (y_norm - y_pred))
    db = (-2 / n) * np.sum(y_norm - y_pred)

    w -= lr * dw
    b -= lr * db

# Convert back to original scale
w_orig = w * salary.std() / experience.std()
b_orig = salary.mean() + salary.std() * b - w_orig * experience.mean()

print(f"Gradient descent result: salary = {w_orig:.0f} × experience + {b_orig:.0f}")
print(f"Expected:               salary ≈ 5000 × experience + 30000")
print(f"MSE history: {history[0]:.4f} → {history[49]:.4f} → {history[99]:.4f}")

▶ Output

Gradient descent result: salary = 5027 × experience + 29716
Expected:               salary ≈ 5000 × experience + 30000
MSE history: 1.0000 → 0.0613 → 0.0613

What happened here: Starting from a flat line (w and b both at 0), gradient descent walked downhill for 100 rounds. The MSE history tells the story: it began at 1.0, dropped to 0.0613 by step 50, and barely moved after that, which means the line had already settled. After converting back from the scaled values, we recovered salary = 5027 × experience + 29716, almost dead-on the true 5000 and 30000 we baked into the data. We never solved an equation. We just kept stepping toward less error. That same idea, scaled up, is how neural networks learn.

With scikit-learn in 3 Lines

You have now done it the long way twice. Here is the version you will actually ship. scikit-learn wraps the whole closed-form solve behind three calls: create the model, fit it, predict with it. Same data, same idea, a fraction of the code.

📄 linear_regression_sklearn.py: the same fit in three lines

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, r2_score

rng = np.random.default_rng(42)

# Same data
n = 100
experience = rng.uniform(1, 15, n).reshape(-1, 1)
salary = 5000 * experience.ravel() + 30000 + rng.normal(0, 5000, n)

# Split
X_train, X_test, y_train, y_test = train_test_split(
    experience, salary, test_size=0.2, random_state=42
)

# Fit. Under the hood sklearn solves the normal equation, not gradient descent.
model = LinearRegression()
model.fit(X_train, y_train)

print(f"Coefficient (slope): {model.coef_[0]:.0f}")
print(f"Intercept: {model.intercept_:.0f}")
print(f"Equation: salary = {model.coef_[0]:.0f} × experience + {model.intercept_:.0f}")

# Evaluate
y_pred = model.predict(X_test)
print(f"\nTrain R²: {model.score(X_train, y_train):.4f}")
print(f"Test R²:  {model.score(X_test, y_test):.4f}")
print(f"Test RMSE: ₹{np.sqrt(mean_squared_error(y_test, y_pred)):,.0f}")

# Predict new values
new_exp = np.array([[5], [10], [15]])
predictions = model.predict(new_exp)
for exp, pred in zip(new_exp.ravel(), predictions):
    print(f"  {exp} years experience → ₹{pred:,.0f} predicted salary")

▶ Output

Coefficient (slope): 5088
Intercept: 29656
Equation: salary = 5088 × experience + 29656

Train R²: 0.9358
Test R²:  0.9454
Test RMSE: ₹4,633
  5 years experience → ₹55,096 predicted salary
  10 years experience → ₹80,536 predicted salary
  15 years experience → ₹105,976 predicted salary

What happened here: Three real lines of work (LinearRegression(), fit(), predict()) land us in the same neighborhood as the from-scratch run: a slope near 5000 and an intercept near 30000. They are not identical to the gradient descent numbers, and that is expected, because this version trains on only 80 of the 100 points (the rest are held back for testing) and draws its random noise in a slightly different order. The test R-squared (0.9454) sits right next to the training R-squared (0.9358), which is the calm signal you want: the model is not overfitting.

RMSE of ₹4,633 says a typical prediction is off by roughly ₹4,600, sensible given the ₹5,000 of noise we added on purpose.

Multiple Regression: Many Features

Think of splitting a restaurant bill by what each person actually ordered: the total is the sum of separate line items, each with its own price. Multiple regression works the same way. Real predictions rarely hang on a single number. A house price depends on size, bedrooms, age, and how far it sits from the center, all at once. Multiple regression handles this by giving each feature its own slope, then adding them up. Same math, just more of it. Each coefficient answers a clean question: holding everything else steady, how much does this one feature move the price?

📄 multiple_regression.py: one slope per feature

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split

rng = np.random.default_rng(42)
n = 300

# House price data with multiple features
df = pd.DataFrame({
    "area_sqft": rng.integers(600, 3000, n),
    "bedrooms": rng.integers(1, 5, n),
    "age_years": rng.integers(0, 30, n),
    "distance_km": rng.uniform(1, 25, n),
})

# True relationship: price = 50*area + 200000*bedrooms - 15000*age - 8000*distance + noise
df["price"] = (50 * df["area_sqft"]
               + 200000 * df["bedrooms"]
               - 15000 * df["age_years"]
               - 8000 * df["distance_km"]
               + 500000
               + rng.normal(0, 50000, n))

X = df.drop("price", axis=1)
y = df["price"]
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

model = LinearRegression()
model.fit(X_train, y_train)

print("Multiple Regression Coefficients:")
for feat, coef in zip(X.columns, model.coef_):
    print(f"  {feat}: ₹{coef:,.0f} per unit")
print(f"  Intercept: ₹{model.intercept_:,.0f}")
print(f"\nR² (train): {model.score(X_train, y_train):.4f}")
print(f"R² (test):  {model.score(X_test, y_test):.4f}")

# Predict: 1500 sqft, 3 BHK, 5 year old, 10 km from center
house = pd.DataFrame({"area_sqft": [1500], "bedrooms": [3],
                       "age_years": [5], "distance_km": [10]})
pred = model.predict(house)[0]
print(f"\n1500 sqft, 3 BHK, 5yr old, 10 km → ₹{pred:,.0f}")

▶ Output

Multiple Regression Coefficients:
  area_sqft: ₹51 per unit
  bedrooms: ₹198,980 per unit
  age_years: ₹-14,941 per unit
  distance_km: ₹-8,181 per unit
  Intercept: ₹500,648

R² (train): 0.9702
R² (test):  0.9656

1500 sqft, 3 BHK, 5yr old, 10 km → ₹1,018,074

What happened here: The model nearly nailed the hidden recipe we used to build the data: about ₹51 per sqft, ₹199K per bedroom, ₹15K knocked off for each year of age, and ₹8K off for each km from the center. Read each coefficient as a plain English sentence. Add one bedroom and the price climbs roughly ₹2 lakh; let the building age a year and it slips about ₹15K. That readability is the whole appeal. A neural network might match the accuracy, but it cannot hand you a one-line reason for every rupee. The final prediction, ₹1,018,074 (about ₹10.2 lakh), uses Python’s default comma grouping in the output; if you want the Indian lakh format you would format it yourself.

When Linear Regression Fails

Try tracing a winding river with a straight ruler and you will always cut across the bends. A straight line has the same problem. It is honest about its limits: it only works when the real relationship is roughly straight and when no wild points are yanking it around. Two situations break it fast: data that genuinely curves, and outliers. Let us watch both happen, because knowing when a tool fails is more useful than knowing when it works.

📄 assumptions.py: two ways a straight line breaks

import numpy as np
from sklearn.linear_model import LinearRegression

rng = np.random.default_rng(42)

# Case 1: a curved (sine wave) relationship that no line can follow
X = rng.uniform(0, 10, 100).reshape(-1, 1)
y_nonlinear = np.sin(X.ravel()) * 10 + rng.normal(0, 1, 100)

model = LinearRegression().fit(X, y_nonlinear)
print(f"Non-linear data: R² = {model.score(X, y_nonlinear):.3f} (poor, use polynomial or a tree)")

# Case 2: a clean straight relationship, then two stray points dropped in
X_clean = rng.uniform(0, 10, 98).reshape(-1, 1)
y_clean = 3 * X_clean.ravel() + 5 + rng.normal(0, 2, 98)

# Two big apartments whose rent came out near zero (bad data, or a steal)
X_outlier = np.vstack([X_clean, [[9.5], [9.5]]])
y_outlier = np.append(y_clean, [0, 0])

model_clean = LinearRegression().fit(X_clean, y_clean)
model_outlier = LinearRegression().fit(X_outlier, y_outlier)
print(f"\nWithout outliers: slope={model_clean.coef_[0]:.2f}, R²={model_clean.score(X_clean, y_clean):.3f}")
print(f"With 2 outliers:  slope={model_outlier.coef_[0]:.2f}, R²={model_outlier.score(X_outlier, y_outlier):.3f}")
print("Just two stray points dragged the slope down and knocked R² from 0.95 to 0.72.")

▶ Output

Non-linear data: R² = 0.000 (poor, use polynomial or a tree)

Without outliers: slope=2.97, R²=0.950
With 2 outliers:  slope=2.64, R²=0.723
Just two stray points dragged the slope down and knocked R² from 0.95 to 0.72.

What happened here: In case 1 the data rises and falls like a wave, and a single straight line cannot chase that shape, so R-squared collapses to essentially 0. The fix is to give the model some curve: polynomial features or a tree-based model. Case 1 is why you always plot your data before trusting a line. In case 2 the data really is straight, but just two bad points (big apartments listed at almost zero rent) tugged the slope from 2.97 down to 2.64 and dropped R-squared from 0.95 to 0.72. Least squares squares every error, so a couple of far-off points carry enormous weight.

That is the headline weakness: linear regression is not robust to outliers. Clean your data, or reach for a more robust method, before you trust the line.

Common Mistakes

Mistake 1: Judging a model by R-squared alone

R-squared tells you the fraction of variation explained, but it says nothing about the size of your errors in real units. Two models can share almost the same R-squared while being off by ₹40,000 in one case and 0.3 of a point in the other. Always report RMSE next to it.

❌ Mistake: trusting R-squared without RMSE

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score, root_mean_squared_error

rng = np.random.default_rng(0)

# Same data shape, very different scales
X = rng.uniform(0, 10, 200).reshape(-1, 1)

# Model A: salaries in rupees
y_big = 50000 * X.ravel() + 100000 + rng.normal(0, 40000, 200)
mb = LinearRegression().fit(X, y_big)
print(f"Salary model:  R² = {r2_score(y_big, mb.predict(X)):.3f}, "
      f"RMSE = ₹{root_mean_squared_error(y_big, mb.predict(X)):,.0f}")

# Model B: a 0 to 5 rating
y_small = 0.4 * X.ravel() + 1 + rng.normal(0, 0.3, 200)
ms = LinearRegression().fit(X, y_small)
print(f"Rating model:  R² = {r2_score(y_small, ms.predict(X)):.3f}, "
      f"RMSE = {root_mean_squared_error(y_small, ms.predict(X)):.2f}")
print("Similar R², wildly different error. R² alone hides the real miss size.")

▶ Output

Salary model:  R² = 0.929, RMSE = ₹40,834
Rating model:  R² = 0.943, RMSE = 0.30
Similar R², wildly different error. R² alone hides the real miss size.

Why this bites: The two models look almost equally good through the R-squared lens (0.929 and 0.943), yet one is off by tens of thousands of rupees and the other by a third of a point. R-squared is unit-free, so it cannot warn you that a “95% accurate” salary model still misses by ₹40K. RMSE speaks in your data’s own units, so report both, every time.

Mistake 2: Fitting a line without looking at the data first

❌ Wrong

# Curved data? A line will "run" but the fit is meaningless.
model = LinearRegression().fit(X, y_curved)   # no plot, no residual check
print("R² =", model.score(X, y_curved))        # could be ~0 and you would not know why

✅ Correct

# Plot the scatter first, then plot residuals after fitting.
# If the data curves, add polynomial features or switch to a tree model.
# A line is only honest when the relationship is roughly straight.

Why: linear regression always returns a line, even for data shaped like a wave. The code does not raise an error, so the only warning you get is a low R-squared or a residual plot with an obvious pattern. Look at the scatter before you fit, and look at the residuals after. A pattern left in the residuals means the line missed something the data was trying to tell you.

Practice Exercises

1. Exercise 1: Build your own tiny dataset of 6 points, fit LinearRegression, and plot the data plus the best-fit line with matplotlib. Mark each residual as a vertical line.
2. Exercise 2: Take clearly curved data (try y = X**2 plus noise), then use PolynomialFeatures to compare R-squared for degrees 1 through 5. Watch where it stops improving.
3. Exercise 3: Reuse the gradient descent code above, but try learning rates of 0.001, 0.1, and 1.5. Print the MSE history for each and explain what too small and too large look like.

Conclusion

You built linear regression three ways and they all agreed. You saw the two closed-form formulas produce a line by hand, watched gradient descent walk downhill to the same answer, then let scikit-learn do it in three lines. You learned to read a model out loud through its coefficients, to trust R-squared and RMSE together instead of either alone, and to spot the two things that break a straight line: curved data and outliers. That is a complete, honest mental model of the most widely used algorithm in machine learning.

Next up is polynomial and regularized regression (Ridge and Lasso), where you give the line permission to curve without letting it overfit. Everything you just learned carries straight over. If you want the full path from Python basics to deep learning in order, start at the Python + AI/ML tutorial series home.

Frequently Asked Questions

Interview Questions on Linear Regression

These come from real screens and onsites. Practice answering before you read each answer.

Q: What does the slope coefficient in a linear regression mean, in plain words?

It is the amount the prediction changes when that one feature goes up by a single unit, holding every other feature fixed. If a house-price model gives area a coefficient of 51, each extra square foot adds about 51 rupees to the predicted price. The sign matters too: a negative coefficient (like age at -15000) means the target falls as that feature rises.

Q: You fit LinearRegression on a million rows. What algorithm actually computes the weights, and at what point would you switch to SGDRegressor?

A direct closed-form least-squares solve: under the hood scikit-learn calls scipy.linalg.lstsq, an SVD-based solver, so there is no learning rate or iteration count to tune. It avoids the textbook normal equation because forming XᵀX is numerically ill-conditioned. That is why it is fast and reliable on small to medium data. For very large or streaming datasets where the matrix math gets expensive, you switch to SGDRegressor, which does use gradient descent.

Q: Why report RMSE alongside R-squared instead of R-squared alone?

R-squared is unit-free and only tells you the fraction of variation explained, so two models with nearly identical R-squared can miss by tens of thousands of rupees versus a fraction of a point. RMSE speaks in the data’s own units and tells you the typical size of an error. One gives relative quality, the other gives the real-world miss, so you want both.

Q: What does a negative R-squared tell you?

It means the model predicts worse than a flat line that always outputs the mean of y. A mean-only baseline scores exactly 0, so anything below 0 says your fit is actively harmful, usually because the relationship is not linear or the model is applied to the wrong data. It is a red flag to plot the data and rethink, not to tune harder.

Q: You fit a linear model and R-squared on the training data is 0.02, almost zero. What do you check first?

Plot the scatter before anything else. A near-zero R-squared usually means the true relationship is curved (a wave or a U-shape) and a straight line cannot follow it, exactly the sine-wave case in this post. If the shape is curved, add polynomial features or move to a tree-based model. Also confirm you did not accidentally feed a constant or shuffled target.

Q: A colleague adds two new rows to the dataset and the slope suddenly shifts a lot while R-squared drops. What likely happened?

Those two rows are probably outliers. Least squares squares every error, so a couple of far-off points carry huge weight and can drag the whole line toward them, just like the two zero-rent apartments that pulled the slope from 2.97 to 2.64 here. Inspect the extremes, remove them if they are genuine data errors, or switch to a robust method like RANSACRegressor.

Q: Why do you scale features before running gradient descent, but not for the normal equation?

Gradient descent takes steps proportional to each feature’s scale, so if one feature is in thousands and another in single digits, the steps become lopsided and convergence is slow or unstable. Standardizing both to a similar range makes the descent smooth and even, which is why the from-scratch example normalized the columns. The normal equation solves the whole thing algebraically in one step, so it does not care about scale.

Q: When would you pick linear regression over a neural network or gradient-boosted model?

When the relationship is roughly linear and you need to explain every prediction, linear regression is hard to beat: it is fast, needs little data, and each coefficient is a plain sentence a business team can follow. Use it as the baseline first, and only move to a heavier model if it clearly beats the line by enough to justify the lost interpretability and speed. If a boosted model barely edges out the line, the line usually wins.

Series: Python + AI/ML Cookbook, Part 5: Machine Learning

Reference: the complete, always-current details live in scikit-learn documentation.

Related Posts

Previous: ML: The Math Behind Machine Learning, Intuition Before Formulas

Next: ML: Polynomial & Regularization (Ridge, Lasso)

Series Home: Python + AI/ML Tutorial Series