The Ultimate Guide to Linear Regression:
From Theory to Production

What Is Linear Regression — and Why Should You Care?

Linear regression is the bedrock of supervised machine learning. Before you can understand neural networks, gradient boosting, or any modern deep learning architecture, you need to truly grasp what regression is doing mathematically and computationally.

More importantly, linear regression is not a "beginner" tool you graduate away from. In production environments, it remains one of the most frequently deployed models because it's interpretable, fast, and surprisingly robust when features are well-engineered.

At its core, linear regression answers a specific kind of question: given a set of input features X, can we learn a linear function that predicts a continuous output y? For example:

• Predict housing prices from square footage, location, and age
• Forecast quarterly sales revenue from marketing spend and seasonality
• Estimate patient readmission risk from clinical biomarkers

The Mathematical Foundation

The linear regression hypothesis is straightforward. For a dataset with n features, we predict Å· as a weighted sum of inputs:

Å· = β₀ + β₁x₁ + β₂xâ‚‚ + ... + βâ‚xâ‚

Where β₀ is the intercept (bias) and β₁ through β₠are the coefficients (weights) we want to learn. In matrix notation, this simplifies elegantly to:

ŷ = Xβ

Our goal is to find the vector β that minimizes prediction error. The standard metric is Mean Squared Error (MSE):

MSE = (1/n) * Σ(yᵢ - ŷᵢ)²

Implementation from Scratch with NumPy

Rather than jumping straight to Scikit-Learn, let's implement linear regression from scratch using only NumPy. This exposes the mechanics that every ML library is hiding under the hood.

import numpy as np
import pandas as pd
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

# Load dataset
housing = fetch_california_housing(as_frame=True)
X = housing.data.values
y = housing.target.values

# Split and scale
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

class LinearRegression:
    """
    Linear Regression implemented from scratch using
    the Normal Equation and Gradient Descent.
    """
    def __init__(self, lr=0.01, n_iter=1000):
        self.lr = lr
        self.n_iter = n_iter
        self.weights = None
        self.bias = None
        self.history = []

    def fit(self, X, y):
        n_samples, n_features = X.shape
        self.weights = np.zeros(n_features)
        self.bias = 0

        for i in range(self.n_iter):
            y_pred = self._forward(X)
            loss = np.mean((y - y_pred) ** 2)
            self.history.append(loss)

            # Compute gradients
            dw = -(2/n_samples) * X.T @ (y - y_pred)
            db = -(2/n_samples) * np.sum(y - y_pred)

            # Update parameters
            self.weights -= self.lr * dw
            self.bias -= self.lr * db

        return self

    def _forward(self, X):
        return X @ self.weights + self.bias

    def predict(self, X):
        return self._forward(X)

    def score(self, X, y):
        y_pred = self.predict(X)
        ss_res = np.sum((y - y_pred) ** 2)
        ss_tot = np.sum((y - np.mean(y)) ** 2)
        return 1 - (ss_res / ss_tot)  # R²

# Train and evaluate
model = LinearRegression(lr=0.05, n_iter=2000)
model.fit(X_train_scaled, y_train)

print(f"Training R²: {model.score(X_train_scaled, y_train):.4f}")
print(f"Test R²:     {model.score(X_test_scaled, y_test):.4f}")

Running this code on the California Housing dataset, we achieve an R² of approximately 0.61 on the test set — respectable for a vanilla linear model with no feature engineering.

Production-Ready with Scikit-Learn

In most production scenarios, you'll use Scikit-Learn's optimized implementation. It uses a more numerically stable algorithm (SVD decomposition) and integrates cleanly with the broader ML ecosystem.

from sklearn.linear_model import LinearRegression, Ridge, Lasso
from sklearn.metrics import mean_squared_error, r2_score
import numpy as np

# Standard OLS
lr = LinearRegression()
lr.fit(X_train_scaled, y_train)

y_pred = lr.predict(X_test_scaled)
rmse = np.sqrt(mean_squared_error(y_test, y_pred))
r2 = r2_score(y_test, y_pred)

print(f"RMSE: {rmse:.4f}")
print(f"R²:   {r2:.4f}")
print(f"Coefficients: {lr.coef_}")
print(f"Intercept:    {lr.intercept_:.4f}")

When to Use Ridge vs Lasso vs Plain OLS

Pure OLS minimizes MSE with no constraints on coefficient size. This works well when features are few and informative. But in practice, you often face:

• Multicollinearity — features that are highly correlated, making coefficient estimates unstable
• High dimensionality — many features relative to samples, leading to overfitting
• Irrelevant features — noise variables that inflate test error

Ridge (L2 regularization) adds a penalty on the sum of squared coefficients, shrinking all weights toward zero but keeping all features. Use Ridge when you believe most features are relevant but their relationships might be noisy.

Lasso (L1 regularization) adds a penalty on the absolute sum of coefficients and performs implicit feature selection by forcing some weights to exactly zero. Use Lasso for sparse feature selection in high-dimensional data.

Deploying Your Model with FastAPI

A model that lives in a notebook is not a data science solution — it's an experiment. Here's a minimal but production-quality deployment using FastAPI and joblib:

import joblib
from fastapi import FastAPI
from pydantic import BaseModel
import numpy as np

# Save model
joblib.dump({'model': lr, 'scaler': scaler}, 'linear_model.pkl')

# FastAPI app
app = FastAPI(title="Housing Price Predictor", version="1.0")

class HousingFeatures(BaseModel):
    MedInc: float; HouseAge: float; AveRooms: float
    AveBedrms: float; Population: float; AveOccup: float
    Latitude: float; Longitude: float

artifacts = joblib.load('linear_model.pkl')

@app.post("/predict", summary="Predict housing price")
def predict(features: HousingFeatures):
    X = np.array([[
        features.MedInc, features.HouseAge, features.AveRooms,
        features.AveBedrms, features.Population, features.AveOccup,
        features.Latitude, features.Longitude
    ]])
    X_scaled = artifacts['scaler'].transform(X)
    prediction = artifacts['model'].predict(X_scaled)[0]
    return {
        "predicted_price_usd": round(prediction * 100_000, 2),
        "model": "OLS Linear Regression",
        "version": "1.0"
    }

With this, you have a fully functional REST API serving real-time housing price predictions. Run it with uvicorn app:app --reload and test it at /docs.